Higher-form symmetry breaking at Ising transitions

In recent years, new phases of matter that are beyond the Landau paradigm of symmetry breaking are mountaining, and to catch up with this fast development, new notions of global symmetry are introduced. Among them, the higher-form symmetry, whose symmetry charges are spatially extended, can be used to describe topologically ordered phases as the spontaneous breaking of the symmetry, and consequently unify the unconventional and conventional phases under the same conceptual framework. However, such conceptual tools have not been put into quantitative test except for certain solvable models, therefore limiting its usage in the more generic quantum manybody systems. In this work, we study Z2 higher-form symmetry in a quantum Ising model, which is dual to the global (0-form) Ising symmetry. We compute the expectation value of the Ising disorder operator, which is a non-local order parameter for the higher-form symmetry, analytically in free scalar theories and through unbiased quantum Monte Carlo simulations for the interacting fixed point in (2+1)d. From the scaling form of this extended object, we confirm that the higher-form symmetry is indeed spontaneously broken inside the paramagnetic, or quantum disordered phase (in the Landau sense), but remains symmetric in the ferromagnetic/ordered phase. At the Ising critical point, we find that the higher-form symmetry is also spontaneously broken, even though the 0-form symmetry is preserved. We discuss examples where both the global 0-form symmetry and the dual higher-form symmetry are preserved, in systems with a codimension-1 manifold of gapless points in momentum space. These results provide non-trivial working examples of higher-form symmetry operators, including the first computation of one-form order parameter in an interacting conformal field theory, and open the avenue for their generic implementation in quantum many-body systems.


I. INTRODUCTION
Global symmetries are instrumental in organizing our understanding of phases of matter. The celebrated Landau paradiagm classifies phases according to broken symmetries, which also determines the universality classes of transitions between phases. Symmetry principles become even more powerful from the point of view of long wavelength, low-energy physics, as the renormalization group fixed points (i.e. IR) often embody more symmetries than the microscopic lattice model (i.e. UV), which is the phenomenon of emergent symmetry [1][2][3][4][5] . A common example is the emergence of continuous spacetime symmetries in the field-theoretical description of a continuous phase transition 6 . It is even plausible that a critical point is determined up to finite choices by its full emergent symmetry, which is the basic philosophy (or educated guess) behind the conformal bootstrap program 7 .
Modern developments in quantum many-body physics have significantly broadened the scope of quantum phases beyond the Landau classification 8 . For these exotic phases, more general notions of global symmetry are called for to completely characterize the phases and the associated phase transitions. Intuitively, these "beyond Landau" phases do not have local order parameters. Instead, non-local observables are often needed to characterize them. For a well-known example, confined and deconfined phases of a gauge theory are distinguished by the behavior of the expectation value of Wilson loop oper-ators 9,10 . To incorporate such extended observables into the symmetry framework, higher-form symmetries [11][12][13] , and more generally algebraic symmetries 14,15 have been introduced. These are symmetries whose charged objects are spatially extended, e.g. strings and membranes. In other words, their symmetry transformations only act nontrivially on extended objects. Most notably, spontaneous breaking of such higher symmetries can lead to highly entangled phases, such as topological order 13 . Therefore, even though topologically ordered phases are often said to be beyond the Landau paradiagm, they can actually be understood within a similar conceptual framework once higher symmetries are included. In addition, just as the usual global symmetries, higher-form symmetries can have quantum anomalies 13 , which lead to strong non-perturbative constraints on low-energy dynamics 16 .
In this work, we make use of the prototypical continuous quantum phase transition, the Ising transition, to elucidate the functionality of the higher-form symmetry. The motivation to re-examine the well-understood Ising transition is the following: in addition to the defining 0form Z 2 symmetry, the topological requirement that Z 2 domain walls must be closed (in the absence of spatial boundary) can be equivalently formulated as having an unbreakable Z 2 (D − 1)-form symmetry, where D is the spatial dimension. Gapped phase on either side of the transition spontaneously breaks one and only one of the two symmetries. Therefore to correctly determine the full emergent internal symmetry in the Ising CFT, the Z 2 higher-form symmetry should be taken into account. For D = 2, the 1-form symmetry manifests more clearly in the dual formulation 17 , namely as the confinementdeconfinement transition of a Z 2 gauge theory, which will shed light on higher-form symmetry breaking transitions in a concrete setting.
A basic question about a global symmetry is whether it is broken spontaneously or not in the ground state. For clarity, let us focus on the D = 2 case. It is wellknown that the Ising symmetric, or "quantum disordered" phase, spontaneously breaks the higher-form symmetry, and the opposite in the Ising symmetry-breaking phase. The fate at the critical point remains unclear to date. To diagonose higher-form symmetry breaking, we compute the ground state expectation value of the "order parameter" for the higher-form symmetry -commonly known as the disorder operator in literature [18][19][20][21][22] , which creates a domain wall in the Ising system. Spontaneous breaking of the Z 2 1-form symmetry is signified by the perimeter law for the disorder operator. In the dual formulation, the corresponding object is the Wilson loop operator. Through large-scale QMC simulations, we find numerically that at the transition, the disorder operator defined on a rectangular region scales as l s e −a1l , where l is the perimeter of the region, and s > 0 is a universal constant. We thus conclude that the 1-form symmetry is spontaneously broken at the (2+1)d Ising transition, and it remains so in the disordered phase of the model. This is in stark contrast with the D = 1 case, where the disorder operator has a power-law decay.
To corroborate the numerical results, we consider generally disorder operator corresponding to a 0-form Z 2 symmetry in a free scalar theory in D dimensions, which is a stable fixed point for D ≥ 3. We show that for the kind of Z 2 symmetry in this case, the disorder operator can be related to the 2nd Renyi entropy. Therefore, the disorder operator also obeys a "perimeter" (i.e. volume of the boundary) scaling, with possibly multiplicative power-law correction. Whether the higher-form symmetry is broken or not is determined by the subleading power-law corrections. We also discuss other free theories, such as a Fermi liquid, where the decay of the disorder operator is in between the "perimeter" and the "area" laws, and therefore no higher-form symmetry breaking.
The rest of the paper is organized as follows. In Sec. II we review higher-form symmetry and its spontaneous breaking, and its relevancy in conventional phases. We also consider higher-form symmetry breaking in free and interacting conformal field theories. In Sec. III we specialize to the setting of quantum Ising model in (2+1)d and define the disorder operator. Sec. IV presents the main numerical results from quantum Monte Carlo simulations, which reveal the key evidence of the 1-form symmetry breaking at the (2 + 1)d Ising transition. Sec. V outlines a few immediate directions about the higherform symmetry breaking and their measurements in unbiased numerical treatments in other quantum many-body systems.

II. GENERALIZED GLOBAL SYMMETRY
Consider a quantum many-body system in D spatial dimensions. Global symmetries are unitary transformations which commute with the Hamiltonian. Typically the symmetry transformation is defined over the entire system, and charges of the global symmetry are carried by particle-like objects.
An important generalization of global symmetry is the higher-form symmetry 13 . For an integer p ≥ 0, pform symmetry transformations act nontrivially on pdimensional objects. In other words, "charges" of pform symmetry are carried by extended objects. In this language, the usual global symmetry is 0-form as the particle-like object is of 0-dimension. p-form symmetry transformations themselves are unitary operators supported on each codimension-p (i.e. spatial dimension (D − p)) closed submanifold M D−p . In particular, it means that there are infinitely many symmetry transformations in the thermodynamic limit. In this work we will only consider discrete, Abelian higher-form symmetry, so for each submanifold M D−p the associated unitary operators form a finite Abelian group G. Physically, higherform symmetry means that the certain p-dimensional objects are charged under the group G, and the quantum numbers they carry constrain the processes of creation, annihilation and splitting etc. In particular, these extended objects are "unbreakable", i.e. they are always closed and can not end on (p − 1)-dimensional objects.
For a concrete example, let us consider (2+1)d Z 2 gauge theory definend on a square lattice. Each edge of the lattice is associated with a Z 2 gauge field (i.e. a qubit), subject to the Gauss's law at each site v: Here e runs over edges ending on v.
The divergence-free condition implies that there are no electric charges in the gauge theory. In other words, all Z 2 electric field lines must form loops. An electric loop can be created by applying the following operator along any closed path γ on the lattice: The corresponding Z 2 1-form symmetry operator is defined as for any closed path γ on the dual lattice. Here the subscript m in W m indicates that this is actually the string operator for Z 2 flux excitations. In field theory parlance, W e is the Wilson operator of the Z 2 gauge theory, and W m is the corresponding Gukov-Witten operator 23 . We notice that the W m (γ ) operator is in fact the product of Gauss's law term v∈e τ x e for all v in the region enclosed by γ . In other words, the smallest possible γ is a loop around one vertex v, and the fact tht W m (γ ) is conserved by the dynamics means that the gauge charge at site v must be conserved (mod 2) as well. Therefore, the Z 2 gauge theory with electric 1-form symmetry is one with completely static charges, including the case with no charges at all. For applications in relativistic quantum field theories, it is usually further required that the 1-form symmetry transformation is "topological", i.e. not affected by local deformation of the loop γ , which is equivalent to the absence of gauge charge as given in Eq. (1).
It is instructive to consider how the 1-form charge of an electric loop can be measured. This is most clearly done in space-time: to measure a p-dimensional charge, one "wraps" around the charge by a (D − p)-dimensional symmetry operator. Appying the symmetry transformation is equivalent to shrinking the symmetry operator, and in (D + 1) spacetime, because of the linking the two must collide, and the non-commutativity (e.g. between W e and W m ) measures the charge value. We illustrate the process for p = 0 ( Fig. 1(a)) and p = 1 ( Fig. 1(b)), in three-dimensional space-time. Now consider the following Hamiltonian of Ising gauge theory: where J, K > 0. When J K, the ground state is in the deconfined phase, which can be viewed as an equal-weight superposition of all closed Z 2 electric loops. In this phase, the Z 2 1-form symmetry is spontaneously broken. When J K, the ground state is a product state with τ x e = 1 everywhere, and the 1-form symmetry is preserved. This is the confined phase. Similar to the usual boson conden-sation, the expectation value of the electric loop creation operator W e (γ) can be used to characterize the 1-form symmetry breaking phase, which obeys perimeter law in the deconfined phase.
This example shows that higher-forms symmetry naturally arises in gauge theories. In condensed matter applications, gauge theories are usually emergent 3,24 , which means that dynamical gauge charges are inevitably present and the electric 1-form symmetry is explicitly broken. Even under such circumstances, at energy scales well below the electric charge gap, the theory still has an emergent 1-form symmetry 25 .
Let us now discuss more generally the spontaneous breaking of higher-form symmetry 13,26,27 . We will assume that the symmetry group is discrete. For a p-form symmetry, a charged object is created by an extended operator W (C) defined on a p-dimensional manifold C. When the symmetry is unbroken, we have where Area(C) is the volume of a minimal (p + 1)dimensional manifold whose boundary is C. t p+1 can be understood as the "tension" of the (p + 1)-dimensional manifold. This generalizes the exponential decay of charged local operator for the 0-form case. On the other hand, when the symmetry is spontaneously broken, where Perimeter(C) denotes the "volume" of C itself. Importantly the expectation value only depends locally on C, which is the analog of the factorization of the correlation function of local order parameter for 0-form symmetry. One can then redefine the operator W (C) to remove the perimeter scaling and in that case W (C) would approach a constant in the limit of large C 28 . At critical point, however, subleading corrections become important, which will be examined below.
The Z 2 gauge theory is famously dual to a quantum Ising model 29 . In fact, more generally, there is a duality transformation which relates a system with global Z 2 0-form symmetry (in the Z 2 even sector) to one with global Z 2 (D − 1)-form symmetry, a generalization of the Kramers-Wannier duality in (1+1)d.
Let us now review the duality in (2+1)d. The dual Ising spins are defined on plaquettes, whose centers form the dual lattice. For a given edge e of the original lattice, denote the two adjacent plaquettes by p and q, as shown in the figure below: The duality map is defined as follows: Note that the expression automatically ensures p σ x p = 1 in a closed system, so the dual spin system has a Z 2 0-form symmetry generated by S = p σ x p , and the map can only be done in the Z 2 even sector with S = 1 30 . Conversely, the mapping also implies v∈e τ e x = 1, and in fact W m (γ ) = 1 for any γ * , i.e. the Z 2 1-form symmetry is strictly enforced.
In the dual model, the electric field line of the Z 2 gauge theory becomes the domain walls separating regions with opposite Ising magnetizations. Therefore, a Wilson loop W e (γ) maps to where ∂M = γ, i.e. M is the region enclosed by γ. Physically X M flips all the Ising spins in the region M , thus creating a domain wall along the boundary γ. It is called the disorder operator for the Ising system, which will be the focus of our study below. Under the duality map, the Hamiltonian becomes The phases of the gauge theory can be readily understood in the dual representation. For K J, the Z 2 gauge theory is in the deconfined phase, which means that the ground state contains arbitrarily large electric loops. For the dual Ising model, the ground state is disordered, with all σ x p = 1. If we work in the σ z eigenbasis (which is natural to discuss symmetry breaking), the ground state wavefunction is given by Namely we pick any basis state and apply the ground state projector. Expanding out the projector, one can see that the wavefunction is an equal superposition of all domain wall configurations, i.e. a condensation of domain walls. Since the domain walls carry Z 2 1-form charges, the condensation breaks the 1-form symmetry spontaneously, much like the Bose condensation spontaneously breaks the conservation of particle numbers In the other limit K J, the gauge theory is confined. Correspondingly, the dual Ising model is in the ferromagnetically ordered phase: there are two degenerate ground states |↑ · · · ↑ and |↓ · · · ↓ . There are no domain walls at all in the limit K → 0. When a small but finite K/J is turned on, quantum fluctuations create domain walls on top of the fully polarized ground states, but these domain walls are small and sparse.
A. Non-invertible anomaly and gapless states A notable feature of the duality map is that on either side, only one of two symmetries, the Z 2 0-form and the Z 2 1-form symmetries, is faithfully represented (in the sense that the symmetry transformation is implemented by a nontrivial operator, even though the duality is supposed to work only in the symmetric sector). The other symmetry transformation is mapped to the identity at the operator level. Physically, only one of them is an explicit global symmetry, while the other one appears as a global constraint (e.g. on the Ising side, domain walls of the 0-form global symmetry are codimension-1 closed manifolds, which is the manifestation that they are charged under a (D − 1)-form symmetry).
A closely related fact is that the ordered phase for one symmetry is necessarily the disordered phase of the other, and any non-degenerate gapped phase must break one and only one of the two symmetries. This has been proven rigorously in one spatial dimension 31 , and is believed to hold in general dimensions as well.
It is clear from these results that these two symmetries can not be considered as completely independent. Recently, Ref. [32] proposed that the precise relation between the two dual symmetries is captured by the notion of a non-invertible quantum anomaly. Intuively, the meaning of the non-invertible anomaly in the context of the Z 2 Ising model can be understood as follows: the charge of the Z 2 0-form symmetry is an Ising spin flip, while the charge of the Z 2 1-form symmetry is an Ising domain wall. These two objects have nontrivial mutual "braiding", in the sense that when an Ising charge is moved across a domain wall, it picks up a minus sign due to the Ising symmetry transformation applied to one side of the domain wall. In other words, the charge of the 1-form symmetry is actually a flux loop of the 0form symmetry. Ref. [32] suggested that two symmetries whose charged objects braid nontrivially with each other can not be realized faithfully in a local Hilbert space. If locality is insisted, then the only option is to realize the D spatial dimensional system as the boundary of a Z 2 toric code model in (D +1) spatial dimension. In this case, the charged objects are in fact bulk topological excitations brought to the boundary. The nontrivial braiding statistics between the two kinds of charges reflects the topological order in the bulk. Such an anomaly is fundamentally different from more familiar 't Hooft anomaly realized on the boundary of a symmetry-protected topological phase (which is an invertible state). We refer to Ref. [32] for more thorough discussions of the non-invertible anomaly.
Since any gapped state must break one of the two symmetries, it is a very natural question to ask whether there are gapless states that preserve both symmetries. An obvious candidate for such a gapless state is the symmetrybreaking continuous transition. At the transition, the two-point correlation function of the Ising order parameter decays algebraically with the distance, implying that the Z 2 0-form symmetry is indeed unbroken. For the dual (D − 1)-form symmetry, the Kramers-Wannier duality maps the disorder operator, which is a string operator in the Ising basis, to the two-point correlator of the Ising order parameter. Therefore the expectation value of the disorder operator also exhibits power-law correlation, and the dual 0-form symmetry is preserved. Therefore the Ising conformal field theory in (1+1)d indeed provides an example of symmetric gapless state with non-invertible anomaly 32 . But for the case of D > 1, the situation is far from clear and that is what we will address in this paper. First we analyze the expectation value of the disorder operator in a free field theory.

B. Scaling of disorder operator in field theory
We now discuss the scaling form of the disorder operator at or near the critical point from a field-theoretical point of view. The natural starting point is the Gaussian fixed point, i.e. a free scalar theory, described by the following Hamiltonian The real scalar φ can be thought of as the coarse grained Ising order parameter, and π is the conjugate momentum of the real scalar φ. The Z 2 symmetry acts as φ → −φ.
The disorder operator X M is basically defined as the continuum version of Eq. (8), where the Z 2 symmetry is applied to a finite region M . Interestingly, for the free theory the expectation value of the disorder operator can be related to another wellstudied quantity, the 2nd Renyi entanglement entropy S 2 . More precisely, for a region M , we have Here S 2 (M ) is the 2nd Renyi entropy of the region M . To see why this is the case, recall that the 2nd Renyi entropy S 2 for a region M of a quantum state |Ψ is given by where ρ M is the reduced density matrix for the region M , obtained from tracing out the degrees of freedom in the complement M : ρ M = Tr M |Ψ Ψ|. In the following we denote the ground wave functional of the state |Ψ by Ψ(φ): The Renyi entropy can be calculated with a replica trick, which we now review in the Hamiltonian formalism. Consider two identical copies of the system, in the state |Ψ ⊗ |Ψ . In the field theory example, the fields in the two copies are denoted by φ (1) and φ (2) , respectively. We denote the basis state with a given field configuration φ (i) in the i-th copy by |φ M is the field configuration restricted to M and similarly φ (i) M for the complement of M . Since the two copies are completely identical, there is a swap symmetry R acting between the two copies R : φ (1) ↔ φ (2) . R M then swaps the field configurations only within the region M : M .
(15) The expectation of R M on the replicated ground state |Ψ ⊗ |Ψ is then given by Therefore the Renyi entropy is the expectation value of the disorder operator for the replica symmetry. For a free theory, we rotate the basis to φ ± = 1 √ 2 (φ (1) ± φ (2) ). In the new basis, the swap symmetry operator becomes: It is straightforward to check that the Hamiltonian of the replica takes essentially the same form in the new basis: The ground state again is factorized: |Ψ ⊗ |Ψ = |Ψ + ⊗ |Ψ − , where |Ψ ± is the state of the φ ± field, with the same wave functional as φ: φ ± |Ψ ± = Ψ(φ ± ) as defined in Eq. (14).
We can now compute the expectation value of R M : where we used the fact that R acts as the identity on φ + . For φ − , R M is nothing but the disorder operator X M . The 2nd Renyi entropy of a free scalar has been wellstudied 33-39 and we summarize the results below.
It is important to distinguish the case where the boundary is smooth and those with sharp corners on the boundary.
First consider a smooth boundary. For a sphere of radius R, in D = 1, 2, 3 we have: Here is a short-distance cutoff, e.g. the lattice spacing, a 1 , a 2 non-universal coefficients and γ a universal constant. For a more general smooth entangling boundary, in 2D the same form holds although the constant correction γ depends on shape of the region. In 3D, it is known that the coefficient of the logarithmic divergent part of the Renyi entropy can be determined entirely from the local geometric data (e.g. curvature) of the surface in a general CFT 40,41 .
If the boundary has sharp corners then there are additional divergent terms in the entropy. The prototypical case is D = 2 when the entangling region has sharp corners. In that case where l is the perimeter of the entangling region and s is an universal function that only depends on the opening angles of the corners. For real free scalar, the coefficient of the logarithmic correction is s ≈ 0.0260 for a square region (so four π/2 corners, as those in Fig. 2) 34,37 . Qualitatively, it is important that for D = 2, 3 the leading term in S 2 always obeys an "perimeter" law, i.e. it only depends on the "area" (length in 2D) of the entangling boundary. If instead we view S 2 as the disorder operator for the Z 2 replica symmetry, the non-universal, cutoff-dependent perimeter term can be removed by redefining the disorder operator locally along the boundary, and the remaining term is universal. For D = 2, the subleading term is either a negative constant when the boundary is smooth, or a ln l correction with a negative coefficient. So according to the relation Eq. (12), the disorder parameter X M , after renormalizing away the perimeter term, does not decrease with the size of M , and therefore the corresponding (D − 1)-form symmetry is spontaneously broken. This is consistent with the fact that the replica symmetry itself must be preserved as there is no coupling between the two copies.
Although the free Gaussian theory is unstable against quartic interactions below the upper critical dimension, and the actual critical theory is the interacting Wilson-Fisher fixed point, results from the free theory can still provide useful insights. It is well-known that for D = 1, for M an interval of length R the disorder operator X M ∼ R −1/4 , the same power-law decay as that of the Ising order parameter due to Kramers-Wannier duality. For D = 2, we will resort to numerical simulations below to address the question.
Notice that the relation between X and S 2 essentially holds for all free theories, including free fermions. For example, the disorder operator associated with the fermion parity symmetry is also equal to S 2 . Interestingly, for a Fermi liquid, it is well-known that ln X = −S 2 ∼ −l D−1 ln l 42,43 , where here l is the linear size of the region. This is an example of a gapless state where the (D − 1)-form symmetry is preserved. Similar results hold for non-interacting bosonic systems with "Bose surface" 44 , an example of which in 2D is given by the exciton Bose liquid 45,46 : In other words, to preserve both the 0-form symmetry and the dual (D − 1)-form symmetry, it is necessary to have a surface of gapless modes in the momentum space. While analytical results discussed in this work are limited to free theories, we conjecture that similar scaling relations hold for interacting CFTs as well. To see why this is plausible, we notice that the entanglement Hamiltonian of a CFT is algebraically "localized" near the boundary of the subsystem 47 , which suggests that even for a nonlocal observable, such as the disorder operator, the major contribution is expected to come from the boundary, and hence a perimeter law scaling. We leave a more systematic study along these lines for future work. In Sec. IV we numerically confirm our conjecture for the Ising CFT in (2+1)d.
We now briefly discuss what happens if a small mass is turned on in Eq. (11). Suppose we are in a gapped phase, and denote by ξ the correlation length. In general, we expect that S 2 obeys an perimeter scaling in the gapped phase, namely the leading term in S 2 is given by a R . In 2D for a disk entangling region of radius R, we have 48 Here a c is the value of a at the critical point (which was denoted by a 1 in Eq. (20)). The function f (x) satisfies Here r is an universal constant (once the definition of ξ is fixed). Suppose the transition is tuned by an external parameter g and the critical point is reached at g c . Since ξ ∼ (g−g c ) −ν where ν is the correlation length exponent, one finds that

III. ORDER AND DISORDER IN ISING SPIN MODELS
In the following we study 1-form symmetry breaking in the transverse field Ising (TFI) model which gives rise to the (2 + 1)d Ising transition. We have reviewed the connection with the Z 2 gauge thory in Sec. II, as well as the 1-form symmetry in the Ising spin system. We will now focus more on the quantitative aspects of the TFI model. Even though the TFI model and the Z 2 lattice gauge theory are equivalent by the duality map, we choose to work with the TFI model here because the numerical simulation is more straightforward.
We will now consider a square lattice with one Ising spin per site, and the global Ising symmetry is generated by S = r σ x r . There are, generally speaking, two phases: a "disordered" phase, where the Ising symmetry is preserved by the ground state 49 , and an ordered phase where the ground states spontaneously break the symmetry. They are separated by a quantum phase transition, described by a conformal field theory with Z 2 symmetry. It is well-understood how to characterize the Ising symmetry breaking (and its absence) in the three cases: consider the two-point correlation function of the order parameter σ z r . The asympotic forms of the correlation function σ z r σ z r for large |r − r | distinguish the three cases: In both the disordered phase and the quantum critical point, the Ising symmetry is preserved because of the absence of long-range order. The prototypical lattice model that displays all these features is the TFI model defined on a square lattice: Note that this is the same as Eq. (9) We will be interested in the disorder operator: where M is a rectangle region in the lattice, illustrated in Fig. 2. In Ref. 14 this operator is called the patch symmetry operator. When X M is applied to e.g. |↑ · · · ↑ , a domain wall is created along the boundary of the region M . These operators are charged under the dual Z 2 1-form symmetry. One can easily see that ψ h=∞ | X M |ψ h=∞ = 1, and ψ h=0 | X M |ψ h=0 = 0. More generally, when M is sufficiently large compared to the correlation length. Here l is the perimeter of the boundary of M , and A is the area of M . The coefficients a and b can be computed perturbatively in the limit of large and small h. In 2D, take M to be a square of perimeter l, so Perimeter(M ) = l and the Area(M ) = l 2 /16. We can find that for large l:

IV. NUMERICAL SIMULATIONS
In this section we study the disorder operator in the (2+1)d TFI model. We employ the Stochastic Series Expansion (SSE) quantum Monte Carlo method 52-55 to simulate the Hamiltonian in Eq. (27). In particular, to be able to directly access the disorder operator in Eq. (28), instead of implementing the algorithm in the conventional σ z basis we choose to work in the σ x basis and construct the highly efficient directed loop algorithm therein 53 . The implementation details of the SSE-QMC algorithm are given in Appendix A.
In our numerical simulations, we choose M to be a rectangular region of size R 1 ×R 2 (i.e. the region contains R 1 R 2 sites), and denote the perimeter l = 2(R 1 + R 2 ). As shown in Fig. 2 (a) and (b), for finite-size studies, we fix the aspect ratio R 2 /R 1 = 1 of square shape and 2 of rectangle shape. The linear system size of the lattice is L and at the critical point we scale the inverse temperature β = 1/T ∼ L to access the thermodynamic limit. O OQ ; 3. − ln( X ) versus l at L = 32 for different h in disordered phases. We use the straight line to fit the data of different external fields and put the obtained slopes in the inset, one sees that as h hc, the fitted slopes (blue circles) approach the predicted relation y = 1 8h 2 (red line). The fitting errors are negligible compared to the circle size.

A. Disordered phase h > hc
First we present results in the disordered phase h > h c . As shown in Eq. (29), we expect that the disorder operator obeys a perimeter law scaling, and for h h c the coefficient is given in Eq. (30). Fig. 3 shows the QMC-obtained ln X M as a function of l for different values of h. The temperature is taken to be β = 10, and we have checked that the results already converge for this value of β. We observe a clear linear scaling, and the inset shows that for large field h h c , the slopes of the ln X M are indeed given by 1/8h 2 asympototically. Now we consider the other limit, when h is approaching the critical point h c from the disordered side. To test the scaling given in Eq. (25), we measure the disorder operator and find the slope a by a linear fit. Fig. 4 shows a c −a as a funtion of h−h c in a log-log plot. A clear power law manifests in the data, and the exponent is found to be ν = 0.63 (2). Considering the finite-size effect, the result agrees very well with the 3D Ising correlation length exponent.

B. Critical point h = hc
The central question to be addressed is whether the Z 2 1-form symmetry is spontanously broken at the critical point. To this end, we measure the disorder operator X at h = h c and scale the inverse temperature β = L in these simulations. We have also checked that finite-β effect is negligible in our calculations. As we have explained, the boundary of M generally contributes to the disorder operator a term proportional to the perimeter. To detect 1-form symmetry breaking, we need to check whether X depends on the area or not. For this purpose, we consider rectangular regions with different aspect ratios: one with 1 : 1 (Fig. 2 (a)) and the other with 1 : 2 ( Fig. 2 (b)), and present the results of X at the h = h c together in Fig. 6. It can be seen that the two sets of data basically fall on the same curve, indicating that the disorder parameter only depends on the perimeter. Given the relation between X M and the Renyi entropy in the free theory, let us examine possible corner contributions to X M , which is parameterized in the coefficient s of Eq. (21). We fit the data points in Fig. 6 to Eq. (21), which yields s = 0.0272 ± 0.004, close to the free value. We perform the same fit for data points with aspect ratio 1 : 2 and obtain essentially the same results (s = 0.0279 ± 0.003). The agreement between the fitting results for regions with different aspect ratios again lends strong support for the perimeter dependence of X M even beyond the leading order, and consequently the 1-form symmetry breaking at the (2 + 1)d Ising CFT.
The convergence of the coefficients a 1 , s and a 0 versus the linear system size L is given in Fig. 8 in Appendix A 3.

C. Ordered phase h < hc
For h < h c where Ising spins order ferromagnetically, our algorithm becomes inefficient because we choose to work in the σ x basis to facilitate the computation of the disorder operator. Nevertheless, simulations indeed find that the disorder parameter decays much more rapidly with the linear size of the region, consistent with the area law in Eqs. (5) and (29). − ln X M as a function of l 2 is shown in Fig. 7  ordered phase, the slope b increases as expected and the data points converge to a straight line for large l 2 . For h = 3.0 very close to the critical point, one can observe that for relatively small values of l 2 the data points do not scale linearly, which can be attributed to a subleading perimeter dependence.

V. CONCLUSION AND DISCUSSION
As discussed in the beginning of the paper, in recent years, new types of quantum phases and phase transitions that are "beyond Laudau" are flourishing, exhibiting topological order, emergent gauge field and fractionalization. Higher-form symmetries and their spontaneous breaking are new conceptual tools introduced to provide an unified framework for both conventional and exotic phases. A quantum phase, gapped or gapless, is fundamentally characterized by its emergent symmetry and the associated anomaly. While the philosophy went back to the Landau classification of phase transitions, the power of this perspective has only begun to unfold recently with the introduction of generalized global symmetries.
Here we re-examine the familiar Ising symmetrybreaking transition, arguably the simplest conformal field theory, from the emergent symmetry perspective. A Ddimensional Ising system has a "hidden" Z 2 (D − 1)form symmetry, whose charges are Ising domain walls. Gapped phases in this system are associated to the spontaneous breaking of the enlarged symmetry (0-form and (D − 1)-form symmetries). It is then of great interest to determine the symmetry breaking pattern at the critical point, to complete our understanding of the global phase diagram from the emergent symmetry point of view.
In this work we determine the scaling form of the disorder operator in Ising CFTs when D > 1. The most challenging case is D = 2 where the transition is described by the interacting Wilson-Fisher fixed point, and we exploit large-scale quantum Monte Carlo simulations. We use the disorder operator of the Ising system to probe the breaking of the dual higher-form symmetry. We find numerically that at the critical point of the 2D quantum Ising model, the one-form disorder operator exhibits sponatenous symmetry breaking as in the disordered phase, whereas in the ordered phase, the one-form symmetry is intact.
The disorder operator is intimately related to a line defect (also called a twist operator) in a Ising CFT, around which the spin operator sees an anti-periodic boundary condition. In fact, a line defect is nothing but the boundary of a disorder operator. It is believed that in general such a line defect can flow to a conformal one at low energy, which is indeed consistent with a perimeter law scaling for the expectation value of the disorder operator 56 . Local properties of disorder line defects have been previously investigated in Ref. [57] and [58]. It will be interesting to understand the relation between the local properties with the universal corner contributions to the disorder operator 59 .
Our findings, besides elucidating the physics of quantum Ising systems from a new angle, provides a working example of higher-form symmetry at practical use. Similar physical systems can be studied, for example, the disordered operator constructed in this work is readily generalized to the (2 + 1)d XY transition and can be measured with unbiased QMC simulations. Another important direction is to study other higher-form symmetry breaking transitions, such as 1-form symmetry breaking transition in 3D systems. It would also be interesting to investigate the ultility of the disorder operator in the topological Ising paramagnetic phase. More applications in quantum lattice models are awaiting to be explored, and will certainly lead to new insight for a new framework that unifies our understanding of the exotic quantum phases and transitions going beyond the Landau paradigm and those within.
Note added.-We would like to draw the reader's attention to few closely related recent works by X.-C. Wu, C.-M. Jian and C. Xu 20,22 and by some of the present author on scaling of disorder oeprator at (2 + 1)d U(1) quantum criticality 21 . zero matrix elements for site operators and bond operators are The updating scheme 52 includes the diagonal update which either inserts or removes a diagonal operator between two states with probabilities regulated by the detailed balanced condition, and the cluster update which flips all the spins on the cluster with the Swendsen-Wang Scheme. The configurations of the updating scheme are shown in Fig. 8. We describe the updating scheme in the following steps:

Diagonal update
We go through the operator strings and either remove or insert a diagonal operator according to the following procedures.
(a) For a diagonal operator (H 0,a or H 1,a ), we removed it with probability, where N denotes the number of lattice sites, and N b denotes the number of bonds. (b) For a null operator (H 0,0 ) , we substitute it with a diagonal operator H 1,a or H 0,a by the procedures below. i. Firstly we make the decision of which kind of diagonal operators to insert. We choose the type of H 1,a with probability or the type H 0,a with probability 1−P (h).
ii. After the decision is made, accept the insertion of an operator with probability and after that we choose an random and appropriate site or bond to insert the operator. If the chosen bond to insert a bond operator has an anti-parallel configurations, then the insertion of a bond operator at this place is prohibited.
(c) For an off-diagonal operator, we ignore it and go to the next operator in the operator strings.

Cluster update
(a) We generally follow two rules to construct the clusters: (1) clusters are terminated on siteoperators H −1,a or H 0,a ; and (2) the four legs of a bond operator H 1,a belong to one cluster. Carry out this procedure until all the clusters are bulit, and a configuration of clusters are shown in Fig. 8.
(b) Clusters identified from the above rules are then flipped with probability 1/2 (which is the Swendsen Wang cluster updating scheme).
Since the disorder operator is a product of σ x , i.e., X M = r∈M σ x r , it is a measurement of an offdiagonal operator in the {σ z } basis. In the σ z basis, the off-diagonal operator can be measured if the operator is a product of operators in the Hamiltonian. It is proved in Ref. 53 that would grow to a very large value as m increases, thus N (k 1 , ..., k m ) would be too small to measure within the limited computing power. So the measurement of X in the {σ z } basis seems hopeless. To solev this problem, we need to change the basis to make σ x diagonal.

SSE on σ x basis
Since we need to measure the disorder operator which is defined as the non-local product of off-diagonal operators, it is extremely hard to measure it in the traditional σ z basis, we then turn to σ x basis as the complete set of basis of the system, and we can use directed loop algorithms 53,54 to simulate this model.