Introduction

Since the discovery of the fractional quantum Hall effect1,2, topological phases of matter have been intensively explored. Exactly solvable models have been constructed3,4,5, and a mathematical framework for classifying topological phases of matter has been developed6,7. Recently, various topologically ordered states have been experimentally investigated with quantum computers and simulators as well8,9,10. Unlike conventional phases of matter, topological phases cannot be characterized by simple local order parameters. Instead non-local string order parameters are required, for example, as proposed by Fredenhagen and Marcu (FM) in the context of the \({{\mathbb{Z}}}_{2}\) gauge-Higgs model11,12,13,14. Recently, this FM string order parameter has also been used to detect \({{\mathbb{Z}}}_{2}\) topological order in a quantum dimer model realized by the Rydberg quantum simulator9,15. The construction of the FM order parameter is based on the following idea. In the presence of an exact Wilson loop symmetry, which is an exact 1-form symmetry16,17,18,19,20,21, the charge condensation transition can be detected using a string order parameter carrying charges at the two ends. However, in the absence of the exact 1-form Wilson loop symmetry, the expectation value of the string order parameter decays exponentially to zero with the string length. Therefore, it cannot be used to detect phase transitions anymore. This problem is circumvented by the FM string order parameter in which the expectation value of the string operator is divided by the square root of the expectation value of a Wilson loop operator of doubled length, such that the exponential decay caused by explicitly violating the exact 1-form Wilson loop symmetry is canceled, see Fig. 1a.

Fig. 1: FM string order parameter and phase diagram.
figure 1

a Schematic of the FM string order parameter associated with charge excitations, where black lines indicate the underlying lattice on a torus, red lines are a string and a loop of Z operators, respectively, and 1 is the trivial anyon flux penetrating through the torus. b Phase diagram of the (2+1)D toric code model in a field, including a toric code phase (blue region), a duality symmetry breaking phase (green dotted line along MK) and a trivial phase (white+red+gray regions). The trivial phase is separated into three regions by emergent 1-form symmetries: the charge condensation region (white) with the emergent Wilson loop symmetry and the flux condensation region (gray) with the emergent 't Hooft loop symmetry, as well as a region without any emergent 1-form symmetry (red). Boundaries of regions in the trivial phase are shown schematically. Inset: definition of the vertex and the plaquette operators on a square lattice.

Full size image

In contrast to the usual definition of the order parameters for spontaneous symmetry-breaking phases and symmetry-protected topological phases, the FM order parameter is not defined based on the exact symmetry of the model, and the FM string order parameter is a non-linear function of the ground states. Therefore, it is not guaranteed that the FM string order parameter is quantitatively well-behaved, i.e., one cannot ensure that the FM string order parameter is smooth in a gapped phase and continuous when crossing a second-order phase transition. Even more so, it is unclear whether the FM string order parameter exhibits critical behavior, as usual order parameters do, in the vicinity of a quantum critical point.

In this work, we study properties of the FM string order parameter of the toric code model in a field and elucidate the role of emergent higher-form symmetries. We evaluate the FM string order parameter associated with charge excitations in the limit of an infinitely long string using transfer matrices of infinite Projected Entangled Pair States (iPEPS). Technically, we approximate ground states of the model using variational iPEPS optimization22,23, in which energy gradients are calculated by automatic differentiation24. We find that near the topological transition between the toric code and the Higgs region, Fig. 1b, the FM string order parameter exhibits a critical exponent of the (2 + 1)D Ising* universality class25. In the presence of the electric-magnetic (EM) duality symmetry, the phase transition has been previously studied with a local order parameter26,27,28. We show that the FM string order parameter exhibits critical scaling in the vicinity of this phase transition but with a new critical exponent, which is different from that of the local order parameter. The universality of this multi-critical point has been recently under debate26,27,28,29 because although the numerical results from the scaling of local operators are consistent with the multi-critical point belonging to the XY* universality class, the semionic statistics between charge and flux is not included in the XY* field theory in an obvious way. We find the critical exponent from the FM string order parameter to be compatible with the (2 + 1)D XY* universality class as well. However, due to the comparatively large numerical uncertainty this does not rule out other field theoretic descriptions of the multi-critical point. In general, our results show that despite being a nonlinear function of the ground state, the FM string order parameter exhibits universal scaling behavior near critical points, suggesting that it is a quantitatively well-behaved order parameter for topological phase transitions accompanied by charge condensation.

We emphasize that the FM string order parameter associated with charge excitations can only be used to reliably detect topological phase transitions when the ground states possess an emergent 1-form Wilson loop symmetry in the infrared (blue and white shaded regions in Fig. 1b). Thus, some prior knowledge of the underlying higher-form symmetries are needed when applying the FM string order parameter to a model with an unknown phase diagram. For instance, in the flux condensation region (gray region in 1b) in which there is no emergent 1-form Wilson loop symmetry, the FM order parameter can be a discontinuous function (i.e., it jumps from zero to a finite value) of ground states even in the absence of phase transitions. Moreover, in the absence of an emergent 1-form symmetry we find the FM order parameter to be numerically unstable for variational iPEPS, see Supplemental Materials30 (similar observations were found in Monte Carlo simulations12,31). In particular, a dual FM string order parameter associated with magnetic flux is required to detect the topological flux condensation transition where an emergent 1-form ’t Hooft loop symmetry exists [blue and gray shaded region in Fig. 1b].

Results

Toric code model and the FM string order parameter

We consider the toric code model on a square lattice in a field:

$${H}_{{\rm{TC}}}({h}_{x},{h}_{z})=-\sum _{v}{A}_{v}-\sum _{p}{B}_{p}-{h}_{x}\sum _{e}{X}_{e}-{h}_{z}\sum _{e}{Z}_{e},$$
(1)

where Av = ∏e∈vXe and Bp = ∏e∈pZe are the vertex and plaquette operators, as shown in the inset of Fig. 1b, and Xe and Ze are Pauli matrices defined on the edges e of the lattice. The phase diagram of the model in Fig. 1b consists of a toric code phase at weak fields and a trivial phase at strong fields32,33,34.

Recent work has introduced higher-from symmetries to describe topological order16,17,18,19,20,21. For certain choices of the magnetic fields, the toric code model in Eq. (1) has exact 1-form symmetries that commute with the Hamiltonian: (i) for hx = 0, the 1-form Wilson loop symmetry commutes with HTC(0, hz), [∏e∈LZe, HTC(0, hz)] = 0, where L is a closed loop on the primal lattice and (ii) for hz = 0 the 1-form ’t Hooft loop symmetry commutes with \({H}_{{\rm{TC}}}({h}_{x},0),[{\prod }_{e\in \hat{L}}{X}_{e},{H}_{{\rm{TC}}}({h}_{x},0)]=0\), where \(\hat{L}\) is a closed loop on the dual lattice. (On a torus, the 1-form symmetries include both contractible and non-contractible loop operators, while on a sphere the 1-form symmetries only consist of contractible loop operators.) Crucially, even when the exact symmetries are explicitly broken, higher-form symmetries can still be emergent at low energies, capturing the robustness of topological order. Concretely, both of the 1-form symmetries are emergent in the toric code phase, see Fig. 1b. In this formalism, topological order can be interpreted as a spontaneous breaking of such emergent symmetries, drawing analogies to conventional symmetry breaking phases19,20,21,35,36. Furthermore, the trivial phase can be separated into several regions by emergent 1-form symmetries26,35. The charge condensation region (white area in Fig. 1b), also known as the Higgs region, only has the emergent 1-form Wilson loop symmetry. The flux condensation region (gray area in Fig. 1b), also known as the confined region, only has the emergent 1-form ’t Hooft loop symmetry. The red area in Fig. 1b does not possess any of the two emergent 1-form symmetries.

As the toric code model in Eq. (1) has the exact 1-form Wilson loop symmetry at hx = 0, a string operator \({\prod }_{e\in {L}_{1/2}}{Z}_{e}\) creates a pair of charges at its two ends, where L1/2 is an open string along the primal lattice. Thus, a string order parameter \({\tilde{O}}_{Z}\) can be constructed to detect charge condensation37:

$${\tilde{O}}_{Z}=\mathop{\lim }\limits_{| {L}_{1/2}| \to \infty }\sqrt{| {\tilde{C}}_{Z}(| {L}_{1/2}| )| },\quad \quad {\tilde{C}}_{Z}(| {L}_{1/2}| )=\langle \Psi | \prod _{e\in {L}_{1/2}}{Z}_{e}| \Psi \rangle ,$$
(2)

where \(\left\vert \Psi \right\rangle\) is a normalized ground state of the toric code model, ∣L1/2∣ is the distance between two ends of L1/2. The bulk of the string order parameter \({\tilde{O}}_{Z}\) commutes with the overlapping local Hamiltonian terms but not its ends. In the toric code phase, the string order operator thus creates two anyonic charge excitations that are orthogonal to the ground state, leading to a vanishing string order parameter in the infinite string limit, \({\tilde{O}}_{Z}=0\). By contrast, in the Higgs phase, the string order parameter can be nonzero because the hz field induces charge fluctuations such that charges condense in the Higgs phase. The string order parameter can also be interpreted as the “disorder parameter” of the exact 1-form Wilson loop symmetry. In the deconfined topological phase, where the ground states spontaneously break the exact 1-form Wilson loop symmetry, the string order parameter vanishes. Conversely, the Higgs regime is disordered under the exact 1-form Wilson loop symmetry and the string order parameter is nonzero. An alternative way of interpreting this string operator is by mapping HTC(0, hz) to the (2 + 1)D transverse field Ising model38,39. This transforms \({\tilde{O}}_{Z}\) and \({\tilde{C}}_{Z}\) to the Ising order parameter and its correlation function, respectively. From that also follows directly that the critical point I [see Fig. 1b] at \({h}_{z}={h}_{zc}^{(I)}=0.328474(3)\)40 belongs to the (2 + 1)D Ising* universality class, where the “*" indicates that in the effective Ginzburg-Landau-Wilson theory the order parameter field ϕ can only be created in pairs and ϕ and − ϕ are physically indistinguishable25. Near the critical point I: \({\tilde{O}}_{Z} \sim {({h}_{z}-{h}_{zc}^{(I)})}^{\beta }\) and \(\xi \sim {({h}_{z}-{h}_{zc}^{(I)})}^{-\nu }\), where β = 0.326418(2) is the critical exponent of the order parameter41, ξ is the correlation length defined via \(| {\tilde{C}}_{Z}(| {L}_{1/2}| )-{\tilde{O}}_{Z}^{2}| \sim {e}^{-| {L}_{1/2}| /\xi }\), and ν = 0.629970(4) is the critical exponents of the correlation length41.

When hx ≠ 0, the Wilson loop operator ∏e∈LZe is no longer an exact 1-form symmetry of the toric code model and the bulk of the string \({\prod }_{e\in {L}_{1/2}}{Z}_{e}\) cannot deform freely. Therefore \({\tilde{C}}_{Z}(| {L}_{1/2}| )\) vanishes on either side of the topological phase transition exponentially with the length of the string ∣L1/2∣. However, when hx is small, the toric code model has an emergent 1-form Wilson loop symmetry17,26,35. In the limit of large hx, the 1-form Wilson loop symmetry cannot emerge, which we indicate by the gray and red areas in Fig. 1b. In the presence of an emergent 1-form Wilson loop symmetry, one can in principle conceive to construct a dressed string operator with an extended width42,43. This is however a challenging task in practice. This problem can be circumvented, by dividing out the “bulk” contribution of the string order paramter, as proposed by Fredenhagen and Marcu, leading to the FM string order parameter11,12:

$${O}_{Z}=\mathop{\lim }\limits_{r\to \infty }\sqrt{| {C}_{Z}(r)| },\,\,\quad {C}_{Z}(r)=\frac{\langle \Psi | {\prod }_{e\in {L}_{1/2}}{Z}_{e}| \Psi \rangle }{\sqrt{\left\langle \Psi \right\vert {\prod }_{e\in L}{Z}_{e}\left\vert \Psi \right\rangle }},$$
(3)

where r = ∣L1/2∣ is the length of the string L1/2. Notice that ∣L1/2∣ should be understood as the length of L1/2 instead of the distance between its two ends because the bulk of the Z-string cannot be deformed freely when hx ≠ 0. In addition, for the \({{\mathbb{Z}}}_{2}\) gauge Higgs model, we should replace the string operator \({\prod }_{e\in {L}_{1/2}}{Z}_{e}\) with \({Z}_{v}\left({\prod }_{e\in {L}_{1/2}}{Z}_{e}\right){Z}_{{v}^{{\prime} }}\), where v and \({v}^{{\prime} }\) are two vertices at the ends of the string L1/2. L is a loop whose length is twice the length of the string L1/2, see Fig. 1a. Compared to Eq. (2), the FM string order parameter in Eq. (3) contains a square root of the expectation value of the Wilson loop operator in the denominator. When hx ≠ 0 the bare Wilson loop operator is no longer a symmetry and its expectation values decays in the presence of an emergent 1-form symmetry of charge excitations with a perimeter law (except in the flux condensation phase with exact 1-form ’t Hooft loop symmetry (hz = 0), where the expectation value of the bare Wilson loop operator decays with an area law): \(\left\langle \Psi \right\vert {\prod }_{e\in L}{Z}_{e}\left\vert \Psi \right\rangle \sim \exp (-{\alpha }_{Z}| L| )\), where αZ is the perimeter law coefficient and ∣L∣ is the length of the loop L. The denominator in the FM string order parameter compensates the perimeter law decay of the numerator, such that only the contribution from the endpoints of the string L1/2 is taken into account, where charges are created. We also separately summarize the behaviors of the numerator and denominator of the FM string order parameter in Table 1, from which we immediately see that the FM string order parameter is always 0 when hz = 0 and it cannot detect the phase transition between the deconfined phase and the confined phase. In addition, we emphasize that the FM string order parameter is a nonlinear function of the ground state. Therefore, one needs to carefully analyze whether it can serve as a bona fide order parameter for topological phase transitions.

Table 1 Z string and Z loop operators
Full size table

We now use iPEPS algorithms to directly evaluate the FM string order parameter in the limit of infinitely long strings on a torus, see Methods. When L is a non-contractible loop on the torus, which is convenient for tensor network methods, care has to be taken when evaluating the denominator in Eq. (3), as it could vanish for certain linear combinations of the degenerate ground states. Nonetheless, we are able to show that contractible and non-contractible loops are equivalent when the ground state is a minimally entangled state44 with a trivial anyon flux penetrating through the torus; Fig. 1a. Therefore, we use this choice for the ground state in the numerical evaluation of the FM order parameter; technical details are discussed in the Methods section.

FM string order parameter of the variational wave function

We will now analyze the general properties of the FM string order parameter across a topological phase transition with charge condensation. In particular, we want to analyze whether the FM string order parameter is continuous upon crossing a second-order phase transition, and if so, we will analyze its scaling behavior. To this end, we first consider a cut through the phase diagram at hx = 0.3, that crosses the transition line IM in Fig. 1b at a point J. The criticality of the topological transition along the line IM is expected to be described by the (2 + 1)D Ising* universality class. There are two important parameters that systematically control the error of the approximation when numerically evaluating the FM string order parameter: the bond dimension D of the iPEPS itself and the bond dimension χ of the environment of the iPEPS; the larger bond dimensions provide better approximations. Therefore, we evaluate FM string string order parameter numerically along hx = 0.3 from iPEPS with various bond dimensions, see Fig. 2a. From this we find the critical point J on the hx = 0.3 line is at \({h}_{z}={h}_{zc}^{(J)}=0.335(1)\), consistence with \({h}_{zc}^{(J)}=0.333(1)\) from previous quantum Monte Carlo simulation40. In the toric code phase, the FM string order parameter vanishes while it is finite in the Higgs phase.

Fig. 2: FM string order parameter of the variational iPEPS along hx = 0.3.
figure 2

a The FM string order parameter from iPEPS with various bond dimensions (D, χ), see legend, is continuous across the topological phase transition. b Data collapse of the FM string order parameter with ν = 0.629970(4) and β = 0.326418(2). c Double-log plot used for extracting the critical exponent βFM from the scaling of OZ in the vicinity of the critical point, where we we used \({h}_{zc}^{(J)}=0.335(1)\). Data is obtained by linearly extrapolating OZ in the inverse iPEPS bond dimension 1/D.

Full size image

Crucially, the FM string order parameter is continuous across a second order phase transition, such that we can extract the critical exponent βFM defined via \({O}_{Z} \sim {({h}_{z}-{h}_{zc}^{(J)})}^{{\beta }_{{\rm{FM}}}}\), Fig. 2c, where we obtain OZ by performing linear extrapolation in 1/D and ignore the small χ dependence. The extracted critical exponents βFM = 0.33(5) is consistent with the exponent β = 0.326418(2) of the (2 + 1)D Ising* universality class41, see Table 2 for a summary of the exponents. We furthermore collapse the data from different bond dimensions, as shown in Fig. 2b, based on the theory of finite entanglement scaling45,46,47. The data for D > 2 indeed collapses on a single curve near criticality. These results show that the FM string order parameter exhibits the correct critical behavior controlled by the (2 + 1)D Ising* universality class near the critical point J.

Table 2 Critical exponents of the Ising* field theory in (2+1)D and (2+0)D
Full size table

Previous order parameters for topological phase transitions were defined on the virtual legs of the iPEPS48,49,50,51,52, because there exist the virtual symmetries in terms of matrix product operators53,54,55. However, without the iPEPS representation, these virtual order parameters cannot be obtained, i.e., they are not physical. This should be contrasted with the FM string order parameter that is defined on the physical level and thus also does not depend on the iPEPS gauge.

Next, we consider the FM string order parameter of the variational wave function along the self-dual line hx = hz, where the toric code model in Eq. (1) has a global electric-magnetic duality symmetry, which exchanges the primal lattice and the dual lattice, as well as X and Z. There is a gapped electric-magnetic duality symmetry breaking phase MK along the self-dual line26,32,33. The transition between the toric code phase and the duality symmetry breaking phase is a multi-critical point M shown in Fig. 1b. From our iPEPS simulation we find that M is located at \({h}_{x}={h}_{z}={h}_{zc}^{(M)}=0.3397(2)\), which is again close to 0.340(2) from the quantum Monte Carlo simulations32 and to 0.3406(4) obtained from the higher-order perturbation expansion33. The phase transition crossing the multi-critical point M can be characterized by a local symmetry-breaking order parameter ∣〈X − Z〉∣, which is the difference between expectation values of of Xe and Ze on a single edge e, as shown in Fig. 3a, from which we extract a critical exponent βlocal = 0.83(5) defined via \(| \langle X-Z\rangle | \sim {({h}_{z}-{h}_{zc}^{(M)})}^{{\beta }_{{\rm{local}}}}\), see Fig. 3c. Our result is consistent with recent Monte Carlo simulation26. We also evaluate the FM string order parameter along the self-dual line crossing the multi-critical point M; Fig. 3d. For the FM string order parameter, we extract the critical exponent βFM = 0.34(4) defined via \({O}_{Z} \sim {({h}_{z}-{h}_{zc}^{(M)})}^{{\beta }_{{\rm{FM}}}}\), Fig. 3f, which is quite different from the critical exponent of the local symmetry-breaking order parameter βlocal = 0.83(5).

Fig. 3: Local duality symmetry-breaking and FM string order parameters of the variational iPEPS along the self-dual line hx = hz.
figure 3

a Local order parameter ∣〈X − Z〉∣ from iPEPS with various bond dimensions (D, χ), see legend. b Data collapse of the local order parameter with ν+ = 0.67175(1) and β− = 0.83048(2). c Critical exponent βlocal, where ∣〈X − Z〉∣ is obtained by linear extrapolation in 1/D and \({h}_{zc}^{(M)}=0.3397(2)\). d FM string order parameter from iPEPS with various bond dimensions, see legend. e Data collapse of the FM string order parameter with βϕ = 0.34870(7). f Critical exponent βFM, where OZ is again obtained by linear extrapolation in 1/D.

Full size image

How can we understand the distinct critical exponents βlocal and βFM? Since the phase boundary IM belongs to the Ising* universality class and the FM string order parameter exhibits an Ising critical exponent, we can assume that the FM string order parameter OZ corresponds to a field ϕz of an effective Ginzburg-Landau-Wilson theory describing the Ising* transition. When two Ising* transition lines IM and GM in Fig. 1b meet at the multi-critical point M, once could conceive that the effective field theory of the multi-critical point M is the XY* model with a Lagrangian \({\mathcal{L}}={\left(\partial \phi \right)}^{2}/2+{m}^{2}{\phi }^{2}/2+g{\phi }^{4}/(4!)\)26,27, up to some irrelevant terms, which possesses O(2) symmetry. Here, ϕ = (ϕx, ϕz) is a two-component vector field. We summarize the scaling dimensions of the order parameter fields and relevant critical exponents in Table 3. We find that βFM is consistent with βϕ and βlocal with β− of the XY* field theory. This can be understood by identifying ϕz with the FM string order parameter Oz, and \({\phi }_{x}^{2}-{\phi }_{z}^{2}\) with X − Z26,27,28. We can separately check consistency with the XY* field theory by performing data collapse of the FM string order parameter and the local order parameter using the critical exponents of the field theory; Fig. 3b and e, where we find a reasonable collapse for iPEPS dimension D > 3.

Table 3 Critical exponents of the XY* field theory in (2+1)D and (2+0)D
Full size table

Some comments are in order. First, the two degenerate ground states in the duality symmetry-breaking phase correspond to the predominant condensation of charges (fluxes) satisfying 〈X − Z〉 < 0 (〈X − Z〉 > 0). Because our FM string order parameter is defined with a Z string, it detects the charge condensation, and thus, we should use the charge condensation-dominated state to evaluate the FM string order parameter. In contrast, the dual FM string order parameter associated with the flux excitations is not well-behaved when applied to the same ground state; see Supplemental materials30. Second, ref. 26 argues that the multi-critical point M may not belong to the XY* universality class because the mutual semionic statistics between charges and fluxes is not included in the XY* field theory in an obvious way. Although our numerical data is consistent with the XY* field theory, we cannot rule out other field theoretic descriptions of the transition due to the comparatively large uncertainties in the critical exponents.

FM string order parameter for the deformed toric code state

Instead of variationally solving the ground state Hamiltonian in Eq. (1), one can analytically construct a deformed toric code state, which shares similar physics with the toric code Hamiltonian in Eq. (1). The advantage of the deformed toric code state is that the wavefunction is exact, so many analytical results can be derived. The disadvantage is that the dimensionality of the universality class of the quantum critical points is fine-tuned and reduced by 1 compared to that of the generic quantum critical points56,57. The deformed toric code state is defined as58,59,60:

$$\left\vert \psi ({g}_{x},{g}_{z})\right\rangle =\prod _{e}(1+{g}_{x}{X}_{e}+{g}_{z}{Z}_{e})\left\vert {\rm{TC}}\right\rangle ,$$
(4)

where \(\left\vert {\rm{TC}}\right\rangle\) is a ground state of the fixed point toric code Hamiltonian, and gx and gz are tuning parameters satisfying \({g}_{x}^{2}+{g}_{z}^{2}\le 1\). The phase diagram of the deformed toric code state is similar to that of the toric code Hamiltonian in Eq. (1), as shown in Fig. 4a. Besides the toric code phase, there are two trivial phases with either charge or flux condensation, and they are separated by a Berezinskii-Kosterlitz-Thouless (BKT) transition line. The phase transition lines from the toric code phase to the trivial phase belong to the (2+0)D Ising* universality class because of the dimensionality reduction at the fine-tuned critical points of the wave function. Because the deformed toric code state is mapped to the partition function of the 2D classical Ashkin-Teller model60, the multi-critical point \({M}^{{\prime} }=(1-1/\sqrt{2},1-1/\sqrt{2})\) at the self-dual line gx = gz is described by the free boson conformal field theory (CFT) compactified on orbifolds with a radius \(R=2\sqrt{2}\), which describes the critical point of the 2D XY* theory61. Thus, we can say the multicritical point \({M}^{{\prime} }\) of the deformed toric code state belongs to the (2+0)D XY* universality class.

Fig. 4: FM string order parameter for the deformed toric code wave function.
figure 4

a Phase diagram of the deformed toric code state of Eq. 4. We also plot a color-map of the FM order parameter OZ in the phase diagram. b The FM string order parameter along gx = 0.14 [vertical dashed line in a], calculated using iMPS with various bond dimensions χ. c Data collapse of the FM string order parameter with gzc = 0.2227(3), β = 1/8 and ν = 1. d Critical exponent βFM of the FM string order parameter, obtained by extrapolating OZ for different χ. e Double-log plot of the FM string order parameter OZ as a function of the correlation length ξ of the boundary MPS at the point (g, θ) = (0.65, π/4) in the BKT phase. Oz vanishes in the limit of infinite bond dimension χ → ∞. f Finite and infinite length FM string order parameters along \({g}_{x}^{2}+{g}_{z}^{2}=0.6{5}^{2}\) [dashed quarter circle in a parametrized with \(\theta =\arctan ({g}_{x}/{g}_{z})\)] from the iMPS with χ = 20. The legend shows the size of the area surrounded by the loop operator.

Full size image

Akin to the variational case, we first calculate the FM string order parameter along a line gx = 0.14 crossing the toric code phase and the Higgs phase by contracting the D = 2 exact iPEPS of the deformed toric code state, using the boundary infinite matrix product states (iMPS) with various bond dimensions χ, see Fig. 4b. We extract the critical exponent βFM of the FM string order parameter according to \({O}_{Z} \sim {({g}_{z}-{g}_{zc})}^{{\beta }_{{\rm{FM}}}}\), for which we numerically determine the critical point gzc = 0.2227(3) and extrapolate OZ using different iMPS bond dimensions χ. The extracted βFM = 0.117(5) is close to the β = 1/8 = 0.125 from the 2D Ising* universality class, see Fig. 4d and Table 2. We perform the data collapse to the FM string order parameter in Fig. 4c. These results indicate that the FM string order parameter is a well-behaved order parameter also for the deformed wavefunction. In contrast to the variatinoal case, along the self-dual line the phase transition from the toric code phase to the BKT phase along the self-dual line is not a charge condensation transition, and the FM string order parameter is zero in both the toric code phase and the BKT phase, see Fig. 4e. Therefore, the FM string order parameter cannot be directly used to detect the BKT transition.

The FM string order parameter evaluation of the deformed toric code is numerically stable, since the wave function can be exactly expressed in terms of an iPEPS without the need of variational optimization. We evaluate the FM string order parameter in the entire phase diagram of the deformed toric code state, as shown in Fig. 4a as a color plot, which surprisingly exhibits a sharp transition in the flux condensation phase even in the absence of a quantum phase transition. To emphasize this behavior, we show the FM string order parameter along a path \({g}_{x}^{2}+{g}_{z}^{2}=0.6{5}^{2}\) in Fig. 4f. It is discontinuous at a certain angle \(\theta =\arctan ({g}_{x}/{g}_{z})\approx 0.4\pi\) even in the absence of a bulk phase transition. In order to exclude the possibility that there are some artifacts of our method that cause the discontinuity, we also evaluate the FM string order parameter with finite string length r = ∣L1/2∣ in Fig. 4f. As the string length increases, the results from the thermodynamic limit are obtained, which indeed implies that the FM string order parameter becomes discontinuous for infinite string lengths. Moreover, there is a small drop at the self-dual line in Fig. 4f, which slowly approaches 0 with increasing bond dimension; Fig. 4e.

Although surprising, the discontinuity in the flux condensed phase is possible because the FM string order parameter is a non-linear function of the ground state. The underlying reason for this abrupt change of the FM order parameter is that the parity of the dominant eigenvectors of the transfer matrices, whose overlap determines the FM order parameter, change at the discontinuity30. In the Supplemental Materials30, we show that for the Hamiltonian in Eq. (1) a similar singular behavior of the FM string order parameter can be found when calculating the FM order parameter using an iPEPS constructed via perturbing the infinite-field limit product states in the confined region62. These unexpected results indicate that in the flux condensation region, the FM string order parameter that creates charges at its end cannot be used as an order parameter anymore, and the FM string order parameter does not have a physical significance. We discuss the consequences of this behavior in the next section.

Discussion

We have evaluated the FM string order parameter in the infinitely long string limit using the iPEPS simulation and found that it exhibits universal scaling controlled by the underlying critical points of charge condensation transitions.

Our results indicate that the FM string order parameter can be discontinuous in the flux condensation region, which does not possess an emergent 1-form Wilson loop symmetry, see Fig. 4a and f. Hence, we argue that only in the presence of an emergent 1-from Wilson loop symmetry, the associated FM string order parameter can be a quantitatively well-behaved order parameter for topological phase transitions. Our argument can be formulated using the idea of quantum error correction63,64. A quantum state \(\left\vert \Psi \right\rangle\) having an emergent 1-form Wilson loop symmetry implies that we can apply a recovery map constructed from quantum error correction to \(\left\vert \Psi \right\rangle\) to get an exact 1-form Wilson loop symmetric state \(\left\vert {\Psi }_{0}\right\rangle\) within the same phase of \(\left\vert \Psi \right\rangle\). Applying this recovery map to both the numerator and denominator of the FM order parameter, we have \(\left\langle \Psi \right\vert {\prod }_{e\in {L}_{1/2}}{Z}_{e}\left\vert \Psi \right\rangle \sim {e}^{-\alpha | {L}_{1/2}| }\left\langle {\Psi }_{0}\right\vert {\prod }_{e\in {L}_{1/2}}{Z}_{e}\left\vert {\Psi }_{0}\right\rangle\) and \(\left\langle \Psi \right\vert {\prod }_{e\in {L}_{1/2}}{Z}_{e}\left\vert \Psi \right\rangle \sim {e}^{-\alpha | L| }\)65, where α is a decay coefficient that depends on \(\left\vert \Psi \right\rangle\) and the recovery map. The FM order parameter evaluated with \(\left\vert \Psi \right\rangle\) using Eq. (3) then reduces to the string order parameter evaluated with \(\left\vert {\Psi }_{0}\right\rangle\) using Eq. (2). As a consequence, the FM order parameter is a well-behaved order parameter in the presence of emergent 1-form Wilson loop symmetry. By contrast, if there is no emergent 1-form Wilson loop symmetry, it will not be guaranteed that a recovery map exists. In the confined regime, the magnetic 1-form ’t Hooft loop symmetry is emergent, and one has to construct an FM order parameter of X operators along the dual lattice to detect the phase transition. Therefore, when detecting a topological phase transition using FM string order parameters, either with numerical simulations or experiments, knowledge about the underlying emergent 1-from symmetries is required.

Our work opens several questions and research directions. First, measuring the FM string order parameter in quantum simulation experiments can provide insights into topological phase transitions with charge condensation. Second, in ref. 65, we argue that 1-form symmetries can be detected with quantum error correction. It will be interesting to precisely analyze the boundaries of the three regions of different emergent 1-form symmetries in the trivial phase of the toric code model in Eq. (1), shown schematically in Fig. 1b, and investigate the nature of the transitions between presence and absence of the 1-form symmetries. Third, the FM string order parameter can be applied to different lattice gauge theories66,67,68,69. It is an interesting direction to define the FM string order parameters for Kitaev’s quantum double models3 and Levin-Wen string-net models5 and apply them to study various topological phase transitions driven by anyon condensation50,51,70. Moreover, since the Higgs phase can be interpreted as a phase protected simultaneously by both a 1-form symmetry and a global symmetry71, one could use the FM-type string order parameters to detect such kind of symmetry-protected topological order, especially when the protecting 1-form symmetry becomes an emergent symmetry72.

Methods

iPEPS optimization

In this section, we show the technical details of optimizing the ground states of the toric code model using iPEPS. We approximate a ground state using the iPEPS ansatz proposed in ref. 73. The iPEPS has a 2 × 2 unit cell, and it is parameterized by a rank-5 tensor A (B is obtained from A by a π/2 rotation) with the virtual bond dimension D and the physical dimension d = 2, as shown in Fig. 5a and b.

Fig. 5: iPEPS ansatz and the CTMRG.
figure 5

a The iPEPS ansatz for the ground state of the toric code model with a bond dimension D. b The iPEPS tensors A and B are related by rotation. The iPEPS tensor A is invariant under two reflections, and we can also impose the virtual \({{\mathbb{Z}}}_{2}\) symmetry when necessary, where the red dots are matrices ZD. c A double tensor and another double tensor sandwiching a Z matrix. d The environment of the iPEPS is approximated by the corner tensors (squares) and the edge tensors (rectangles). e The hermitian eigenvalue decomposition (eigh) of the top-left corner of the tensor networks in d, the isometries can be obtained from eigenvectors corresponding to the χ largest eigenvalues (in absolute value). f CTMRG procedures updating the corner and edge tensors using the isometries.

Full size image

We impose the square lattice symmetry onto the tensor A such that the iPEPS tensor is invariant under two reflections Rv and Rh, see Fig. 5b. Because of the symmetry, the number of independent variational parameters is less than 2D4. We can parameterize such a symmetric tensor A using the following method. We first construct the 2D4 × 2D4 matrix representations of Rh and Rv applying on the tensor A, and a projector \({P}_{R}=({\mathbb{1}}+{R}_{h})({\mathbb{1}}+{R}_{v})/4\). A subspace spanned by the eigenvectors of PR with an eigenvalue 1 is \(\{\left\vert {v}_{i}\right\rangle | {P}_{R}\left\vert {v}_{i}\right\rangle =\left\vert {v}_{i}\right\rangle \}\). If 〈vi∣vj〉 ≠ δij, we can orthonormalize them using the QR decomposition. Reshaping \(\left\vert {v}_{i}\right\rangle\) to the tensors with the dimensions D × D × D × D × 2, we can parameterize the iPEPS tensor as A = ∑iλivi, where λ = {λi} are variational parameters. If we also want to impose the virtual \({{\mathbb{Z}}}_{2}\) symmetry to the tensor A as shown in Fig. 5a, we need another projector \({P}_{Z}=({{\mathbb{1}}}_{D}^{\otimes 4}\otimes {{\mathbb{1}}}_{d}+{Z}_{D}^{\otimes 4}\otimes {{\mathbb{1}}}_{2})/2\), where ZD is a D × D matrix representation of the non-trivial element in \({{\mathbb{Z}}}_{2}\), i.e., \({Z}_{D}^{2}=1\); and we consider the projector P = PZPR, where [PZ, PR] = 0. Using the orthonormal basis of the subspace spanned by the eigenvectors of P with eigenvalue 1, we can parameterize the tensor A using λ.

The iPEPS \(\left\vert \Psi (\lambda )\right\rangle\) can be constructed from A(λ). The energy expectation value \(E=\left\langle \Psi (\lambda )\right\vert {H}_{{\rm{TC}}}({h}_{z},{h}_{x})\left\vert \Psi (\lambda )\right\rangle /\langle \Psi (\lambda )| \Psi (\lambda )\rangle\) can be evaluated by contracting iPEPS \(\left\vert \Psi (\lambda )\right\rangle\). We contract (the squared norm of) the iPEPS using the corner transfer matrix renormalization group (CTMRG) algorithm. As shown in Fig. 5d, we approximate the environment of the double tensors (shown in Fig. 5c) in a 2 × 2 unit cell using corner (rectangles) and edge tensors (squares) with a bond dimension χ. Since we impose the square lattice symmetry to the tensor A, we can contract the iPEPS using the symmetric CTMRG73. The bond dimensions of the corner and edge tensors grow to D2χ after absorbing the double tensors, and we should truncate the bond dimension back to χ. Figure 5e shows that using hermitian eigenvalue decomposition we obtain the isometries (triangles), which can be used to truncate the bond dimensions, and we just use eigenvectors corresponding to the χ largest eigenvalues (in absolute value) to construct the isometry. With the isometries, we can update corner and edge tensors, as shown in Fig. 5f.

In order to optimize the iPEPS, one has to provide the energy gradient ∂E/∂λ. The best way to calculate the energy gradient is using automatic differentiation (AD)24, which calculates the gradient through a backward propagation along the computational graph based on the chain rule in calculus. One problem of applying AD to calculate ∂E/∂λ is that the gradient can be infinite when eigenvalues are degenerate. Although one can add a small perturbation to lift the degeneracy, numerical instability can still happen with a small probability. When we get an infinity gradient, we can detach the isometries from the computation graph and get an approximate gradient; a trade-off between the stability and the accuracy. A possibly better solution is to use the approaches shown in ref. 74.

Given the energy expectation value and its gradient, we use the BFGS (Broyden-Fletcher-Goldfarb-Shanno) algorithm to minimize the energy expectation value. When the optimization is converged, we have an iPEPS \(\left\vert \Psi (\lambda )\right\rangle\) approximating a ground state of the toric code model. For this work, the iPEPS optimization was performed by PyTorch on the NVIDIA A100 80 GB GPU cards. We use the checkpoint function of PyTorch to reduce the huge memory cost of backward AD calculation. Moreover, when getting the isometries, we should use the hermitian eigenvalue decomposition rather than singular eigenvalue decomposition, because the former is about ten times faster than the latter on the GPU. It takes about two weeks (3 days) to calculate a single curve containing more than 100 data points with bond dimension D = 5 (D = 4) on a single A100 GPU card. Near critical points, CTMRG needs about 200 to 300 iterations to converge for a single energy evaluation.

iPEPS for topologically degenerate ground states

Here, we discuss two kinds of iPEPS representations of the toric code ground states. Understanding them is useful for initializing the iPEPS optimization and evaluating the FM string order parameter. When hx = hz = 0, the toric code Hamiltonian defined on a torus commutes with two Wilson loop operators and two ’t Hooft loop operators:

$${W}_{x}^{Z}=\prod _{e\in {L}_{x}}{Z}_{e},\quad {W}_{y}^{Z}=\prod _{e\in {L}_{y}}{Z}_{e},\quad {W}_{x}^{X}=\prod _{e\in {\hat{L}}_{x}}{X}_{e},\quad {W}_{y}^{X}=\prod _{e\in {\hat{L}}_{y}}{X}_{e},$$
(5)

where Lx (Ly) is a non-contractible loop along x(y) direction on the primal lattice, and \({\hat{L}}_{x}\) (\({\hat{L}}_{y}\)) is a non-contractible loop along x(y) direction on the dual lattice. They satisfy

$$\begin{array}{lll}\left[{W}_{x}^{Z},{W}_{y}^{Z}\right]\,=\,\left[{W}_{x}^{X},{W}_{y}^{X}\right]=\left[{W}_{x}^{Z},{W}_{x}^{X}\right]=\left[{W}_{y}^{Z},{W}_{y}^{X}\right]=0,\\ \left\{{W}_{x}^{Z},{W}_{y}^{X}\right\}\,=\left\{{W}_{x}^{X},{W}_{y}^{Z}\right\}=0.\end{array}$$
(6)

The four ground states can be labeled by eigenvalues of a pair of two commuting loop operators, i.e., common eigenstates of \({W}_{x}^{X}\) and \({W}_{y}^{X}\):

$$\left\vert {+}_{x}{+}_{y}\right\rangle ,\quad \left\vert {+}_{x}{-}_{y}\right\rangle ,\quad \left\vert {-}_{x}{+}_{y}\right\rangle ,\quad \left\vert {-}_{x}{-}_{y}\right\rangle ,$$
(7)

or common eigenstates of \({W}_{x}^{Z}\) and \({W}_{y}^{Z}\):

$$\left\vert {0}_{x}{0}_{y}\right\rangle ,\quad \left\vert {0}_{x}{1}_{y}\right\rangle ,\quad \left\vert {1}_{x}{0}_{y}\right\rangle ,\quad \left\vert {1}_{x}{1}_{y}\right\rangle .$$
(8)

In particular, the common eigenstates of \({W}_{x}^{Z}\) and \({W}_{x}^{X}\) (alternatively one can use \({W}_{y}^{Z}\) and \({W}_{y}^{X}\)) are minimally entangled states44:

$$\left| {0}_{x}{+}_{x}\right\rangle \equiv \left| {\mathbb{1}}\right\rangle ,\quad \left| {1}_{x}{+}_{x}\right\rangle \equiv \left| {\boldsymbol{m}}\right\rangle ,\quad \left| {0}_{x}{-}_{x}\right\rangle \equiv \left| {\boldsymbol{e}}\right\rangle ,\quad \left| {1}_{x}{-}_{x}\right\rangle \equiv \left| {\boldsymbol{f}}\right\rangle .$$
(9)

Let us focus on the ground states \(\left| {0}_{x}{0}_{y}\right\rangle\) and \(\left| {+}_{x}{+}_{y}\right\rangle\), which can be exactly expressed in terms of the so-called “single-line" and “double-line" iPEPS with the toroidal boundary condition75. They are included in the 2 × 2 unit cell iPEPS ansatz in Fig. 5a. As shown in Fig. 6a, by defining two of rank-3 tensors, we can obtain the A tensor for the iPEPS \(\left| {+}_{x}{+}_{y}\right\rangle\) shown in Fig. 6b or the iPEPS \(\left| {0}_{x}{0}_{y}\right\rangle\) shown in Fig. 6c. Because the two kinds of nonequivalent of iPEPS representations of the fixed point toric code states are included in the 2 × 2 unit cell iPEPS ansatz, we expect that it performs better than other kinds of iPEPS ansatz for the toric code model, especially along the self-dual line.

Fig. 6: The exact iPEPS at the fixed point hx = hz = 0 of the toric code phase.
figure 6

a The two rank-3 three tensors are defined to construct the tensor A. The tensor A for \(\left\vert {+}_{x}{+}_{y}\right\rangle\) and \(\left\vert {0}_{x}{0}_{y}\right\rangle\) has different virtual \({{\mathbb{Z}}}_{2}\) symmetry. iPEPS for b \(\left\vert {+}_{x}{+}_{y}\right\rangle\) and c \(\left\vert {0}_{x}{0}_{y}\right\rangle\) on a torus. The gray lines indicate the primal lattice.

Full size image

When evaluating the FM string order parameter on a torus geometry for a non-contractible loop L, which as we have seen is convenient for our tensor network methods, the denominator in Eq. (3) depends on the choice of the topologically degenerate ground states in the toric code phase. This can be easily understood at the fixed point of the toric code phase. When using the ground state \(\left| {+}_{x}{+}_{y}\right\rangle\) to evaluate the FM string order parameter OZ, we find that it is not well-defined because \(\left| {+}_{x}{+}_{y}\right\rangle\) spontaneously breaks the \({W}_{x}^{Z}\) symmetry, i.e., \({W}_{x}^{Z}\left| {+}_{x}{+}_{y}\right\rangle =\left| {+}_{x}{-}_{y}\right\rangle\), such that the denominator \(\left\langle {+}_{x}{+}_{y}\right| {W}_{x}^{Z}\left| {+}_{x}{+}_{y}\right\rangle =0\). This problem can be circumvented by considering either the ground state \(\left| {0}_{x}{0}_{y}\right\rangle\) or the minimally entangled state of the trivial topological sector \(\left\vert 1\right\rangle\), such that the denominator of the FM string order parameter OZ is 1. However, if we want to correctly evaluated FM string order parameters OX and OZ simultaneously, we have to consider the minimally entangled state \(\left| {\mathbb{1}}\right\rangle\). The FM string order parameters OX and OZ evaluated using the minimally entangled state in the trivial topological sector which mimics the one defined on a contractible loop.

Away from the fixed point of the toric code model, the exact Wilson and ’t Hooft loop operators in Eq. (5) are not symmetries the Hamiltonian anymore. Since we mainly focus on the FM string order parameter OZ, we just need to find the ground state \(\left| {0}_{x}{0}_{y}\right\rangle\), which is the simultaneous eigenstate of emergent non-contractible Wilson loop operators along x and y directions. Because the optimized iPEPS usually converge to \(\left| {0}_{x}{0}_{y}\right\rangle\) or \(\left| {+}_{x}{+}_{y}\right\rangle\), we can control which ground state it converges to by initializing the iPEPS tensor as A + ϵR, where A is given in Fig. 6a and R is a random tensor satisfying the required lattice symmetry, and ϵ is a small number. We emphasis that we do not imposing the virtual \({{\mathbb{Z}}}_{2}\) symmetry shown in Fig. 5b to the tensor during the optimization.

Evaluation of the FM string order parameter using iPEPS

With the optimized iPEPS tensor, we can evaluate the FM string order parameter efficiently in the limit of an infinitely long string using transfer matrices of iPEPS. To this end, we construct from the iPEPS of the ground state two transfer matrices; one is the usual transfer matrix \({\mathbb{T}}\) and the other is \({{\mathbb{T}}}_{Z}\) containing Z operators; see Fig. 7a for graphical notations. The numerator of the FM string order parameter consists of r transfer matrices \({\mathbb{T}}\) followed by r transfer matrices \({{\mathbb{T}}}_{Z}\) containing a Z string and the denominator consists of 2r transfer matrices \({{\mathbb{T}}}_{Z}\) containing a Z loop:

$$\begin{array}{rcl}{O}_{Z}&=&\mathop{\lim }\limits_{r\to \infty }{\left[\frac{{\rm{Tr}}({{\mathbb{T}}}^{r}{{\mathbb{T}}}_{Z}^{r})/{\rm{Tr}}({{\mathbb{T}}}^{2r})}{\sqrt{{\rm{Tr}}({{\mathbb{T}}}_{Z}^{2r})/{\rm{Tr}}({{\mathbb{T}}}^{2r})}}\right]}^{1/2}\\ &=&{\left[\mathop{\lim }\limits_{r\to \infty }\frac{{\rm{Tr}}({{\mathcal{T}}}^{r}{{\mathcal{T}}}_{Z}^{r})/{\rm{Tr}}({{\mathcal{T}}}^{2r})}{\sqrt{{\rm{Tr}}({{\mathcal{T}}}_{Z}^{2r})/{\rm{Tr}}({{\mathcal{T}}}^{2r})}}\right]}^{\frac{1}{2}}\end{array}$$
(10)
$$=\mathop{\lim }\limits_{r\to \infty }\frac{{t}^{\frac{r}{2}}{t}_{Z}^{\frac{r}{2}}/{t}^{r}}{{t}_{Z}^{\frac{r}{2}}/{t}^{\frac{r}{2}}}| \langle V| {V}_{Z}\rangle | =| \langle V| {V}_{Z}\rangle | ,$$
(11)

where we compress the transfer matrices \({\mathbb{T}}\) and \({{\mathbb{T}}}_{Z}\) to the transfer matrices \({\mathcal{T}}\) and \({{\mathcal{T}}}_{Z}\) with a dimension D2χ2 using the edge tensors from the CTMRG, see Fig. 7b, and t (tZ), V (VZ) are dominant eigenvalue and eigenvector of \({\mathcal{T}}\) (\({{\mathcal{T}}}_{Z}\)). Here, we assume that the dominant eigenvectors are non-degenerate and discuss the degenerate case later. In the limit of r → ∞, the action of the transfer matrices is set by the dominating eigenvalue and eigenvector, \({{\mathcal{T}}}_{Z}^{\infty }={t}_{Z}^{\infty }\left\vert {V}_{Z}\right\rangle \left\langle {V}_{Z}\right\vert\) and \({{\mathcal{T}}}^{\infty }={t}^{\infty }\left\vert V\right\rangle \left\langle V\right\vert\), as shown in Fig. 7c, so Eq. (10) can be simplified to Eq. (11). Moreover, the perimeter law coefficient can be obtained from the dominant eigenvalues of the transfer matrices: \({\alpha }_{Z}=-\log ({t}_{Z}/t)\).

Fig. 7: Evaluating the FM string order parameter using iPEPS.
figure 7

a Using the double tensors in Fig. 5c, two iPEPS transfer matrices are defined to evaluate the FM string order parameter. b Compressing the iPEPS transfer matrices using the edge tensors from the CTMRG. c The infinite power transfer matrices are given by their fixed points.

Full size image

Next, we consider degenerate fixed points of \({\mathcal{T}}\) and \({{\mathcal{T}}}_{Z}\). This could happen when a ground state in the toric code phase is chosen as the minimally entangled state. In the limit \(r\to \infty ,{{\mathcal{T}}}^{r}\) (\({{\mathcal{T}}}_{Z}^{r}\)) can be expressed as:

$${{\mathcal{T}}}^{\infty }=\mathop{\sum }\limits_{\alpha =1}^{d}{t}^{\infty }\left\vert {V}_{\alpha }\right\rangle \left\langle {V}_{\alpha }\right\vert ,\quad {{\mathcal{T}}}_{Z}^{\infty }=\mathop{\sum }\limits_{{\alpha }_{Z}=1}^{{d}_{Z}}{t}_{Z}^{\infty }\left\vert {V}_{Z,{\alpha }_{Z}}\right\rangle \left\langle {V}_{Z,{\alpha }_{Z}}\right\vert ,$$
(12)

where d (dZ) denotes the number of the dominant eigenvectors of \({\mathcal{T}}\) and \({{\mathcal{T}}}_{Z}\), and α (αZ) specifies the degenerate dominant vectors. Substituting Eq. (12) in Eq. (10), we can calculate the FM string order parameter even when the transfer matrix fixed points are degenerate:

$${O}_{Z}={\left[\frac{1}{\sqrt{d{d}_{Z}}}{\rm{Tr}}\left(\mathop{\sum }\limits_{\alpha = 1,\beta = 1}^{d,{d}_{Z}}\left| {V}_{1,\alpha }\right\rangle \left\langle {V}_{1,\alpha }| {V}_{Z,\beta }\right\rangle \left\langle {V}_{Z,\beta }\right| \right)\right]}^{\frac{1}{2}}.$$
(13)

Note Added

While finalizing the manuscript we became aware of related work on the stability of the FM order parameter76.