Finite-temperature critical behavior of long-range quantum Ising models

We study the phase diagram and critical properties of quantum Ising chains with long-range ferromagnetic interactions decaying in a power-law fashion with exponent $\alpha$, in regimes of direct interest for current trapped ion experiments. Using large-scale path integral Monte Carlo simulations, we investigate both the ground-state and the nonzero-temperature regimes. We identify the phase boundary of the ferromagnetic phase and obtain accurate estimates for the ferromagnetic-paramagnetic transition temperatures. We further determine the critical exponents of the respective transitions. Our results are in agreement with existing predictions for interaction exponents $\alpha>1$ up to small deviations in some critical exponents. We also address the elusive regime $\alpha<1$, where we find that the universality class of both the ground-state and nonzero-temperature transition is consistent with the mean-field limit at $\alpha = 0$. Our work not only contributes to the understanding of the equilibrium properties of long-range interacting quantum Ising models, but can also be important for addressing fundamental dynamical aspects, such as issues concerning the open question of thermalization in such models.


Introduction
Systems featuring long-range interactions are central in condensed matter and statistical physics, due to both their widespread presence in nature and the wide range of characteristic physical phenomena they display, the latter often being at odds with well-known predictions and results concerning short-range models (see, e.g, [1] for a review). Within the last decade, the interest in quantum long-range interacting models has further surged due to the progress in manipulating and controlling these systems at an unprecedented level [2][3][4][5][6]. Specifically, these experimental platforms naturally realize long-range quantum Ising or Heisenberg models, with the possibility to engineer many-body interaction potentials decaying proportionally to d −α as a function of distance d, ranging from van-der-Waals-like (α = 6) and dipolar interactions (α = 3) in the context of Rydberg atoms [3,6], to Coulomb (α = 1) and infinite-range (α = 0) potentials for trapped ions [2,5].
Recent experiments in such long-range interacting models have mostly centered on the investigation of inherent dynamical phenomena, such as many-body localization [7], discrete time crystals [8,9], prethermalization [10], Kibble-Zurek mechanism [11,12], or dynamical quantum phase transitions [13,14]. Despite of recent progress [15,16] one key question has, however, remained open: especially in the limit of small interaction exponents, it is not known whether these long-range systems follow the fundamental principle of thermalization as expected for generic short-range models. In the first place, this obviously requires a thorough understanding of the thermal properties of the system of interest, which have only been partially explored even in paradigmatic Hamiltonians such as the one-dimensional long-range quantum Ising model.
In particular, the ground-state properties of the latter in the case of ferromagnetic (FM) interactions have been the focus of investigation via analytical and renormalization group (RG) techniques [17][18][19], as well as linked-cluster expansions [20], tensor network approaches and/or density matrix RG [21,22], Monte Carlo methods [23] and, very recently, Stochastic Series Expansion (SSE) Monte Carlo [24] investigation in the α > 1 region, demonstrating, e.g., that the critical behavior of the model belongs to the mean-field and short-range universality class (UC) for 1 < α < 5/3 and α ≥ 3, respectively. The antiferromagnetic case has also been intensely studied via the use of several approaches [20,[24][25][26][27], with notable results including, among others, the demonstration that the half-chain entanglement entropy displays area-law violations in the intermediate regime 1 < α < 2 [25]. Similarly, p-wave superconductors with long-range pairing [28] have been shown to display exotic critical behavior, even if, due to the presence of Jordan-Wigner strings, those models do in general differ from Ising chains with similarly decaying interactions. Considerable effort has also been dedicated to the theoretical investigation of the dynamical properties of this type of model [29][30][31][32][33][34][35].
Oppositely with respect to the zero-temperature case, the finite-temperature regime is still poorly understood. Indeed, the latter has been predicted by general theoretical arguments [36] to belong to the universality class of the corresponding classical long-range Ising model, with quantum effects not changing this description at the qualitative level. While this picture has been essentially confirmed for the case α = 2 by SSE studies [37], the latter demonstrated, in the proximity of the ground-state critical point, the presence of considerable finite-size effects induced by strong quantum fluctuations, which all but prevent observation of the expected classical regime even at very large system sizes.
In the light of the experimental realizations of these models discussed above, investigating the thermal critical behavior of these Hamiltonians remains therefore of great importance, in order to determine the role and strength of the quantum effects in perturbing the predicted classical picture. Furthermore, (numerically) exact analysis of the finite-temperature regime is essential to determine non-universal details such as, e.g., the position of thermal critical points, which are influenced in a key way by quantum effects, and whose knowledge is crucial for laboratory realizations. Such a study is of especially great interest in the extremely longranged regime 0 < α < 1, which, to our knowledge, has not been the object of this kind of investigation, and (as mentioned above) is directly realizable in trapped-ions setups.
In this work, we study both the ground-state and finite-temperature phase diagram of the long-range FM quantum Ising model in one spatial dimension, by means of numerically exact, large-scale Path Integral Monte Carlo simulations. We perform our calculations for two representative values of α: namely, we choose α = 0.05 and α = 1.50, within the extremely long-range region α < 1 and intermediate region 1 < α < 2, respectively. We employ a wide variety of well-known finite-size scaling techniques to determine the position (i.e., the critical points) and critical exponents of both the ground-state and finite-temperature paramagnetic (PM)-FM transitions displayed by the model, obtaining the phase diagram displayed in Fig. 1.
We determine the critical points and critical exponents for the ground-state FM-PM transition. Our results for critical point positions and correlation length critical exponents are in agreement with existing predictions in the literature where the latter are available (i.e., α = 1.50), while we encounter relatively small (∼ 7%) deviations with respect to our estimate for the magnetization critical exponent. We then obtain accurate results for the position of the critical points in the finite-temperature regime for several values of the interaction strength. Concomitantly, our estimated correlation length critical exponents at α = 1.50 essentially confirm the theoretical prediction of no qualitative deviations from the classical universality class due to quantum fluctuations, while discrepancies (up to 10% in the strongly interacting region) appear in the susceptibility critical exponent.
The structure of the paper is the following. Sec. 2 introduces the Hamiltonian, the numerical technique employed for its study, and the finite-size scaling approaches we employed to analyze its critical behavior. Sec. 3 discusses our obtained results on the critical behavior of the model. Finally, in Sec. 4 we outline the conclusions of our work and offer an outlook for future direction of research.   (1), displaying the ground-state and finite-temperature phase boundary and critical exponents obtained using finite-size scaling techniques. Panels (a) and (b) correspond to α = 0.05 and α = 1.50, respectively. Here, T is the system temperature in units of the Boltzmann constant, and V is the interaction strength in units of the transverse field (see below). The displayed results for the effective thermal exponent and its product with the magnetization and susceptibility critical exponent are those obtained via data collapse (see below).

Hamiltonian and known results
The model analyzed in this work is described by the Hamiltonian where V > 0 is the interaction strength, i, j run over the sites 1, . . . , L of a one-dimensional lattice with periodic boundary conditions, r ij is the distance between sites i and j, S z i (S x i ) is the component along z (x) of the spin-1/2 operator acting on site i, and K(L) ≡ (L − 1) −1 i =j r −α ij is the Kać renormalization factor. The latter ensures the existence of a proper thermodynamic limit in the regime α ≤ 1, while for α > 1 it amounts to a rescaling of the interaction strength, and does not change the universal features of the critical behavior of the model. We remark that the presence of this renormalization factor is directly related to how interactions with α < 3 are engineered in trapped ions experiments. The latter exploit coupling between the ions and collective modes of the ion chain (phonons), mediated via a single laser shined over the full sample. Increasing the number of ions while keeping the lattice spacing constant naturally leads to a reduced coupling strength, that translates into the fact that the energy of the full system is still extensive -as reflected by Kać normalization. In the following, periodic boundary conditions are taken into account following the minimum-image convention, and h = 1 will be taken as unit of energy.
For very small interaction strength V , the ground state of the system in the thermodynamic limit is a paramagnet, characterized by a vanishing value of the magnetization along the z direction |m z | ≡ L −1 | i S z i |. On the contrary, for V ≫ 1 the system is in a FM phase, displaying a finite |m z |. The existence of a finite-V phase transition connecting these two states can be proven via analytical arguments (see, e.g., [17]); its UC depends strongly on the value of the decay parameter α. Indeed, the α = 0 case, also referred to as Lipkin-Meshkov-Glick model [38], can be described in an exact fashion at the mean-field level [39], and the PM-FM transition has been proven to belong to the Gaussian UC in the 1 < α < 5/3 region. In contrast, in the regime α ≥ 3, the critical point belongs to the short-range UC (i.e., the one of the FM-PM transition in the nearest-neighbor limit α → ∞).
In the finite-temperature regime, generic scaling arguments [36] predict that the model should display the same critical behavior as its classical (i.e., h = 0) counterpart, due to the finiteness of the system size in the imaginary time dimension (see below). The critical behavior of the classical model has been studied via both analytical (see, e.g., [40]), RG (see, e.g., [41]) and numerical techniques (see, e.g., [42]) in the α > 1 regime. Here, the system displays a second-order FM-PM thermal phase transition for 1 < α < 2, with the region 1 < α < 3/2 belonging to the mean-field regime, while in the point α = 2 the model undergoes a finite-temperature transition of the BKT type, and the short-range regime is reached (i.e., no finite-temperature transition takes place) for α > 2.

Numerical techniques and finite-size scaling
We perform our investigation of the Hamiltonian in eq. (1) via Path Integral Monte Carlo (PIMC) [43], a numerically exact technique for the study of unfrustrated systems of bosons and quantum spins. In this approach, one maps the features of a quantum model of interest to those of an equivalent, higher-dimensional classical one, which is then studied via Metropolis Monte Carlo simulations. The quantum-to-classical mapping described above maps the partition function of the extended transverse-field Ising model in eq. (1) into the one of an anisotropic extended Ising model on a rectangular lattice, via a procedure known as Suzuki-Trotter breakup. Here, in addition to the original spatial dimension, one also considers a discretized and periodic one, known as imaginary time, which extends in the interval [0, β], where β = 1/T is the inverse system temperature in units of the Boltzmann constant. The number of sites M along this direction (also known as slices) is a free parameter which affects the accuracy of the mapping: indeed, the latter is exact up to O(β/M ) corrections, which vanish in the limit M → ∞.
In the spatial direction, the extended Ising model resulting from the mapping displays the same FM long-range interactions present in the spin-spin term of the model in eq. (1), while spin-spin couplings are nearest-neighbor in the imaginary time direction. Our PIMC algorithm combines conventional Wolff cluster updates [44] in imaginary time with efficient long-range cluster updates [42] in the spatial direction. The choice of these two state-of-theart techniques allow to accurately analyze large system sizes (up to L = 8192 sites) at low enough temperatures (down to β = 1024) to reach the ground state regime. The Suzuki-Trotter corrections mentioned above are kept into account by performing simulations with increasing number of slices (up to M = 65536), until a value M = M * is found such that the corresponding values of the observables of interest were determined to be identical, within statistical error, to those obtained for M = 2M * . The same protocol (with β in the place of M ) is adopted to ensure the T → 0 limit is reached in the investigation of the ground state regime.
The PIMC algorithm gives us direct access to observables commuting with the S z i operators, including the integer powers of |m z |. This allows us to compute quantities such as the Binder cumulant where . . . stands for statistical averaging, which is expected to converge to 1 (0) in a FM (PM) phase [45]. We also compute the "classical" susceptibility which, in proximity of a finite-temperature critical point of a quantum model, approximates well the exact functional form of the magnetic susceptibility [37].
In order to extract reliable information on the critical behavior of the model in the thermodynamic limit, we exploit the well known finite-size scaling (FSS) theory [45]. In this framework, scaling relations of various quantities in terms of the correlation length ξ, which diverges when approaching a critical point, are exploited to obtain finite-size information by noting that in a finite system ξ will saturate to a value O(L), where L is the system size. Features such as the position of the critical point or the critical exponents, on which the original scaling relations depended, can then be directly extracted via numerical fits as a function of L. In the following section, when discussing the fitting procedures to obtain such quantities, we will offer detailed formulae regarding FSS predictions for observables such as U and χ.

Results
We investigate the critical properties of the model in eq. (1) in the ground-state and finitetemperature regime for α = 0.05 and α = 1.50.

Ground-state critical behavior
The first step in our analysis is the determination of the PM-FM critical point V c in the ground-state regime, which we accomplish by fitting to our numerical data for the Binder cumulant U its expected FSS behavior. The Binder cumulant curves U (V ) for system sizes L and, e.g., 2L are expected to cross at size-dependent points V = V U (L), which will follow (to the leading order) the FSS scaling [24,46] where V c is the critical point, and the effective thermal exponent θ t is linked to the correlation length critical exponent ν.
In the ground-state regime ν −1 = θ t outside of the mean-field region; conversely, when the latter is entered, corrections to the leading scaling behavior can be taken into account [24] via the generalized expression ν −1 = (d uc (α)/d) θ t , where d is the dimensionality and d uc (α) = 3(α − 1)/2 is the upper critical dimension for the value of α of interest.
Comparison of eq. (4) with the predicted leading-order FSS behavior for the value of the Binder cumulant at the V U (L)s, allows us to obtain estimates for V c and θ t , by fitting our computed results for the crossing features [see Fig. 2(a)] with the functional forms above. Fig. 2(b-c) display examples of the FSS fitting procedures mentioned above; the obtained values of the critical point and of the effective thermal exponent θ t are listed in Table 1. In order to gain more insight into the ground-state critical behavior of the model, we perform a data collapse analysis by directly exploiting the FSS predictions for the behavior of the squared magnetization close to a critical point [24,45], where β m is the magnetization critical exponent, up to corrections of higher order in 1/L. This scaling law implies that the rescaled magnetization curves y m L ≡ m 2 z (L)L +2βmθt for different system sizes should coincide if plotted as a function of for a wide range of candidate values of V c , θ t and β m , choosing as our final estimates for these quantities the values which resulted in the fit with the lowest chi-square value. While it is hard to assign a rigorous error bar to the results of a data collapse analysis, we estimate the order of magnitude of the error on our results by performing the same fits in a considerably larger (i.e, containing of the order of double the number of points) window around the critical point, and taking the difference between the optimal values of V c , θ t , and β m for the two windows as the order of their numerical uncertainty.
Our collapsed data is displayed in Fig. 3(a-b); the obtained estimates for V c , θ t and β m are listed in Table 1. We note that the data collapse behavior takes place over a fairly wide range of values of the rescaled order parameter x V L , despite relatively narrow fitting windows for the scaling behavior in eq. (6) (the intervals between dashed lines in Fig. 3). This highlights the faithfulness of the data collapse scaling description of our numerical data, which translates to highly reliable estimates of the critical properties of the system. Examination of our results points out i) the remarkable agreement of the critical point estimates obtained via the Binder cumulant FSS and the data collapse, and ii) conversely, the incompatibility between the two estimates for the effective thermal exponent θ t . Due to the arguments mentioned above, we believe the data collapse estimates for the critical features to be more reliable in this regard.
For α = 1.50, we find agreement for θ t and deviations of the order of 7% for 2β m θ t from the independent SSE predictions in Ref. [24] which, in our notation, are θ t ≃ 2β m θ t ≃ 0.667. We also find good agreement with the estimate V c ≃ 0.42 (in our notation) given in [24] for the position of the ground-state critical point, by performing a data collapse where the rescaled

Finite-temperature critical behavior
Once the boundary of the ground-state FM phase is determined, we investigate whether or not FM order survives for T > 0, and more in general the details of the critical behavior of the model in this regime. To this end, we perform finite-temperature calculations for fixed values of V belonging to the FM phase in the ground state regime. We apply the FSS framework to quantities such as the Binder cumulant and the susceptibility, computed as a function of T , to estimate features of the temperature-driven critical behavior. Indeed, our results for the Binder cumulant as a function of β at fixed V and different system sizes immediately confirm the presence of a finite-temperature phase transition, as pointed out by the appearance of the crossing behavior discussed above [see Fig. 4(a)] at sizedependent points β U (L, V ). We determine the V -dependent critical temperatures β c (V ) and the associated θ t (V ) via fitting of the FSS relations in eqs. (4)-(5) to our computed crossing features, with the thermal critical points β c and β taking the role of V c and V , respectively. If the hypothesis of essentially classical critical behavior for the finite-temperature quantum model holds (as we argue below) one may link [47]   Examples of this analysis are displayed in Fig. 4(b-c): the obtained critical parameters are listed in Table 2. We remark here that our application of this approach encountered in some cases strong difficulties due to significant finite-size effects in proximity of the β c (V, L). In particular, the relatively large numerical uncertainties on the values of the Binder cumulant in this region led to the necessity to perform conservative estimates of the finite-size crossing points. In turn, this prevented us in some cases from obtaining meaningful (i.e., with small enough error bars) estimates for θ t .
In order to obtain an independent estimation of our quantities of interest, we investigate the finite-temperature behavior of the magnetic susceptibility for the same values of V selected in our Binder cumulant analysis. At finite system size and fixed interaction strength, χ is expected to display peaks at size-dependent temperatures β χ (L, V ); the FSS framework predicts for the latter [24,45] the leading scaling behavior as a function of the system size.
Our numerical data confirm the expected behavior of χ [see Fig. 5(a)]. Fitting the FSS functional form in eq. (7) to the computed peak positions [see Fig. 5(b) for an example] allows us to directly estimate the critical temperatures and effective thermal exponents as a function of the interaction strength (see Table 2 for a list of results).
While also requiring conservative estimates (and therefore large error bars) for the peak positions, due to strong finite-size effects, we found the susceptibility-based approach to be much less sensitive to this issue than the Binder cumulant FSS discussed above. In particular, we encountered problematic results only for V = 2.5, for both values of α considered in this work, where our estimates were strongly dependent on the set of system sizes considered in the fitting procedure (the reported results correspond to the fits with all sizes considered). We finally analyze the critical properties of the model by performing a data collapse analysis for the behavior of the magnetic susceptibility close to the finite-temperature critical points [24,42,45], where γ is the susceptibility critical exponent, up to corrections of higher order in 1/L. The analysis follows the same protocol outlined in our discussion of the ground-state regime, with the rescaled dependent and independent variables here being y χ L ≡ χ(L)L −γθt and  Table 2: Summary of the computed estimates for β c , θ t , and γθ t (see text) for the finitetemperature transitions at our investigated values of α and V . Our results are categorized according to the methodology employed to derive them: namely, FSS of the Binder cumulant crossings (U ), FSS of the magnetic susceptibility peak position (χ), and data collapse of the susceptibility (χ dc ). Estimates marked with an asterisk ( * ) did not converge with respect to the choice of minimum size to be included in the fitting procedure.
(β c − β) L θt , respectively. Fig. 6 displays our collapsed data for all the values of α and V investigated in this work; the corresponding optimal (in the sense discussed above) results for β c , θ t and γ are displayed in Table 2. As in the ground-state regime, we observe that the parameter range in which the data collapse scaling ansatz is respected noticeably exceeds our fitting window (and vastly so, in most cases), highlighting the accuracy of this approach in describing the critical behavior of the model. Furthermore, this protocol does not require the estimation of size-dependent features, sush as the curve crossings for the Binder cumulant, or the peak position for the susceptibility, allowing us to obtain much more reliable and systematics-free results. We also note that high degree of accuracy with which the scaling law in eq. (7) can be applied to describe the behavior of the "classical" susceptibility in eq. 3 is a strong indication of the goodness of the latter as an approximation for the complete functional form of the magnetic susceptibility.
A direct analysis of the results for the critical exponents listed in Table 2 shows that our estimates obtained via FSS of the Binder cumulant crossings, where meaningful in the sense discussed above, are consistent within error bar with the ones obtained via susceptibility data collapse. Concomitantly, in some points we observe differences (which remain consistently small, except for the point α = 1.50, V = 5.00) between the latter and the results of the susceptibility peak position FSS for the values of V in which the latter have converged with respect to the system sizes employed in the fitting procedure. In the points where this did not happen, the θ t result from the susceptibility peak position fit decreased, shifting towards the data-collapse results, when smaller sizes were discarded.
According to the arguments mentioned in Sec. 2, the universality class of the T > 0 FM-PM transition should be the same of the corresponding transition in the classical counterpart of model eq. (1). For α = 1.50, the classical Hamiltonian is in the mean-field regime, and RG predictions, confirmed by classical Monte Carlo calculations [42], yield the estimates θ t = γθ t = 1/2. Direct comparison with our most representative and reliable results in Table 2 (i.e., the one obtained via data collapse of the magnetic susceptibility) shows that our estimates for θ t are in essential agreement with the classical prediction (with deviations outside of the estimated order of magnitude of the error only appearing for V = 5.0). Compatibility between our estimate and the theoretical predictions, even for V = 5.0, is confirmed by the results obtained via FSS of the Binder cumulant, while the susceptibility FSS estimates, where converged, show appreciable deviations only for V = 5.0. Conversely, our estimates for γθ t show relatively consistent deviations (up to the order of 10%), which increase with the interaction strength.
These differences with the predicted results may be in principle due to several causes, including i) the "classical" approximation employed for the study of the susceptibility in our analysis, or ii) genuine quantum effects which introduce deviations with respect to the predicted classical behavior. However, we find it unlikely that either (i) and/or (ii) may be the dominant physical mechanism underlying the observed deviations, since both effects are essentially quantum in nature, and are expected to become weaker for larger values of V , where in contrast our results are more at odds with the classically predicted values. Indeed, for higher interaction strengths quantum effects are expected to weaken, due to both the larger value of V (in comparison to the transverse field h) and the higher temperature at which the critical region is located. This consideration leads us to the conclusion that despite these deviations (which may be caused by finite-size effects, or by higher-order corrections) the critical behavior of the model in this regime follows the classical UC.
As in the ground-state case, we find essential compatibility with the (classical) mean-field exponents at α = 0; in particular, we match the predicted values [39] θ t = γθ t = 1/2 up to relatively small deviations (of up to 2.5%) for the latter quantity, which also become larger in the strongly interacting regime, and are therefore likely not due to genuine quantum effects as argued above.

Conclusions and outlook
We study the ground-state and finite-temperature phase diagram and critical behavior of the long-range quantum Ising model in one spatial dimension, for values of the interaction exponent parameter of direct interest for current experiments in trapped ion setups. We perform numerically exact, large-scale PIMC simulations within both the extremely longrange region and intermediate long-range regime, respectively, employing a wide variety of finite-size scaling techniques to determine the location (i.e., the critical points) and critical exponents of both the ground-state and finite-temperature phase transitions displayed by the model.
We determine transition points and critical exponents for the ground-state FM-PM transition. We find essential agreement with existing predictions for these quantities, where available (up to small deviations for the value of the magnetization critical exponent), and compatibility of our extremely-long-range results with the fully-connected universal properties. We then accurately estimate the position of the critical points in the finite-temperature regime for several values of the interaction strength. Here, our estimated critical exponents in the intermediatelong-range region essentially confirm the theoretical prediction of classical universality. In particular, in the intermediate long-range regime our estimated correlation length critical exponent is fully consistent with the classical predictions, while the susceptibility exponent displays deviations at most up to the order of 10%. Similarly, in the extremely long-range region we find compatibility with the (classical) mean-field universality class up to deviations of the order of 2.5% in the value of the correlation length critical exponent. For future works, it would be interesting to verify if some of these findings also apply to long-range p-wave superconductors [28], that, while described by free theories, could still display some of the phenomenology we discuss.
Beyond exploring the equilibrium phase diagram and the nature of critical points, our work is also directly relevant for another open question appearing in the context of quantum Hamiltonians with long-ranged interactions. This concerns quantum thermalization and equilibration during coherent quantum dynamics without coupling to an environment, which appears all but settled. In the infinitely-connected limit of α = 0 it is already well known that thermalization does not occur [48,49]. Furthermore, numerical works close to this infinitelyconnected limit have already observed indications that thermalization could be prevented at least on the achievable time scales [50,51]. In order to settle this fundamental question, the understanding of the thermal equilibrium phases and properties, to which this work contributes, represents a first key step. While thermalization corresponds to ensemble equivalence of the thermal ensemble with the diagonal ensemble, capturing the long-time steady states during dynamics [52], it is also not known to which extent such long-range models exhibit ensemble equivalence on a general level. This concerns for instance the equivalence of the thermal and microcanonical ensemble, which is of central importance from the statistical physics point of view.