This is an interesting question that popped through my mind. Some of us should know what is meant by „gauge transformations”, „gauge invariance/symmetry” and are used to seeing these terms whenever lectures on quantum field theory are read. But the electromagnetic field in vacuum (described in a specially relativistic fashion by the tensor ##F_{\mu\nu}##) is a classical one (i.e. it exists also in a non quantum setting), so one has the right to ask. Given E,B or ##F##, what does it mean to say about the "classical electromagnetism" to be U(1) invariant?

You really should be working in terms of the 4-potential and not the field strengths as this is your dynamical field. The gauge invariance you are looking for is the same as that typically introduced in classical electromagnetism. Without a gauge fixing condition your 4-potential will not be uniquely determined (the gauge symmetry being a local U(1) symmetry). This should be described in detail (not only for EM but also for general non-commutative gauge fields) in any textbook covering classical Yang-Mills theory.

A very good book is vol. III of Scheck's textbook on theoretical physics, where he treats electromagnetism right away from the modern point of view, which is, of course, that it is describing a massless spin-1 field, which is necessarily a gauge field (if you don't want unobserved continuous intrinsic spin-like degrees of freedom).

Classcial electrodynamics is gauge invariant also under the presence of sources (electric charge-current distributions).

As Orodruin said, it's described in terms of the potentials of the field, which is introduced using the homogeneous Maxwell equations
$$\vec{\nabla} \cdot \vec{B}=0, \quad \vec{\nabla} \times \vec{E} + \frac{1}{c} \partial_t \vec{B}=0.$$
From the first equation you introduce the vector potential
$$\vec{B} = \vec{\nabla} \times \vec{A}.$$
Plugging this into the 2nd equation yields
$$\vec{\nabla} \times \left (\vec{E}+\frac{1}{c} \partial_t \vec{A} \right)=0,$$
which means that the expression in the parantheses is a gradient:
$$\vec{E}=-\vec{\nabla} \Phi - \frac{1}{c} \partial_t \vec{A}.$$
Now, obviously for given ##\vec{E}## and ##\vec{B}## you can add to ##\vec{A}## an arbitrary gradient,
$$\vec{A}'=\vec{A} - \vec{\nabla} \chi,$$
and then in order to get also the same ##\vec{E}## you must have
$$\vec{E}=-\nabla \Phi'-\frac{1}{c} \partial_t \vec{A}'=-\nabla (\Phi'- \frac{1}{c} \partial_t \chi') - \frac{1}{c} \partial_t \vec{A} \; \Rightarrow \; \Phi'=\Phi+ \frac{1}{c} \partial_t \chi.$$
Thus everything must be unchanged under the gauge transformation
$$\vec{A}'=\vec{A}-\vec{\nabla} \chi, \quad \Phi'=\Phi + \frac{1}{c} \partial_t \chi.$$
The inhomogeneous Maxwell equations are
$$\vec{\nabla} \times \vec{B} - \frac{1}{c} \partial_t \vec{E}=\frac{1}{c} \vec{j}, \quad \vec{\nabla} \cdot \vec{E}=\rho,$$
or
$$\vec{\nabla} \times (\vec{\nabla} \times \vec{A})+\frac{1}{c} \left (\vec{\nabla} \Phi +\frac{1}{c} \partial_t \vec{A} \right) = \frac{1}{c} \vec{j}, \quad -\Delta \Phi - \frac{1}{c} \partial_t \vec{\nabla} \cdot \vec{A}=\rho.$$
Using
$$\vec{\nabla} \times (\vec{\nabla} \times \vec{A})=\vec{\nabla} (\vec{\nabla} \cdot \vec{A})-\Delta \vec{A},$$
you can see that it is advantegeous to make use of your gauge freedom by partially fixing the gauge field ##\chi## introducing the Lorenz-gauge condition (note it's Lorenz here, not Lorentz for historical justice!):
$$\frac{1}{c} \partial_t \Phi + \vec{\nabla} \cdot \vec{A}=0,$$
because then the four vector fields decouple into
$$\Box \Phi=\rho, \quad \Box \vec{A}=\frac{1}{c} \vec{j},$$
where the d'Alembert operator is devined as
$$\Box=\frac{1}{c^2} \partial_t^2 - \Delta.$$
Note that for this to be consistent without any reference to matter-field (or classical point-mechanics) equation of motion, you must have the continuity equation for electric charge as an integrability condition,
$$\partial_t \rho + \vec{\nabla} \cdot \vec{j}=0.$$

Another important point of gauge invariance is that the variation of the action for the motion of charged particles,
$$S=\int \mathrm{d} t [-m c^2 \sqrt{1-\dot{\vec{x}}^2/c^2} + q \Phi - q \dot{\vec{x}} \cdot \vec{A}]$$
is gauge invariant.

It's also very easy to see that verything can be made manifestly covariant by making
$$(A^{\mu})=(\Phi,\vec{A}), \quad (j^{\mu})=(c \rho,\vec{j})$$
to four-vectors, leading to the covariant form of Maxwell's equations
$$F_{\mu \nu} = \partial_{\mu} A_{\nu}-\partial_{\nu} A_{\mu}, \quad \partial_{\mu} F^{\mu \nu}=\frac{1}{c} j^{\nu}.$$
The obvious convenience of the Lorenz gauge comes not the least from the fact that it's a Lorentz-invariant constraint,
$$\partial_{\mu} A^{\mu} =0.$$

The full beauty of gauge invariance comes into view in quantum theory, where you can heuristically argue for the introduction of vector fields by making global intrinsic symmetries like the invariance of the Schrödinger, Klein-Gordon, or Dirac equation under changes of the phase of the fields (building the symmetry group U(1)) local.

The electromagnetic four-potential is the affine connection of local U(1) symmetry. I think, this becomes obvious in a natural way only in quantum theory.

I would say it is also rather evident in the structure of classical Yang-Mills theory as it is appearing in precisely the appropriate manner in the covariant derivative. Specialising that to a U(1) symmertry of course directly gives you classical electromagnetism.

Because from the point of view of classical physics, there's no physical reason to consider several "species" of electromagnetic potentials.

That it does not describe something we observe does not necessarily make it uninteresting. In addition, classical EM is a (Abelian) classical Yang-Mills theory. Calling YM theory uninteresting classically is calling EM uninteresting classically.

using only the electromagnetic fields and potentials, how do I know that the gauge symmetry is U(1)?

I don't think the detour into Yang-Mills was helpful. So let's back up a bit.

The central idea of group theory in physics is that the same equations have the same solutions, and the encountered symmetries are reflected in these solutions. Thus, it makes sense to use these symmetries as labels in categorizing these equations (and their solutions).

Ultimately, a human being needs to be able to identify the symmetries - just as it takes a human being to notice that an equilateral triangle is symmetric about 120 degree rotation, and symmetric about exchange of any two sides. Some human being has to look at the equations of motion and say "the properties of X mean it has the Y symmetry". There is no "operator" that you feed equations of motion and out pops "U(1)".

These effects are not subtle.

Take E&M, and add the light-by-light scattering term E dot B into the classical Lagrangian. All of a sudden ε0 and μ0 stop being numbers and start being tensors (and problems quickly become unsolvable).If I go to an SU(2) theory, I not only have that, but my fields themselves become charged, and charge is no longer expressible as a number - it needs to be a matrix. The equations of motion will be very different than the Maxwell equations, and an expert can recognize them as "aha - that's what SU(2) looks like.:"

Since the (EM) gauge transformations are parameterised by the arbitrary real number [itex]\Lambda[/itex], the underlying group (if any) must be a Lie group, call it [itex]G[/itex]. Denote its Lie algebra by [itex]\mathcal{L}(G) \cong T_{e}(G)[/itex]. Now, the infinitesimal transformation [itex]\delta_{\Lambda}A_{\mu} = \partial_{\mu}\Lambda[/itex] implies that [tex][\delta_{\Lambda} , \delta_{\Omega}] A_{\mu} = 0 .[/tex] This means that [itex]\mathcal{L}(G)[/itex] is an Abelian Lie algebra. Since [itex]\Lambda : \mathbb{R}^{(1,3)} \to \mathbb{R}[/itex], then [itex](i \Lambda)[/itex] must map an open subset of space-time into the complex line [itex]i \mathbb{R}[/itex]. But, we know that [itex]\mathcal{L}\left(U(1)\right) \cong i \mathbb{R}[/itex]. We therefore conclude that [itex]i \Lambda \in \mathcal{L}\left(U(1)\right)[/itex]: indeed, topologically [itex]U(1)[/itex] may be identified with circle [itex]S^{1}[/itex], and the tangent space at the identity is obtained by differentiating the curves [itex]t \to e^{i t \theta}[/itex] at [itex]t = 0[/itex] for all real numbers [itex]\theta (x)[/itex], giving as the Lie algebra the complex line [itex]i \mathbb{R}[/itex]; thus [itex]\mathcal{L}\left(U(1)\right) \cong i \mathbb{R}[/itex]. Now, exponentiation (of the Lie algebra element [itex]i \Lambda[/itex]) leads to [tex]g(x) \equiv e^{i \Lambda (x)} \in U(1) .[/tex] This allows us to rewrite the EM-gauge transformation as [tex]A_{\mu} \to A_{\mu} + i g(x) \partial_{\mu} g^{-1}(x), [/tex] where [itex] g \partial_{\mu} g^{-1} \in \mathcal{L}\left(U(1)\right)[/itex].

Of course a proper (but very lengthy) way to answer your question, is by deriving the Maxwell’s equations from the action of [itex]U(1)[/itex](considered as the set of all sections which are vector bundle isomorphisms) on the fibre bundle over [itex]\mathcal{M}^{(1,3)}[/itex].

[...]
Of course a proper (but very lengthy) way to answer your question, is by deriving the Maxwell’s equations from the action of [itex]U(1)[/itex](considered as the set of all sections which are vector bundle isomorphisms) on the fibre bundle over [itex]\mathcal{M}^{(1,3)}[/itex].

Do you, perhaps, have a reference for this (book, article)? I'm interested in this calculation, too.
Thank you very much for your effort.

Do you, perhaps, have a reference for this (book, article)? I'm interested in this calculation, too.
Thank you very much for your effort.

I believe, many textbooks on differential geometry do explain the [itex]U(1)[/itex]-connection and its curvature on [itex]\mathcal{M}^{(1,3)}[/itex]. The most elementary textbook which treats the subject is

But, your question has a simple one-sentence-answer: If the current is conserved, i.e., [itex]\mbox{d}j = 0[/itex], then the following Maxwell’s action [tex]S[A] = \int_{\mathbb{R}^{(1,3)}} \left( \frac{1}{2} F \wedge ~^*\!F + j \wedge A \right) ,[/tex] is invariant under [tex]A \to A + e^{\Gamma} \mbox{d} e^{-\Gamma} ,[/tex] where [itex]e^{\Gamma}[/itex] is a function on spacetime with values in the group [itex]U(1)[/itex]. Clearly, the above answer is a lot easier than showing that the [itex]U(1)_{p \in \mathcal{M}}[/itex] transformation of the fibre [itex]E_{p}[/itex] does not depend on the local trivialization for the vector bundle at [itex]p \in \mathcal{M}[/itex].

The fiber bundle formulation of classical field theory has always been a subject that gave me headaches which could not be treated by quite a number of books. For example, the simplest question I still have as „unanswered” is: OK, even if F,j,A „live” in a (vector, principle, associated) bundle, their functional dependence is still on ##x^{\mu}##, that is a parameter set for the base space (4D Minkowski spacetime). Therefore, when you write ##dj = 0##, is that d the exterior differential in the gauge bundle, or the exterior differential in the de Rham bundle over space-time?

what does it mean to say about the "classical electromagnetism" to be U(1) invariant?

Here is my try. Classical electromagnetism (without matter) is invariant under the transformation
$$A_{\mu}\rightarrow A_{\mu}+\partial_{\mu}\lambda.$$
Since ##\lambda## is a single continuous function, this is a symmetry under a one-parameter continuous transformation. This means that the group of transformation is a one-parameter continuous group. Every one-parameter continuous group is locally isomorphic to U(1). If, in addition, we require that the group must be compact, I think that U(1) is the only possibility even globally.