Spontaneously Broken Subsystem Symmetries

We investigate the spontaneous breaking of subsystem symmetries directly in the context of continuum field theories by calculating the correlation function of charged operators. Our methods confirm the lack of spontaneous symmetry breaking in some of the existing continuum field theories with subsystem symmetries, as had previously been established based on a careful analysis of the spectrum. We present some novel continuum field theory constructions that do exhibit spontaneous symmetry breaking whenever allowed by general principles. These interesting patterns of symmetry breaking occur despite the fact that all the theories we study are non-interacting.


Introduction
Fracton lattice models [1][2][3][4][5] have forced a rethinking of some of the longest held believes about quantum field theories. Standard lore would posit that, at very long distances, any quantum field theory should either flow to a gapped system, potentially described by a non-trivial topological field theory (TFT), or a gapless system described by a conformal field theory (CFT). Fractonic lattice models do not fit into either of these frameworks; their limited mobility excitations as well as their system size dependent groundstate degeneracy is not describably by any standard TFT or CFT.
A genuine continuum field theory framework capable of describing fractonic physics has been laid out in [6][7][8]. The major new ingredient in these field theories are subsystem symmetries: global symmetries that do not act on the entire system but only on degrees of freedom localized along certain sub-manifolds. Such subsystem symmetries are often present in fractonic lattice models; for the exotic field theories constructed in [6][7][8] they become a definining feature. While some of these new field theories don't even exhibit genuine fractonic excitations, they all have some kind of subsystem symmetry.
At face value these exotic field theories look somewhat conventional. The simplest examples are free field theories; in fact for many examples the symmetries prevent any relevant or marginal interactions. On the surface the only new features are unusual kinetic terms. For example, in the simplest 2+1 dimensional model of [7] based on a single real scalar, the standard (∂ x φ) 2 + (∂ y φ) 2 kinetic term gets replaced with (∂ x ∂ y φ) 2 . This is the lowest dimension kinetic term consistent with a subsystem shift symmetry demanding the action be invariant under shifts of φ that can depend on either x or y, δφ = f x (x) + f y (y). This kinetic term, however, leads to dramatic consequences. Field configurations with arbitrary large discontinuities in the x direction have low energy, and are hence unsuppressed in the path integral, as long as they are smooth in y, and vice versa. The presence of these discontinuous field configurations gives rise to a surprising interplay between UV and IR physics which was originally observed in [7,8], and later expanded on in [9][10][11]. Most notably, many quantities in the continuum theory remain sensitive to the lattice spacing of the underlying discrete model, that is to UV scale physics.
One place where this UV/IR mixing becomes important is when trying to understand whether the subsystem symmetries are spontaneously broken. Already in the original examples of [6][7][8] the issue of symmetry breaking turned out to be quite subtle. The theory has a large number of classical zero energy modes, carrying a large momentum in one direction but no momentum in the other. One could identify these modes with Goldstone modes of a putatively broken subsystem symmetry. However, quantum mechanically these modes get pushed to energies of order the UV cutoff. The symmetry is in fact restored. Closely related is the fact that the lightest states charged under the subsystem symmetry that could detect any breaking get pushed to UV scale energies as well.
A very good probe of symmetry breaking are correlation functions of charged operators. One thing we do in this work is to calculate these correlation functions directly in the continuum theories of [6][7][8] and verify that the subsystem symmetries do, in fact, remain unbroken as suggested by the spectrum. This same conclusion has also been recently obtained in [11], where the same correlation functions were obtained by carefully taking the limit of the discretized theory.
Given these results, one may wonder whether and when spontaneous breaking of subsystem symmetries in a continuum field theory can happen. As in ordinary quantum field theories, where the Mermin-Wagner-Coleman theorem forbids breaking of a continuous symmetry in 1+1 dimensions and below, spontaneous symmetry breaking should not always be allowed. A theorem due to Batista and Nussinov [12] shows that the criteria of whether a subsystem symmetry acting on a n dimensional submanifold can be broken or not are the same as for an ordinary symmetry for a field theory living in an n dimensional spacetime. For sufficiently high n there is nothing that forbids breaking and so one may wonder how to realize this phenomenon in exotic field theories.
One option to analyze the spontaneous breaking of subsystem symmetries is to work directly with the lattice models. Indeed, one of the early studies of fractons showed that the spontaneous breaking of a discrete subsystem symmetry accounts for the large ground state degeneracy in the X-cube model [13]. More recently progress has been made in studying the phase transition between the broken and unbroken phase of discrete subsystem symmetries in lattice models [10,[14][15][16][17][18]. The breaking of the continuous subsystem symmetry of [7,8] can also be achieved by taking a different limit of the lattice model, namely an infinte volume limit at finite lattice spacing [11]. While this limit does give rise to a broken phase, it does not seem to be amendable to continuum field theory language.
In contrast to these lattice constructions, in this work we present some novel continuum field theories that do exhibit subsystem symmetry breaking. Interestingly enough, all these theories are free theories, just like the original work of [6][7][8], and so our results are exact. Additionally, we find that whenever the theorem of [12] does allow symmetry breaking, our models indeed do lead to spontaneously broken subsystem symmetries.
The outline of this work is as follows. In section 2 we will review why we chose correlations of charged operators to detect symmetry breaking and show that our methods reproduce the known results for the theories of [6][7][8], consistent with the recent paper [11]. In section 3 we introduce the classical XY-plaquette model as one example of a field theory that does in fact lead to a spontaneously broken symmetry in dimensions where this is not prohibited by the theorem of [12]. In section 4 we generalize our model to a large class of closely related theories that all follow the same pattern. We end with a discussion and future directions in section 5.

Diagnosing the spontaneous breaking of subsystem symmetries
Let us consider a field theory living in d spacetime dimensions with a subsystem symmetry acting on n dimensional submanifolds. According to the modified Mermin-Wagner-Coleman theorem [12], the rules for whether spontaneous breaking of these subsystem symmetries can or can not occur are equivalent to the standard considerations for a quantum field theory living in n dimensions. In particular, spontaneous breaking of a continuous global subsystem symmetry can only happen if n ≥ 2. However this theorem in and of itself does little to predict if a specific subsystem symmetry is spontaneously broken, or even how to diagnose if such a breaking occurs.
The classical approach to diagnosing a spontaneously broken symmetry is to define a local order parameter that is charged under the symmetry, and then calculate if this order parameter acquires an expectation value in the vacuum state. This computation is delicate as the expectation value of a charged operator vanishes for any finite size system which obeys the symmetry. Hence the correct procedure is to calculate the expectation value of the order parameter in the presence of a symmetry breaking background field, which is then taken to zero. 1 The canonical example of this procedure is to consider a complex field in d dimensions, φ, and a Lagrangian with a Mexican hat potential This theory has a continuous U (1) symmetry which shifts φ by a phase, so φ is a natural order parameter for this theory. Notice that the minimum of the potential is not invariant under this symmetry, so classically we would expect that the vacuum expectation of φ to be | φ | = v. However this may not be the case in the quantum field theory. To compute the true expectation value of φ we need to introduce a symmetry breaking term to the action and then compute the expectation value of φ in the presence of this term. The true expectation value of φ will be given by the limit Though this approach was used to prove the Mermin-Wagner-Coleman theorem for subsystem symmetries in a lattice theory [12], it not clear how to define an order parameter for a subsystem symmetry in a continuum field theory, or how to introduce a symmetry breaking term that only breaks the desired subsystem symmetry.
A different approach to diagnosing spontaneous symmetry breaking is to compute the long range correlation functions of charged operators, and see if they vanish. This approach is used in the context of higher form symmetries to diagnose if a p-form symmetry is spontaneously broken [20,21]. For the above example of the complex scalar field φ, one can compute the two point function φ † (0)φ(r) .

(2.4)
If the field theory is local then the expectation values of a product of distant operators should factorize, and thus lim Hence if lim r→∞ φ † (0)φ(r) = 0 then the symmetry is not spontaneously broken, while if lim r→∞ φ † (0)φ(r) = 0 then φ acquires an expectation value and the symmetry is spontaneously broken. We note that the local charged operator φ may need to be renormalized to give a finite answer. In such a case the UV renormalization needs to be distinguished from the IR behavior of the 2-point function. This same method can be used to identify spontaneous symmetry breaking of subsystem symmetries, though we must ensure that the long distance correlation functions used to diagnose symmetry breaking is uncharged under the symmetry. Thus, for a subsystem symmetry, the correlation functions used to identify symmetry breaking will not be 2-point functions, but rather higher point correlation functions. We present an example of such a correlation function in the following subsection.

Symmetry breaking in the XY-plaquette model
As simple example of identifying spontaneous symmetry breaking of subsystem symmetries, we shall look at continuum field description of the XY plaquette model [1,7,8]. In [7,8] it was claimed that no spontaneous symmetry breaking happens, even in 3 + 1 dimensions where spontaneous symmetry breaking is permitted. We would like to verify this result using our diagnostic for spontaneous symmetry breaking.
The original XY-plaquette model in d + 1 dimensions has a continuum field theory formulation consisting of a compact scalar field φ, and the action (in Euclidean signature) This model has the global subsystem shift symmetry which acts on the primary 1 + (d − 1) dimensional planes and is given by where f i (x i ) is an arbitrary function of the i'th coordinate.
To diagnose if this symmetry is spontaneously broken we should compute a long-range correlator of charged operators, which for this theory is the operator e iφ (as φ is not a good operator in the continuum field theory [7]). Furthermore, we need the correlator we compute to be uncharged under the global subsystem shift symmetry. The simplest correlator that meets this criteria is the equal time correlator of four e ±iφ operators inserted on one of the spatial planes where without loss of generality we took the operators to lie in x 1 -x 2 plane, and where the remaining spatial coordinates are equal. As this is a Gaussian theory, we can take the functional integral directly, and we find that Starting with the case d = 2, we have that This integral is UV divergent, and so must be regulated. As the integral factorizes, it is natural to introduce independent cutoffs in the x and y directions so that this integral becomes 2 (2.12) This UV divergence cannot be regulated by a renormalization of the operators (or equivalently by local counter term) due to the presence of the terms proportional to log(x) log(Λ y ) and log(y) log(Λ x ). Not only do these terms make the result unregularizable, their presence indicates a UV-IR mixing in the charged states of this model, as observed in [7,8]. This same correlation function was computed in the lattice theory in [1] and [11], and both found the same result. 3 We can do a similar calculation in d = 3, where we can take the k 3 integral by introducing a cutoff, so (2.13) As before, this integral will include mixed terms of the form f (x, y) log(Λ z ) which cannot be regularized by a local counter term.
In both these instances this diagnostic of spontaneous symmetry breaking agrees with the results in [7,8,11], namely that the symmetry is not spontaneously broken. Furthermore, these results confirm that the charged states have energies of order the lattice spacing and so correlation functions involving charged operators cannot be regularized to give a consistent IR result.
We would like to construct theories with a spontaneously broken continuous subsystem symmetry, as the IR physics of such theories would be described by the Goldstone modes of the broken symmetry. Finding such theories would be a first step to studying the phase transition between the broken and unbroken phases, a new type of critical phenomena that has yet to be considered in the literature. In the next two section we will analyze a family of statistical models inspired by the XY-plaquette model which do indeed have spontaneously broken subsystem symmetries. 2 We note that the form of (2.12) should not depend on the particular regularization scheme. For example on can compute K2(x, y) using the standard spherical momentum cutoff of k 2 1 + k 2 2 ≤ Λ 2 and find that K2(x, y) =

The classical XY-plaquette model
To realize a theory with spontaneously broken subsystem symmetries, we can try to start from the lattice model related to the XY Plaquette model in d-dimensions [1]. Consider a square lattice in d dimensions with an angular variable φ i on each lattice site. In the original XY model we consider the classical statistical mechanics model where the energy of each configuration of spins is given by where the sum runs over all sites i, j which are nearest neighbors on the lattice. The partition function of this model is given by the path integral As the partition function is given by a Euclidean path integral, we can consider this model as a quantum mechanical model with a Euclidean action S E = E, while the d − 1 dimensional Hamiltonian will be related to this action by a Legendre transform.
We can follow the exact same procedure to construct an XY-plaquette model by taking the interactions to be on plaquettes rather than between nearest neighbor sites. This will result in the energy where the sum runs over all square plaquettes S containing the sites {s 1 , s 2 , s 3 , s 4 } in cyclic order, ∆ S φ = φ s 1 − φ s 2 + φ s 3 − φ s 4 , and J is the interaction strength. The partition function of the statistical theory is described by the Euclidean path integral This model is similar to the XY-plaquette model considered in [1] and [7], as the plaquette interaction term is the same, though in our case E XY-plaquette (φ) is the Euclidean action of the quantum model rather than the Hamiltonian, and it does not contain a kinetic term.
This XY-plaquette model has the subsystem symmetry of shifting all φ i 's on any d − 1 dimensional hyper-plane in the lattice by a constant angle. Similar to the original XY model, we would expect this model to be in a broken phase at low energy if the dimension allows for such a phase to exist. The continuum limit of this broken phase should correspond to taking the continuum limit of the of the Euclidean action on the lattice. This was done carefully in [7], and the resulting model is a compact scalar field φ in d dimensions with the Lagrangian This Lagrangian has the global dipole symmetry of shifting φ by any arbitrary function of a single variable As φ is 2π periodic, these functions are defined up to an integer multiple of 2π. Thus neither φ nor ∂ i φ are good operators in the continuum limit, while ∂ i ∂ j φ is a good operator so long as i = j, as is e iφ , similar to the analysis of [7]. Furthermore, we will only focus on field configurations that have a finite action in the continuum limit. This excludes any winding modes that have infinite action (due to a computation that is identical to the computation done in section 4.2 of [7]). 4 The spontaneous breaking of this dipole symmetry is what we are interested in analyzing. As this subsystem symmetry acts on d − 1 dimensional hyper-planes, the minimal dimension we can hope to have a spontaneously broken symmetry is d = 4, and indeed we will see that this is the case. However before we analyze the breaking of the dipole symmetry we present an analysis of the spectrum of the model. We note that a further study the correlation functions of uncharged operators is presented in appendix A.

Spectrum
To understand the spectrum of the model, we start from the equation of motion for φ 1≤i<j≤d Then we can expand φ in terms of Fourier modes, which leads to the dispersion relation If we consider the theory in Lorentzian signature, that is take one of the d spatial coordinates and Wick rotate it, then the dispersion relation becomes We note that this dispersion relation is continuous around k = 0, and indeed ω → 0 as |k| → 0 from any direction. Next we quantize the Lorenzian theory on a d − 1-torus of equal length , so that the momenta are quantized k i = 2π n i . In terms of the momentum modes the Lagrangian becomes Each of these modes for generic values on n behaves as a simple harmonic oscillator with ground state energy where ω given by the dispersion relation (3.9). The rest of the states are built via the generic Fock space construction on these ground states. If, however, all but one of the integers in n are zero then ω n = 0 and these fields have no harmonic potential. We may need to quantize these modes more carefully, as they can have large momenta in the remaining direction.
To quantize the Lagrangian for the modes where all but one of the n's are zero, we follow [8] and write these modes in position space by writing the field as Then the Lagrangian for these modes is (3.14) We see that unlike in [8] these modes are decoupled, and the "gauge" symmetry that connects these modes is just a subset of the global symmetry of the Lagrangian φ → φ + f (t). Overall these modes behave just like the zero modes of a standard periodic free scalar field, and do not acquire a mass of the order of the lattice spacing, which is very different than the theories considered in [7,8]. This analysis shows that the even for the modes where all but one of the n's are zero the Lagrangian in (3.11) is the correct Lagrangian, and the Hamiltonian for these modes in the i'th direction is The conjugate momenta of these modes, denoted by π i,n , are the generators of the global symmetry φ → φ + f i (x i ). These conjugate momenta do indeed commute with the Hamiltonian, as expected for a global symmetry, though their spectrum is continuous. Choosing a ground state for the system is the same as choosing the eigenvalues of these modes, and so the π i,n 's shift us from one ground state to another, in analogy to typical Goldstone modes. Thus this picture looks like a spontaneously broken symmetry, with the φ being an analogous Goldstone boson.
There is also the special mode n = 0, for which the Lagrangian vanishes and so is not a physical mode. This is very different than in [8] where this mode is the unique ground state of the system. The reason this mode disappears from the theory is the time translation symmetry φ → φ + f (t), which is generated by this zero mode.
Finally we note is that the resulting Hamiltonian is non-local in space for all of the modes. Indeed for the continuum theory the conjugate variable to φ is so the Hamiltonian of the system is which has a non-local kinetic term. This Hamiltonian still has a fine expansion in terms of Fourier modes, as seen above, and so can be quantized in the standard way. It is interesting that a local Lagrangian gave rise to a non-local Hamiltonian, which raises the question which is more fundamental. From an effective field theory prospective, one writes down the most general local Lagrangian that obeys the desired symmetries. Then one keeps only the renormalizable terms, which are found by computing the scaling dimensions. Based on this philosophy a local Lagrangian is more fundamental to understanding the IR physics, as it contains the relevant information about the symmetries and scaling dimensions. The model we presented fits nicely into this picture, as it is the simplest Lagrangian that obeys the desired subsystem symmetry, and indeed it seems that the resulting quantum field theory is local and well defined.
There is however the alternate bottom up approach where one starts from a local UV Hamiltonian, say on a lattice, and then tries to infer the physics at large distances. In this case one assumes that the Hamiltonian is local, and so it does not seem that the classical XY model can describe the IR physics that emerge from such a construction. This raises many interesting questions relating to the existence of phases of matter that fit into an effective field theory picture, but that cannot be described by a local lattice Hamiltonian.
There is also the possibility that certain phases of matter may arise from a local lattice Hamiltonian that has a non-local Lagrangian. For example one can introduce a nearest neighbor interaction of the conjugate momenta variables to a lattice Hamiltonian, which would keep the Hamiltonian local but result in a non-local Lagrangian. Such systems have been considered in the condensed matter literature, for example [22] studied similar Hamiltonians with subsystem symmetries, though it is not clear to us how such models fit into the canonical picture of Wilsonian quantum field theory.

Spontaneous subsystem symmetry breaking
To identify if the subsystem symmetry (3.6) is spontaneously broken, we will look at long range correlation functions of the charged local operator e iφ . 5 The 1, 2, and 3-point functions of these operators must vanish by the global dipole symmetry. Hence the simplest nonzero correlation function of such operators is the 4-point rectangle function, similar to what was considered in (2.8). This correlator is where we have chosen the rectangle to lie in the (x 1 , x 2 ) plane, and suppressed the dependence on the remaining d − 2 coordinates. Taking the Gaussian functional integral we can write this expectation value as where (3.20) Notice that this integral is well defined in the IR, as around k 1 , k 2 1 this integral behave as (k 2 1 k 2 2 )/( i =j k 2 i k 2 j ), which is finite around zero for a fixed angle. For d = 2 we can take this integral analytically, and we find This grows as we take x and y to be large, hence the full correlator (3.18) vanishes, and thus the subsystem symmetry is unbroken. For d = 3 this integral evaluates to where α ≡ |x| |y| . Again the integral grows as x or y (or both) are taken to be large, implying that the subsystem symmetry is still unbroken. This is expected as the minimal dimension allowed for spontaneous symmetry breaking in this model is d = 4.
For d = 4 we can write the rectangular function as where K(m) is the complete elliptic integral of the first kind. Notice that I 4 (x, y) is only a function of the ratio x/y, and so we expect it to be a constant for large rectangles with the ratio x/y fixed. However this integral is logarithmically divergent in the UV, and this UV divergence neads to be regularized. Moving to polar coordinates and introducing a UV cuttoff, we can take the radial integral to find (3.24) Notice that the angular integral multiplying the log Λ term is finite, hence we can renormalize the operator by a local counterterm to cancel the logarithmic UV divergence, leaving the remaining integral finite and cutoff independent. Canceling the UV cutoff does however make the integral logarithmically diverge as the area of the rectangles is taken to be large while fixing the ratio of the two sides. However this term approaches a constant if only one side of the rectangle is taken to be large while the other is kept fixed (notice that this is not the case for I 2 and I 3 which both diverge even in this limit, though slower than the divergence for fixed ratio). This implies that the subsystem symmetry is indeed spontaneously broken for d = 4.
We expect a similar behavior in dimensions d > 4, where an operator renormalization can take care of the UV divergence, while no IR divergent will exist at all. This is in alignment with the previous discussions of spontaneously breaking subsystem symmetries based on [12]. For completeness we present plots of I 3 (x, y) and I 4 (x, y) in figures 1 and 2.

Other classical XY models with subsystem symmetries
The classical XY-plaquette model can be generalized to a family of models in d dimensions that has a subsystem symmetry which acts on d − m + 1 dimensional hyperplanes. As before we start with a square lattice in d dimensions with an angular variable φ i on each lattice site, though now the energy of each configuration is given by The next class of models is the case m = 3, which is a classical version of the XYcube model considered in [15,23]. We present a brief analysis of this model below, before considering the general m case.

The classical XY-cube model
The classical XY-cube model has the subsystem symmetry of shifting all φ i 's on any d − 2 dimensional hyper-plane in the lattice by a constant angle. Similar to the previous model, we would expect this model to be in a broken phase at low energy if the dimension allows for such a phase to exist, which for this model is d > 4. The continuum limit of this broken phase should correspond to taking the continuum limit of the of the Euclidean action on the lattice. This was done carefully in [23], and the resulting model is a compact scalar field φ in d dimensions with the Lagrangian This Lagrangian has the global dipole symmetry of shifting φ by any arbitrary function of any two variables As φ is 2π periodic, these functions are defined up to an integer multiple of 2π. Thus neither φ nor ∂ i φ or ∂ i ∂ j φ are good operators in the continuum limit, while ∂ i ∂ j ∂ k φ is a good operator so long as i = j = k, as is e iφ , similar to the analysis of [23]. As before, we will only focus on field configurations that have a finite action in the continuum limit, excluding any winding modes that have infinite action.
To find the spectrum of the theory, we can expand φ in terms of Fourier modes, which leads to the dispersion relation If we consider the theory in Lorentzian signature, that is take one of the d spatial coordinates and Wick rotate it, then the dispersion relation becomes (4.5) As before, this dispersion relation is continuous around k = 0. We can quantize the Lorenzian theory on a d − 1-torus of equal length , so that the momenta are quantized k i = 2π n i . In terms of the momentum modes given in (3.10), the Lagrangian becomes Each of these modes for generic values on n behaves as a simple harmonic oscillator with ground state energy with ω given by the dispersion relation (4.5). The rest of the states are built via the generic Fock space construction on these ground states. If, however, all but two of the integers in n are zero then ω = 0 and these modes have no harmonic potential. As before, these modes behave just like the zero modes of a standard periodic free scalar field, and do not acquire a mass of the order of the lattice spacing, which is very different than the theories considered in [7,8,23]. This analysis shows that, even for the modes where all but two of the n's are zero, the Lagrangian in (4.6) is correct. The momenta conjugate to these modes, denoted by π i,j,n,m , are the generators of the global symmetries φ → φ + f ij (x i , x j ). These conjugate momenta do indeed commute with the Hamiltonian, as expected for a global symmetry. Choosing a ground state for the system is the same as choosing the eigenvalues of these operators, and so the π i,j,n,m 's shift us from one ground state to another, just as the zero modes for the usual Goldstone boson. The modes corresponding to where all but one of the n's are zero now have vanishing action, and correspond to the symmetries involving the time direction φ → φ + f 0j (x 0 , x j ). Finally we note that the resulting Hamiltonian is non-local in space, and the discussion regarding this non-locality from section 3.1 applies to this model as well.
As in the XY-Plaquette model, to check for which dimensions the symmetry is spontaneously broken we need to compute correlation functions of charged operators. In this model the simplest invariant correlation function is the cubic correlator of eight e iφ given by As this theory is Gaussian, we can take the path integral directly resulting in (4.10) For d = 3, the lowest dimension for which the model makes sense, we have that As before, this implies that the symmetry is not broken as the long range correlation functions grows. For d = 4 this integral is This integral is completely convergent for any finite values of x, y, z, though a closed form expression for the integrand is beyond our reach. We can however take various large distance limits of this function. Just via scaling arguments it is clear that if we take x, y and z to be large while keeping their ratios fixed then (4.14) If we take only x and y to be large while keeping z and the ratio x/y fixed then Finally, taking x → ∞ while keeping y and z fixed we see that the integral scales as so even in this limit the cube correlation function diverges. As this symmetry acts on d − 2 dimensional hyper-planes, it cannot be spontaneously broken in dimensions d ≤ 4, so this analysis does in fact make sense. For this symmetry to be spontaneously broken we need to be in at least d = 5, where the integral evaluates to (1 − cos(k 1 x))(1 − cos(k 2 y))(1 − cos(k 3 z)) (4.17) This integral is still UV finite, and converges for any values of x, y and z, though it has no known closed form expression. We can however do the same limit analysis to how it behaves when we take the cube to be large. If we take x → ∞ while keeping y and z fixed then the integral approximates as 6 , (4.18) which is a constant value independent of x. Hence this symmetry is indeed broken in d = 5, as expected.

Classical XY-hypercube models
We now turn our attention back to the general m case, the classical XY-hypercube model of degree m, with the energy given by (4.1). This model has a subsystem symmetry which acts on d − m + 1 dimensional hyperplanes. The continuum field theory of these models consists of a periodic scalar field φ with the Lagrangian This model has the subsystem shift symmetry given by which acts on d − m + 1 dimensional hyperplanes.
The spectrum of the model with arbitrary m has many of the same features as the spectrum of the plaquette and cubic models. The dispersion relation in Lorentzian signature will be given by The quntization of the theory on a torus follows the same procedure as was done for the plaquette and cubic models. Generic modes will behave like Harmonic oscillators with ground state energy E = 1 2 ω, where ω is given by the dispersion relation (4.21). Modes where all but m − 1 of the k's have zero harmonic potential but a non-zero kinetic term, and so are the zero modes of the theory. These modes generate the subset of shift symmetries in (4.20) that do not involve the time direction. Similarly, the modes where all but m − 2 of the k's are zero have a vanishing Lagrangian and generate the symmetries that do involve the time direction.
The dimensions for which the subsystem symmetry (4.20) can be spontaneously broken is d > 1 + m (by the generalized Merman-Wagner-Coleman theorem [12]). We have shown that the symmetry is indeed spontaneously broken when d = 2 + m for m ≤ 3 7 , and expect this trend to continue for arbitrary m.
To diagnose whether the subsystem symmetry (4.20) is spontaneously broken, one would need to compute a correlation function of e ±iφ operators placed on the corners of a mdimensional hypercube, as this is the simplest correlation function of charged operators that is uncharged under the subsystem symmetry. Calling this correlation function F m (x 1 , x 2 , . . . , x m ), one can preform the Gaussian path integral to find that where (4.23) The behavior of I m,d (x 1 , x 2 , . . . , x m ) at large separations determines the spontaneous breaking of the subsystem symmetry; if I m,d grows at large separation then the symmetry is unbroken while if it approaches a constant then the symmetry is spontaneously broken. We have shown that I 2,4 and I 3,5 approach a constant when one of the coordinates is taken to be large, and expect I m,m+2 to behave the same way for arbitrary m. We can actually compute this limit directly from (4.23). Notice that I m,m factorizes, and thus we can directly compute Hence the symmetry is not spontaneously broken for d = m. For d = m + 1 we can take the limit x 1 → ∞ while keeping the rest of the variables fixed, 8 and we find that (4.25) again implying the symmetry is not spontaneously broken. For d = m + 2 the integral I m,m+2 is UV finite assuming m > 2, so we can take the same limit x 1 → ∞ while keeping the rest of the variables fixed, resulting in the finite value 9 (4.26) As in this limit I m,m+2 approaches a constant (i.e. is independent of x 1 ) the symmetry is indeed spontaneously broken. Hence the subsystem symmetry (4.20) is spontaneously broken when d = m + 2 for arbitrary m, but not when d ≤ m + 1. This is in agreement with the expectations based on the generalized Merman-Wagner-Coleman theorem [12].

Discussion and future directions
Our main goal in this paper was to discuss spontaneous breaking of subsystem symmetries, and to construct well behaved continuum quantum field theories that model the spontaneously broken phase. Though we successfully built field theories that exhibit a spontaneous broken subsystem symmetry, there are still many open questions relating to these models, and to quantum field theories with subsystem symmetries in general. In this section we discuss in further detail two general areas of study that we find most interesting. The first is to study the 8 The simplest way to do this calculation is to first take both x1 and x2 to be large, when the limit reduces to an integral similar to (4.15). Then one can take the limit x1 → ∞ while keeping x2 fixed and obtain (4.25). 9 This limit is obtained in the same way as (4.18).
unbroken phase and the phase transition in the classical XY-plaquette model we constructed in section 3. The second is to analyze generalizations and interacting theories built upon the classical XY-plaquette model. There are however many other interesting aspects to keep investigating with regards to the models we constructed. One natural question to ask is whether the subsystem symmetries are still spontaneously broken at finite temperature, assuming such breaking is allowed by the number of dimensions. As our models have no free parameters, the symmetry is either broken or unbroken at all finite temperatures. To understand which is the case, we need to compute the same 4-point correlation functions of exponential operators we considered when taking one of the spatial dimensions to be compact. Though one can evaluate the sum over the Matsubara frequencies in the compact direction, the resulting integrals over the momenta in the remaining directions are unwieldy, and must be UV regularized for the classical XY-plaquette model. Hence we decided to leave this problem for a future endeavour.
Another interesting observation is that the Lagrangian (3.5) in d = 2 is similar to the Lagrangian of a zero velocity chiral boson, L = α∂ 0 φ∂ 1 φ. The chiral boson also possesses the same subsystem symmetry (3.6), which is not spontaneously broken in d = 2. One can try to generalize the chiral Lagrangian to higher dimensions, namely consider the Lagrangian L = α i =j ∂ i φ∂ j φ which also possesses the same subsystem symmetry. However this Lagrangian breaks the reflection symmetry of the lattice, and also has a Hamiltonian that is unbounded from below. One can overcome some of the difficulties in quantizing the d = 2 theory using the conformal symmetry, and perhaps in higher dimensions there is a different conformal symmetry which emerges, similar to the ones studied in [24], which may be worth studying.

The unbroken phase and phase transition
Up to now we have focused on the symmetry broken phases of the classical XY models we considered. However, just like the standard XY model, the lattice models have a phase where the symmetry is unbroken, and a phase transition between the two phases. To model the phase transition and the unbroken phase, we would like to construct a Lagrangian for a complex scalar field Φ that has the desired U (1) subsystem symmetry. Then tuning the potential of Φ in this Lagrangian would (at least classically) give Φ a vacuum expectation value, and thus take us between the spontaneously broken and unbroken phases.
This Lagrangian can be built in a similar manner to the Lagrangian's with global dipole symmetries constructed in [25,26]. Focusing on the classical XY-plaquette model, the lowest order terms in this effective Lagrangian are In the symmetry broken phase we would choose the potential in (5.1) to be a Mexican hat. We can then expand the Lagrangian around one of the minima of the potential, and to leading order in the fields we would arrive at the continuum Lagrangian for the XY-plaquette model given in equation (3.5). This is a simple check that the Lagrangian we wrote down has the desired symmetry broken phase that we analysed in section 3.
To access the unbroken phase one would need to choose a generic potential with a minimum at |Φ| 2 = 0. Around this minimum the Lagrangian to lowest orders reads Unfortunately this effective Lagrangian is non-Gaussian even at lowest order due to the kinetic term containing 4 Φ fields. Hence an analytic or even perturbative description of the dynamics in the unbroken phase is beyond our reach (as was also the case in [25].) This phase can still be studied numerically, either on the lattice using (3.3) or in the continuum field theory (5.2). This phase may also be accessible analytically using non-perturbative techniques, like perhaps taking some sort of large N limit. It would also be of great interest to study the phase transition between broken and unbroken phases. However, as the effective Lagrangian (5.2) is non-Gaussian and strongly coupled in the vicinity of the phase transition, even a mean field approximation of the transition is not possible. It is not even clear if the phase transition is continuous or first order, and the literature on phase transitions in similar models is uncertain.
For example the plaquette Ising model [27], which is the Z 2 version of the classical XY-plaquette model (3.3), has been shown to undergo a weakly first order phase transition [28,29]. On the other hand a recent study of the phase transition in models related to the Z N versions of the X-cube model indicated that the transition is continuous for N > 4 [17].
As the plaquette Ising model is dual to the Z 2 X-cube model [13], these two results are in agreement. One may then think of the XY-plaquette model as a N → ∞ limit of Z N models, in which case the phase transition is likely continuous, though more numerical or analytical evidence is needed to confirm this conjecture.

Generalizations and interacting theories
Aside from modelling the unbroken phase and the phase transition, it would also be interesting to construct interacting continuum field theories based on the classical XY-plaquette model. We present two general systematic constructions of such theories: the first method is to gauge the U (1) subsystem symmetry to a tensor gauge theory in the spirit of [25], while the second idea is to promote a global U (1) symmetry to a subdimentional symmetry.
The most straightforward generalization of the classical XY-plaquette model is to gauge the symmetry (3.6). This can be done by starting with the generalized Lagrangian 5.1, and then replacing the derivative term Φ∂ i ∂ j Φ−∂ i Φ∂ j Φ with a "covariant" derivative term similar to [25], Here A ij is a symmetric U (1) tensor gauge field with the transformation properties A ij → A ij + ∂ i ∂ j α, while Φ → e iα Φ under a gauge transformation.
We would also need to add a kinetic term for this gauge field, and in general we one could investigate this tensor gauge theory in its own right. The gauge invariant field strength of this theory would be [8] F  [7,8,25,30].
A different approach to systematically construct interacting theories from the classical XY-plaquette model is by promoting a global U (1) symmetry to a subsytem symmetry. This procedure is reminiscent of the gauge principle, where one promote a global symmetry to a local one.
For example we can start out with a field theory involving a single complex scalar field ψ which is charged under a global U (1) symmetry, and a Lagrangian Then we would like to partially gauge the global symmetry into a subsystem symmetry. This can be done by introducing a "covariant" derivative where φ acts as the analog of a gauge field. Under the subsystem symmetry φ transforms as (3.6), while ψ → e ig µ fµ(x µ ) ψ. The full action of this theory would then become The resulting theory is an interacting field theory with a subsystem symmetry. In essence one can partially gauge any theory with a conserved U (1) charge, resulting in a whole new family of theories to play with. It would be interesting to investigate how these partially gauged theories behave, how to generalize them to non-Abelian symmetries, and perhaps to understand if they have analogous Coulomb and Higgs phases. Finally, we note that this partial gauging procedure can also be used to construct interacting theories coupled to the XY-plaquette models of [7,8], which would also be interesting to investigate.

A Correlation functions in the Classical XY-plaquette model
The structure of correlation functions of simple uncharged operators is very non-standard in models with subsystem symmetry. In particular these correlation functions tend to diverge even when the operators are far apart, but one of the spatial variables is close [11]. Such UV divergences in IR correlation functions has been considered a hallmark of UV/IR mixing. As such, we are interested in computing the correlation functions of the simple operators ∂ i ∂ j φ, and to observe any patterns of UV/IR mixing in the classical XY-plaquette model. We will focus our analysis on the theory in two and three dimensions for simplicity, and we think this also provides a general qualitative idea of the structure of the correlation functions even in higher dimensions.

dimensions:
For d = 2 we can compute the correlation function directly as This correlation function is somewhat simplistic, but does contain some interesting features. Notice that the correlation function is indeterminate when x 1 or x 2 are zero, which resembles the UV/IR mixing of correlation functions in [11]. However, unlike in the XY plaquette model of [7,11], this divergent behavior cannot be regularized by an IR regulator (say by putting the system on a torus,) but rather requires a UV regulator. A similar picture appears in 3 dimensions, though there the correlation functions have more structure.

dimensions:
In three dimension we can compute the 2-point correlation function of ∂ 1 ∂ 2 φ to be = g 8π 2 dθdrr 2 | sin θ cos θ|e ir(cos θx 1 +sin θx 2 )−r| sin θ cos θ||x 3 | = g 4π 2 2π 0 dθ x 3 sin 2 θ cos 2 θ x 2 3 sin 2 θ cos 2 θ − (x 1 cos θ + x 2 sin θ) 2 x 2 3 sin 2 θ cos 2 θ + (x 1 cos θ + x 2 sin θ) 2 3 (A.2) This result is finite if x i = 0 for i = 1, 2, 3, but diverges when the two operators share one of the same coordinates. We can expand this integral around these limits to see how it diverges. If we take x 3 → 0 then we see that On the other hand if we take x 1 → 0 then and taking x 2 → 0 results in a similar expression with x 2 replacing x 1 . We can also work in cylindrical coordinates, defining x 1 = ρ cos(δ), x 2 = ρ sin(δ), in which case the 2-point function becomes where the integral is only a function of ρ/x 3 . In the limit x 3 ρ we get while in the limit ρ x 3 the integral evaluates to In all these limits the 2-point function diverges even though the operators are separated. This is similar to the divergence of correlation functions as two operators are brought close together, only in this case the operators stay separated and only the separation in one of their coordinates is taken to be close. This behavior is similar to the one obsercved in [11].
As this divergence cannot be regulated by IR regulator, it is unclear to us how these correlation functions fit into the picture of UV/IR mixing. Furthermore, one can define a re-normalized operator O 12 ∼ Λ 3/2 ∂ x ∂ y φ, where Λ is some UV scale. Then this renormalized operator would have a finite 2-point function of the form O 12 (x 1 , x 2 , x 3 )O 12 (y 1 , y 2 , y 3 ) = C, x 3 − y 3 = x 1 − y 1 = 0, or x 3 − y 3 = x 2 − y 2 = 0, 0, otherwise, (A.10) where C is some constant that depends on the re-normalization scheme. This renormalized operator has restricted mobility as it can only correlate along lines of fixed x 1 and x 2 in the x 1 -x 2 plane, and in terms of these renormalized operators the theory would have no knowledge of the UV physics.
However we are unconvinced that this is the correct way to describe the IR theory, and are skeptical if this divergence of the 2-point function should be regulated at all. Rather, we think this divergence is contains important information about the IR physics, and the causal structure of the theory. An analogous situation arises in relativistic quantum field theory, where it is common for propagators diverge on the light cone. This light cone singularity is a genuine property of the IR physics and the structure of space-time, and not something which is needs to be regulated away. In any case it would be interesting to further understand the consequences of these divergences, and to fit them more squarely into the picture of UV/IR mixing.