First- and second-order quantum phase transitions in the long-range unfrustrated antiferromagnetic Ising chain
Víctor Herráiz-López Sebastián Roca-Jerat Manuel Gallego Ramón Ferrández Jesús Carrete David Zueco Juan Román-Roche





These two authors contributed equally
Departamento de Física Aplicada, Universidad de Zaragoza, Zaragoza 50009, Spain


These two authors contributed equally
Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física de la Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain




Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física Teórica, Universidad de Zaragoza, Zaragoza 50009, Spain


Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física de la Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain


Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física de la Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain


Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física de la Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain


Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain
Departamento de Física de la Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain






  


§ ABSTRACT 
    We study the ground-state phase diagram of an unfrustrated antiferromagnetic Ising chain with longitudinal and transverse fields in the full range of interactions: from all-to-all to nearest-neighbors. First, we solve the model analytically in the strong long-range regime, confirming in the process that a mean-field treatment is exact for this model. We compute the order parameter and the correlations and show that the model exhibits a tricritical point where the phase transition changes from first to second order. This is in contrast with the nearest-neighbor limit where the phase transition is known to be second order. To understand how the order of the phase transition changes from one limit to the other, we tackle the analytically-intractable interaction ranges numerically, using a variational quantum Monte Carlo method with a neural-network-based ansatz, the visual transformer. We show how the first-order phase transition shrinks with decreasing interaction range and establish approximate boundaries in the interaction range for which the first-order phase transition is present. Finally, we establish that the key ingredient to stabilize a first-order phase transition and a tricritical point is the presence of ferromagnetic interactions between spins of the same sublattice on top of antiferromagnetic interactions between spins of different sublattices. Tunable-range unfrustrated antiferromagnetic interactions are just one way to implement such staggered interactions.



First- and second-order quantum phase transitions in the long-range unfrustrated antiferromagnetic Ising chain
    Juan Román-Roche
    September 9, 2024
==============================================================================================================





§ INTRODUCTION



Recent advances in cold-atom simulators have led to renewed interest in systems with long-range interactions <cit.>. Long-range interactions decay as r^-α, with r the distance between interacting degrees of freedom <cit.>. Despite being ubiquitous in nature, e.g. dipolar or gravitational interactions, these systems have been less studied than their nearest-neighbor counterparts because long-range interactions complicate numerical and analytical treatments. Nevertheless, some results exist that showcase remarkable differences in behavior between long- and short-range models. Some examples are the spontaneous breaking of continuous symmetries in one dimension <cit.>, the existence of an area law of entanglement <cit.>, the existence of Majorana modes <cit.>, the spreading of correlations <cit.> and topological properties <cit.>. 
Notably, long-range interactions have been shown to induce first-order phase transitions in classical dipolar gases <cit.> and (quantum) Bose-Hubbard <cit.> and XX <cit.> models.


Quantum phase transitions constitute one of the fundamental phenomena of condensed matter physics. They manifest as nonanaliticities in the ground state energy as some critical parameter is varied <cit.>. Leaving aside the more exotic topological phase transitions <cit.>, quantum phase transitions, like thermal phase transitions, can be first or second order, depending on whether the order parameter or its derivatives are discontinuous at the critical point.

Just as the Ising model is the paradigmatic model in classical statistical mechanics, its quantum analogue, the (nearest-neighbour) transverse field Ising chain (TFIC), is the paradigmatic example of a solvable model featuring a quantum phase transition <cit.>. The TFIC is also a starting point to devise more sophisticated models.




To elucidate the effect of long-range interactions on the order of phase transitions,  it is convenient to construct a minimal model. The trivial long-range generalization of the TFIC is known to feature only a second-order phase transition <cit.>. In fact, there is no phase transition at all for antiferromagnetic interactions at extremely long ranges. At the same time, it has been shown that the Ising model with antiferromagnetic nearest-neighbor interactions and ferromagnetic next-nearest-neighbor interactions presents a tricritical point (TP) where the critical line changes from first to second order <cit.>.
To shed more light on this issue, we have generalized the antiferromagnetic Ising chain to feature tunable-range unfrustrated antiferromagnetic interactions and studied its ground-state phase diagram analytically in the strong long-range regime using the technique described in Ref. . We find that the phase diagrams of the nearest-neighbor and strong long-range models differ significantly. The nearest-neighbour model presents a second-order critical line between an antiferromagnetic and a paramagnetic phase, with a singular point -for vanishing transverse field, where the model becomes classical- where the phase transition becomes first order <cit.>. In contrast, our results show that in the strong long-range model, a significant portion of the critical line between the antiferromagnetic and paramagnetic phases becomes first order. The critical line changes order at a tricritical point that occurs at non-zero transverse field. 
To understand how the critical line morphs from one limit to the other, we employ a variational quantum Monte Carlo method with a neural-network-based ansatz, the visual transformer <cit.>, in the full range of interactions 0 < α < ∞, with α→∞ corresponding to the nearest-neighbor limit. These numerical results show how the position of the tricritical point smoothly moves from zero to a finite transverse field. Remarkably, the first-order phase transition survives well into the weak long-range regime (α > d with d the dimension of the lattice, and specifically d=1 for the chain). Finally, we check whether the first-order phase transition is present if the ferromagnetic intrasublattice interactions are removed. In this case, we find that the full critical line is second order. This indicates that the key to stabilize a first-order phase transition is the simultaneous presence of antiferromagnetic intersublattice and ferromagnetic intrasublattice interactions. Tunable-range antiferromagnetic interactions are just one way to implement and tune these staggered interactions.

The rest of the paper is organized as follows. In Sec. <ref> we present the Hamiltonian for the tunable-range unfrustrated antiferromagnetic Ising chain. Section <ref> is dedicated to the exact solution of the model in the strong long-range regime, with a characterization of the ground state phase diagram through the order parameter and the correlations. The numerical results are presented in Sec. <ref>. We end the paper with a discussion of the results and the source of the first-order phase transition in Sec. <ref> and provide technical details and complementary results in the appendices.



§ THE MODEL



We consider a one-dimensional spin chain (d=1) with tunable-range Ising interactions and subject to transverse and longitudinal fields. The Hamiltonian reads


    H = -ω_z ∑_i=1^N S_i^z - ω_x ∑_i=1^N S_i^x - ∑_i,j=1^N J_ijS_i^xS_j^x   ,


where S_i^x,z are spin-s operators acting on site i. The interactions are J_ij = (-1)^i+jΓJ̃(𝐫_ij)/Ñ, with


    J̃(𝐫_ij)= 
    	b     if   𝐫_ij=0 
    
    	|𝐫_i j|^-α    otherwise


where the distance 𝐫_ij is given by the nearest image convention using periodic boundary conditions (PBC). Γ >0 is the interaction strength, b is a parameter that can be tuned to shift the spectrum of J, Ñ = ∑_i J̃_ij is Kac's renormalization factor and α is the coefficient that tunes the range of interaction. The alternating sign makes the interactions antiferromagnetic when spins are separated by an odd number of lattice parameters and ferromagnetic when spins are separated by an even number of parameters, see Fig. <ref>. According to the classification of long-range interactions valid for both classical and quantum models, the strong long-range regime corresponds to α < d=1 <cit.>. This regime is characterized by the loss of additivity, although extensivity is preserved by Kac's renormalization factor, Ñ, ensuring a well-defined thermodynamic limit. In the limit α→∞, the Hamiltonian [Eq. (<ref>)] corresponds to the nearest-neighbor Ising chain with transverse and longitudinal fields.
The tunable-range staggered interactions that we consider here are a generalization of the antiferromagnetic nearest-neighbor interactions in unfrustrated lattices. The alternating sign prevents frustration and allows for the formation of two sublattices, as sketched in Fig. <ref>. For vanishing fields, the ground state is the antiferromagnetic configuration, with the spins fully polarized along the x axis in alternating directions for even and odd spins.  Staggered tunable-range interactions have been proposed previously for the Heisenberg model <cit.>







§ EXACT SOLUTION IN THE STRONG LONG-RANGE REGIME





 §.§ Canonical partition function


Reference  presents an exact analytical solution for quantum strong long-range models in the canonical ensemble and the thermodynamic limit (N →∞), extending previous classical results <cit.>. We sketch its application to Hamiltonian (<ref>) here; the details can be found in App. <ref>.

The Hamiltonian is divided into free and interacting parts


    H = H_0 - ∑_i,jJ_ijS_i^x S_j^x   ,
 

First, the interaction matrix is diagonalized


    J_ij = 1/N∑_k=0^N-1λ_ik D_k λ_jk  ,


where D_k are its eigenvalues and λ_ik/√(N) is an orthogonal matrix because J_ij is symmetric. The smallest eigenvalue is set to zero by tuning the on-site interaction parameter b (<ref>). For s≠ 1/2, setting b≠0 introduces new terms in the Hamiltonian, but due to Kac's renormalization factor, Ñ, these are negligible in the thermodynamic limit.
Then, the Hamiltonian can be mapped to a generalized Dicke model <cit.>, replacing the long-range interactions by effective interactions between each particle and a set of M real fields, { u_k }, where M is the number of non-zero eigenvalues of J.
The mapping is only valid if lim_N→∞ M/N = 0, restricting the applicability of the method to models in which the number of non-zero eigenvalues of the interaction matrix is a negligible fraction of the total.
This is the reason why the solution is restricted to the strong long-range regime, where α<d ensures that this condition is met.
Following Wang's procedure, we can obtain the canonical partition function of the corresponding Dicke model <cit.>, thus solving our original model. This yields


    Z = exp(-N β f[u̅_k]) ∏_p=0^M-1√(2  D_p)  ,


where f[u̅_k] = min_{u_k} f [u_k] and f [u_k] = f(u_0, …, u_M-1). By noticing that in the thermodynamic limit f[u̅_k] is the free energy per particle: f[u̅_k] = lim_N →∞ - (Nβ)^-1ln Z, it becomes apparent that computing the partition function is reduced to a minimization of the variational free energy per particle


    f[u_k] = ∑_k=0^M-1u_k^2/D_k + f_ m[u_k]


with respect to the auxiliary fields {u_k}. Here f_ m = F_ m / N, with


    F_ m[u_k] = -1/β∑_i=1^N ln[  sinh( (2s+1)βε_i[u_k] )/sinh(βε_i[u_k] )]  ,


and



    2ε_i [u_k] = √(ω_z^2 + ( ω_x + 2∑_k=0^M-1λ_ik u_k  )^2)  .


Although computing the canonical partition function has been reduced to solving a minimization problem, finding the global minima of a multivariate function is not a simple task, and success is not guaranteed in most cases. However, for α = 0 (homogeneous all-to-all couplings which are ferromagnetic or antiferromagnetic depending on the distance between the spins) the problem naturally becomes univariate, as there is only one non-zero eigenvalue of the interaction matrix J, so M=1. This is equivalent to setting u_k ≠ 0 = 0, with u_0 corresponding to the largest eigenvalue: D_0 = max{D_k} = Γ and λ_i0 = (-1)^i. Additionally, we find that this solution with u_k ≠ 0 = 0 is the global minimum also for 0 ≠α < 1. This can be shown analytically for ω_x = 0 (See. Appendix <ref>) and has been verified numerically otherwise. With this, computing the canonical partition function has been reduced to a univariate minimization of f(u) ≡ f(u, 0, …, 0) in terms of u ≡ u_0 for all α < 1, which can be tackled analytically or numerically. In the remainder of this section, we do so to obtain the ground state phase diagram of the model and the susceptibilities. This also implies that the equilibrium properties of the model are universal for all α < 1. In fact, f(u) coincides with the free energy obtained with a mean-field approach, proving that mean field is exact also for strong long-range unfrustrated antiferromagnetic models. This is consistent with the fact that unfrustrated antiferromagnetic models can be mapped to ferromagnetic models in a staggered field and a mean-field description of strong long-range ferromagnetic models has been shown to be exact <cit.>.



 §.§ Ground State Phase Diagram



In order to study the ground state of the model, we define the variational ground-state energy as the zero-temperature limit of the variational free energy e_0(u) = lim_β→∞ f(u), such that


    e_0(u) = u^2/Γ
        - s ( ε_+(u) +  ε_-(u))  ,


with


    2 ε_±(u) = √(ω_z^2 + (ω_x ± 2u)^2) .


It can be shown that the global minimum, u̅, is proportional to the staggered magnetization (See App. <ref>) 


    m̅_ s = 1/N∑_i=1^N (-1)^i ⟨ S_i^x ⟩ = u̅/Γ  .


The staggered magnetization is the order parameter of unfrustrated antiferromagnetic models such as Hamiltonian (<ref>). It is zero in the paramagnetic phase and it measures how close the ground state is to a perfect antiferromagnetic configuration in the antiferromagnetic phase. Figure <ref>(a) shows the phase diagram of the model. The staggered magnetization is obtained from u̅, by minimizing e_0, and plotted as a function of the longitudinal and transverse fields. The model exhibits a quantum phase transition (QPT) between an antiferromagnetic phase and a paramagnetic phase. The nature of the phase transition changes along the critical line from a second- to a first-order QPT. The change can be described analytically by applying the Landau theory of phase transitions to a series expansion of e_0(m_ s), with m_ s = u/Γ the variational staggered magnetization.  By studying the change of sign of the expansion coefficients we obtain the tricritical point


    ω_x,  tp = sΓ8/5√(5)  , 
    ω_z,  tp = sΓ16/5√(5)  ,


and the equation for the second-order portion of the critical line


    4s^2Γ^2 ω_z^4 = ( ω_z^2 + ω_x^2 )^3   ,


when |ω_z | > 2|ω_x |.
The first-order critical condition cannot be obtained from this analysis because the global minima corresponding to the antiferromagnetic configurations fall outside of the radius of analyticity of the series expansion of e_0(m_ s). Nevertheless, having identified the second-order portion, the rest of the critical line is first order by exclusion. We verify this graphically in Fig. <ref>(c) where we show the landscape of minima of e_0 for different values of ω_x when ω_z / (sΓ) = 0.2. For ω_x / (sΓ)=0.8 there exists a local minimum at m_ s=0, corresponding to the paramagnetic state, and two degenerate global minima at m_ s≈± 1, corresponding to the symmetric antiferromagnetic ground states. For ω_x / (s Γ)=1.2, the minimum at m_ s = 0 has become the global minimum,  following a first-order phase transition between the antiferromagnetic and paramagnetic phases. The same behavior is observed for any ω_z < ω_z,  tp.

We have thus found the existence of a first-order QPT on a finite portion of the critical line. We emphasize that this is valid for all strong long-range models, α < 1. This constitutes a qualitative difference with respect to the nearest-neighbor limit, α→∞, of the model, in which the QPT is second-order along the full critical line except for the point where ω_z = 0, when the model is classical <cit.>. This is sketched in Fig. <ref>(b) for comparison.  It is worth noting that there is some evidence that the phase transition could be mixed order <cit.> around the classical point in the nearest-neighbor limit, but not purely first order like we have just described for the strong long-range regime <cit.>.  


 





For α = 0 it is possible to do a classical analysis of the model that predicts the same phase diagram that we just described, with distinct finite regions of first- and second-order phase transitions. See App. <ref> for more details. This illustrates that a classical analysis can serve as an exploratory tool for spin models with all-to-all interactions. It offers an alternative perspective on the phase transition based on the magnetizations of each sublattice rather than the order parameter (the staggered magnetization in this case). In any case, it is not a replacement for an exact solution due to the fact that its application is limited to α = 0 and that it does not allow for the computation of susceptibilities. 



 §.§ Analysis of correlations


Introducing a perturbative field to the Hamiltonian (<ref>), H → H -∑_i=1^N h_i S_i^x, we can compute the susceptibilities as


    χ_ij = 1/βlim_{h_n}→ 0∂^2 ln Z[h_n]/∂ h_j ∂ h_i  .



The susceptibility is proportional to the Kubo correlator <cit.>[Chap. 4], and hence a measure of correlations between spins.
For a translation invariant model the susceptibility can be computed analytically <cit.>[Chap. 6] (see details in App. <ref>). At zero temperature it can be written as


    χ_ij = (A^-1)_ij Y_j  ,


where the matrix is defined by A_ij = δ_ij - 2Y_iJ_ij and Y_i = s ω_z^2 / (8 ε̅_i) with ε̅_i = ε_(-1)^i(u̅) as defined in Eq. (<ref>).









In Figs. <ref>(a) and (b) we show the behavior of the susceptibility matrix, χ_ij, for α=0 in the two phases of the model. Although Eq. (<ref>) has been obtained in the thermodynamic limit (N →∞), evaluating it requires that we fix a finite value of N. We set N=8 and mask out the diagonal elements to better display the structure of the correlation matrix. For α = 0, in the antiferromagnetic phase the correlations show an alternating pattern with intrasublattice correlations being positive and intersublattice correlations being negative. In addition, the intrasublattice correlations of the sublattice that is aligned with the longitudinal field (even sites) are weaker than the intrasublattice correlations of the sublattice that is aligned against the field (odd sites). This effect accentuates with increasing longitudinal field, to the point that the correlations almost vanish for the sublattice aligned with the field for large longitudinal field. The spins of the sublattice that is aligned with the longitudinal field behave increasingly as paramagnetic free spins despite the model still being in an antiferromagnetic state. In the paramagnetic phase the two sublattices become equally magnetized along the combined external field and the correlation matrix recovers a perfectly alternating pattern that matches the staggered interactions, with intra- and intersublattice correlations being of equal magnitude and opposite sign.

Setting α≠ 0 introduces a spatial dependence in the correlations in both phases. The structure of the correlation matrix remains the same but the elements are modulated by the distance between the corresponding spins. 
The correlations exhibit a power-law decay with distance within each of the possible families: intersublattice, even intrasublattice and odd intrasublattice. This is the case regardless of the proximity to the critical point. This behaviour is typical of strong long-range systems and contrasts with weak long-range and short range systems where power-law decay is only present at the critical point <cit.>. Since the model is translation invariant, we can define χ_r, 01≡χ_0  2r+1, χ_r, 00≡χ_0   2r and χ_r, 11≡χ_1   2r+1 for the intersublattice, even intrasublattice and odd intrasublattice correlations; they all follow a power-law decay: χ_r ∝ r^-α_χ with the same exponent, α_χ. 
In the paramagnetic phase 
χ_r, 00 = χ_r, 11. The rate of decay of correlations depends linearly on the rate of decay of interactions, i.e. α_χ = a α + b. Fig. <ref> (d) shows the slope, a, across the phase diagram. The numerical fit of the relation also shows that b ≈ 0 in all cases, which is consistent with the fact that the model must become distance independent for α = 0. Close to the second-order critical line, α_χ becomes independent of α, with a = 0. This is consistent with the behavior described in Ref.  for the strong long-range ferromagnetic Ising model. In contrast, the first-order phase transition is marked by a discontinuity in the slope of the linear dependence between two non-zero values.

The effect of the first- and second-order phase transitions is also apparent in the behavior of single correlation matrix elements across the phase diagram. Figure <ref>(c) shows the value of a correlation matrix element as a function of the longitudinal and transverse fields, for α = 0. They exhibit a divergence at the second-order phase transition and a finite discontinuity at the first-order phase transition. The behavior is analogous for any matrix element and any value of α. 


 






§ NUMERICAL SOLUTION FOR ARBITRARY INTERACTION RANGE



Having established a significant difference between the phase diagram of the strong long-range and nearest-neighbors models, it is natural to wonder how this difference is interpolated for arbitrary values of the range of interactions, α. Since only the strong long-range regime is analytically tractable, in this section we resort to zero-temperature variational quantum Monte Carlo (qMC) <cit.> simulations with a visual transformer (ViT) <cit.> ansatz for finite system sizes to study the model in the full range of interactions. We show the results obtained for spins with s = 1/2 using this architecture, given its recent success in describing spin-1/2 chains with long-range interactions <cit.>. Details on the particular implementation used here can be found in App. <ref>.


Using this technique we obtain the ground state of finite-size chains with periodic boundary conditions along the whole parameter space (ω_x, ω_z, α). In order to study how the phase transition evolves with α, we choose as order parameter the squared staggered magnetization m_s^2 [cf. Eq. (<ref>)]. This choice is justified by the fact that, in finite-size chains, the symmetry of the ground state is not broken and the staggered magnetization is zero throughout the whole parameter space. The squared magnetization, on the contrary, exhibits all the features related to first- and second-order phase transitions and allows us to carry out this numerical characterization.

In Fig. <ref> we observe how ⟨ m_s^2⟩ behaves for different slices of the phase diagram (cf. Fig. <ref>) as a function of the different values of α for a fixed chain size of N = 50. In panels (a), (b) and (c), which correspond to vertical slices (fixed transverse field), we see that for certain values of α the nature of the phase transition evolves from first to second order as the transverse field, ω_z, increases. Panel (a), corresponding to ω_z/(sΓ) = 0.2, merits special attention. It shows that the discontinuity characteristic of first-order phase transitions is clearly present up to α = 2.5 (See App. <ref> for further confirmation). In panels (b) and (c) we see how this transition becomes second order, with a continuous change of the order parameter at the critical point. If instead of comparing a given value of α across panels we focus on any given panel, (a), (b) or (c), we observe that the transition develops a discontinuity as α is lowered. In all cases the transition is clearly discontinuous for α≤ 1, in agreement with the analytical predictions. 




It is interesting to note that analytical results of Sec. <ref> predict that for any α≤ 1, the critical behavior should be universal. In contrast, the numerical results do exhibit a dependence of the critical point and the value of the order parameter on α. We attribute this to finite-size effects, which are expected to be particularly important in long-range systems. The relatively small size considered here allows us to witness dependencies on α that are washed away in the thermodynamic limit. Additionally, we observe that the ansatz encounters difficulties in obtaining the correct ground state near the critical point for α≤ 1. One would expect the critical point to recede toward smaller values of ω_x as α increases. However, for α≤ 1 the states corresponding to the ordered and paramagnetic phases are practically degenerate in energy and the ViT fails to systematically converge to the correct ground state, so this trend can no longer be numerically confirmed.

In Fig. <ref>(d) we show a horizontal slice of the phase diagram (fixed longitudinal field). The phase transition is second order for all values of α, as expected from the analytical results of Sec. <ref>. Following Figs. <ref>(a) and (b) we expected a decrease of the second-order critical point as α increases. This is confirmed by the numerical results.
The critical point tends to the predicted values of ω_z/(sΓ) = 2 for α≤ 1 (strong long-range regime) and ω_z/(sΓ) = 0.5 in the limit of α→∞ (nearest-neighbors regime). We observe that from α≳ 7 the curves collapse and thus it can be considered as the numerical limit from which the model operates in the nearest-neighbors regime. On the other hand, the numerical instabilities now appear at a point where the model is classical (no transverse field).

To better certify the order of the different phase transitions, beyond a visual analysis of the discontinuities (or lack thereof) of the order parameter, we perform a finite-size scaling analysis. In Fig. <ref> we show this analysis for α = 2 (See App. <ref> for other values of α). As it can be seen from panel (a), all the curves corresponding to N = 50, 70, 100 coincide, showing an independence with size characteristic of first-order transitions, in addition to the discontinuity in the order parameter. In contrast, in panel (c), the curves exhibit the typical size dependence of second-order phase transitions. The continuous change of the order parameter becomes increasingly non-analytical as the size increases. The crossover point marks the second-order critical point expected in the thermodynamic limit.





Figure <ref> summarizes the numerical results discussed in this section, along with additional results reserved for App. <ref>. We sketch the critical line for different values of α, highlighting the regions where it is first and second order. We have established that the first-order phase transition at finite transverse field is present for α≤ 2.5. The portion of the critical line that is first order decreases progressively with increasing α until for α≳ 3 the full critical line is second order. This is remarkable because it indicates that the first-order phase transition is present not only in the strong long-range regime but also beyond the mean-field threshold and presumably (see the next paragraph) in the full weak long-range regime. Although the precise boundaries, α_ MF and α^*, between the regime in which the model exhibits mean-field critical exponents and the weak long-range regime and between the weak long-range and short-range regimes, have not been established for the antiferromagnetic model under consideration, it is known that for the analogous ferromagnetic model they lie at α_ MF = 5/3 ≈ 1.66 and α^* = 3 <cit.>.
Additionally, we have witnessed the shrinking of the antiferromagnetic phase as the model approaches the short range regime, driven by a decrease of the critical transverse field at zero longitudinal field from ω_z/(sΓ) = 2 for α > 1 to ω_z/(sΓ) = 0.5 for α≳ 7. 

It is worth noting that the first-order phase transition is very sensitive to finite-size effects. This is evident from the curve corresponding to α = 0.5. We know from our analytical results that the tricritical point where the transition changes from first to second order occurs at ω_z / (sΓ) = 16/(5 √(5)) ≈ 1.43. However, from our numerical results with N=100 we are only able to certify a first-order phase transition up to ω_z / (sΓ) = 0.8. This leads us to believe that our numerical results are also underestimating the region where the critical line is first order for all the other values of α, and therefore also for how low a value of α the first-order phase transition survives. 




§ DISCUSSION



In this paper we have generalized the Ising chain in transverse field to feature tunable-range antiferromagnetic interactions. We have solved the model analytically in the strong long-range regime, showing that it presents a tricritical point in the phase diagram where the critical line changes from first to second order. This is in contrast with the nearest-neighbours limit, where the critical line was known to be always second order but for the point of vanishing transverse field, where the model is classical. To understand the transition from one limit to the other, we have studied the model numerically in the full range of interactions. We have found that the first-order phase transition is present beyond the strong long-range regime, apparently in the full weak long-range regime, only disappearing when the model enters the short-range regime beyond α≈ 3. This confirms that the range of interactions can influence the nature of phase transitions in a model. 

The fact that similar behaviour has been observed in a model with only antiferromagnetic nearest-neighbor and ferromagnetic next-nearest-neighbor interactions <cit.> suggests that the key ingredient that stabilizes a first-order phase transition is the presence of intrasublattice ferromagnetic interactions, in contrast with models with only intersublattice antiferromagnetic interactions. To verify this hypothesis we have also considered Hamiltonian (<ref>) but without intrasublattice ferromagnetic interactions, i.e. with J_ij = 0 for i+j even. This model is not tractable with the analytical method used in Sec. <ref>, so an exact solution in the strong long-range regime is not possible. Nevertheless, a classical analysis of the α = 0 limit, analogous to the one described in App. <ref>, shows that the phase transition is always second order, except for the point of vanishing transverse field, ω_z = 0, where the model is classical. Thus, we can conclude that intrasublattice ferromagnetic interactions (on top of intersublattice antiferromagnetic interactions) are the key ingredient to generate a first-order phase transition. In that sense, in the tunable-range unfrustrated antiferromagnetic Ising model of Hamiltonian (<ref>), α acts as a tuning knob to change the ratio between inter- and intrasublattice interactions. For low enough α, intrasublattice ferromagnetic interactions become large enough to stabilize a first-order phase transition. The same mechanism might explain the first-order phase transitions reported in other quantum long-range models with staggered interactions <cit.>.

Finally, the appearance of a first-order phase transition at zero temperature suggests that the strong long-range unfrustrated antiferromagnetic Ising chain may exhibit ensemble inequivalence at finite temperature. Ensemble inequivalence has been reported to appear in the strong long-range regime in both classical <cit.> and quantum models featuring a first-order phase transition <cit.>.


 






§ ACKNOWLEDGEMENTS

The authors acknowledge funding from the grant
TED2021-131447B-C21 funded by MCIN/AEI/10.13039/501100011033 and the EU
‘NextGenerationEU’/PRTR, grant CEX2023-001286-S financed by MICIU/AEI /10.13039/501100011033, the Gobierno de Aragón (Grant E09-17R Q-MAD), Quantum Spain and the CSIC Quantum
Technologies Platform PTI-001. This research project was made possible through the access granted by the Galician Supercomputing Center (CESGA) to its supercomputing infrastructure. The supercomputer FinisTerrae III and its permanent data storage system have been funded by the Spanish Ministry of Science and Innovation, the Galician Government and the European Regional Development Fund (ERDF). J. R-R acknowledges support from the Ministry of Universities of the Spanish Government through the grant FPU2020-07231. S. R-J.  acknowledges financial support from Gobierno de Aragón through a doctoral fellowship.




§ SOLVING THE STAGGERED ANTIFERROMAGNETIC ISING CHAIN



After diagonalizing the interaction matrix, J, as in Eq. (<ref>), Hamiltonian (<ref>) can be mapped onto a generalized Dicke model <cit.> given by

    H_D = H_0 + ∑_k=0^M-11/D_k a_k^† a_k - ∑_k=0^M-1∑_i=1^N (a_k + a_k^†) λ_ik/√(N) S_i^x   ,

where H_0 = -ω_z∑_i^N S_i^z - ω_x ∑_i^N S_i^x and a_k and a_k^'^†, with [a_k,a_k^'^†] = δ_k,k^', are the creation and annihilation operators of M bosonic modes.
The mapping is valid if lim_N→∞ M/N=0, where M is the number of non-zero eigenvalues of the interaction matrix, J, from the Hamiltonian (<ref>) (the smallest eigenvalue is always set to zero by tuning the parameter b from Eq. (<ref>)).
We have thus replaced the interaction between the spins by an interaction between each spin and a set of M bosonic modes.

Following Wang's procedure <cit.>, the canonical partition function of that generalized Dicke model can be obtained by first computing the partial trace over the matter degrees of freedom.
Then, replacing the photonic degrees of freedom by a collection of M real Gaussian integrals, we obtain

    Z = ∫∏_k=0^M-1√(ND_k/π)du_k  exp{-Nβ f[u_k]}  ,

where the variational free energy per particle, f[u_k], is given by Eq. (<ref>).

The variational free energy has two contributions: one accounting only for the real fields associated with the bosonic modes, {u_k}, and the function F_ m[u_k] [Eq. (<ref>)] which accounts for the spins and their interaction with the bosonic modes. It is computed as

    F_ m[u_k] = -1/βln(Z_ m)   ,

with

    Z_ m[u_k] =  Tr_m[ exp{ -β( H_0 - ∑_i=1^N ∑_k=0^M-1 2λ_ik u_k S_i^x ) }]   .

Z_ m[u_k] represents the partial trace of the generalized Dicke model over its matter degrees of freedom. Using the fact that it factorizes over the N different spins we obtain

    Z_ m[u_k] = ∏_i=1^N ( ∑_m=-s^+sexp{-2mβεu_k}) 
    
        = ∏_i=1^N sinh((2s+1)βε_i)/sinh(βε_i)  ,

where ε_i[u_k] is given by Eq.(<ref>).

Hence, using Wang's procedure the partition function of the original model can be expressed as the integral over a set of M real fields, Eq. (<ref>).
The explicit linear dependence on N of the exponent allows using the saddle-point method, which is exact in the limit N→∞.
Hence, the integral can be replaced by the value of its integrand at its maximum, as done in Eq. (<ref>), which corresponds to finding the global minimum of the variational free energy per particle.





§ FINDING THE GLOBAL MINIMUM OF THE VARIATIONAL FREE ENERGY PER PARTICLE IN ABSENCE OF LONGITUDINAL FIELD



We want to find the global minimum of the variational free energy per particle, f[ u_k ] from Eq. (<ref>), for all α < 1 and ω_x = 0.
Defining for each site on the lattice the variable


    μ_i = ∑_k=0^M-1λ_ik u_k   ,


allows writing Eq, (<ref>) as


    2ε_i[u_k] = √(ω_z^2 + 4 μ_i^2)  .


The relation from Eq. (<ref>) can be inverted to obtain u_k = ∑_i=1^N μ_i λ_ik/N.

The vanishing gradient condition of f[u_k] in its local minima, {u̅_k }, can be written as


    u̅_k/D_k = s/N∑_i=1^N λ_ik∑_l=0^M-1λ_ilu̅_l/ε̅_i B_s(2sβε̅_i )   ,


with B_s(x) the Brillouin function (B_1/2(x) = tanh(x)) <cit.>[Chap. 11] and ε̅_i = ε_i[u̅_k ].
A linear combination of this set of M equations allows us to find a relation satisfied in the stationary points of the variational free energy,

    μ̅_j = s ∑_i=1^N J_ijμ̅_i/ε̅_i B_s (2sβε̅_i)   .

This relation allows us to particularize the expression of f[u̅_k], i.e. the expression of the variational free energy per particle in its stationary points, as

    f[u̅_k] 
       = 1/N∑_j=1^N [ sβμ̅_j^2/ε̅_j B_s(2sβε̅_j) - ln( sinh((2s+1)βε̅_j)/sinh(βε̅_j))   ]   .

Hence, in the stationary points the variational free energy can be written as a sum of a certain function in each site of the lattice

    f[u̅_k] = 1/N∑_i=1^N ξ(μ̅_i^2)   ,

as ε̅_i does not depend on μ̅_i itself but on μ̅_i^2 (see Eq. (<ref>)).

Equation (<ref>) means that among all the possible local minima, the global one must satisfy that |μ̅_i | = μ̅ for each lattice site, in order to minimize each of the N terms ξ(μ̅_i^2).
Due to the structure of the matrix λ_ik, the only possible configuration that meets |μ̅_i | = μ̅ is

    μ̅_i = (-1)^i μ̅  .

It can be translated into a condition for the fields {u_k} inverting the relation from Eq. (<ref>). This condition must be satisfied in the global minimum of the variational free energy per particle,

    u̅_k≠0=μ̅ δ_k0  .

Hence, it analytically shows that the multivariate minimization problem can be turned into a single-variable one when ω_x = 0 for any value of α.
Then, when α=0 the problem is naturally single-variable as M=1 for any value of ω_x.
For any other α<1 it has been numerically checked that the problem remains single-variable, as was analytically proved in the strong long-range ferromagnetic model <cit.>.





§ ANALYTICAL CALCULATION OF THE MAGNETIZATION AND THE SUSCEPTIBILITY




As explained in the main text, introducing a perturbative logitudinal field in Hamiltonian (<ref>) allows us to obtain the magnetization of the model and the susceptibilities between spins. In this appendix we will present the analytical exact expression for both magnitudes in the strong long-range regime.

By introducing a perturbative magnetic field, Hamiltonian (<ref>) is replaced by H → H - ∑_i=1^N h_s S_i^x, which leads to a partition function which depends on the perturbative field: Z[h_n] = Z(h_1,…,h_N), defined by Eq. (<ref>) with the replacement ω_x →ω_x + h_i in Eq. (<ref>).
To compute the perturbative-field-dependent partition function, the variational free energy will depend on {h_n} through Eq. (<ref>), modifying the minimization problem we need to solve.

First, we can define field-dependent magnetization and susceptibilities as


    ⟨ S_i^x ⟩ [h_n] = 1/β∂ln Z[h_n]/∂ h_i  , 
    χ_ij[h_n]  = 1/β∂^2 ln Z[h_n]/∂ h_j ∂ h_i  ,


which allow us to compute the magnetization and the susceptibilities from Hamiltonian (<ref>) as ⟨ S_i^x ⟩ = lim_{h_n}→ 0⟨ S_i^x ⟩ [h_n] and χ_ij = lim_{h_n}→ 0χ_ij [h_n].
For a translation invariant model, as considered, we will be able to compute both magnitudes analytically <cit.>[Chap. 6].

A direct computation of the derivative (<ref>) allows us to write


    ⟨ S_i^x  ⟩[ h_n ] = s B_s(2sβε̅_i ) ω_x + h_i + 2 ∑_k=0^M-1u_k λ_ik/2 ε̅_i .


Note that here the solution to the minimization problem will depend on the perturbative fields, i.e. u̅_k ≡u̅_k[h_n].

Using Eq. (<ref>) the vanishing gradient condition of the variational free energy per particle in its global minimum can be writen as


    1/D_ku̅_k = 1/N∑_i=1^N λ_ik⟨ S_i^x ⟩ [h_n]   .


In the antiferromagnetic model λ_i0=(-1)^i and D_0=Γ, which means that u̅_0/Γ = m̅_ s, proving Eq. (<ref>) in the limit {h_n}→ 0 as a consequence of the null gradient condition, valid for α<1.

Computing the derivative of Eq. (<ref>) leads to


    1/D_k∂u_k/∂ h_j = 1/N∑_i=1^N λ_ikχ_ij[h_n]   ,


where susceptibilities are defined by Eq. (<ref>). Merging both expressions, we can obtain a relation between susceptibilities, which taking the limit {h_n}→0 can be written as


    χ_ij = Y_i ( δ_ij + 2 ∑_r=1^N J_irχ_rj)   ,


with


    Y_i = s/2ε̅_i^2[  ε̅_i B_s(2sβε̅_i) + (ω_x + 2∑_l=0^M-1λ_ilu̅_l )^2 (sβ/2 B_s^' (2sβε̅_i) - 1/4ε̅_iB_s(2sβε̅_i))  ]   .




We can manipulate Eq. (<ref>) to write a matrix equation,


    ∑_r=1^N ( δ_ir - 2Y_iJ_ir) χ_rj = Y_i δ_ij  ,


which allow us to obtain the susceptibilities by inverting the matrix A_ij = δ_ij - 2Y_iJ_ij, writing


    χ_ij = A^-1_ij Y_j   ,


as particularized in Eq. (<ref>) for β→∞.



§ CLASSICAL ANALYSIS OF THE PHASE TRANSITION FOR Α = 0



For α=0 the interaction is homogeneous and alternating in sign and we can express Hamiltonian (<ref>) in terms of total spin operators for each of the sublattices, J_Λ^γ = ∑_i ∈ℒ(Λ) S_i^γ, with Λ = A, B the sublattice index and ℒ(Λ) the corresponding sublattice,


    H = -ω_z (J_A^z + J_B^z) -ω_x (J_A^x + J_B^x ) - Γ/N(J_A^x - J_B^x)^2  .


This Hamiltonian commutes with the total spin operators, J_Λ^2 = (J_Λ^x)^2 + (J_Λ^y)^2 + (J_Λ^z)^2. Consequently, it connects only states with the same total spins J_A^2 and  J_B^2. Using exact diagonalization for finite sizes, we show in Appendix <ref> that the ground state always lies on the subspace of maximum total spins. This implies that the ground state properties of the model are perfectly captured by a model of two spins of sizes, j_A = j_B = sN/2, interacting according to Eq. (<ref>). In the thermodynamic limit, N →∞, the description is further simplified by the fact that these spins can be described exactly in their classical limit <cit.>. We can therefore study the energy resulting from replacing the spins with classical magnetizations J_Λ^γ→ s N/2 m_Λ^γ. The magnetizations are unit vectors that we can assume to lie on the x, z plane m_Λ = (cosθ_Λ, 0, sinθ_Λ). With this, the energy per site reads


    e(θ_A, θ_B)=     -sω_z/2(sinθ_A + sinθ_B) 
       -sω_x/2(cosθ_A + cosθ_B) 
        - s^2Γ/4(cosθ_A - cosθ_B)^2  .


The ground state staggered magnetization is given by m_s = s (cosθ̅_A - cosθ̅_B) / 2, with θ̅_A and θ̅_B the minimizers of E(θ_A, θ_B).

In Fig. <ref> we show the energy landscape as a function of θ_A and θ_B for different values of ω_x and ω_z. For ω_x / (s Γ) ≪ 1 and ω_z / (s Γ) ≪ 2 there are two global minima at (θ̅_A, θ̅_B) = (0, π) and (θ̅_A, θ̅_B) = (π, 0) corresponding to the two symmetric antiferromagnetic configurations with m_s = ± s. If ω_x is kept fixed and ω_z is increased until ω_z / (s Γ) > 2, the two global minima merge into a single one at (θ̅_A, θ̅_B) = (π/2, π/2) corresponding to a paramagnetic configuration with m_s = 0. The coalescence of minima is indicative of a second-order phase transition. Contrarily, if ω_z is kept fixed and ω_x is increased a new local minimum progressively develops at (θ_A, θ_B) = (0, 0) until at ω_x / (s Γ) ≥ 1 it becomes the global minimum, corresponding to a paramagnetic configuration with m_s = 0. The former global minima become local minima. The formation of a new local minimum that grows until becoming the global minimum is indicative of a first-order phase transition. In summary, Fig. <ref> shows that the energy landscape encoded in E(θ_A, θ_B) displays the typical behavior associated with first- and second-order phase transitions.




We find that the phase diagram predicted by this simple classical model coincides with the phase diagram computed exactly and shown in Fig. <ref>. It is therefore an equivalent formulation of the minimization problem presented in Eq. (<ref>), in terms of two variables θ_A, θ_B instead of the single minimization variable u. In this sense, we can understand that Eq. (<ref>) provides the most concise formulation of the minimization problem, as a univariate function of the order parameter, that lends itself particularly well to the analysis in terms of the Landau theory of phase transitions that we performed in Sec. <ref>.



§ VERIFICATION USING EXACT DIAGONALIZATION



For α = 0 there are two equivalent descriptions of the model in terms of either individual spins at each site, as in Hamiltonian (<ref>), or in terms of total spin operators of each sublattice, as in Hamiltonian (<ref>). As explained in App. <ref>, the total spin of each sublattice, J_Λ^2, is a conserved quantity and the Hamiltonian is block diagonal with blocks of fixed total spins. Here we diagonalize the subspace of maximum total spins, j_A = j_B = sN/2 and compare its low-energy spectrum to the low-energy spectrum of the full Hamiltonian (<ref>). In Fig. <ref> we show that the ground-state energies of the two Hamiltonians coincide, with differences appearing in the multiplicity of excited states. We show this for N = 10 and s = 1/2 fixing ω_x / (sΓ) = 0.2 and varying ω_z / (s Γ) and vice versa, but the same behavior is observed for any accessible system size, N, and spin size, s, and for any values of ω_x/ (s Γ) and ω_z/ (s Γ) of the phase diagram. This implies that the low energy physics are well captured by the subspace of maximum total spins, j_A = j_B = sN/2. This fact is exploited in App. <ref> to provide a classical description of the ground state phase diagram. 






§ NUMERICAL PARAMETERS

The visual transformer (ViT), the variational ansatz used to obtain the numerical results discussed in the main text, is a type of neural network composed of several blocks. A detailed block diagram with all the operations involved, as well as the hyperparameters that define the architecture, can be found in Ref. . For this work the only modification that has been made is to add more attention blocks, the so-called core blocks in said reference, which have been fundamental to capture all the long-range correlations presented by the model. Specifically, we have used three of those blocks. The new training procedure is described in the following:

The optimization process is carried out by Stochastic Gradient Descent (SGD) with custom schedules for the learning rate, λ, combined with the Stochastic Reconfiguration (SR)  method <cit.>, characterized by a stabilizing parameter called diagonal shift that we denote as Δ_sr. The training protocol for λ consists of a linear warm-up followed by an exponential decay. This protocol is defined by an initial value of the learning rate, λ_0, a maximum value, λ_max that it attains after n_warm iterations corresponding to the linear warm-up and a ratio γ for the exponential decay. To obtain the results presented throughout this manuscript, a maximum of 500 iterations were used for training, along with the following parameters: λ_0 = 0.1, λ_max = 6.0, n_warm = 150, γ = 0.999. The stabilizer parameter present in the SR method was kept constant with the value of Δ_sr = e-4.




§ ADDITIONAL NUMERICAL RESULTS


Figure <ref> includes further numerical results that reinforce the points discussed in section <ref> of the main text. There we show different vertical slices of the phase diagram for different values of α. These slices allow us to monitor how the tricritical point evolves along the parameter space. These results show that the first-order transition extends to non-negligible transverse fields such as ω_z/(sΓ) = 0.4 beyond the strong long-range regime, for α = 1.5. The first-order phase transition is present up to α = 2.5 and a non-zero transverse field of ω_z/(sΓ) = 0.2. In the cases of α = 0.5 and α = 1.0 we see how the numerical instabilities are alleviated by increasing the transverse field strength.







