First order phase transitions in the square lattice "easy-plane" J-Q model
Nisheeta Desai Ribhu K. Kaul





Department of Physics & Astronomy, 
University of Kentucky, 
Lexington, KY 40506.

Department of Physics & Astronomy, 
University of Kentucky, 
Lexington, KY 40506.



We study the quantum phase transition between the superfluid and valence bond solid in “easy-plane" J-Q models on the square lattice. The Hamiltonian we study is a linear combination of two model Hamiltonians: (1) an SU(2) symmetric model, which is the well known J-Q model that does not show any direct signs of a discontinuous transition on the largest lattices and is presumed continuous, and (2) an easy plane version of the J-Q model, which shows clear evidence for a first order transition even on rather small lattices of size L≈16. A parameter 0≤λ≤ 1 (λ=0 being the easy-plane model and λ=1 being the SU(2) symmetric J-Q model) allows us to smoothly interpolate between these two limiting models.  We use stochastic series expansion (SSE) quantum Monte Carlo (QMC) to investigate the nature of this transition as λ is varied - here we present studies for λ=0,0.5,0.75,0.85,0.95 and 1. While we find that the first order transition weakens as λ is increased from 0 to 1, we find no evidence that the transition becomes continuous until the SU(2) symmetric point, λ=1. We thus conclude that the square lattice superfluid-VBS transition in the two-component easy-plane model is generically first order. 



First order phase transitions in the square lattice “easy-plane” J-Q model
    Ribhu K. Kaul
    December 30, 2023
===========================================================================





§ INTRODUCTION



The quantum transition from Néel or superfluid to a valence bond solid (VBS) has been proposed to be described by the deconfined criticality scenario.<cit.> In this scenario it is generically possible to have a direct continuous Néel-VBS transition. A number of field theoretic formulations that describe this putative critical point at long distances have been put forward and interesting connections between different representations have been conjectured via duality arguments. <cit.> Establishing these fascinating connections non-perturbatively by lattice simulations is an exciting field of current research. In the original study two kinds of symmetries were highlighted for their possibility as platforms for deconfined criticality, an SU(2) symmetric system and a U(1)×Z_2 symmetric system. Physically, the SU(2) field theory could be a description for a rotationally symmetric S=1/2 anti-ferromagnet and its transition to a valence bond solid. The U(1)×Z_2  system can be thought of as a model for the same Néel-VBS transition in  magnet with easy-plane anisotropy or alternatively as a model for a superfluid to Mott transition. [We will use the terms magnetic and superfluid interchangeably throughout this manuscript] 

In the years since the original proposal, it has been demonstrated that the Néel-VBS transition and many of its variants can be studied in sign problem free quantum spin Hamiltonian models on large lattices.<cit.> Through extensive numerical simulations in the SU(2) symmetric models many aspects of the proposal have been borne out and no direct evidence for a first order transition has been observed. <cit.> Numerical studies of classical statistical mechanics models of tightly packed loops and dimer models in three dimensions that have been argued to realize the same universal physics as the SU(2) Néel-VBS transition are also consistent with the deconfined criticality scenario. <cit.> Despite this large body of evidence for the deconfined criticality scenario, numerical studies have observed scaling violations whose origin is currently unclear.<cit.>



In  contrast in the easy-plane case where SU(2) is broken to U(1)× Z_2 a number of numerical studies have concluded that the transition is first order.<cit.> 
Recently however it has been claimed that a continuous transition has been found in a square lattice model with somewhat weaker easy-plane anisotropy,<cit.> suggesting that perhaps a large easy-plane anisotropy could result in a first order transition, and the first order and second order regime are separated by a multicritical point.[We note for completeness that a model of hard core bosons at 1/3 filling on the Kagome lattice has been argued to be described by a similar field theory and host a putative easy-plane deconfined critical point.<cit.> Although it was originally determined to have a first order transition,<cit.> it has been claimed in recent work to host a continuous transition.<cit.>]
Motivated by this study, we address the issue of how the easy-plane transition is connected to the symmetric one,  by studying a model that interpolates between these two limiting cases on the square lattice. For the symmetric model we use the popular J-Q model which shows no direct evidence for first order behavior even on the largest studied lattice sizes. For the easy plane case we use an easy-plane J-Q that was shown to have a first order transition already visible on L≈ 16. The interpolating model introduced in detail below is slightly different from the one studied in Ref. meng2017:dualities,ma2018:prb where the easy-plane anisotropy was introduced only in the J-term; both models are believed to have the same universal features however. In this work we present studies on larger lattices and a more thorough analysis. Contrary to the previous study, we find no evidence for new continuous easy-plane criticality. Instead we find a first order transition for 0≤λ <1 that weakens as λ is increased and we approach the symmetric point (λ=1) at which all our direct signals of a first order transition vanish and the transition is presumed continuous. This is the primary result of our paper and is summarized in Fig. <ref>. Although no numerical study can rule out that the transition becomes continuous for a very small but finite window close to λ=1 (with finite easy-plane anisotropy), we find this rather unlikely given our results below. We thus conclude that the easy-plane Néel-VBS transition is generically first order on the square lattice.



§ THE MODEL





The Hamiltonian studied here is an S=1/2 system on an L × L square lattice,

    H^JQ_λ=λ H^JQ_s + (1-λ)H^JQ_ep,

 and is a linear combination of two parts, H^JQ_s is the SU(2) symmetric part and H^JQ_ep is the easy plane part that explicitly breaks the SU(2) symmetry. λ is an anisotropy 
parameter that allows us to smoothly interpolate between the easy plane limit (λ=0) and the SU(2) symmetric limit (λ=1). We define the 
singlet projection operator on a bond between two sites i and j as,

    P_ij= 1/4- S⃗_⃗i⃗.S⃗_⃗j⃗.

S^μ_i=1/2σ^μ_i, are standard spin-1/2 operators where σ^μ_i are Pauli matrices.

Then H^JQ_s, which is the well known J-Q model, <cit.> can be written as, 

    H^JQ_s = -J∑_⟨ ij ⟩ P_ij - Q ∑_⟨ ijkl ⟩ P_ij P_kl

The second term in the above equation is a sum over all elementary plaquettes ijkl. H^JQ_s has full SU(2) symmetry inherited from P_ij. Similarly if we define 

    P̃_ij= S^x_iS^x_j+S^y_iS^y_j

the easy plane Hamiltonian, H^JQ_ep can be written as,<cit.>

    H^JQ_ep = J∑_⟨ ij ⟩P̃_ij  - Q ∑_⟨ ijkl ⟩P̃_ijP̃_kl

P̃_ij has a symmetry of U(1) × Z_2, which corresponds to U(1) rotations about the ẑ-axis and the Z_2 operation of a π rotation about the x̂-axis. 


We study the quantum phase transition from the magnetic phase to the valence-bond solid (VBS) phase as g ≡ Q/J is varied for a fixed λ. While in the easy plane limit, i.e. λ=0, this transition has been found to be first order, <cit.> it has been argued to be to continuous in the SU(2) symmetric limit, λ=1. <cit.> In this work we interpolate between the two limiting models with the aim of elucidating the evolution of the nature of the quantum transition and in particular to investigate whether the transition is continuous for any λ<1. 




§ NUMERICAL SIMULATIONS




The numerical results presented below have been obtained using the stochastic series expansion (SSE) quantum Monte Carlo method. <cit.> We use the directed loop algorithm <cit.> to carry out global loop updates on our Monte Carlo configurations (see Appendix <ref>).





















Fig. <ref> shows a phase diagram obtained from numerical simulations as a function of the coupling g=Q/J and the anisotropy parameter λ. For a given λ, on increasing g we find a quantum phase transition from the magnetic to the VBS phase. Here, we work in units where √(J^2+Q^2)=1 and at an inverse temperature β=L for an L× L lattice. All data presented has been tested to be in the T=0 limit as demonstrated in Appendix <ref>. 




 §.§ Measurements



When λ < 1, the presence of a small amount of anisotropy makes the spins preferentially align in the XY plane. Therefore as we vary g in our simulations, we look for a phase transition between the XY order (superfluid) and VBS phases. We define the following quantities to detect magnetic order,


    S_m_⊥^2(k⃗) = ∑_r⃗  e^i k⃗.r⃗⟨ S^x_0⃗S^x_r⃗+S^y_0⃗S^y_r⃗⟩



    S_m_∥^2(k⃗) = ∑_r⃗  e^i k⃗.r⃗⟨ S^z_0⃗S^z_r⃗⟩


The square of the superfluid order parameter is ⟨ m_⊥^2⟩ =S_m_⊥^2(π,π). The quantity ⟨ S^x_0⃗S^x_r⃗+S^y_0⃗S^y_r⃗⟩ is measured during the loop update by keeping track of the distance between the head and the tail of the loop when they are at the same time slice. <cit.> We also define a Néel order parameter square as ⟨ m_∥^2 ⟩=S_m_∥^2(π,π). When λ=1.0, S^x, S^y and S^z are equivalent and therefore ⟨ m^2_⊥⟩ and ⟨ m^2_∥⟩ are equal upto a normalization. 

VBS order is one where each spin forms a singlet with its neighbour and the pattern of singlets forms columnar order. We construct the VBS order parameter square from the following quantity,

    S_ϕ^2_x(k⃗) = ∑_r⃗ e^i k⃗.r⃗⟨𝒬_x(0⃗)𝒬_x(r⃗) ⟩.

Here 𝒬_x(r⃗) is a plaquette operator which equals the sum of all the operators in the Hamiltonian acting on the plaquette at r⃗ as described in Appendix <ref>. 


The spin stiffness ρ_s is defined as,

    ρ_s=∂^2 E(ϕ)/∂ϕ^2|_ϕ=0

Here E(ϕ) is the energy of the system with a twist of ϕ in the boundary condition in either the x or the y direction. In the QMC, this quantity is related to the winding number of loops in the direction that the twist has been added, <cit.>

    ρ_s= ⟨ W^2 ⟩/β

 where β is the inverse temperature. ρ_s goes to a finite value in the magnetically ordered phase but goes to 0 otherwise. The quantity Lρ_s is expected to show a crossing for different values of L at the coupling at which magnetic order is destroyed.
 
In order to detect the ordered phase we make use of ratios defined as, 

    ℛ_o_p=1-|S_o_p(k'_o)|/|S_o_p(k_o)|.
 

Here o_p=m_⊥^2,m_∥^2,ϕ^2; k_o and k'_o are the ordering momentum and momentum closest to the ordering momentum respectively. In the ordered phase ℛ_o_p goes to 1 and in the disordered phase it goes to 0 on increasing system size, therefore they are expected to cross for different system sizes at the critical point. 



 §.§ Numerical Results





 
In this work, we focus on four values of the anisotropy parameter, λ=0.5, λ=0.75, λ=0.85 and λ=0.95. We have included a comparison with the symmetric case λ=1 when appropriate.



  §.§.§ Crossing Analysis


Fig. <ref> (a) and (b) show ratios ℛ_m_⊥^2 and ℛ_ϕ_x^2  (defined above) crossing for different L for λ=0.5. This indicates a transition from the magnetic to VBS phase. Fig. <ref>(c) shows crossing of these ratios for the same L.  As shown in <ref>(d), the crossing analysis from <ref>(c) yields the transition point to be at g^*_c=12.111(3), which is close to the value at which the couplings at the crossing points, g_c(L), converge. This extrapolation has been done only using small system sizes, L≤ 64. We notice that smaller system sizes can be seen to smoothly converge to g^*_c≈12.1, bigger system sizes start deviating from this trend. This is because of the double peaked structure that starts to develop in the order parameter estimators, making it difficult to reliably extrapolate g_c(L) using bigger lattices. Fig. <ref> shows crossings of the quantity  Lρ_s for both λ=0.5 and λ=0.75, which also indicates transition out of the magnetic phase.


To investigate the nature of the transition, we study the extrapolation of observables with system size at the critical point. For a continuous transition, all the observables described above (ρ_s, ⟨ m^2 ⟩, ⟨ϕ^2 ⟩), should go to zero at the critical point as L→∞. Figs <ref> and <ref> shows values of ⟨ m^2_⊥⟩ and ⟨ϕ_x^2 ⟩ at the crossing points of ℛ_m_⊥^2 and ℛ_ϕ_x^2 at L extrapolated to the infinite system size limit using two different fitting forms (as described in the caption). The extrapolated values of ⟨ m_⊥^2 ⟩ and ⟨ϕ_x^2 ⟩ are clearly finite for λ=0.5, λ=0.75 and λ=0.85.  Small positive extrapolated values of these quantities using both the fit forms can also be seen for λ=0.95. Fig. <ref> shows how the order parameter values at the transition point get progressively smaller on increasing λ, indicating a weakening of the first order nature. They tend to approach zero as λ→1. Therefore we can argue that the first order transition continuously evolves on varying λ to a second order transition at λ=1.0. 

The stiffness extracted from crossings of Lρ_s for L and L/2 in Fig <ref>(a),(b) is plotted as a function of 1/L in <ref>(c). ρ_s clearly extrapolates to a finite value for 1/L → 0 for both λ=0.5 and λ=0.75. Fig. <ref> shows the same analysis for λ=0.85 and λ=0.95. This points to a first order transition for λ=0.5,0.75,0.85 and 0.95. For λ=1, on the other hand it is apparent from our data that is hard to argue for a finite order parameter for ⟨ m_⊥^2 ⟩ and ⟨ϕ_x^2 ⟩ from the data we have. A more thorough analysis of the λ=1 is available in Ref. Shao2016:science. We note that these extrapolations become hard to do on very large system sizes because of ergodicity issues that we discuss below and that we argue arise fundamentally at first order transitions. 



  §.§.§ Histograms


 
To further elucidate the nature of the transition we carefully study the histograms of observables near the critical point. Fig. <ref> shows the probability distributions of the QMC estimators for ρ_s, m^2_⊥ and ϕ_x^2 at the transition for λ=0.5. There are clearly two peaks in the histograms of ρ_s for L=48, 64, 96, one at 0 and the other at a finite value. This double peak feature is clearly noticeable in m^2_⊥ and ϕ_x^2 only for L=96. The double peak gets more pronounced with system size which indicates that the first order behavior survives in the thermodynamic limit. The time series data shows switching between the two orders: one order parameter is finite when the other goes to 0, thus one order is present when the other is not. This is characteristic of a first order transition. This system exhibits clear first order behavior only for L>64, therefore we conclude that this transition is a weak first order transition. The first order nature of the transition is even weaker for λ=0.75, the double peak in the histograms of stiffness appears for L>96 as shown in Fig.  <ref>. We find no evidence of double peaked histograms for λ=0.85 and 0.95 for the largest system sizes studied here. However, as explained before, we do find other evidence of first order behaviour in these two cases. There is also no evidence of double peaked histograms for λ=1.0 for the largest system size we have studied. This is consistent with the numerical findings in the past. <cit.>  Therefore we conclude that the transition is first order for λ=0.5,0.75,0.85 and 0.95. The first order behaviour gets progressively weaker as λ gets closer to 1, eventually disappearing at λ=1.
 
 
 


§ CONCLUSIONS




We studied an interpolation of two previously known and well studied models, the J-Q model <cit.> which hosts a continuous Néel-VBS transition and the easy-plane J-Q model <cit.> which hosts a first order superfluid-VBS transition. By studying the phase transition as a function of the parameter λ that interpolates between the two limits, we found the phase diagram shown in Fig. <ref>. Our main conclusion is that whenever the easy-plane anisotropy is present the transition is first order. All signs of discontinuity vanish  only at the symmetric point λ=1. This indicates that the easy-plane anisotropy is a relevant perturabation at the SU(2) symmetric deconfined critical point and results in a runaway flow to a first order transition.



We acknowledge helpful discussions with A. Sandvik. Partial financial support was received through NSF DMR-1611161 and Keith B Macadam Graduate Excellence Fellowship. Computing resources were obtained through NSF's XSEDE award TG-DMR-140061 and the DLX computer at the University of Kentucky.






§ DETAILS OF NUMERICAL SIMULATIONS AND CHECKS







 §.§ Lattice Hamiltonian


The Hamiltonian defined by Eq. <ref> is sign problem free on a bipartite lattice [The following unitary transformation: S_x → -S_x and S_y → -S_y, on one of the sublattices of the bipartition, yields all negative off-diagonal elements for the Heisenberg exchange which is the sign problem free condition] and therefore we use the SSE QMC algorithm with directed loop updates to simulate it. The easy plane limit of this model (λ=0) has no diagonal terms. Hence, to make the model easier to simulate in this limit, we add a constant to the Hamiltonian to generate diagonal matrix elements.<cit.> The easy plane part of the model defined in Eq. <ref> then becomes: 

    H^JQ_ns=J∑_⟨ ij ⟩ (P̃_ij+1_ij)  - Q ∑_⟨ ijkl ⟩P̃_ijP̃_kl

To make the loop update more convenient we treat all bonds as plaquettes by multiplying an identity to the adjacent bond, for e.g. the P_ij operator in Eq. <ref> gets replaced in the following way:

    P_ij=1/N^b_plaq∑_klP_ij.1_kl

Here 1_kl is an identity operator, the sum in this equation is over all four site plaquettes ijkl such that kl is adjacent and parallel to ij. N^b_plaq is number of plaquettes each bond is a part of, which is 2 in the square lattice case. After making these substitutions the full Hamiltonian described by Eq. <ref> becomes:


    H=λ{J/2∑_ijkl(P_ij.1_kl+1_ij.P_kl)+Q  P_ij.P_kl)}+ 
     (1-λ){J/2∑_ijkl(P̃_ij.1_kl+1_ij.P̃_kl+1_ij.1_kl)+Q P̃_ij.P̃_kl}
 









 
 


 §.§ Plaquette Operator


The plaquette operator, 𝒬_x(r⃗) in Eq. <ref> is the sum of all operators in the Hamiltonian acting on the plaquette at r⃗. Let r⃗_ijkl be the position vector of the lower left site of the plaquette ijkl 


    𝒬_x(r⃗_ijkl)=J/2{λ (P_ij.1_kl+1_ij.P_kl)
    +(1-λ) (P̃_ij.1_kl+1_ij.P̃_kl+1_ij.1_kl)}
    +Q {λ  P_ij.P_kl+(1-λ) P̃_ij.P̃_kl}
 






 §.§ QMC vs ED

Tables <ref> and <ref> show comparison of the groundstate energy per unit site (e), spin stiffness (ρ_s) and square of the order parameters, ⟨ m_⊥^2 ⟩ and ⟨ϕ_x^2⟩, got from QMC and from exact diagonalization for λ=0.5 and λ=0.75 on 4×4 lattices. ρ_s, m^2_⊥ and ϕ_x^2 are as defined in Sec. <ref>.



 §.§ Convergence to T=0


Fig. <ref> shows the behaviour of the observables we have measured for L× L square lattices as a function of inverse temperature β. The measurements have been done close to the critical points (g=12 for λ=0.5 and g=16 for λ=0.75). These quantities can be seen to saturate to the T=0 value on increasing the value of β. The β at which this saturation occurs depends on the system size L. As we increase the system size these values saturate to the value at β=L faster, therefore we pick β=L for our simulations.




apsrev



