Finite-temperature critical behaviors in 2D long-range quantum Heisenberg model

The Mermin-Wagner theorem states that spontaneous continuous symmetry breaking is prohibited in systems with short-range interactions at spatial dimension $D\le 2$. For long-range interactions with a power-law form ($1/r^{\alpha}$), the theorem further forbids ferromagnetic or antiferromagnetic order at finite temperature when $\alpha\ge 2D$. However, the situation for $\alpha \in (2,4)$ at $D=2$ is not covered by the theorem. To address this, we conduct large-scale quantum Monte Carlo simulations and field theoretical analysis. Our findings show spontaneous breaking of $SU(2)$ symmetry in the ferromagnetic Heisenberg model with $1/r^{\alpha}$-form long-range interactions at $D=2$. We determine critical exponents through finite-size analysis for $\alpha<3$ (above the upper critical dimension with Gaussian fixed point) and $3\le\alpha<4$ (below the upper critical dimension with non-Gaussian fixed point). These results reveal new critical behaviors in 2D long-range Heisenberg models, encouraging further experimental studies of quantum materials with long-range interactions beyond the Mermin-Wagner theorem's scope.

The Mermin-Wagner theorem states that spontaneous continuous symmetry breaking is prohibited in systems with short-range interactions at spatial dimension D ≤ 2. For long-range interactions with a power-law form (1/r α ), the theorem further forbids ferromagnetic or antiferromagnetic order at finite temperature when α ≥ 2D.However, the situation for α ∈ (2, 4) at D = 2 is not covered by the theorem.To address this, we conduct large-scale quantum Monte Carlo simulations and field theoretical analysis.Our findings show spontaneous breaking of SU (2) symmetry in the ferromagnetic Heisenberg model with 1/r α -form long-range interactions at D = 2.We determine critical exponents through finite-size analysis for α < 3 (above the upper critical dimension with Gaussian fixed point) and 3 ≤ α < 4 (below the upper critical dimension with non-Gaussian fixed point).These results reveal new critical behaviors in 2D long-range Heisenberg models, encouraging further experimental studies of quantum materials with long-range interactions beyond the Mermin-Wagner theorem's scope.
In recent years, the importance of the studies on longrange(LR) lattice models have been gradually noticed, due to the fact that they exhibit intrinsically different properties from their short-ranged(SR) counterparts.For example, LR Heisenberg models at spatial dimension D = 2 acquires anomalous magnon dispersion different from the linear and quadratic spin-waves in the SR antiferromagnetic and ferromagnetic models [1,2].In addition, the violation of Mermin-Wagner theorem and unconventional critical properties in LR systems also attracted much attention in investigations of both quantum spin models and interacting fermionic models .
These phenomena also have immediate experimental relevance.Due to the fast development in the Rydberg atom arrays [24][25][26][27][28], the magic angle twisted bilayer Graphene and other 2D quantum moiré materials  and the programmable quantum simulators [67,68] such as quantum gases coupled to optical cavities [69].LR interactions in the forms of van der Waals, dipole-dipole and Coulomb have given rise to a plethora of correlated topological and quantum phases of matter beyond the semi-classical or mean-field type descriptions, and new theoretical paradigm that could cope with these fast emergent experimental facts are critically called for.
One particularly interesting direction is to explore the critical properties of phase transitions with continuous symmetry breaking, outside the realm of the established Mermin-Wagner theorem.For 1D LR antiferromagnetic Heisenberg chain [23] and Heisenberg ladders [14] with 1/r α -form LR interactions, the phase diagram as well as the critical exponents have been addressed and it has been found that there is a upper critical value α c above which there is no phase transitions for these systems.Below α c , the transition exists and the critical exponents are dependent on α, as identified by both field theory analysis and numerical evidence.However, for 2D LR Heisenberg models with finite-temperature transitions, it was only known that, for D = 2 Heisenberg model with ferromagnetic LR interaction 1/r α , a finite-temperature ferromagnetic phase will not exist when α ≥ 4 which has been proved analytically in Ref. [70], and for α ≤ 2 the system is gapped due to the generalized Higgs mechanism [1,2] and the finite-temperature ferromagnetic order should be allowed.However, the situation in α ∈ (2, 4) is not well understood.Although there are classical field theory predictions and renormalization group analysis on this issue [3,4,7], which state there is a Gaussian fixpoint for 2 < α < 3 and a non-Gaussian fixed-point for 3 ≤ α < 4, a thorough numerical treatment on the 2D quantum Heisenberg model has not been performed to date.Such unbiased numerical analysis of this model is crucial not only because the field-theory scenario needs to be impartially examined on the realistic lattice models, but also due to the fact that the Heisenberg model is one of the most central toy models in condensed matter and statistic physics and a complete clarification of the critical properties of this model will serve as the cornerstone of further studies on LR quantum many body systems.
Here we bridge these gaps by large-scale QMC simulations and field theory analysis.We find clear evidence of the breakdown of the Mermin-Wagner theorem with finite-temperature phase transitions in α ∈ (2, 4), as shown in Fig. 1.By performing the state-of-the-art finite-size scaling analysis, as illustrated in Fig. 2 As the temperature is reduced, the system undergoes a continuous phase transition from paramagnetic phase to ferromagnetic phase in entire region of α ∈ (2, 4).The black dots are the critical points determined from QMC simulations, as exemplified in Fig. 2 and Fig. 3.The standard error of the mean (SEM) is used when estimating the errors of the physical quantities.
obtain the accurate critical exponents of the phase transition as a function of α as shown in Fig. 3, and demonstrate these results nicely satisfy the field-theory predictions both for α < 3 where the system is above the upper critical dimension with Gaussian fixed point and for 3 ≤ α < 4 where the system is below the upper critical dimension with non-Gaussian fixed point.Our results explicitly show the critical behaviors for α ∈ (2, 4) in LR Heisenberg model at D = 2 and will intrigue further theoretical and experimental physics and even mathematics studies of systems with LR interactions beyond the realm of the Mermin-Wagner theorem [3][4][5][6][7][8]71].

Results.
Model.The Hamiltonian of the LR ferromagnetic Heisenberg model is where denotes the LR coupling and r ij is the nearest distance between site i and site j under the periodic boundary condition.In order to alleviate the strong finite-size effects in systems with LR interactions arising from the cut-off of LR interactions under the periodic boundary condition, we replace J ij with the Ewaldcorrected coupling Jij [19,72] which takes the form of This modified coupling parameter Jij counts all the possible distances between two sites under the periodic boundary condition ,so that the effect of cutting off the tail of LR interactions is minimized, and this trick has been shown to be very useful in the simulation of many LR systems [14,18,19,72].For 2D there is no closed form for Eq. ( 2), so we truncate the summation at |m| max , |n| max = 1000 for α < 3 which is large enough to have the well-converged finite-size scaling behavior, as shown in Fig. 2. For α ≥ 3 the finite-size effects are mainly from crossovers to SR case, and we find the original coupling J ij is fine to obtain converged results.When α ≥ 2D the system reduces to the SR case where there is no spontaneously continuous symmetry breaking phase at finite-temperature.When α ≤ D, the Hamiltonian is no longer extensive and there is no well-defined thermodynamic limit.Between α ∈ (2, 4) we carry out the QMC simulations [73][74][75] up to the linear system size of L = 256, as shown in Fig. 2, to determine the precise phase boundary as well as the critical exponents ν, β and η.Note that because of strong finite-size effects, we only compute the region of α ∈ [2.3, 3.7] where our QMC simulations can obtain well-converged results.The origins of finite-size effects as α approaches the two boundaries, α = 2 and α = 4, exhibit inherent distinctions.When α → 2, the finite-size effect arises from the escalating intensity of LR (long-range) interactions, which fundamentally reduces the efficiency of the Ewald-corrected scheme.Conversely, as α → 4, the system approaches the regime where finite-temperature phase transitions do not exist.Consequently, near this boundary, the convergence of data points becomes exceedingly slow to be overcome.The results are shown in Figs. 1 and 3 and will be discussed in the critical exponents section.The QMC implementation is explained in the Supplementary Note 1.
Note that when α ≤ D, the Hamiltonian defined in Eq. ( 1) can actually be Kac-normalized [10,76] to be extensive with the addition of a factor N −1 i<j Jij to the Hamiltonian.Although this is not the focus of our paper, we examine the Kac-normalized Hamiltonian and the results are shown in the Supplementary Note 2.
Critical exponents.Fig. 2 shows our results at α = 2.5.We first use the crossing points of the Binder ratios to locate the critical temperature T c .The crossing points of U (T, L) with U (T, 2L) are denoted as T * (L), and through fitting to Eq. ( 9) the precise value of T c can be obtained.We then use the value of T c to perform data collapse according to Eq. (10) and Eq. ( 13) separately for 3 ≤ α < 4 and α < 3, to obtain the critical exponents ν ′ and β.To obtain the anomalous dimension η Q , we measure the correlation function G(L/2) at the obtained critical temperature T c and obtain the anomalous dimension separately by fitting to Eq. ( 11) for 3 ≤ α < 4 and Eq. ( 14) for α < 3.
According to the conventions defined in Eq. ( 16) and field theory results of the mean-field critical exponents in Eq. ( 7), we can extract the expression for the three critical exponents in the Gaussian region which are ν ′ = 1, β = 1 2 and η Q = 1.Outside the Gaussian region, we have η = 4 − α and γ defined in Eq. ( 8), and the value of β and ν can be obtained via solving the scaling relations between the critical exponents with ν = γ 2−η and β = γη 2(2−η) .The critical exponents we have obtained are shown in Fig. 3.We find that within the region we simulated, our QMC-obtained critical exponents ν ′ (α) , β(α), and η Q (α) match nicely with the prediction of both LR Gaussian theory (for α < 3) and the two-loop perturbative RG (for 3 ≤ α < 4) , although there is a sign of deviating from two-loop RG predictions when α approaches 4. The possible deviation might be explained by the increasing finite-size effects near the boundary or the inefficiency of two-loop perturbative RG predictions when α is away from α = 3.The results can be further improved by either considering higher-order RG corrections or by pushing the QMC simulations to larger system sizes.No-tably, the predicted form of anomalous dimension η receive no corrections at any α ∈ (2, 4) [3] and our results confirm this argument with η matching with η = 4 − α well in the whole region.

Discussions.
Our investigation reveals a finite-temperature phase transition point in the 2D LR Heisenberg model, occurring for values of α within the range of α ∈ (2, 4), which separates the ferromagnetic phase from the paramagnetic phase.We observe that the phase transition point exhibits distinct behaviors: a Gaussian fixed point characterizes the transition for α ≤ 3, while a non-Gaussian fixed point emerges for 3 < α < 4. Similar phenomena have been observed in various LR systems [6-10, 14, 18, 19, 23].However, it is important to note that LR Ising-like systems differ intrinsically from LR Heisenberg-like systems.The former does not adhere to the Mermin-Wagner theorem, guaranteeing a finite-temperature transition for all α > 0, while the latter exhibits an upper critical value α c beyond which the Mermin-Wagner theorem precludes the existence of phase transitions.In conclusion, our results clearly point out the LR quantum many-body system exhibit unconventional critical properties beyond the realm of the Mermin-Wagner theorem, which are also worthwhile to pursue in future experimentalrealizations ,such as the quantum simulators.

Field theory analysis.
We review here the field theory description of the model at the thermodynamic limit dating back to Ref. [3].The action can be written as to match the lattice model, we need α = d + σ.Under the scaling symmetry the kinetic term remains unchanged when ∆ ϕ = D−σ 2 .The coupling constant of ϕ 4 interaction, on the other hand, scale as λ → s 2α−3D λ. ( When α < 3D 2 , the coupling constant decays at larger length scale, which means the λϕ 4 term is an irrelevant operator.The Gaussian fixed point at λ = 0 is a stable fixed point.Notice when λ = 0, the action is in a purely quadratic form, hence named "Gaussian" fixed point.This was established mathematically in Ref. [6].When α > 3D 2 , the λϕ 4 term becomes relevant, which triggers a renormalization group towards a different non-Gaussian fixed point [3].One can perform standard renormalization technique to calculate the scaling dimension of various operators, by evaluating Feynman diagrams with non-conventional propagators.Such a calculation was first performed in [3].Since the kinetic term in Eq. ( 3) is no-local, which can not receive corrections from any local counter terms, the scaling dimension of ϕ will not be renormalized (This can be easily seen by analyzing the the Callan-Symanzik equation for the two point function ⟨ϕ(x)ϕ(y)⟩, see for example, Ref. [77]).Equivalently, we have η = 2∆ ϕ − D + 2. Our numerical result clearly confirms such a theoretical prediction.For a fixed σ in Eq. ( 3), we can define the upper critical dimension as the space-time dimension at which the ϕ 4 term is marginal.The ∆ ϕ 4 = 4∆ ϕ = D uc gives us We now focus on the D = 2 case.When α < 3, the critical behavior is controlled by the λ = 0 Gaussian fixed point.The critical behavior is similar to the usual Ising model at D > 4, due to the effect of dangerously irrelevant operators [78], the critical exponents are given by For example, the β = 1/2 exponent can be seen from the following argument.Deform the action (3) by a mass term dx D tϕ(x) 2 with negative t and minimize the potential, we get ⟨ϕ⟩ ∝ (−t/λ) β , with β = 1/2.The other exponents can be calculated by similar mean field theory analysis.The critical exponent η controls the two point function ⟨ϕ(x)ϕ(y)⟩ only at the strict thermodynamic limit.At finite sizes, the power law behaviour will be modified to (14), which follows from analysing the effect of dangerously irrelevant operators carefully [79].
When α > 3, on the other hand, the second term in Eq. (3) becomes relevant, and renormalization group flows towards a different non-Gaussian fixed point [3].The critical exponent η will remain at its mean field theory value [3] as in Eq. (7).The other exponents, on the other hand receives correction at O (α − 3) 2 .The twoloop perturbation results for γ is where ψ(z) is the logarithmic derivative of the gamma function.The other critical exponents can be obtained by scaling relations between them.
When α > 4, the long-range model becomes equivalent to short-range models, due to the Mermin-Wagner theorem [80][81][82], the system will be gapped at finitetemperature.In the field-theory language, the value of α at which such a long-range to short-range crossover happens when the scaling dimension of ϕ equals to the scaling dimension of the short range model.In two dimensions, this gives α = 4 [3,4].
Finite-size scaling analysis.To identify the phase transitions and obtain the critical exponents, we compute the square magnetization ⟨m 2 ⟩, the correlation function G(r), and the Binder ratio in the QMC simulation.The crossing point of U (T, L) with U (T, 2L) is denoted as T * (L) and it is expected to converge to the thermodynamic limit critical temperature T c following the scaling relation: Given the values of T * (L) with sufficiently small errors and large enough system sizes L , the critical point T c can be precisely located as shown in Fig. 2. To obtain the critical exponents ν, β and η, when D ≤ D uc , the standard finite-size scaling behavior (FSS) [83,84] allows us to perform a data collapse near the critical points with the relation The anomalous dimension can also be obtained by fitting to the correlation function at the critical point T c However, when D > D uc , which is our case when α < 3, the system enters the mean-field region where the hyperscaling relation breaks down, famously due to the effect of dangerously irrelevant operator [79,85,86].The scaling of the correlation length in this region shall follow the relation ξ L ∼ L Duc D instead of ξ L ∼ L [10,18,19,79,[85][86][87], and this leads to the modification of hyperscaling relation with where ν ′ = Duc D ν and α H is the critical exponent associated with the specific heat.For our system Eq.( 1), the upper critical dimension is D uc = 2(α − D), which we will explain later in the field theory analysis section.Accordingly, Eq. ( 10) also needs to be modified and the correct relation for data collapse in mean field region is [10,18,19,79,85] The scaling of correlation function for α < 3 is also modified with where By fitting to Eq. ( 14), the modified anomalous dimension η Q as well as η can be obtained.
To unify the conventions, we define and Then ν ′ , β and η Q will be obtained with the same scaling functions for both α < 3 and 3 ≤ α < 4.
Rong conducted the field theory analysis .All authors contributed to the analysis of the results and the preparation and revision of the draft.
= − a(i<j) H 1,a(ij) + H 2,a(ij) and a(ij) is the bond index.Take the eigenstates of σ z as basis, the non-zero matrix elements in Eq. ( 18) are In the SSE QMC simulation, the loop update scheme is maintained the same with the ferromagnetic Heisenberg model with nearest-neighbor couplings.However, to efficiently carry out the diagonal update scheme, we choose the candidate bonds for inserting diagonal operators with an importance sampling procedure with P choose ∝ J ij .
The diagonal update scheme is thus 1.If a diagonal operator (H 1,a ) is visited, remove it with probability 2. If an identity operator (H 0,a ) is visited, insert a diagonal operator according to: • First choose a candidate bond a to make the insertion with probability • Then accept the insertion of a diagonal operator at this position with probability To generate a set of random bond index according to the probability defined in Eq. ( 21), we use the naive Walker's method [88] with complexity of O(N 2 ).Despite there is optimization of this method which reduces the complexity to O(N ) [72], we find for the system size we simulate the original method is sufficient and easy to implement.The above procedure ensures that diagonal operators with higher matrix elements have higher probability to be inserted and compared with randomly choosing candidate bonds this strategy certainly has better efficiency.The system also undergoes a continuous phase transition from FM to PM as temperature is increased.The black dots are the critical points determined from QMC simulations.The standard error of the mean (SEM) is used when estimating the errors of the physical quantities.

𝑈(𝑇
phase transition at α ∈ (2, 4).However, when α ≤ 2 the system is no longer extensive and it is still unclear whether such a system still hold a phase transition point.To make the system extensive, a Kac normalization factor can be added to the Hamiltonian and the normalized Hamiltonian is In this case, the energy density is ⟨↑ • • • ↑↑ |H/N | ↑↑ • • • ↑⟩ = −(N − 1)/N which will be a constant for any value of α.For the Hamiltonian defined in Eq. ( 24), we perform the QMC simulations and find that there is a continuous phase transition for α = 1.8 and the critical exponents also satisfies the prediction of mean-field theory, as shown in Supplementary Figure 4.The similar phenomena has also been observed in 1D LR quantum Ising models [10], where at α = 0.05 which is below the system dimension, the critical exponents are still consistent with mean-field predictions.
In addition, we also examine this system at other values of α and we find that the phase transition point T c remains the same for all the α ≤ 2 we consider, as indicated in Supplementary Figure 5.There is an intuitive understanding for this finding: the Kac normalization [10,76] makes the energy scale to be the same for all α ≤ 2 which somehow suppresses the effect of different decaying exponent α and the phase transition temperature T c is thus scaled to be the same for all α.This point should be further examined by more robust analysis.

FIG. 1 .
FIG. 1. Phase diagram of the 2D LR ferromagneticHeisenberg model.As the temperature is reduced, the system undergoes a continuous phase transition from paramagnetic phase to ferromagnetic phase in entire region of α ∈ (2, 4).The black dots are the critical points determined from QMC simulations, as exemplified in Fig.2and Fig.3.The standard error of the mean (SEM) is used when estimating the errors of the physical quantities.

FIG. 2 .
FIG. 2. The determination of the critical point and exponents at α = 2.5.(a) Binder ratio U (T, L) versus temperature T for different system sizes.(b) Crossing points of Binder ratios T * (L) versus 1/L.The solid line represents a fitting of the data points with Eq. (9).The fitted curve is T * (L) = −2.935L−1.491 + 3.5776.(c) Data collapse of the order parameter ⟨m 2 ⟩ near the critical point Tc.Notice here we replace the correlation length exponent ν with ν ′ as in Eq. (13).(d) ln[G(L/2)] versus ln(L) for different system sizes L = 16, 24, 36, 54, 80, 120, 180.The data is fitted with a straight line as in Eq. (14) and the fitted result is ln[G(L/2)] = −0.999(1)ln(L).The errors of ln[G(L/2)] are smaller than the symbol sizes and SEM is used when estimating the errors of the physical quantities.

FIG. 3 .
FIG. 3. Critical exponents ν ′ , β and ηQ in the region of α ∈ [2.3, 3.7] obtained from data collapse and from fitting to the correlation function G(r).The black and blue solid lines in (a), (b) and (c) are the predictions of LR Gaussian theory (α < 3) and two-loop perturbative RG predictions (3 ≤ α < 4) for Gaussian and interacting (non-Gaussian) fixed points.[3, 4, 6].SEM is used when estimating the errors of the physical quantities.

FIG. 4 .FIG. 5 .
FIG. 4. The determination of the critical point and exponents at α = 1.8.(a) Binder ratio U (T, L) versus temperature T for different system sizes.(b) Data collapse of the order parameter ⟨m 2 ⟩ near the critical point Tc.The obtained results are ν ′ = 1.00(5) and β = 0.505(8).The standard error of the mean (SEM) is used when estimating the errors of the physical quantities.