Boundary Modes in the Chamon Model

We study the fracton phase described by the Chamon model in a manifold with a boundary. The new processes and excitations emerging at the boundary can be understood by means of a diagrammatic framework. From a continuum perspective, the boundary theory is described by a set of scalar fields in similarity with the standard $K$-matrix Chern-Simons theory. The continuum theory recovers the gapped boundaries of the lattice model once we include sufficiently strong interactions that break charge conservation. The analysis of the perturbative relevance of the leading interactions reveals a regime in which the Chamon model can have a stable gapless fractonic phase at its boundary.


Introduction
Fracton phases  constitute an intriguing and novel type of quantum matter, whose understanding sits at a confluence of different research fields.Fractonic phases are often classified in terms of the restricted mobility presented by their emergent quasiparticles, which can either move under mild restrictions (type-I) or are completely immobile (type-II) due to the proliferation of new excitations at every step of their motion.Three-dimensional gapped fracton phases also display a robust ground state degeneracy.However, unlike conventional two-dimensional (2D) topological phases, their ground state degeneracy depends on some geometric data, typically growing exponentially with the linear size of the system.
Boundary theories of fracton phases may present peculiar properties due to the geometric sensitivity natural to these systems.In the context of lattice theories, gapped boundaries of the 2 X-cube model [3,8] have been analyzed in Refs.[43,44].It has been pointed out that braiding of fractons in the bulk is geometry dependent and insufficient to classify the different types of boundaries [43].More recently, the authors of Ref. [45] approached the X-cube model from a continuum perspective and showed that its boundary theory can be viewed as a generalized Kmatrix Chern-Simons theory resembling the situation in standard topological orders [46][47][48].
In this work, we study the boundary theory for the 2 Chamon model [1] from both perspectives, lattice and continuum.First we show that the excitations in the lattice model can be represented by a combination of a few fundamental diagrammatic structures, closely related to cage nets in other fracton models [49,50].At the boundary, additional diagrams are allowed and account for new processes that violate fracton parity constraints and affect the mobility of the excitations.
On the continuum side, we extend on previous work [20] and show that the Chamon model can be described by a Chern-Simons-like theory with higher derivatives and two gauge fields labeled by a layer index.In the presence of a boundary, the effective field theory gives rise to physical boundary degrees of freedom described by a K-matrix theory.The boundary action is equivalent to the one recently derived for the X-cube model [45], despite the different bulk excitation spectrum of these models.However, unlike the X-cube model, the description of a (001) boundary in the Chamon model requires imposing a boundary condition for the normal derivative of the fields.As one of our main results, we show that our prescription leads to a consistent boundary theory with the correct number of degrees of freedom and encodes the information that line operators that terminate at the boundary create a single fracton at the endpoint.
To recover the gapped boundary spectrum of the lattice model, we investigate the effects of < l a t e x i t s h a 1 _ b a s e 6 4 = " + U t 9 s q v j A W v 9 B 9 y (1).(b) fcc lattice of octahedra.We distinguish between two types of (001) planes, in which the octahedra belong to either A and B or to C and D sublattices.
perturbations allowed by the discrete symmetry and the compactification of the bosonic fields.
While the boundary phases of the lattice model can be identified with the strong coupling regime, we also discuss the perturbative relevance of the interactions in the weak coupling regime.Another main result is that we predict the existence of a stable gapless boundary phase where all perturbations that break charge conservation become irrelevant.Such a phase has emergent continuous subsystem symmetries and is unexpected from the point of view of the microscopic model.The paper is organized as follows.In Sec. 2 we present the 2 Chamon model both in the bulk and in the presence of a boundary.In Sec. 3 we review the formalism of Ref. [20] and examine the symmetries of the effective field theory.Section 4 contains our results on the continuum description of the boundary modes.In Sec. 5 we offer some final remarks and point out possible routes for future work.Finally, the diagonalization of the boundary Hamiltonian and the calculation of correlation functions of charged operators are detailed in the appendices.

Lattice Model
The Chamon code [1,2] is an exactly solvable spin model defined on the fcc lattice which displays type-I fractonic behavior.Consider a cubic lattice Ω with sublattices Ω even and Ω odd .We say that a given site r = a s (i, j, k), with i, j, k ∈ and a s the lattice spacing, belongs to Ω even (Ω odd ) if i + j + k is even (odd).We assume that both Ω even and Ω odd contain an even number of sites.We place a spin 1/2 with Pauli operators σ I r , where I ∈ {x, y, z} ≡ {1, 2, 3}, at every site r ∈ Ω even .For every r ∈ Ω odd we define the stabilizer operator acting on the vertices of the octahedron centered at r [see Fig. 1(a)]: where ê1 = a s (1, 0, 0), ê2 = a s (0, 1, 0), and ê3 = a s (0, 0, 1).The Hamiltonian of the Chamon model is given by Since all stabilizers commute among themselves and square to the identity, O 2 r = 1, the model is exactly solvable with eigenstates labeled by the eigenvalues ±1 of the stabilizers.A ground state |Ψ〉 of the Hamiltonian obeys the condition ( The ground state degeneracy of the Chamon model on a 3-torus was calculated in Ref. [2].The ground states can be labeled by the eigenvalues of rigid line operators that wind around the lattice.If all lattice dimensions are even, with lengths L x = 2p x , L y = 2p y , and L z = 2p z , the number of ground states is 2 4 gcd(p x ,p y ,p z ) where gcd(p, q, r) stands for the greatest common divisor of the integers p, q, r.
We refer to a single defect stabilizer with O r = −1 as a fracton.In a lattice with periodic boundary conditions in all three directions, the model has the constraints where A, B, C, D denote the four sublattices of Ω odd represented in Fig. 1(a) as red, green, blue, and yellow dots, respectively.More specifically, we label the sublattices as follows: A : r = a s (2l + 1, 2m, 2n), B : r = a s (2l, 2m + 1, 2n), C : r = a s (2l, 2m, 2n + 1), D : r = a s (2l + 1, 2m + 1, 2n + 1), (5) where l, m, n ∈ .Equation (4) implies that defects can only appear in pairs on each of the sublattices.In fact, the action of a Pauli operator on a ground state creates four defects that belong to two sublattices.For instance, the operator σ z r anticommutes with four stabilizers contained in a (001) plane, namely O r±ê 1 and O r±ê 2 .Applying σ z r on the ground state, we create four defects belonging to A and B sublattices if r • ê3 is even or to C and D sublattices if r • ê3 is odd; see Fig. 1(b).Importantly, the independent 2 constraints in Eq. ( 4) allow us to distinguish between two sets of (001) planes (and equivalently for the other spatial directions) and to label the defects by the pair of sublattices on which they are created.
The model admits a diagrammatic description similar to that of cage-net fracton models [49].We start by representing the action of a local spin operator by strings connecting the four flipped stabilizers on a given plane, see Fig. 2(a).Due to the 2 character of the model, the superposition of two strings on a link between two stabilizers is equivalent to no string.Thus, applying spin operators along a line we can create an enlarged ribbon structure that accounts for the motion of pairs of defects as shown in Fig. 2(b).In this representation, the defects can be identified as the corners at which the strings form 90 • angles.Following the nomenclature of Ref. [2], we will refer to the pairs of defects at the ends of the strings as dipoles.Note that in this 2 model the dipole is made up of two defects with the same 2 charge, but these defects cannot annihilate each other due to the constraints in Eq. ( 4).Importantly, these dipoles move along rigid lines aligned with face diagonals of the lattice.By contrast, a fracton (or monopole) cannot move without creating additional excitations, whereas there are no mobility constraints for four defects (a quadrupole) which can be created by the action of local operators on the ground state.This nomeclature < l a t e x i t s h a 1 _ b a s e 6 4 = " + U t 9 s q v j A W v 9 B 9 y is further motivated by the higher-moment conservation laws in U(1)-symmetric effective field theories for fracton models [51,52].We will discuss the gauge theory for the Chamon model in Sec. 3.
We can now ask about all the basic diagrams representing operators which commute with the Hamiltonian.By construction, the action of such operators on a ground state must preserve the condition in Eq. (3).The allowed configurations are shown in Fig. 3(a).The first type of diagram consists of a vertex with four strings emanating from the same site and forming 90 • corners in two different planes.The second type is simply a straight line passing through a stabilizer site.From these basic elements we can build cage nets without any defects at their endpoints, as shown in Fig. 3(b).This type of cage-net tetrahedron in the lattice model is consistent with the form of gauge-invariant operators in the Chern-Simons theory as discussed in Ref. [9].
We now turn to the model with a boundary.We consider a semi-infinite volume with an open boundary placed at z = 0, with the system occupying the region z ≤ 0. At the boundary we define "broken" stabilizers given by five-site operators for r = a s (i, j, 0) ∈ Ω odd .We refer to a single defect of a five-site stabilizer as a boundary fracton.Physically, we can understand this choice of the boundary stabilizers by considering that we start from the Chamon model in the thermodynamic limit and apply a strong magnetic field in the z direction, H Z = −h r,z>0 σ z r , to freeze out the spins in the half space z > 0. The effective Hamiltonian is then obtained by projecting the stabilizers in the z ≤ 0 region onto the sector where all spins in the z > 0 region are polarized with σ z r = +1.On the boundary plane, this projection reduces the standard octahedral stabilizers in Eq. ( 1) to the five-site stabilizers in Eq. ( 6).
The boundary stabilizers violate two out of the four constraints in Eq. ( 4).For the boundary at z = 0, the parities associated with C and D sublattices are no longer conserved.Similarly, shifting the boundary to z = −1 would break the parities associated with A and B sublattices.Note that, unlike the X-cube model [43,44], removing a layer in the Chamon model does not alternate between "smooth" and "rough" boundaries in analogy with the terminology for the toric code [38].In contrast, the boundaries at z even and z odd are geometrically similar, but they can be classified by the pair of constraints that they break.
From the action of spin operators on stabilizers close to the boundary, we obtain new string configurations in the diagrammatic description.The new basic diagrams take the form of triangles in xz and yz planes with the longest edge at the boundary, as depicted in Fig. 4. Importantly, the creation of an odd number of defects out of the vacuum is now possible due to the broken parities for two sublattices.
The triangular diagrams are associated with distinctive processes involving the excitations that reach the boundary.The simplest example is the transmutation of a dipole into a single defect.We first bring the dipole to the boundary by successive application of spin operators along a line, as depicted in Fig. 5(a).We then glue the appropriate triangle, which has the effect of converting the Another important process is the reflection of a dipole at the boundary.After converting the dipole into a boundary fracton, we can again compose the diagram with another triangle that transforms the fracton back into a dipole, but now propagating in a perpendicular direction in the same or in a different plane.The diagram for reflection in the same plane is shown in Fig. 5(c).Alternatively, we can view the diagram in Fig. 5(c) as a process in which two dipoles meet at the boundary and annihilate.Importantly, there is no process which allows for a single dipole to condense at a (001) boundary of the Chamon model.In contrast, pairs of e or m particles can condense in the rough or smooth boundaries of the X-cube model [43,44].
The possibility of converting a bulk dipole into a boundary fracton also modifies the mobility of dipoles along the surface.When an (a, b) dipole in the z = 0 plane meets a bulk dipole that impinges on the boundary, the bulk dipole can be reflected back while the dipole that moves on the boundary plane changes its direction of propagation.Thus, a boundary dipole can lift its mobility restrictions by undergoing an elastic collision with a bulk dipole.
Finally, we can associate boundary half cages (BHCs) [43] with the (001) boundaries of the Chamon model.By definition, a BHC operator does not create any excitations in the bulk but creates isolated excitations on the boundary.We construct such an operator by bringing a cagenet tetrahedron to the boundary and terminating four ribbons with boundary fractons.Once again, the pair of constraints respected by the boundary determines the sublattice indices for the isolated monopoles that can appear in the BHCs.A particular example is illustrated in Fig. 6.

Effective Field Theory in the Bulk
The physics of the Chamon model can be understood from a effective theory perspective by means of a Chern-Simons-like theory.A first attempt of describing the model in this way through a topdown approach can be found in [9] and later on a description of the model via a bottom-up approach and the generalization of it in higher dimensions was obtained in [20].In this section SciPost Physics Submission < l a t e x i t s h a 1 _ b a s e 6 4 = " R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f 6 q u N B g = = < / l a t e x i t > y  we review the framework of Ref. [20] and argue about the necessity of a layer index, that is explicit in the microscopic description, that has an important effect in the continuum theory of the bulk, leading to a BF-like theory instead of the previously Chern-Simons-like descriptions.

Effective Action
We start by representing the Pauli operators via the exponential map where we assume an implicit sum over repeated indices m, n ∈ {1, 2}.Here K is the 2 × 2 antisymmetric matrix The Pauli operators are expressed in terms of two lattice fields θ n (r).The identity σ x r σ y r σ z r = i1 implies a neutrality condition for the t-vectors, 3 I=1 t I m = 0. We can single out the z direction and fix the t-vectors simply as The algebra of the Pauli operators translates into the new variables through the relations provided that we impose the commutation relations The lattice fields are compact since the shift θ n (r) → θ n (r) + 2πm n with m n ∈ leaves the Pauli operators invariant.
The representation with the t-vectors in Eq. ( 9) is spatially anisotropic.Given the special role of the z direction, we label the fields at each point by whether r ∈ Ω odd belongs to the set L 1 of (001) planes with A and B sublattices (r • ê3 even) or to the set L 2 of planes with C and D sublattices (r • ê3 odd).We introduce a layer index l and write θ n (r) → θ (l)  n (r), with l = 1 for AB planes and l = 2 for CD planes.We also define the forward and backward lattice derivatives along with the second derivative The lattice Hamiltonian in Eq. ( 2) can then be written as where we have symmetrized the complex exponentials to render the Hamiltonian explicitly Hermitian.
We can write down an effective lattice action describing the ground state subspace of the Hamiltonian in Eq. ( 14).We achieve this by pinning the argument of the cosine to its minimum.This procedure leads to the action where θ (l) 0 can be regarded as a Lagrange multiplier implementing the ground state projection.Note that the lattice action has a well defined scaling if the microscopic fields θ (l)  m are dimensionless in the standard sense, but the Lagrange multiplier θ (l) 0 must have dimension of inverse time.The action is invariant under the gauge transformation with arbitrary dimensionless functions ζ (l) (r, t).To take the continuum limit, we define the rescaled fields where a s is the lattice spacing.Note that the dimension of A (l) 0 has a contribution from the dimension of θ (l) 0 .If we consider a short-time cutoff a t ∼ a [t]  s with dynamic exponent [t] = 2, then A (l) 0 has the same spacetime dimension as A (l) n .We also replace 1 As a result, the continuum action reads with the derivative operators Equation ( 18) is equivalent to the Chern-Simons-like theory proposed in Refs.[9,20,21] with an additional copy.The two copies with l = 1, 2 are related by the lattice translation r → r + a s (1, 1, 1), which exchanges the sublattices AB ↔ CD.This property is reminiscent of Wen's plaquette model on the square lattice [53], where two types of excitations are related by a lattice translation.We note that the lattice spacing does not completely disappear from the continuum theory.For instance, it shows up in the gauge transformation This is a manifestation of the UV-IR mixing inherent to fracton theories.This microscopic parameter will eventually need to be restored to make sense of the gauge structure as well as some physical properties computed in the field theory.
To keep track of the spatial derivatives, we find it convenient to switch to the notation

SciPost Physics
Submission We then define the magnetic fields We can use the antisymmetric differential operator J for general values of the indices I, J ∈ {1, 2, 3}.Note that, in particular, D x y = D xz − D yz .Along the same lines, we consider the antisymmetric field A I J (r) = −A J I (r) with the x y component defined as A x y (r) ≡ A xz (r) − A yz (r).With this convention, we define the electric fields as so that E (l)  x y (r) = E (l) xz (r) − E (l) yz (r).Note that the electric and magnetic fields are gauge-invariant operators.We can then work with the action in Eq. ( 18) in the following equivalent form:

Symmetries and Line Operators
From the equation of motion of A (l) I J , we can read from Eq. ( 25) Integrating Eq. ( 26) over the J K plane, we obtain the conservation law where we introduced the coordinates ξ I = |ε I J K | 2 (x J + x K ) and ξI = ε I J K 2 (x J − x K ), corresponding to the directions of motion of dipoles in the lattice model.In this notation, we have D J K = 2∂ ξ I ∂ ξI .We refer to Eq. ( 27) as the electric symmetry, since it emerges from the vanishing electric field.
We can derive additional subsystem conserved quantities if we restrict the integration to be along lines instead of planes.Assuming that the theory is defined on a finite volume with periodic boundary conditions in all spatial directions, we obtain where the integrations are along lines that wind around the system in the direction of ξ I and ξI .There are also conserved quantities associated with ribbons of arbitrary width w: The generators of the above symmetry, exp can be interpreted as dipoles with charges ±q separated by a distance w.
The model also exhibits a magnetic symmetry with a conserved charge given by the integral of the magnetic field over the entire system: We can also define conserved quantities on manifolds of codimension one, of the form with u ∈ {ξ x , ξx } and v ∈ {ξ y , ξy }.

Boundary Theory
In this section we derive the the continuum description of the Chamon model with a boundary.We point out some analogies with K-matrix theory and discuss the constraints of subsystem symmetries on correlations of charged operators.We then analyze the effects of perturbations and how they relate to possible gapped and gapless boundary phases.

Boundary Action and Symmetries
The action in Eq. ( 25) is gauge invariant up to boundary terms.Let us consider the theory defined on a manifold M = × V, with representing the time direction and V the spatial volume with a boundary at z = 0, i.e., V = x × y × (−∞, 0].The variation of the action under a general transformation of the fields reads For the gauge transformation in Eq. ( 20), we obtain As in standard Chern-Simons theories, gauge invariance may be restored by restricting the gauge transformations in the presence of a boundary.In this process, the gauge fields become physical degrees of freedom at the boundary.One possibility is to impose 0 vanish as well.A more general gauge-fixing condition that yields δS = 0 is with the convention of summing over the repeated index l .Here we allow for a linear combination of the fields with real coefficients κ l l obeying κ l l = κ ll .We can satisfy the zero-flux condition B (l) = 0 in the ground state sector by writing the gauge fields in the form where ϕ l are dimensionless scalar fields to be associated with the boundary degrees of freedom.Substituting Eq. ( 36) into the bulk action Eq. ( 25), we find that the integrand is a total derivative that integrates to the boundary action independently of the choice of coefficients κ ll in Eq. (35).Note that the boundary action involves the derivative of the scalar fields with respect to the direction perpendicular to the boundary.This derivative is related to the difference between fields in adjacent layers.To match the correct number of degrees of freedom between bulk and boundary, we demand that the derivative be a linear combination of the fields: with real coefficients c l l .These coefficients may depend on microscopic details of the boundary interactions in the generic theory, beyond the exactly solvable lattice model.Using integration by parts, we realize that only the antisymmetric part of c ll contributes to the action, and we obtain Assuming c 21 > c 12 , we can rescale the fields ϕ l → (c 21 − c 12 ) −1/2 ϕ l to cast the boundary action in the form with the same K matrix that appears in the bulk theory in Eq. ( 8).This action implies that upon quantization the pair of boundary fields must obey the equal-time commutation relation We stress that the number of independent fields in the boundary theory is associated with two sets of planes in the bulk theory.Had we not labeled the two gauge fields from the start, we would have to drop the index l in Eq. ( 37), but we would be forced to treat ∂ z ϕ as independent from ϕ.
In any case, we would reach the same conclusion about the doubling of the modes and would end up with a boundary action formally equivalent to Eq. ( 40).We will return to the interpretation of Eq. ( 38) in Sec.4.3.
The action in Eq. ( 40) does not give rise to any dynamics since the corresponding Hamiltonian vanishes identically.This result is analogous to the derivation of the action for chiral edge modes from the Chern-Simons theory for quantum Hall states [54,55].To determine the boundary dynamics, we note that the action is invariant under the shift symmetries where f l (ξ z ) and fl ( ξz ) are arbitrary functions of ξ z = 1 2 (x + y) and ξz = 1 2 (x − y).The usual global U(1) symmetry corresponds to choosing f l and fl to be constants, but the general form of Eq. ( 42) is connected with the subsystem symmetries of the fractonic theory.We then add to the action a term that only involves spatial derivatives and respects the above symmetries: where M is a symmetric matrix, as imposed by the even number of spatial derivatives in the last term.The shift symmetries are generated by the currents which obey the continuity equations The latter are simply the equations of motion derived from the action in Eq. ( 43).Solving these equations in terms of Fourier modes with momentum p = (p x , p y ), we find the dispersion relation where we define Thus, the stability of the boundary theory requires µ > 0. In analogy with the velocity matrix of chiral edge states in the quantum Hall effect [54], hereafter we assume that M is a positivedefinite matrix.Note that µ is the single parameter governing the dispersion of the elementary excitations.It is also interesting to note that, using Φ = (ϕ 1 , ϕ 2 ), we can define a duality transformation Φ → σ 1 Φ, where σ 1 is the Pauli matrix acting on the two-component scalar field.This transformation maps K → −K and M → σ 1 M σ 1 without affecting the dispersion parameter µ.
The sign change of the K matrix preserves the algebra of the physical operators, and can be used to implement the proper re-scaling of the fields above Eq.( 40) in case c 21 < c 12 .
In addition to the conserved currents in Eq. ( 44), the boundary action has a dipole winding symmetry generated by the currents The corresponding charges of the winding symmetry are dipole configurations that wind around the system in either ξ z or ξz directions: The latter are reminiscent of the bulk theory, cf.Eq. ( 30).
The boundary Hamiltonian derived from Eq. ( 43) reads To diagonalize the Hamiltonian, we assume periodic boundary conditions in the x y plane and employ the mode expansion where L and L are the lengths in the directions of ξ z and ξz , respectively, and p n = 2πn/L and pn = 2πn/ L are the discrete momenta.Defining ρ l,n,n = 2p n pn ϕ l,n,n , we find that the normalmode operators obey a generalized U(1) Kac-Moody algebra In terms of ρ l,n,n , the Hamiltonian becomes Using a Bogoliubov transformation (see Appendix A for details), we obtain where b † n,n and b n,n are bosonic creation and annihilation operators, respectively, associated with momentum p = 1 2 (p n + pn ), 1 2 (p n − pn ) and energy ω n,n = 2µ|p n pn |.In the thermodynamic limit L, L → ∞, the Hamiltonian in Eq. ( 50) predicts gapless boundary modes with an anisotropic quadratic dispersion.In fact, the excitation energy vanishes along the lines p y = ±p x in momentum space for arbitrary values of |p|.As a consequence, shortwavelength modes can contribute to the low-energy physics, another hallmark of the UV-IR mixing.However, we still need to analyze the effect of symmetry-allowed interactions that may gap out the boundary spectrum, similarly to what happens to edge modes of topological phases without charge-conservation symmetries [39][40][41].To understand this point, we must first identify the operators that create charged excitations at the boundary.

Charged Operators
From Eq. ( 44), we can define the charge operators The operators charged under Q l are the exponentials where we label the fields by the charge vector q = (q 1 , q 2 ), so that Ψ −q (r) = Ψ † q (r).Using Eq. ( 41), we obtain Expressing the position r ∈ ∂ V in terms of the coordinates ξ z and ξz , we can write the commutation relation for the ϕ l fields as where sgn(x) is the sign function.As a consequence, the charged operators obey the algebra In particular, the operators with "elementary" charges q = (1, 0) and q = (0, 1) do not commute with each other.We shall think of Ψ † (1,0) and Ψ † (0,1) as creating monopoles at the boundary of the system.We note that from the field theory alone it is not entirely clear how the charges and the entries of the K matrix should be quantized, since the lattice spacing a s appears in the gauge structure of the model.Here we rely on information from the lattice model, where the elements of the K matrix are quantized according to Eq. ( 10), and impose an analogous quantization condition for the boundary fields.Remarkably, our construction suggests a generalization to fracton models described by higher-level K matrices [20], which can harbor excitations with fractional charge.
Here we focus on the model described by the K matrix in Eq. ( 8) corresponding to the 2 Chamon model.
The symmetries of the model strongly constrain the correlations of charged operators.Let us consider n-point functions of the form The shift symmetries in Eq. ( 42) impose that the correlations are non-vanishing only if with arbitrary functions f l ξ z α and fl ξz α .Beyond the usual charge neutrality condition, Eq. ( 61) implies that the charge of an excitation created by Ψ q α (r α ) is effectively position dependent [56,57].This condition gives us information about the mobility of excitations at the boundary.The immobility of fractons is manifested in the vanishing of the correlation for a set of isolated positions {r α }.On the other hand, correlations of dipoles made up of opposite charges separated by a short distance w ∼ a s in the ξ z direction have the form These correlations have to obey a less restricted rule where f l (x) is the derivative of f l (x).In this case, we obtain a nonzero correlation for two dipoles at positions corresponding to the same value of ξ z but arbitrary values of ξz , as expected because these dipoles can propagate along lines in the ξz direction.Similarly, dipoles composed of charges separated by w ∼ a s in the ξz direction are correlated along the ξ z direction.

Boundary Line Operators
The line operators described in Sec.3.2 must be modified when one of their endpoints lies at the boundary.Analyzing these operators, we would like to verify that the effective field theory is able to recover the boundary processes discussed for the lattice model in Sec. 2. Consider a dipole with charges separated by a small distance w ∼ a s along the ξx direction in a yz plane; cf. the lattice representation in Fig. 5.We bring the dipole from infinity to the boundary on the z = 0 plane along a line Γ in the positive ξ x direction.Let r 0 = (x 0 , y 0 , 0) be the endpoint of Γ on the boundary.The motion of the dipole is described by the operator where ξ x 0 = ξx 0 = y 0 / 2. Using Eq. ( 36) and D yz = 2∂ ξ x ∂ ξx , we can perform the integrals and express the line operators in terms of the ϕ l fields as where r 0 = r 0 + 1 2 (0, w, −w).Expanding for small w, we obtain Using Eq. ( 38) to eliminate the derivative in the z direction, we obtain an expression well defined within the boundary theory: Comparing the above expression with the charged operator in Eq. ( 56), we recognize that the factor e −i 2w a s c ll ϕ l creates a fracton with charge q = 2w a s (c l1 , c l2 ).Charge quantization imposes ( 2w/a s )c l l ∈ .This result provides further interpretation for the boundary condition in Eq. ( 38): The coefficients c l l encode the fusion of a dipoles into a boundary fracton at the endpoint of the line operator.Recall that there are two types of boundaries, distinguished by which sublattice parities are broken and by the types of fractons that can appear at the termination of BHCs.For instance, if dipoles can be converted into fractons with charge q = (1, 0), but not with charge q = (0, 1), we must have c l2 = 0.Moreover, the factor e i 2w∂ y ϕ l (r 0 ) in Eq. (67) can be interpreted in terms of the creation of additional dipoles that propagate along the boundary.Note that we can decompose where the derivatives with respect to ξ z and ξz correspond to the natural orientation of the dipoles in the boundary plane.Similar combinations appear when we consider operators associated with dipoles that propagate to the boundary along a line in an xz plane.

Perturbations and gapping conditions
The lattice model discussed in Sec. 2 has gapped boundary modes since, like in the bulk, defects of the boundary stabilizers cost finite energy.To recover this result within the continuum theory, we must consider perturbations that break the U(1) symmetries of the Hamiltonian in Eq. ( 50) by creating charged excitations at the boundary.Defining such gapping terms can also be motivated by an analogy with the standard theory of K matrix Chern-Simons.In (2+1) dimensions, a bulk M described by an odd number of copies of Chern-Simons theories will always have protected gapless modes at the edge, ∂ M, which are described by copies of chiral boson theories [54].This is a consequence of the gravitational anomaly (non-zero thermal Hall conductance) [58] that is standard in such theories when defined in a manifold with a boundary.For such theories the presence of the gravitational anomaly makes it impossible to fully gap the edge modes.This is no longer the case when the bulk is described by an even number of Chern-Simons theories.In this situation there can be anomaly cancelation, due to the even number of chiral bosons at the edge, resulting in a vanishing thermal Hall conductance.Therefore, there is no obvious obstruction to gapping the edge modes.In fact, these edge theories can be gapped given that the gapping terms obey a list of criterias [39,41].
Inspired by the lower dimensional case we follow a similar route.The boundary theory ( 50) is a generalization of the standard chiral boson theories.The even number of fields at the boundary suggests that there can be an anomaly cancelation and therefore one can look for appropriate interactions that fully gap the boundary modes.These gapping terms are not arbitrary and must obey some criterias in order for the gapped spectrum to be stable, in analogy to the standard cases [39,41].
We start with a perturbation of the form with g q > 0. This interaction corresponds to a process that creates or annihilates defects with charge q = (q 1 , q 2 ).If we assume that the charges are quantized with q l ∈ , the interaction breaks the continuous symmetries in Eq. ( 42) down to the discrete subsystem symmetries with q l m l (ξ z ), q l ml ( ξz ) ∈ .
Let us analyze the effects of δH q in the strong coupling limit g q → ∞.In a semiclassical picture, the cosine potential pins the fields to one of its minima.The expansion of the cosine around ϕ l = 0 generates a mass term δH q ∼ 1 2 g q (q l ϕ l ) 2 .For small momentum, we obtain the dispersion with a gap given by ∆ = (g q M l l q l q l ) 1/2 .Note that ∆ > 0 for any q since M is positive-definite.More generally, the boundary action may contain multiple cosine terms with different charges.It is possible to consistently minimize two cosine potentials δH q and δH q if they commute with one another, which is ensured by the "null statistics" condition [see Eq. ( 59)] q l (K −1 ) ll q l ∈ 16 . (71) For q = q, this condition is trivially satisfied because K −1 is antisymmetric.Interaction terms with q = q that obey Eq. ( 71) are called compatible [41].In this case, the interactions generate a stable gap in the dispersion as in Eq. ( 70).Otherwise, if the interactions are incompatible, the attempt to pin the potentials simultaneously gives rise to singular terms such as ω 2 ∼ (p 2 x − p 2 y ) −1 , which invalidates the expansion of the cosines.
For the Chamon model, the null statistics condition implies that the cosine terms associated with elementary charges q = (1, 0) and q = (0, 1) are incompatible.To generate a gap, we can pin either ϕ 1 = 0 or ϕ 2 = 0, depending on which perturbation has a larger coupling constant.We can interpret this result in terms of the two types of (001) boundaries, distinguished by the types of fractons with broken parity constraints.These gapping conditions can be straightforwardly generalized to boundary theories where the scalar field ϕ l contains more than two components [20].
So far we have considered the limit of large g q , but another important question concerns the relevance of the boundary interactions at weak coupling.Here we should note that the action in Eq. ( 43) is scale invariant with dynamic exponent [t] = 2.However, unlike the quantum Lifshitz theory in 2 + 1 dimensions [59,60], our fractonic theory lacks continuous rotation invariance in the boundary plane and the vanishing of the dispersion along lines in momentum space poses challenges for a standard renormalization group analysis [61].Nevertheless, we may explore a possible effective scaling dimension of the operator in Eq. ( 68) by calculating its correlation function in the unperturbed model.We proceed in analogy with the computation of correlators of vertex operators in conformal field theory [62].We find that the equal-time correlation of charged operators vanishes exactly (see Appendix B): This result can be traced back to the subsystem symmetries and the generalized neutrality condition in Eq. ( 61), which is not satisfied for any r = 0. Thus, the correlation of charged operators is effectively short-ranged in all spatial directions.Similar behavior appears in effective theories for the Bose metal phase in two dimensions [63] and the classical plaquette-dimer model in three dimensions [61,64].The short-range correlations indicate that δH q is irrelevant at weak coupling.This conclusion is supported by a perturbative renormalization group analysis [64], adapted to incorporate the anisotropic dispersion and the large-momentum contributions to the low-energy physics.
However, the discussion in Section 4.2 suggests that operators that create dipoles can have non-vanishing correlations.For example, consider the perturbation where we set the same coupling constant for both terms assuming that the boundary preserves the C 4 rotation symmetry that takes ξ z → ξz , ξz → −ξ z .Taking two points along a line with fixed ξ z , we obtain the correlation (see Appendix B) where α s ∼ a s is a short-distance cutoff and the exponent is given by On general grounds, we expect a relevant operator to be associated with slowly decaying correlations, i.e., a small value of η q [61].This reasoning can be justified more rigorously by an RG analysis as discussed in Ref. [64].Note that the exponent η q is tied to the ratio between nonuniversal microscopic parameters.Nevertheless, we can make some qualitative considerations assuming that α s is of the same order as the dipole length w.
Clearly, the effective scaling dimension of the dipole operator increases with the charges q l ∈ .For the operators with charges (n, 0) and (0, n), we have Similarly to the quantum Lifshitz theory [60], the relevance of the perturbation depends on parameters that govern the dispersion relation of the excitations, see Eq. (46).To understand the dependence on the matrix elements of M , let us first discuss the case where M is diagonal.In this case, we have µ = M 11 M 22 .If we fix w/α s ≈ 1, the exponents in Eq. ( 76) only depend on the ratio M 11 /M 22 .Decreasing M 11 /M 22 increases η (0,n) , while making η (n,0) smaller.Thus, we cannot make both exponents arbitrarily large at the same time.We take this as a sign that, if M is diagonal in the "flavor" basis in which dipoles are created, the theory contains at least one relevant operator at weak coupling.In contrast, for non-diagonal M we have M 12 as an additional parameter.In this case, for fixed M 11 and M 22 , we can get η q → ∞ for all charges by taking |M 12 | → M 11 M 22 , so that µ → 0. In this regime, all dipole perturbations become irrelevant, regardless of the choice of w/α s .Note that this condition can happen close to the edge of stability of the boundary theory, since the spectrum would become imaginary for Our effective field theory then suggests the existence of a stable gapless boundary phase at weak coupling and for sufficiently small µ.Provided that all perturbations are irrelevant, the low-energy boundary modes exhibit emergent continuous subsystem symmetries as described by the quadratic Hamiltonian in Eq. ( 50).Even though we only encountered gapped boundaries in the analysis of the lattice model in Sec. 2, the physical conditions for finding gapless boundary modes can be illuminated by an analogy with the quantum Hall effect.From the perspective of a hydrodynamic approach [54], the velocity of the chiral edge mode in a quantum Hall fluid is proportional to the electric field of the confining potential near the edge.As a consequence, electrons propagate with a lower velocity when the confining potential varies smoothly in the direction perpendicular to the edge.Along the same lines, we expect that in our theory the dispersion parameter µ could be smaller for a smoothly varying boundary between the fractonic phase and a trivial phase.The spatial dependence can be controlled by considering, for instance, the model defined in infinite volume including a Zeeman term H Z = − r h(r) • σ r with an inhomogeneous magnetic field h(r) = λ(z − z 0 ) 2 Θ(z − z 0 )ẑ, where Θ(x) is the Heaviside step function.The solvable lattice model of Sec. 2 corresponds to a sharp boundary in the limit λ → ∞.In contrast, a smooth boundary characterized by a small value of λ introduces a small energy scale in the problem and should give rise to slower boundary modes.Alternatively, one may try to tune the dispersion parameter µ by adding suitable boundary perturbations on a sharp boundary so as to enhance the flavor mixing contained in the off-diagonal matrix element M 12 .
Starting from the fixed point of gapless boundary modes, the system can undergo a fractonic Berezinskii-Kosterlitz-Thouless (BKT) transition when the dipole operators become relevant [61,64].This transition is driven by the proliferation of dipole excitations.If we assume that g q flows to strong coupling and pin the cosine potential in Eq. ( 73), the expansion around the minimum generates terms proportional to ∂ ξ z ϕ l 2 and ∂ ξz ϕ l 2 .These terms break the subsystem symmetries in Eq. ( 42), destroying the fracton physics at the boundary and allowing the excitations to recover full mobility.It is worth mentioning that a similar transition in 2+1 dimensions has been discussed in the context of the Xu-Moore model for p + i p superconducting arrays [65,66].In the latter, it has been argued that the subsystem symmetries imply a dimensional reduction in the effective theory for the transition.Note that we can in principle push the boundary of the Chamon model towards the BKT-type transition by decreasing the parameter M 12 .On the other hand, the gapless phase is also destabilized if |M 12 | is large enough so that the dispersion parameter µ becomes imaginary.Therefore, we expect the gapless phase to appear as an intermediate phase as a function of a control parameter that tunes M 12 .We leave a detailed study of boundary phase transitions in the Chamon model for future work.

Discussion
We have examined the boundary theory of the Chamon model from the lattice and continuum perspectives.We focused our attention on a specific type of boundary on a (001) plane where five-site boundary stabilizers violate two out of four parity constraints.The effective field theory for the boundary was formulated in terms of a two-component bosonic field and resembles the usual K-matrix Chern-Simons descriptions.While our boundary action is compatible with a recent study of the X-cube model [45], the boundary condition in the Chamon model involves the normal derivative of the boundary field.Physically, this boundary condition represents processes that convert dipoles into fractons.We believe that our formulation can be useful to describe other fracton boundaries, for example boundaries of the X-cube model along planes other than the (001) plane.
We also analyzed the effects of perturbations to the quadratic boundary Hamiltonian.At strong coupling, we found that charged operators can open a gap in the boundary spectrum.The theory allows for (at least) two types of gapped boundaries, obtained by pinning either component of the boundary field.We expect the exactly solvable lattice model with a sharp boundary to be in the strong-coupling regime since the model does not contain any small parameter.On the other hand, the weak-coupling regime might be realized in a generic model where perturbations introduce quantum fluctuations and the boundary can be smooth on the scale of the lattice parameter.Analyzing the weak-coupling regime, we found that the charged operators that create monopoles are always irrelevant because their correlations are short ranged in all spatial directions.In this regime, the leading perturbation is related to the creation of dipoles, and it can also be irrelevant depending on the model parameters.As a result, the effective field theory suggests the existence of a stable gapless boundary phase which is beyond the lattice model with a sharp boundary.
Our results raise several questions for future work.One question is whether other types of boundaries of fracton phases can also be described by a generalization of the K-matrix theory discussed here.For instance, it would be interesting to investigate n models that may harbor fractional excitations at the boundary, in analogy with the chiral edge modes of fractional quantum Hall fluids.In addition, the boundary processes are geometry dependent, and different boundaries can be obtained by simply changing the direction of the surface plane.A more ambitious goal would be to map out boundary phase diagrams of fractonic models.For topological phases in 2+1 dimensions, it is possible to tune boundary interactions to drive phase transitions without closing the bulk gap, and the universality classes are described by conformal field theories familiar from symmetry-breaking transitions in 1+1 dimensions [42].One may then wonder whether some exotic 2+1 critical theories might be realized at the boundary of fractonic systems.

B Correlation Functions
We now want to examine the correlation functions of the model.To start, we notice that the mode expansion for the ϕ l fields can be written as Here we have introduce a small-momentum cutoff p c to regularize the infrared divergence of the integral.We must take p c → 0 at the end of the calculation of correlation functions of local operators.We can also calculate the correlation of the derivative of the fields, which are associated with dipoles.In particular, we have We obtain a similar expression for 〈∂ ξz ϕ l (r)∂ ξz ϕ l (0)〉 by exchanging ξ z ↔ ξz .Quadrupoles are associated with the second derivative D x y ϕ l .In this case, we obtain a p c -independent term in the correlator that decays as a power law at large distances.For |x 2 − y 2 | α 2 s , we obtain Along the directions y = ±x, the correlator changes sign and decays more slowly: (96) We are now in position to discuss the correlations of charged operators.To recover the neutrality condition within a direct calculation of correlators of vertex operators, one must be careful with the dependence on the infrared cutoff p c in the field propagators [62].We are interested in correlators of the form e iq l ϕ l (r) e −iq l ϕ l (0) = e q l ql 〈ϕl(r)ϕ l (0)〉− q l q l 2 〈ϕl(r)ϕ l (r)〉− ql ql 2 〈ϕl(0)ϕ l (0)〉 .(97) Using the result in Eq. (93), we obtain e iq l ϕ l (r) e −iq l ϕ l (0) = exp − C ll 32π (q l − ql )(q l − ql ) ln We can now see that the double logarithmic behavior is responsible for the vanishing of the correlator for monopoles.Even if we impose the neutrality condition ql = q l , the argument of the exponential on the right-hand side of Eq. ( 98) retains a dependence on p c that forces the correlator to vanish when we take the limit p c → 0. As a result, the monopole operators are always short-range correlated.As discussed in Sec.4.2, this type of correlation could only be nonzero if the generalized neutrality condition in Eq. ( 61) were satisfied, which is not the case for two monopoles created at different positions.
On the other hand, the correlator of dipoles does not share this problem.The reason is that 〈∂ ξ z ϕ l (r)∂ ξ z ϕ l (0)〉 scales as a simple logarithm, see Eq. (94).Performing the same analysis with regularized propagators, we obtain a power-law decay along the direction in which the dipoles can move: e iwq l ∂ ξ z ϕ l (ξ z =0, ξz ) e −iwq l ∂ ξ z ϕ l (0) = (p c α s ) w 2 C ll (q l −q l )(q l −q l ) Setting ql = q l is a sufficient condition to prevent this correlator from vanishing in the limit p c → 0. For q = q ∈ 2 , the imaginary part in the off-diagonal matrix elements of C does not contribute to the exponent, and we obtain the result in Eq. ( 74) of the main text.

< l a t e x i t s h a 1 _
b a s e 6 4 = " j q y l w c x B a g M h e s C Q 8 V F O W E k S c m o = " > A A A B 6 H i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y N U Y 9 E L x 4 h k U c C G z I 7 9 M L I 7 O x m Z t Z I C F / g x Y P G e P W T v P k 3 D r A H B S v p p F L V n e 6 u I B F c G 9 f 9 d n J r 6 x u b W / n t w s 7 u 3 v 5 B 8 f C o q e N U M W y w WM S q H V C N g k t s G G 4 E t h O F N A o E t o L R 7 c x v P a L S P J b 3 Z p y g H 9 G B 5 C F n 1 F i p / t Q r l t y y O w d Z J V 5 G S p C h 1 i t + d f s x S y O U h g m q d c d z E + N P q D K c C Z w W u q n G h L I R H W D H U k k j 1 P 5 k f u i U n F m l T 8 J Y 2 Z K G z N X f E x M a a T 2 O A t s Z U T P U y 9 5 M / M / r p C a 8 9 i d c J q l B y R a L w l Q Q E 5 P Z 1 6 T P F T I j x p Z Q p r i 9 l b A h V Z Q Z m 0 3 B h u A t v 7 x K m h d l 7 7 J c q V d K 1 Zs s j j y c w C m c g w d X U I U 7 q E E D G C A 8 w y u 8 O Q / O i / P u f C x a c 0 4 2 c w x / 4 H z + A O k n j Q U = < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " + u Q y N R f l h 6 Z f p B t 0 O s l + e 4 s j u B k = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 m k q M e i F 4 8 t 2 F p o Q 9 l s J + 3 a z S b s b o Q S + g u 8 e F D E q z / J m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 I B F c G 9 f 9 d g p r 6 x u b W 8 X t 0 s 7 u 3 v 5 B + f C o r e N U M W y x W M S q E 1 C N g k t s G W 4 E d h K F N A o E P g T j 2 5 n / 8 I R K 8 1 j e m 0 m C f k S H k o e c U W O l 5 q R f r r h V d w 6 y S r y c V C B H o 1 / + 6 g 1 i l k Y o D R N U 6 6 7 n J s b P q DK c C Z y W e q n G h L I x H W L X U k k j 1 H 4 2 P 3 R K z q w y I G G s b E l D 5 u r v i Y x G W k + i w H Z G 1 I z 0 s j c T / / O 6 q Q m v / Y z L J D U o 2 W J R m A p i Y j L 7 m g y 4 Q m b E x B L K F L e 3 E j a i i j J j s y n Z E L z l l 1 d J + 6 L q X V Z r z V q l f p P H U Y Q T O I V z 8 O A K 6 nA H D W g B A 4 R n e I U 3 5 9 F 5 c d 6 d j 0 V r w c l n j u E P n M 8 f 6 q u N B g = = < / l a t e x i t > y < l a t e x i t s h a 1 _ b a s e 6 4 = " V c d u h I m t G 3 1 x t I w C R H / 3 q E w y 9 5 Y = " > A A A B 6 H i c b V D L T g J B E O z F F + I L 9 e h l I j H x R H Y N U Y 9 E L x 4 h k U c C G z I 7 9 M L I 7 O x m Z t Y E C V / g x Y P G e P W T v P k 3 D r A H B S v p p F L V n e 6 u I B F c G 9 f 9 d n J r 6 x u b W / n t w s 7 u 3 v 5 B 8 f C o q e N U M W y w W M S q H V C N g k t s G G 4 E t h O F N A o E t o L R 7 c x v P a L S P J b 3 Z p y g H 9 G B 5 C F n 1

Figure 1 :
Figure 1: Chamon model on the fcc lattice.(a) Black dots represent sites in Ω even where the spin operators act.The positions in Ω odd are divided into A, B, C and D sublattices, represented by the colors red, green, blue, and yellow, respectively.Every position in Ω odd is associated with an octahedron and the corresponding stabilizer defined as in Eq. (1).(b) fcc lattice of octahedra.We distinguish between two types of (001) planes, in which the octahedra belong to either A and B or to C and D sublattices.

Figure 2 :
Figure 2: Diagrammatic representation of excitations.(a) The action of σ x r on the site marked by a black dot flips the sign of four stabilizers in a yz plane with B and C sublattices.The quadrupole is represented by strings connecting the defects.(b) Successive application of spin operators in the same plane can stretch the strings and move the dipoles in the direction of either ê2 + ê3 or ê2 − ê3 .The excitations reside at the corners of the ribbon structure.

Figure 3 :
Figure 3: (a) Basic diagrams representing the string configurations allowed in the ground state sector.(b) Example of a cage-net tetrahedral structure.

Figure 4 :
Figure 4: A (001) boundary that breaks the parities associated with C and D sublattices.The broken symmetries allow for triangular diagrams contained in xz and yz planes.

Figure 5 :
Figure 5: Processes involving dipoles that propagate to the boundary.(a) A dipole can be brought from the bulk to the boundary through successive combinations of the square diagrams in Fig. 2. (b) The dipole can be converted into a boundary fracton using the triangular diagrams in Fig. 4. (c) The dipole can be reflected at the boundary and change its direction of propagation.This process also implies that a pair of dipoles moving in different directions can annihilate on the boundary if the endpoints of the ribbons coincide.

Figure 6 :
Figure 6: A boundary half cage for the (001) boundary at z = 0.The red and green dots mark the positions of the isolated boundary fractons.