This is indeed a very confusing question and I have spent a lot of time parsing the literature looking for an answer. The references I found most useful are section 3 of the paper by Harlow and Ooguri, a set of lectures notes, and a study of SSB in superconductors(1) by van Wezel and van den Brink. Below, I summarize my current understanding; any comments and corrections from gauge theory experts are welcome. Apologies for the long answer, because of all the confusion, I thought it's useful to be careful and detailed. The upshot is that the precise answer depends on a choice of boundary conditions, but I think there is a way in which you can consistently say that there is spontaneous symmetry breaking (SSB) in the gauge theory.
Global vs local symmetries
Firstly, as pointed out in Ruben's answer, a breaking of local gauge symmetry is ruled out both by Elitzur's theorem and by the fact that these are un-physical, 'do-nothing' transformations (we require states in the physical Hilbert space do be invariant under them). The symmetry that is supposed to be broken in the ground state is the global U(1) symmetry. It's important to clarify the relationship between these. Local gauge transformations are defined by the requirement that they tend to the identity at (spatial) infinity. The simplest examples are those that act only on some finite region. The global symmetry operations are those that do not vanish at infinity. More precisely, the global symmetry group is the factor of the group of all gauge transformation by the local ones: $G_\text{global} = G_\text{gauge} / G_\text{local}$. Note that by this definition, $G_\text{global}$ is not a subgroup of $G_\text{local}$
The global symmetry transformations thus defined act entirely on the boundary at spatial infinity (for this reason, they are also called 'asymptotic symmetries'). There is a good physical reason for that. We expect the global symmetry to be the one generated by the total electric charge of the system, $Q_\text{tot}$. But by Gauss's law, $Q_\text{tot}$ can be measured by looking at the outgoing electric flux at infinity. Another visualization is that charged objects always exist at the endpoints of Wilson lines; if $Q_\text{tot} \neq 0$ then the Wilson lines of 'unpaired' charges extend all the way to infinity and you can evaluate the total charge by counting them.
Choices of boundary conditions
Is $G_\text{global}$ a physical symmetry that can be spontaneously broken? In some sense, the answer is a matter of choice. As Harlow and Ooguri discuss (p. 58), when varying the EM action, one picks up a boundary term which is a product of the variation of the gauge field $A$ and the electric flux leaving the system(2). There are therefore two ways of ensuring gauge invariance:
- You can require that there is no outgoing flux, in which case you can freely vary the gauge field. This is equivalent to enforcing that $Q_\text{tot} = 0$ for all physical states.
- Alternatively, if you want to allow configurations with finite charge in your theory, then you have to restrict variations of $A$ to those that vanish on the boundary.
Both choices are still consistent with requiring physical states to be invariant under $G_\text{local}$.
The picture becomes clearer in lattice gauge theory (p. 66-67 in Harlow and Ooguri). One has to make a choice of boundary conditions, whether to enforce the Gauss constraint on the 'boundary sites', i.e. endpoints of edges 'sticking out' of the lattice (see their Fig. 9 for illustration). Doing so is equivalent to option 1 above, because in this case, Wilson lines cannot end at the boundary so $Q_\text{tot} = 0$ by construction. But you can also choose not to enforce the constraint on these boundary sites, in which case you keep a larger Hilbert space that has physical states with $Q_\text{tot} \neq 0$. From these, you could in principle construct the states with fixed phase that would describe spontaneous symmetry breaking. States with different phases would be physically distinct states in the Hilbert space, leading to a degenerate vacuum in the thermodynamic limit just like in a theory with only global $U(1)$ symmetry.
Dirac order parameter
As noted before, Elitzur's theorem precludes any truly local order parameter, since if one existed, it would also be non-invariant under $G_\text{local}$ and hence non-physical. Instead, we have the Dirac order parameter (Eq. 7.5 in the lectures notes)
$\Psi_D(\mathbf{r}) \equiv \Psi(\mathbf{r}) e^{i \int d^d\mathbf{r'} \mathbf{E}_\text{cl}(\mathbf{r}-\mathbf{r'})A(\mathbf{r'})}$
Here $\Psi(\mathbf{r})$ is the local matter (e. g. Cooper pair) field while $\mathbf{E}_\text{cl}(\mathbf{r})$ is a classical function corresponding to the electrostatic field of a point charge, i.e. it satisfies the equation $\nabla \cdot \mathbf{E}(\mathbf{r}) = \delta(\mathbf{r})$. This is clearly non-local; however, there is a gauge choice (Coulomb gauge, $\nabla \cdot \mathbf{A} = 0$) where the expression in the exponential vanishes and $\Psi_D$ just looks like the local fermion field.
Note that $\Psi_D$ is invariant under $G_\text{local}$ but not under $G_\text{global}$ (if you don't require $\mathbf{A}$ to vanish at infinity, you pick up a boundary term in the exponent under gauge transformation), so it is a good order parameter for the putative phase where this latter symmetry is spontaneously broken. Indeed, given the 'asymptotic' nature of this symmetry, it's unsurprising that its order parameter has to be non-local in this way.
A closer look at SSB
This might still seem a bit unsatisfying. There are at least two possible complaints:
- I could have defined my theory on a manifold without a boundary (e.g. by taking periodic boundary conditions), in which case the choice above does not arise and I'm forced to have $Q_\text{tot}=0$.
- Even if I take the case where all charges are allowed, there is a theorem in QFT showing that no physical (=local and gauge invariant) operator can change $Q_\text{tot}$(3). This implies that there is a superselection rule: superpositions of states with different charges are physically meaningless.
A hint that these objections are not fatal is given by the fact that the latter one arises also for superfluids, where there is no gauge redundancy, and in quantum optics, where the photon number is not even conserved.
To resolve the issues, one has to consider more carefully what we mean by spontaneous symmetry breaking. This is done in the latter two references cited at the beginning(4). The main point is that, whenever we say that a system spontaneously breaks a symmetry, we need to have some other, external reference system in mind to which we can compare it. This is true of any symmetry breaking, independently of whether we consider a gauge theory or not. So both for a superconductor and a superfluid, when we assign a phase to it, what we really mean by it is the relative phase with some previously chosen reference system(5). Another way of saying this is that while two decoupled superfluids/superconductors are in fact symmetric under changing their phases independently, an infinitesimal coupling(6) between them is enough to reduce this to the case when the ground state is only invariant under transforming both at the same time.
This gets rid of both complaints above. The more natural thing to consider for SSB is a system that does have a boundary, with the reference system being on the 'outside'; and taking the reference system into account removes the superselection rule. These considerations are the same whether there are gauge fields involved or not. So I think it is fair to say that the $U(1)$ gauge theory has SSB just as much as system with only global $U(1)$ symmetry does (and indeed, in their paper, van Wezel and van den Brink explicitly check that there is a tower of (physical) states that becomes a degeneracy in the thermodynamic limit).
(1) For those less comfortable with condensed matter jargon, feel free to replace 'superconductor' with '$U(1)$ gauge theory' everywhere. In fact, I'm cheating a bit: I have in mind a system where the matter fields that condense carry charge $e$, instead of $2e$ as they wound in an actual superconductor. Near the end, I also mention superfluids for comparison, which you should think of as 'some theory with spontaneously broken global $U(1)$ symmetry (and no gauge fields)'.
(2) Formally, it reads $\delta S_\text{boundary} \propto \int_{\partial M} \delta A \wedge \star F$, where $\partial M$ is the boundary of the manifold $M$ on which the theory is defined, and $\star F$ is the Hodge dual of the EM field tensor $F$. This latter object is the electric flux, as can be shown using Maxwell's equations.
(3) The intuition for this should be clear from the discussion above: to change $Q_\text{tot}$, you need to bring in a charge from infinity.
(4) In particular, there is a nice discussion of this in Sec. 1.5.3 of the lecture notes on SSB.
(5) This relative phase is clearly physical, since it leads to a measurable Josephson current. As pointed out in the lecture notes, this current is proportional to the correlation between the Dirac order parameter, $\langle \Psi_D(\mathbf{r})\Psi_D^\dagger(\mathbf{r'})\rangle$, in the case when the two points, $\mathbf{r}, \mathbf{r'}$, are located in the two different superconductors.
(6) By 'infinitesimal' I mean the same thing as in the definition of SSB by adding an infinitesimal symmetry-breaking field. So the precise statement is that if we first take the size of both superconductors to infinity and then take the coupling between them to zero, the remaining symmetry is still only $U(1)$ and not $U(1) \times U(1)$.