Diagnosing weakly first-order phase transitions by coupling to order parameters

The hunt for exotic quantum phase transitions described by emergent fractionalized degrees of freedom coupled to gauge fields requires a precise determination of the fixed point structure from the field theoretical side, and an extreme sensitivity to weak first-order transitions from the numerical side. Addressing the latter, we revive the classic definition of the order parameter in the limit of a vanishing external field at the transition. We demonstrate that this widely understood, yet so far unused approach provides a diagnostic test for first-order versus continuous behavior that is distinctly more sensitive than current methods. We first apply it to the family of $Q$-state Potts models, where the nature of the transition is continuous for $Q\leq4$ and turns (weakly) first order for $Q>4$, using an infinite system matrix product state implementation. We then employ this new approach to address the unsettled question of deconfined quantum criticality in the $S=1/2$ N\'eel to valence bond solid transition in two dimensions, focusing on the square lattice $J$-$Q$ model. Our quantum Monte Carlo simulations reveal that both order parameters remain finite at the transition, directly confirming a first-order scenario with wide reaching implications in condensed matter and quantum field theory.


Introduction
The theory of phase transitions is an integral component in the understanding of many body phenomena, playing a significant role in fields ranging from statistical mechanics to condensed matter and high energy physics.The most recent advancements in this area have focused on phase transitions beyond the Landau-Ginzburg-Wilson (LGW) paradigm, where models have been proposed and studied extensively both analytically and numerically.In these scenarios the field theory describing the transition is not given in terms of the order parameters, as in the LGW paradigm, but instead is formulated in terms of fractionalized degrees of freedom.This leads to the possibility of a generic, continuous phase transition between two ordered phases whose symmetry groups and broken symmetries are not mutually compatible.
Prominent candidate examples for such exotic transitions are the 'deconfined quantum critical points' (DQCP) between Néel ordered antiferromagnetic and valence bond solid phases in quantum spin systems [1,2].The model believed to epitomize this scenario is the so called J-Q model [3] (to be discussed in detail below), which to date has been the most well studied in this context.Indeed, an impressive body of numerical work has been devoted to the analysis of the nature and the critical exponents of this and related models [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17].While most of the numerical data for the S = 1/2 J-Q model has been interpreted as evidence for a continuous quantum phase transition, albeit with significant corrections to scaling, some authors have however interpreted their data as evidence for a weakly first order transition.The current consensus opinion in the community is that the true nature of the J-Q transition remains to be settled definitively.
On the analytical side, several connections between the original (NCCP 1 ) DQCP theory [1,2] and other field theories of current interest such as the Abelian Higgs model (scalar QED 3 ) [18], the SO (5) nonlinear sigma model with a Wess-Zumino term [19] or fermionic QED 3 [20] have been put forward [21].It was first believed that the NCCP 1 theory is continuous, but Refs.[22][23][24][25] put forward and discussed scenarios of colliding fixed points in theory space, where the annihilation of two real fixed points provides a possible mechanism for pseudo-universal weakly first order behavior through the appearance of complex fixed points.In general, it is an important open question to determine the critical number N c of (bosonic or fermionic) matter field flavors coupled to gauge fields, separating a regime of conformal field theories in the infrared (i.e.continuous transition behaviour) from a regime of pseudo-critical (weakly) first order regime in the infrared.This scenario also applies to the classical Q-state Potts model in two dimensions where the conformal window resides below Q ≤ 4, resulting in (in)famously weak first order transitions near Q > 4 [26][27][28].Additionally, similar ideas might also apply regarding the conformal window of non-abelian gauge theories [23], with possible implications for the hierarchy problem in particle physics.
It is therefore of great interest to refine numerical simulations so that conformal windows can be clearly resolved, which means developing highly sensitive methods for detecting weak first order transitions.However, pseudo-critical behavior implies a finite-but huge-correlation length at the transition.As a result, a conclusive diagnosis seemingly becomes impossible since the general wisdom is that system sizes with a linear extent at least as large as the correlation length are required to resolve the first order nature of the phase transition [29].
In this work we revive the textbook definition of the order parameter from the classical theory of phase transitions, which has all but fallen out of use in the numerical community.We find that applying these ideas in a modern context has the power to resolve weakly first order transitions with exquisite detail.After introducing the basic idea behind our approach in Sec. 2, we demonstrate the power of our method in the Hamiltonian formulation of the well-understood Q-state Potts model in Sec. 3 using an infinite matrix product state (iMPS) implementation, where we find distinct signatures for first order behavior manifesting at correlation lengths of a few lattice spacings, despite the correlation length at the transition itself reaching on the order of a thousand lattice spacings for Q = 5.We then move to two dimensional quantum critical phenomena using a Quantum Monte Carlo implementation, first in Sec. 4, where we show that this method corroborates the established continuous nature of the O(3) Wilson-Fisher CFT quantum critical behavior of a family of explicitly dimerized S = 1/2 quantum magnets.Finally in Sec. 5, we address the Néel to valence bond solid phase transition in the square (and rectangular) lattice S = 1/2 J-Q model and provide direct evidence for a first order scenario, thus resolving a long standing debate in the field, with implications in many directions both in condensed matter and in quantum field theory.

Outline of the Method
In an ordinary symmetry breaking phase transition there are two complementary conceptual approaches to track the order parameter O as a function of a control parameter g (temperature T or a Hamiltonian parameter).We assume that g > g c is the symmetry broken phase, g < g c the paramagnetic phase with g c the transition point.
• In a symmetry preserving setup one measures the square of the order parameter 〈|O| 2 〉 g or functions thereof (such as Binder cumulants or order parameter susceptibilities) for finite systems and then extrapolates to infinite system size using finite-size scaling techniques.
• On the other hand, it has long been known that one can directly measure the order parameter by coupling it to a uniform external field via adding a term λ d d x O(x) to the Hamiltonian (d denotes the space dimension).One then extrapolates 〈O〉 g,λ to infinite size at fixed coupling λ, then takes the limit as λ → 0 + , yielding the order parameter 〈O〉 g in the thermodynamic limit.
In the symmetric setup if the transition occurring at g c is continuous then we expect the standard power law behaviour 〈O〉 g ∼ (g − g c ) β , g → g + c , with β an appropriate critical exponent, while in an external field at the critical point g c : with 1/δ a critical exponent which is controlled by the universality class of the transition at g c1 .
If however the transition at g c is first order then there is a coexistence of the paramagnetic and the symmetry broken phase at g c .The applied field λ then prefers the symmetry broken phase, leading to a finite 〈O〉 g c ≡ m c , the (unique) value of the order parameter discontinuity at the transition.
The central quantity of interest in our work is the following logarithmic derivative at the critical coupling g c : We also refer to this quantity as the "running exponent 1/δ", since we will study this quantity as a function of λ.According to the discussion above we expect [1/δ](λ) to approach 1/δ for λ → 0 + in the continuous case, while the finite value value m c of the order parameter at g c in the first order case drives the logarithmic derivative to zero.
While not of central interest for the present work, one can also track the behavior of [1/δ](λ) for other values of g.In the symmetry broken regime g > g c we expect the running exponent to scale to zero, as in the first order case at g c .The paramagnetic phase g < g c requires some more care.For a unique, gapped ground state in the paramagnetic phase we expect a standard linear response regime, resulting in a running exponent [1/δ](λ) = 1 for small fields λ.This is actually also the expected behaviour for all g in a generic finite size system at very small λ.We will witness this phenomenon for the finite size quantum Monte Carlo simulations presented in Secs. 4 and 5.
While in the early days of Monte Carlo investigations of phase transitions approaches with and without a coupling to the order parameter were pursued [30], the tediousness of the double limit with an external field put this method at a disadvantage compared to symmetric setups that could locate critical points and measure exponents with fewer simulations.Subsequently with the development of powerful cluster algorithms that are tailor made for symmetric models [31,32], the order parameter coupling approach has seemingly fallen into complete disuse in modern numerical simulations of phase transitions.A notable exception is the pinning field method used in Ref. [33] to resolve the nature of the semimetal to Mott insulator transition in the honeycomb Hubbard model.In that work the applied field was confined to a single site, whereas we use a spatially extended coupling to the order parameter in our work.
In our present work we demonstrate that the order parameter coupling approach is a powerful tool for diagnosing weak first order transitions, performing well beyond the abilities of the currently used symmetric approaches.For our purposes we find its apparent drawbacks to be inconsequential in practical simulations, allowing it to be seamlessly integrated into state of the art numerical algorithms, here using infinite matrix product state (iMPS) and finite size quantum Monte Carlo (QMC) algorithms.In fact the presence of an order parameter coupling allows us to devise statistically exact QMC estimators for the running exponents as a function of the external field, eliminating finite-difference errors and facilitating the diagnosis of first order transitions.In the next section we demonstrate for the Q-state Potts model family that the difference in behavior of Eq. ( 1) between continuous and weakly first order transitions is surprisingly stark, and occurs at rather large values of the order parameter coupling λ and correspondingly short correlation lengths.

The Q-state Potts model
We start discussing our method by an application on a challenging problem, the Q-state Potts model in the 1+1D Hamiltonian formulation [28,34,35], which basically corresponds to a spatially anisotropic version of the widely known 2D classical Potts model.We consider the Hamiltonian already fine tuned to the exact value of the quantum phase transition between the ordered and   1), of the first row, highlighting the convergence towards the expected exponents in the left panel, and the existence of distinct maximum in the weakly first order cases in the right panel.In the bottom row we plot the extracted MPS correlation length, showing that the inflection point and the corresponding maximum in the running exponent for the weakly first order instances occur at correlation lengths of only a few lattice spacings.
paramagnetic phase (i.e.working at g c ), and add a symmetry breaking field λ favoring one of the Q ferromagnetically ordered states, here chosen as q = 0: where i runs over the sites of the chain and q over the Q distinct local states q ∈ {0, ..., Q − 1}.The operator M x i has unity matrix elements between any two local spin states on site i.In zero external ("longitudinal") field (λ = 0) this model sits exactly at a continuous quantum phase transition for Q = 2, 3, 4 that becomes discontinuous for (integer) Q > 4 [26].Just above this threshold Q > 4 the first-order discontinuity is exceptionally weak and it is, in fact, notoriously difficult for Monte Carlo simulations to correctly identify the nature of the transition, because of so called pseudo-critical behaviour associated to large remnant correlation lengths at the first order transitions [27,28].This pseudo-critical behaviour is attributed to a fixed point collision in theory space, with complex fixed points arising for Q > 4. We now show that by introducing the external field λ at the transition, a striking qualitative difference can be observed between the continuous and first-order cases at rather large couplings λ and correspondingly short correlation lengths.
In the top row of Fig. 1 we display the zeroth-component magnetization O m = 〈|0〉〈0| i −1/Q〉 λ as a function of the external field λ that are obtained directly in the thermodynamic limit using an infinite system Matrix Product State (iMPS) DMRG algorithm, see App.A for details on the technical aspects.In the left column we display results for the known continuous cases Q = 2, 3, 4, while in the right column results for the weakly first-order cases Q = 5, 6 are shown.The Q = 2, 3, 4 data in the top row exhibits rather clean power law scaling of O m as the coupling λ goes to zero, while in the Q = 5, 6 cases a saturation of O m to a non-zero residual magnetization is observed in the same limit.These limiting values are in very good agreement with the exact results obtained by Baxter [36].In the middle row we calculate numerical finite-difference derivatives of log O m with respect to log λ in order to highlight possible power law behavior O m ∝ λ 1/δ .Indeed at small λ the derivatives for Q = 2, 3, 4 approach the expected values for the known CFTs governing the fixed points: 1/δ = 1/15 for Q = 2, 4 and 1/δ = 1/14 for Q = 3.Note for Q = 4 there is a known logarithmic correction to the power law behavior: m ∝ (λ/ log λ) 1/15 [27], leading to a tiny non-monotonicity for Q = 4.In the weakly first order cases Q = 5, 6 in contrast we observe a pronounced maximum in the derivatives, which has its origin in the inflection point in the original data, c.f. top panel.For our model and Q > 4 there is a saturation in O m both for small and large2 values of λ leading necessarily to (at least) one inflection point at an intermediate value of the external field, which we denote by λ ⋆ (Q).In order to assess at what length scales the pronounced feature of a maximum and the downward drift to zero at smaller λ occurs, we present in the bottom row of Fig. 1 the correlation lengths ξ Q (λ) obtained from the transfer operator of the iMPS wave function for the different values of Q.In the continuous cases Q = 2, 3, 4 we observe again power-law behavior as expected.The most notable observation for the weakly first order cases Q = 5, 6 is that the correlation length ξ Q peak measured at the coupling λ ⋆ (Q) is actually quite small, i.e. about 4 resp. 2 lattice spacings for Q = 5 and Q = 6 respectively.These correlation lengths ξ Q peak are two to three orders of magnitude smaller compared to the huge, albeit finite, correlation lengths at the first order phase transition itself [37,38].
We believe that the pronounced maximum feature of [1/δ](λ) at intermediate values of the coupling λ and the subsequent drift towards zero as λ → 0 + is a robust phenomenon for weakly first order transitions more generally, and it might have its origin in the colliding fixed point scenario advocated for the Q > 4 Potts models.It is notable that the very weak first order transition for Q = 5 has a broader maximum and a relatively slow approach to zero compared to the case Q = 6.It is however striking how different Q = 5 behaves in contrast to the continuous transition The measured running exponents in all three cases monotonically approach the expected value, behaving analogously with the continuous Q = 2, 3, 4 Potts model.We note that the staggered dimer model seems to show an initial fast approach to the O(3) exponent, followed by a slower approach at low fields (right inset).Here we have faded points that we roughly judge by eye to be in the finite size regime.
at Q = 4, despite the presence of a logarithmic correction in the latter case, usually spoiling a clean analysis.These remarkable observations now pave the way to study many open problems in various fields where weakly-first order transitions are hard to discriminate from continuous phase transition with the methods available so far.As an important open question we will address the nature of the phase transition in the J-Q model, which is a candidate for a DQCP.Before tackling this problem, however, let us first validate our approach for a family of 2+1D quantum many body systems with an undisputed continuous quantum phase transition, which we now study using a finite size quantum Monte Carlo method.

Quantum models for the O(3) transition
We consider three models that host a quantum critical point in the 3D O(3) universality class.The first is the well studied square lattice S = 1/2 Heisenberg bilayer system, whose Hamiltonian is given by Here J 1 couples nearest neighbor spins within each square lattice, and J 2 is the coupling between the layers.We have also added an external field λ m that couples to a component of the order parameter, in this case the staggered S z magnetization.
With J 1 , J 2 > 0 this model undergoes a transition from a Néel ordered antiferromagnet for J 2 /J 1 < g c to a dimer singlet phase on the interlayer bonds when J 2 /J 1 > g c , with g c = 2.52181(3) [40].
In order to probe the critical scaling at the transition, we compute the order parameter O m = 〈S z 11 〉 L,λ m on finite size systems of side length L using the stochastic series expansion (SSE) algorithm [41] with directed loops [42].We furthermore have developed a statistically exact Monte Carlo estimator to directly measure the logarithmic derivative ∂ log O m /∂ log λ m that eliminates finite difference errors (see Appendix B), cleanly extracting the running exponent (1/δ) as a function of the field λ m .
In order to paint a more comprehensive picture we also study two other models that host the same O(3) transition but with significant finite-size corrections to critical scaling [43] in one of them.Both models are taken on the square lattice (single layer), where again two different bond strengths (J 2 > J 1 ) are used.Following the nomenclature of [43], we study the columnar dimer model (CDM), consisting of columns of x-oriented strong bonds with a Fourier component (π,0), as well as the staggered dimer model (SDM), consisting of alternating x-oriented strong bonds with a Fourier component (π,π).For the CDM and SDM we use the critical coupling ratios J 2 /J 1 = 1.90951 and J 2 /J 1 = 2.51943, respectively [43].
In Fig. 2 we show all three of the running Néel exponents [1/δ](λ m ) (and the bilayer Néel order parameter in the left inset) at the critical point for the three different models.In the bilayer model we have used J 2 /J 1 = 2.5223, where the difference with the g c quoted above is inconsequential for this plot but shows better agreement with the O(3) exponent in finite size data collapses (see Appendix G).We have used β = 2L for the bilayer model and β = L/2 for the CDM and SDM as in [43], all in units with J 1 = 1 (see Appendix F showing negligible temperature effects).
The bilayer order parameter (left inset) as well as the CDM and SDM order parameters (not shown) all display clean power law scaling at low fields, where the dashed line shows the expected power law for reference 1/δ = 0.2091(1) [39].The measured running exponents (main panels) provide a more fine-grained view of the approach to the expected power law as the field is lowered, where the O(3) exponent is again plotted as a dashed line.In the finite size setup we are using, there is an L-dependent crossover scale for λ m , below which one ultimately observes O m ∝ λ m , a generic result for any finite-size system in the limit λ m → 0. This phenomenon explains why the derivative curves for a given L start to bend upwards at small λ m .As the system size increases, this finite size regime is pushed to smaller values of λ m and a consistent picture representative for finite λ m at L → ∞ emerges.The infinite size (and zero temperature) converged data reveals a monotonic increase of the running exponents, which approaches the 1/δ value expected for a 3D O(3) Wilson-Fisher universality class [39].
Remarkably, even in the SDM, where sizeable corrections to scaling and non-monotonicity of finite-size effective exponents have previously been reported [43,44], we observe a clean monotonic approach to the expected exponent.We note that within our numerical resolution the SDM running exponent seems to contain a regime of fast approach at higher fields, giving way to a much slower approach at lower fields.Although a more careful scaling analysis of the SDM would be desired in this context, we can clearly see the broad picture that is captured by all three of these critical models.While the comparatively slow convergence of [1/δ](λ m ) towards the expected 1/δ value renders our approach less useful to accurately determine 1/δ, we emphasize that the important result here is the absence of a pronounced maximum in [1/δ](λ m ) and the subsequent lack of a drift towards zero as λ m → 0 + .This demonstrates that implementing our method using a finite-size QMC method for a well understood continuous 2+1D quantum phase transition leads to behavior in clear analogy to the continuous phase transitions in the Q = 2, 3, 4 Potts cases studied in 1+1D with iMPS.

O vbs
Figure 3: The Néel and VBS order parameter running exponents for the square and rectangular lattice J-Q models tuned at the transitions J/Q = 0.0447 and J x /Q x = 0.205, respectively.For the square lattice we use β = L/2 (Q = 1) and for the rectangular lattice β = L x /2 (Q x = 1).The running exponents show a local maximum and persistent drift at low fields, behaving as the Q = 5, 6 Potts model.We observe a striking similarity between the known first-order rectangular case and the square lattice case, providing compelling evidence that the transition remains weakly first-order in the square lattice case as well.
5 The J-Q models Finally we turn to a main objective of this work, which is to shed light on the nature of the quantum phase transition between Néel order and valence bond solid (VBS) order in two dimensional S = 1/2 spin systems, thought to be described by the deconfined criticality scenario.An important difference to the previously discussed cases is that the deconfined criticality scenario describes the transition between two ordered phases, therefore we have to track the behaviour of two separate order parameters at the transition point.At a continuous phase transition we expect both running exponents to approach their corresponding values dictated by the universality class in question.In contrast, at a first order phase transition the two order parameters are both expected to be finite, as the two ordered phases coexist at the transition point.
The most well studied model in this context is referred to as the J-Q model [3], written as Here J > 0 is the antiferromagnetic coupling between nearest neighbors S = 1/2 spins on a square lattice, and Q > 0 is a product of two adjacent J terms acting on an elementary square of four spins in both the x and y orientations.When Q = 0 we are left with the Heisenberg model which has Néel order, and conversely when J = 0 the spins form a columnar VBS phase [3].At a small value of the coupling ratio J/Q ≈ 0.04, the transition between these seemingly unrelated orders takes place.Currently, the true nature of the transition in the J-Q model is still under debate.While Refs.[3][4][5][7][8][9] and [6,16] interpret their data as being in favour of a continuous quantum phase transition, it appeared that the extracted critical exponents show pronounced drifts as a function of the maximal system size.However in these simulations no direct evidence for first order behaviour has ever been seen, such as a negative Binder cumulant or multiple peaks in histograms of the energy or order parameters.There are however some papers claiming to observe first order behaviour using the flowgram method [10][11][12].
To probe the nature of this transition we perform two separate studies: one in which a staggered magnetic field coupling to the Néel order parameter is added, and another with a field coupling to the VBS order parameter.As with the O(3) models before, the J-Q model with the Néel field is written with the observable O m = 〈S z 1 〉 L,λ m .We write the model with the VBS field as which preferentially selects one of the four columnar VBS patterns.The VBS order parameter is computed as the expectation value of the difference between even and odd x-bonds Again, we have developed statistically exact QMC estimators for the logarithmic derivatives in both models (see Appendix B).
We furthermore compare the behavior of the J-Q model to a known first-order Néel-VBS transition that is realized by introducing rectangular lattice anisotropy, as was previously studied in [16].Following this methodology, we take spatially anisotropic couplings J y /J x = 0.8 and Q y /Q x = 0.8 on rectangular lattices with L x = 4L y /3.In this model we do not have a prior estimate of the transition, so it was located by scanning the binder ratio of the staggered magnetization (in zero field) for several system sizes (see Appendix E).Here we find a rough estimate of the transition, in this case J x /Q x ≈ 0.205, is more than enough precision for the results we present here.Just as in the square lattice J-Q model, we sit at the transition and introduce separate Néel and VBS fields while measuring the running exponents.
In Fig. 3, we show the running Néel and VBS exponents as in Eq. ( 1) for both the square lattice and rectangular lattice J-Q models tuned to their respective transitions.For the square lattice J-Q model we have used the transition value J/Q = 0.0447 [7].The left panel presents data for each model coupled to the Néel order parameter field, while the right panel similarly presents data for both models coupled to the VBS order parameter field.We clearly observe strong deviations from critical power law scaling for both models with both effective exponents drifting toward zero at small field values, suggesting coexisting Néel and VBS order at the transitions.While our available system sizes do not allow us to track the running exponents all the way to zero, they nevertheless approach closely (surpassing, in the rectangular case) the unitarity bound for scalar operators in a 2+1D CFT ∆ φ ≥ 1/2 [45], yielding a lower bound 1/δ ≥ 1/5.The downward drift of the running exponents is substantial and the contrast to the behavior observed in the continuous O(3) models shown in Fig. 2 on the same vertical scale is striking.We further emphasize the similarity of the running exponents between the known first-order rectangular case and the square lattice case, as well as point out the resemblance to the behavior observed in the Q = 5, 6 Potts model, painting a compelling picture that the square lattice J-Q model is weakly first order.

Discussion and Outlook
Working tangentially to the current symmetry-preserving studies of quantum phase transitions by reintroducing the classic definition of the order parameter in a modern context, we have pushed the sensitivity to diagnose weakly first-order transitions to an unprecedented level.As an important application we have shown that the SU(2) J-Q model on the square lattice does not host a genuine DQCP, but instead a weakly first order transition with coexisting Néel and VBS order at the quantum phase transition.This important result also corroborates recent field theoretical arguments claiming an absence of a genuine DQCP with SO(5) symmetry in 2+1D [24,25,46] and validates early numerical simulations claiming first order behavior based on a flowgram analysis [10][11][12].Furthermore it puts the J-Q model on the same footing as a 3D classical loop model studied in [22], which shows also indications for a weak first order transition, and is expected to realize yet another lattice version of the same NCCP 1 field theory as the J-Q model studied here [1,2].
In view of these results, it is clear that many previously studied models using similar methods need to be revisited [13,16].As an important next step, one can determine where the critical window as a function of N begins in the SU(N ) J-Q models on the square lattice [4,6,15,16].Our results might also have implications for phase transitions out of Dirac spin liquids by virtue of a conjectured duality [21].
Since our approach effectively allows for a controlled study as a function of the correlation length at the transition, it should naturally find applications in 2D tensor networks with applied perturbations (see [47,48] for studies of the 2+1D Ising model and a coupled Heisenberg spin ladder at criticality) to probe the existence of DQCP in frustrated quantum magnets [49,50], where QMC is not applicable due to the sign problem.A closely related important study now within reach is to probe the existence of SU(N ) Dirac Spin Liquids which are conjectured to exist in many frustrated spin models [51][52][53][54][55], and whose field theoretical description is fermionic QED 3 with N f = 2N massless fermion flavors.In the SU(2) Dirac spin liquid context natural perturbations are related to fermion bilinears and monopole operators recently characterized for various lattices in Ref. [56].
The fact that in our approach moderate lattice sizes are typically sufficient to detect weak first order behaviour suggests an immediate applicability for fermionic determinantal QMC methods which typically operate at smaller system sizes compared to the QMC methods used in this work.Exotic quantum phase transitions related to those discussed in the present work have been reported for interacting fermion systems and might warrant an independent confirmation using our technology.
On a more speculative note it will be worthwhile to explore the possibility to transport the ideas developed and demonstrated in this work to lattice field theory simulations of QED, QCD or related theories of importance to high energy physics.
The DMRG and QMC source code as well as the data and plotting scripts for the main figures (Figure 1 to 3) are freely available online [57].

A iMPS for Potts model
In this appendix we present complementary information regarding the iMPS study of the onedimensional Q-state quantum Potts model presented in Sec. 3.
We define the Q-state Potts model on a spin chain of N sites and Q spin states per site denoted as |q i 〉 with q i ∈ {0, 1, . . ., Q − 1}.Writing the Hamiltonian in the eigenbasis of the interaction term we get where K i is a diagonal matrix with its eigenvalues being the Q-th roots of unity, and T being the spin-flip operator: The perturbation strength λ couples to the sum over the projectors onto a single local state, here the |0〉-state.It is customary to also add a perturbation with a transversal field with coupling strength g to the Potts Hamiltonian, for g = 0 the Hamiltonian ( 5) minus an energy density of 3 is equal to the Hamiltonian given in (2).At λ = 0 the Hamiltonian is invariant under the global action of the symmetric group S Q and it exhibits two phases.For g < 0 the system is ordered and the ground state space is Q-fold degenerate with the ground state breaking the S Q symmetry.When g > 0 the system is disordered and the non-degenerate groundstate preserves the symmetry.These phases are separated by a phase transition at λ = 0 and g = 0 which is of second order for Q ≤ 4. For Q = 2 the Potts Hamiltonian reduces to the transverse field Ising model.In the following section we assume that g = 0.
To calculate the properties of its ground state for arbitrary Q and λ in the thermodynamic limit we employ the generalization of the Density Matrix Renormalization Group algorithm to infinite spin chains (iDMRG) [58].It works by decomposing the state into a finite number of rank 3 tensors which are repeated infinitely along the chain and thus form the unit cell.In the present paper only unit cells of size two are used.This set of tensors is called an infinite Matrix Product State (iMPS) which allows us to approximate the state of the system by introducing a cutoff of the tensors' bond dimension χ.The approximation limits the amount of entanglement contained in the state with the entanglement entropy being capped at S = ln(χ).For systems at a quantum critical point the entanglement entropy of the ground state diverges, however since we are perturbing the critical systems with a relevant field the Hamiltonian is gapped, and for large enough bond dimension our iMPS representation is basically numerically exact.The computational difficulty increases with increasing correlation length, i.e. with smaller values of λ.
A necessary requirement for the iDMRG is an efficient Matrix Product Operator (MPO) representation of the Hamiltonian which is easy to achieve for models with only nearest-neighbour interactions.The energy expectation value of an MPS is then expressed as a tensor contraction of MPS and MPO.The algorithm variationally minimizes the energy by sweeping through the system, optimizing the tensors of 2 neighboring sites at a time, until the energy and entanglement entropy converge.In every sweep the eigenvalue problem is projected into these 2 sites and solved using the Lanczos algorithm which is based on calculating the action of the projected Hamiltonian on a wavefunction many times.
In order to speed up the tensor contractions we make use of the Hamiltonian's symmetry.At λ = 0 the Q-state Potts model is invariant under the global action of the non-abelian symmetric group S Q , for λ ̸ = 0 this symmetry is reduced to its subgroup S Q−1 .It is technically much easier to deal with abelian groups thus we only consider the Z Q and Z Q−1 subgroups respectively.In order to exploit the Hamiltonian's symmetry it needs to be decomposed in the right way.To achieve this we apply a global on-site unitary transformation: where the . Note that the zeroth-component is unchanged: | 0〉 = |0〉.In the Z Q case, which is not relevant in our paper, one should make the transformation |ñ i 〉 = Ũ(Q) |n i 〉 instead.In this basis the Hamiltonian can be written as where the projectors are defined as P i = |0〉〈0| i and D i = 1 i − P i .The spin-flip operator is modified to 〉 and R i in this basis is given as We are using TeNPy's implementation of the iDMRG algorithms as well as its methods for optimizing tensor contractions exploiting abelian symmetries [59].
Finally we are calculating the ground state of the Potts model for multiple values of the perturbation strength λ going as close to criticality as numerically feasible.To decrease the compute time we calculate each groundstate for a specific λ by using the previously obtained result of the next higher perturbation strength as initial state to the iDMRG, starting at the product state at λ → ∞.
For each value of λ the expectation value of the projector 〈|0〉〈0| i 〉 λ is evaluated.At λ = 0 the Z Q symmetry of the unperturbed Hamiltonian implies an equal expectation value for all Q components, 〈|0〉〈0| i 〉 λ→0 = 1/Q, thus the order parameter for the Potts model is defined as O m := 〈|0〉〈0| i − 1/Q〉 λ .By numerically computing the logarithmic derivative we get the exponent This is done for all values of Q which are of interest and multiple bond dimensions χ as long as it is numerically feasible, our most sophisticated calculations use χ = 2048 and require up to 3000 sweeps.The results shown in Fig. 1 are converged in χ and have a negligible error in the numerical derivative.
The iMPS description of a state also allows us to easily obtain the correlation length by analyzing the eigenvalue spectrum of the transfer matrix.The correlation lengths shown in Fig. 1 were extrapolated to infinite bond dimension by the method outlined in Ref. [60].

B QMC simulations and exponent estimators
We begin by restating the Hamiltonian in the absence of external fields: The two separate perturbed models are defined by and Here, in order to keep our formulae as simple as possible, we have opted for the notation h ≡ λ m and d ≡ λ vbs .
Beginning with H m JQ , the h−field introduces no sign problem since the perturbation is diagonal.However, the field does prohibit the mapping to a deterministic loop model.We therefore simulate the model using the stochastic series expansion (SSE) algorithm [41] with directed loops [42].On a technical level, we find it convenient to absorb the staggered field into the Heisenberg bond operator, while leaving the Q−term in tact.As a result, loop updates on the Q interactions remain deterministic and non-deterministic decisions only need to be made when updating the bond operators.
Directly measuring the Néel order parameter (staggered magnetization) is straightforward in the presence of the external staggered field, since it aquires a nonzero value.More remarkably, the presence of the field allows us to devise a QMC estimator for the effective critical exponent.
To define this, we first need to be more explicit about the form of the Hamiltonian, referring the reader to [41] for more general information about the SSE framwork.
First note that the external field only effects diagonal matrix elements of the bond operator, which are given by diag(0, J 2 + h B , J 2 − h B , 0) in the basis {↑↑, ↑↓, ↓↑, ↓↓}, and h B is equal to h divided by the coordination number of the lattice.We have also pulled out an overall minus sign.We now shift the bond operators by h B +ε, so the diagonal part becomes diag(h B +ε, J 2 +2h B +ε, J 2 +ε, h B +ε).Here ε has been introduced to lower the bounce probabilities obtained from solving the directed loop equations, and we have used ε = 4h in our simulations.Now that we know the matrix elements, we can see that the weights of the QMC configurations are proportional to where N Q is the number of Q matrix elements, N J is the number of off-diagonal J matrix elements and N D i are the numbers of different diagonal J matrix elements.Differentiating this weight with respect to h B gives We can now compute ∂ O m /∂ h B , which is the main ingredient for the exponent estimator: where 〈mN D i 〉 C ≡ 〈mN D i 〉 − 〈m〉〈N D i 〉 is the "connected" average and m is the staggered S z magnetization per site.Finally we can write the exponent estimator all together as: To make our measurements as precise as possible, we can average over the entire imaginary time history when computing the staggered magnetization.This is facilitated by only considering matrix elements that change the staggered magnetization as one moves through the operator sequence.
We also note that simulations and measurements are nearly identical in the Heisenberg bilayer model, where in that case again the field is incorporated into the nearest-neighbor J 1 term, and updates on J 2 matrix elements are deterministic.
We now describe the measurement of the VBS effective exponent in the H vbs JQ model, which is slightly more complicated but conceptually similar.The Hamiltonian is given by So the even columns of x-bonds have matrix elements J 2 + d 2 , whereas the odd columns of x-bonds and all y-bonds have matrix elements J 2 as normal.This clearly favors one of the four columnar VBS patterns.The associated order parameter is ) is the singlet projector on sites i, j.O vbs is then just the difference in the expectation value of an even x-bond and an odd x-bond.
First note that these expectation values can be measured within the SSE framework as Where N J x e (N J x o ) are the number of even (odd) x-bonds in the operator string, which is why we have divided by N site /2 to get the value on a single bond.It is also necessary to divide by (J + d) and J, since the number operators give averages of the operators appearing in the Hamiltonian, which are multiplied by those factors.We refer the reader to [61] for useful derivations and formulas for bond operator measurements in the SSE.
As before, we now want to compute the derivative with respect to d.We will show how this is done starting with 〈P 1,2 〉: The derivative of 〈P 2,3 〉 is given by: Now in order to compute the derivatives ∂ ∂ d 〈N J x e 〉 and ∂ ∂ d 〈N J x o 〉, we express the QMC weights as previously.This time the configuration weights are proportional to where N Q is the number of Q-operators, N J y is the number of y-oriented J operators, and N J x e (N J x o ) is the number of x-oriented J operators at even (odd) locations from before.The derivative with respect to d is We then have Where α = e, o and again the subscript C means the "connected" part.We can now express Finally the full expression for the exponent estimator is given by In the end, it is just necessary to make measurements of 〈N J x e 〉, 〈N J x o 〉, 〈N J x e N J x e 〉, and 〈N J x e N J x o 〉.
One can then compute the effective exponent using Eq. ( 25) and the statistical error can be computed by bootstrapping the binned data.

C Detuning from phase transition
While the running exponent directly at the phase transition (i.e. at g c ) of the Potts model has been discussed in detail in Sec. 3 it is interesting to investigate what happens if the system is detuned away from the phase transition, for example because g c is not known precisely enough.
By perturbing the Potts model with a transverse field with coupling strength g as introduced in Eq. (5) the system is taken away from criticality, which serves to illustrate the range over which we observe critical scaling as λ → 0, here shown in Fig. 4. We find that, as expected, in the continuous cases when Q = 2, 3 and when g < 0 (favoring the ordered phase), our running exponents drift toward zero as λ → 0. One may worry that if this were the case studying a generic model, one might erroneously conclude a first order transition.However we see that further increasing g there is a wide swath -akin to a critical fan -where the exponent saturates to a consistent finite value before eventually reaching the linear regime of the disordered phase.This is to be contrasted with the case of Q = 5, 6, where the corresponding fan is absent as λ → 0. We conclude that when applying our methodology to generic models, it may be necessary to study several transition point estimates to conclude first order behavior.We note that in the case of the J-Q model, this subtlety does not appear since we separately study both order parameters at the same transition point.

D QMC versus exact diagonalization
In order to confirm the validity of our QMC simulations and exponent estimators, we compare with exact results obtained on small system sizes.Here we focus on the J-Q model on an L = 4 square lattice.Fig. 5 shows both of the exponent estimators compared to exact diagonalization, where finite differences have been used to compute the logarithmic derivatives.For both the staggered magnetization exponent (left panel) and the VBS exponent (right panel) we have used J = 0.0451 and Q = 1.In both cases we observe agreement within the QMC error bars, which can be seen in the insets where the difference between the QMC and ED values are plotted.

E J-Q model with rectangular lattice anisotropy
A simple way of producing a first-order Néel to VBS phase transition is to introduce rectangular lattice anisotropy into the J-Q model, as was previously studied for general SU(N ) spin symmetry in [16].Adopting this same setup, we take spatially anisotropic couplings J y /J x = 0.8 and Q y /Q x = 0.8 on rectangular lattices with L x = 4L y /3.The rectangular lattice anisotropy induces a two-fold degenerate pattern in the VBS phase.We then can estimate the value of the transition based on the binder cumulant of the staggered magnetization, as measured in the pure model without yet introducing the external order parameter fields.The binder cumulant is defined as In Fig. 6 we measure the binder cumulant for different system sizes as function of J x /Q x , taking β = L x /2 in units where Q x = 1.The step in the binder cumulate is an estimate for the transition point, which we roughly estimate to be J x /Q x ≈ 0.205.This is more than enough accuracy than is needed for the system sizes used in the main text (L x ≤ 64).We note that the first-order nature of the transition is strong enough for us to detect conventional symptoms such as double peaked histograms of our binned staggered magnetization measurements, which are shown in the inset.Comparing histograms near the transition shows the peaks becoming more pronounced with increasing system size, indicating a thermodynamic free energy with distinct local minima.Here we focus on the magnetic signal of the transition by measuring the binder cumulant of the staggered magnetization.A rough estimate of the transition at J x /Q x ≈ 0.205 is more than enough precision for the system sizes used in the main text.We further demonstrate the conventional signs of a first-order transition in this model by showing histograms of our staggered magnetization measurements.A clear double peaked structure emerges with increasing system size near the transition, indicating distinct free energy minima.We note that the data only significantly differs in the finite-size region of the curves (too low field values for a fixed system size), whereas the size independent portion is unaffected.The inset shows the best estimate for y h based on pair collapses of system sizes (L, 2L) plotted as a function of 1/L.This procedure was done for both J 2 /J 1 = 2.52205 and 2.5223, where the latter shows better agreement with the exponent estimate from the literature [39], and is shown in the main panel collapse.We note that these field values are significantly lower than the ones used in the main text.

F Zero temperature convergence
We would briefly like to demonstrate the absence of finite temperature effects in the size-independent portion of our QMC data for the running exponents.In all cases we have chosen β ∼ L, with a prefactor larger than the inverse velocity of spin excitations.This ensures that the imaginary time direction grows sufficiently large as a function of L such that only the ground state contributes in the thermodynamic limit.Once the data from different system sizes begins to overlap, we can then be confident that this portion of the the curve is also converged to zero temperature.We demonstrate this for the square lattice J-Q model in Fig. 7, where we have taken J = 0.447 (Q = 1) and two values of the inverse temperature (β = L/2 and β = L).Here we see that the two data sets only differ in the finite-size regions of the curve, whereas the size independent regions are unaffected.

G Additional bilayer data
Here we present supplementary data for the Heisenberg bilayer near the transition.As can be seen in the main paper, plotting the raw effective exponent as a function of the external field does not provide a high precision determination of the order parameter exponent.In an effort to observe fine-grained resolution of the exponent and to further demonstrate unambiguous critical scaling in this model we perform data collapses [62,63] at the transition and as a function of the external field and system size, as can be seen in Fig. 8.Here by plotting the exponent obtained from pair collapses of system sizes (L,2L), we observe high sensitivity with respect to the value of the transition used during data collection (shown for J 2 /J 1 = 2.52205 and 2.5223 in the inset).Of the two values tested, the best agreement with the exponent quoted in the literature [39] is obtained with J 2 /J 1 = 2.5223 and so this value is used in the data presented in the main text.

Figure 1 :
Figure1: Quantum Q-state Potts Chain: infinite system MPS simulations for the quantum critical point with an applied symmetry breaking field λ.The left (right) column displays the exactly known continuous (weakly first order) transitions for Q = 2, 3, 4 (Q = 5, 6).The top row shows the order parameter O m as a function of the perturbing field λ on log-log axes.The inflection points in the right panel are highlighted with circles.The middle row presents the logarithmic derivative [1/δ](λ), i.e. the running exponent defined in Eq. (1), of the first row, highlighting the convergence towards the expected exponents in the left panel, and the existence of distinct maximum in the weakly first order cases in the right panel.In the bottom row we plot the extracted MPS correlation length, showing that the inflection point and the corresponding maximum in the running exponent for the weakly first order instances occur at correlation lengths of only a few lattice spacings.

Figure 2 :
Figure2: The order parameter (staggered magnetization) exponent [1/δ](λ m ), i.e. the running exponent defined in Eq. (1), for three lattice models realizing the O(3) quantum phase transition: the Heisenberg bilayer, the columnar dimer model, and the staggered dimer model.In the left inset we provide the order parameter as a function of the external Néel field on log-log axes, showing clean power law scaling at low fields, where the expected power law 1/δ = 0.2091(1)[39] is shown for reference with the dashed line.The measured running exponents in all three cases monotonically approach the expected value, behaving analogously with the continuous Q = 2, 3, 4 Potts model.We note that the staggered dimer model seems to show an initial fast approach to the O(3) exponent, followed by a slower approach at low fields (right inset).Here we have faded points that we roughly judge by eye to be in the finite size regime.

Figure 4 :
Figure4: Value of the running exponent 1/δ (colorbar) as function of λ and g.Shown are two models with continuous transition (Q = 2, 3 at χ = 64) and two models with a first order transition (Q = 5, 6 at χ = 256).To improve readability the colorbar is adjusted for each model by relating it to values at the phase transition (g = 0), where we used the exactly known running exponent and the peak indicated in Fig.1for Q = 2, 3 and Q = 5, 6 respectively.

∆Figure 5 :
Figure 5: A comparison of the QMC exponent estimators in the J-Q model compared with exact results on an L = 4 system.We have set Q = 1 and J = 0.0451 in both cases.The insets show the difference between the QMC data and the exact values.

Figure 6 :
Figure6: Locating the Néel to two-fold VBS transition in the J-Q model with rectangular lattice anisotropy.Here we focus on the magnetic signal of the transition by measuring the binder cumulant of the staggered magnetization.A rough estimate of the transition at J x /Q x ≈ 0.205 is more than enough precision for the system sizes used in the main text.We further demonstrate the conventional signs of a first-order transition in this model by showing histograms of our staggered magnetization measurements.A clear double peaked structure emerges with increasing system size near the transition, indicating distinct free energy minima.

Figure 7 :
Figure7: Absence of finite temperature effects in the size-independent region of the running exponents.Here we show data from the square lattice J-Q model as a function of the Néel field taken with J = 0.0447 (Q = 1) at two different inverse temperatures.We note that the data only significantly differs in the finite-size region of the curves (too low field values for a fixed system size), whereas the size independent portion is unaffected.

Figure 8 :
Figure8: Here we demonstrate critical finite-size scaling of the staggered magnetization as a function of system size and external staggered field with J 2 /J 1 = 2.5223.Exactly at the critical point, the scaling ansatz is m = L y h −D f (hL y h ), which makes the quantity hL D m only a function of hL y h .Here y h = D/(1/δ + 1).The inset shows the best estimate for y h based on pair collapses of system sizes (L, 2L) plotted as a function of 1/L.This procedure was done for both J 2 /J 1 = 2.52205 and 2.5223, where the latter shows better agreement with the exponent estimate from the literature[39], and is shown in the main panel collapse.We note that these field values are significantly lower than the ones used in the main text.