Phase diagram and quantum criticality of Heisenberg spin chains with Ising-like interchain couplings -- Implication to YbAlO$_3$

Motivated by recent progress on field-induced phase transitions in quasi-one-dimensional quantum antiferromagnets, we study the phase diagram of $S=1/2$ antiferromagnetic Heisenberg chains with Ising anisotropic interchain couplings under a longitudinal magnetic field via large-scale quantum Monte Carlo simulations. The interchain interactions is shown to enhance longitudinal spin correlations to stabilize an incommensurate longitudinal spin density wave order at low temperatures. With increasing field the ground state changes to a canted antiferromagnetic order until the magnetization fully saturates above a quantum critical point controlled by the $(3+2)$D XY universality. Increasing temperature in the quantum critical regime the system experiences a fascinating dimension crossover to a universal Tomonaga-Luttinger liquid. The calculated NMR relaxation rate $1/T_1$ indicates this Luttinger liquid behavior survives a broad field and temperature regime. Our results determine the global phase diagram and quantitative features of quantum criticality of a general model for quasi-one-dimensional spin chain compounds, and thus lay down a concrete ground to the study on these materials.

Motivated by recent progress on field-induced phase transitions in quasi-one-dimensional quantum antiferromagnets, we study the phase diagram of S = 1/2 antiferromagnetic Heisenberg chains with Ising anisotropic interchain couplings under a longitudinal magnetic field via large-scale quantum Monte Carlo simulations. The interchain interactions is shown to enhance longitudinal spin correlations to stabilize an incommensurate longitudinal spin density wave order at low temperatures. With increasing field the ground state changes to a canted antiferromagnetic order until the magnetization fully saturates above a quantum critical point controlled by the (3 + 2)D XY universality. Increasing temperature in the quantum critical regime the system experiences a fascinating dimension crossover to a universal Tomonaga-Luttinger liquid. The calculated NMR relaxation rate 1/T1 indicates this Luttinger liquid behavior survives a broad field and temperature regime. Our results determine the global phase diagram and quantitative features of quantum criticality of a general model for quasi-onedimensional spin chain compounds, and thus lay down a concrete ground to the study on these materials.
Introduction. In low-dimensional correlated electron systems strong quantum fluctuations give rise to quantum phase transitions (QPTs) [1] and a number of exotic quantum phenomena, such as unconventional superconductivity [2,3], non-Fermi liquid behavior [4,5], and quantum spin liquids [6]. In the past decade, tremendous progresses have been made in understanding the nature of QPTs and associated emerging phenomena in quasi-one-dimensional (Q1D) antiferromagnets. These include the E 8 symmetry [7][8][9], manybody string excitations [10,11] and novel quantun criticality [12,13] in transverse field Ising chains, and Bose-Einstein condensation (BEC) and glassy phases in coupled antiferromagnetic (AFM) chains [14,15]. As a paradigmatic model for 1D quantum antiferromagnets, the S = 1/2 Heisenberg chain is well described by a Tomonaga-Luttinger liquid (TLL), where both the longitudinal and transverse spin correlation functions follow algebraic decay. [16] Under a magnetic field, the staggered transverse correlations are always dominant over the longitudinal ones, and a canted AFM order with staggered transverse correlations (denoted as the TAF order) is stabilized when interchain couplings become relevant. In systems with an Ising anisotropy, besides the TAF phase which arises from a spinflop mechanism [17], the peculiar quantum fluctuations in the Ising anisotropic XXZ chain give rise to incommensurate modulation of the longitudinal spin correlations [18] and can stabilize an incommensurate longitudinal spin density wave (LSDW) order [19]. This LSDW state has been recently observed in several Q1D antiferromagnets [20-22, 24, 25].
Recent inelastic neutron scattering (INS) measurements reveal quantum critical TLL behavior of a coupled S = 1/2 chain compound YbAlO 3 with nearly isotropic (Heisenberg) intrachain exchange couplings [26]. A surprising observation is an incommensurate AFM state induced by the applied magnetic field. In this phase, the modulation of the ordering wave vector is proportional to the magnetization, which is a charac-teristic of the LSDW order in coupled Ising anisotropic XXZ chains. This leads to the puzzle on the origin of the incommensurate AFM order in Heisenberg chains. A clue from both experimental observation and theoretical analysis is the relevance of the Ising anisotropic interchain coupling [26,27].
Still a generic open question is that how the interchain Ising anisotropy would affect the phase diagram and low-energy excitations of Heisenberg chains. This poses a major challenge to existed theories based on the interchain mean-field approximation where the interchain fluctuations are neglected [19].
To tackle these issues, in this letter we study the fieldinduced phase diagram of S = 1/2 Heisenberg chains with Ising anisotropic interchain couplings by using large-scale quantum Monte Carlo (QMC) simulations. Our results unambiguously show that the interchain interactions enhance longitudinal spin correlations to stabilize an incommensurate LSDW. With increasing field, the ground state transforms to a TAF state, and is fully polarized above a quantum critical point (QCP) controlled by the (3+2)D XY universality. Increasing temperature from the QCP, the scaling of thermal energy and NMR relaxation rate demonstrate that the system experiences a clear dimension crossover to the universal TLL behavior, exhibiting rich physics and fine structure of the quantum criticality. We then propose NMR measurements as a means to probe the ground states and related low-energy excitations in YbAlO 3 and other Q1D antiferromagnets.
Model and method. We consider a model defined on a three-dimensional (3D) cubic lattice for the S = 1/2 Heisenberg spin chains with weak interchain couplings of the XXZ type under a longitudinal (z-direction) magnetic field. The Hamiltonian reads as Here S i = {S x i , S y i , S z i } is an S = 1/2 spin operator defined on site i. J c and J ab are respectively the intrachain (along c axis) and interchain exchange couplings between the nearest neighbor spins. ε is a parameter characterizing the spin anisotropy of the interchain coupling. g is the gyromagnetic factor, µ B is the Bohr magneton, and H is the applied magnetic field. We take J c as the energy unit and define the reduced temperature t = T /J c and reduced field h = gµ B H/J c . For Ising anisotropic interchain interaction, ε < 1. Here we take ε = 0.25 and J ab = 0.2J c for demonstration. The effects of ε and J ab on the phase diagram of the system will be discussed later. To study the model in Eq. (1) we perform numerically exact quantum Monte Carlo (QMC) simulations based on the stochastic series expansion (SSE) algorithm [28,29]. In the simulations, the largest system size is 32 × 32 × 256 and the lowest temperature accessed is t = 0.003.
Phase diagram and the LSDW phase. Our main results are summarized in the phase diagram of Fig. 1(d). At low temperatures, three ordered phases appear sequently with increasing field, and the spin patterns along a chain in these phases are illustrated in Fig. 1(a). An Ising AFM phase with ordered moments aligned in the z direction is stabilized for h < h 1 ≈ 0.6, and a LSDW state with incommensurate longitudinal spin correlations is stabilized for field regime h 1 < h < h 2 ≈ 0.89. For h > h 2 the ground state becomes a TAF, which is a canted To examine the nature of the ordered states, we calculate the normalized longitudinal and transverse spin structure factors The Ising AFM order is signaled by a peak of S zz (q) at q = (π, π, π). When h > h 1 we find that the peak splits into two located at incommensurate q = (π, π, π ± ∆Q) ( Fig. 1(b)). In this incommensurate AFM phase the ordering wave vector varies with increasing field, satisfying |∆Q| = 2π m z ( Fig. 1(c)), a characteristic reflecting the Q1D TLL physics of the LSDW state [19]. This confirms that the incommensurate order is indeed a LSDW, which in this model arises from the enhancement of longitudinal correlations by interchain Ising anisotropy. For h > h 2 , the peak of S zz is suppressed and the ground state changes to the TAF with a peak of S xy (q) at q = (π, π, π), as shown in Fig.S2 of Supplemental Material (SM) [30]. Quantum criticality. The QPT at h c takes place when the TAF order is suppressed. Since the TAF order breaks the spin U (1) symmetry, the transition can be viewed as a magnetic BEC [14,33] with a dynamical exponent z = 2. The QCP then belongs to the (3 + 2)D XY universality class. To check this, we first study the scaling behavior of critical field h c (t) at low temperatures. As shown in Fig. 2, h c (t) data determined from either the peak of field dependent susceptibility χ zz (h) = ∂m z /∂h or the correlation length ξ ∼ L ( Fig.S3 and Fig.S4 [30]) follow the scaling relation of 3d BEC, h c − h c (t) ∼ t d/2 = t 3/2 . We then study the finitetemperature crossover in the vicinity of the QCP. At each field the temperature dependent susceptibility χ zz (t) develops a broad peak, and the peak position defines the crossover temperature T cr (Fig.S4 [30]). Near a QCP, T cr ∼ |h − h c | νz , where ν is the correlation length exponent. From Fig. 2(a) we find that T cr ∼ |h − h c | on both sides of the QCP, consistent with the (3 + 2)D XY universality z = 2 and ν = 1/2.
Owing to its Q1D structure, the system exhibits finite temperature 1D-3D crossover. This is clearly shown in the scaling of thermal energy, E = H /N , right at the critical field h c . Since dE/dt = C ∼ t d/z , where C is the specific heat, we expect E ∼ t φ E with φ E = d/z + 1 = 5/2. But as shown in Fig. 3(a), this scaling fits only at low temperatures for t 0.06. And for t 0.2, φ E ≈ 3/2, implying an effective dimension d ef f = 1. Careful windowing analysis [34] in Fig. 3(b) finds a gradual increase of φ E from about 3/2 to 5/2 for 0.06 t 0.2, clearly indicating a 1D-3D crossover in this temperature regime. The 1D-3D crossover gives rise to rich quantum scaling behaviors. For example, the genuine 3D nature of the QCP is inherent in the low-temperature scaling of susceptibility data in the disordered phase (Fig. 3(d)), which satisfies χ zz ∼ |h − h c | ν(d+z)−2 X t |h−hc| νz with d = 3, ν = 1/2, and z = 2. In the quantum critical regime it is ex- with the same exponents at low temperatures. But the scaling of QMC data above the dimension crossover temperature in Fig. 3(c) are consistent with d ef f = 1, ν = 1/2, and z = 2, characterizing a quantum critical TLL behavior.
NMR relaxation rate 1/T 1 and the TLL behavior. In the TLL regime of an XXZ chain the spin correlations decay algebraically as S z 0 S z r − (m z ) 2 ∼ cos(2k F r)r −1/η and S x 0 S x r ∼ (−1) r r −η , where k F = π(1/2 − m z ), denoting the Fermi wave number of pseudo-fermions mapped from the spin model by a Jordan-Wigner transformation, and the Luttinger exponent η determines the decay rate. For a Heisenberg hinv is the crossover field where the dominant fluctuation changes from longitudinal (η > 1) to transverse (η < 1) with increasing field. It is found hinv ≈ h2, which separates the LSDW and TAF ground states. chain, η < 1 for all fields and the staggered transverse fluctuations are always dominant. On the other hand, for an Ising anisotropic XXZ chain, there is an η inversion at field h inv , e.g. η > 1 for h < h inv where longitudinal fluctuations dominate and η < 1 for h > h inv where the dominant fluctuations turn to transverse ones [19].
To examine whether the TLL behavior near h c extends to lower fields and to determine the dominant spin fluctuations associated with magnetic ordering we calculate the NMR spin-lattice relaxation rate 1/T 1 , which probes the low-energy spin fluctuations of a magnetic system [31,32]. For simplicity we set the nuclear gyromagnetic ratio and the hyperfine coupling to be unity, and define the longitudinal (zz) and transverse (xy) relaxation rates as 1/T αα 1 = dq χ αα (q, ω 0 )/βω 0 , where β = 1/t, χ αα (q, ω 0 ) is the imaginary part of the dynamical susceptibility at NMR frequency ω 0 → 0, α = x, y, z, and 1/T xy 1 = 1/T xx 1 + 1/T yy 1 . To avoid handling the analytical continuation in QMC simulations, we further adopt an approximation [35,36], where δS α i = S α i − S α i . As shown in Fig.S5 [30], this approximation gives reasonable results for a Heisenberg chain.
The results for the 3D model are shown in Fig. 4(a) and (b). At h = 0.7 where the ground state is a LSDW, the temperature dependent 1/T zz 1 develops a prominent peak at the ordering temperature, signaling enhanced critical fluctuations. But 1/T xy 1 only shows a kink at the transition and is significantly suppressed in the ordered phase. At higher fields where the ground state is the TAF, the peak feature is seen in 1/T xy 1 , and 1/T zz 1 drops rapidly in the ordered phase. Above the transition we find an algebraic temperature dependence of 1/T xy 1 (Fig. 4(c)), signaling a TLL behavior. In a TLL 1/T xy 1 ∼ T η−1 according to bosonization results [37]. Fitting to this function, we extract the η parameter at each field, as shown in Fig. 4(d). Surprisingly η > 1 for h 0.85, indicating dominant longitudinal fluctuations in Heisenberg chains. With increasing field, η decreases monotonically, and an η inversion takes place at h inv ≈ 0.85. Interestingly, h inv ≈ h 2 , which separates the LSDW and TAF ground states. The η inversion and the peak feature of 1/T 1 at transition indicate that the condensation of the dominant fluctuations leads to the corresponding type of magnetic order in this system.
Discussions and Conclusion. Our QMC results provide the first clear evidence of a LSDW phase in a 3D model. The calculated 1/T 1 data unambiguously show that this phase is stabilized by the enhanced incommensurate longitudinal fluctuations in Heisenberg chains. The enhancement of the longitudinal fluctuations originates from the interchain Ising anisotropy and this effect is beyond the conventional meanfield scenario [19] in which the dynamical effects of the interchain couplings are ignored so that the transverse fluctuations always dominate.
When h > h inv the dominant spin fluctuations become transverse and the TAF ground state is correspondingly stabilized. The transverse fluctuations also govern in the quantum critical regime where the system shows quantum critical TLL behavior at intermediate temperatures and converges to the (3 + 2)D XY universality at low temperatures. Such a scenario of quantum criticality generally holds for a broad class of weakly coupled XXZ spin chain systems that have the same symmetry as the model in Eq. (1).
For the model studied here, we find that the phase diagram is sensitive to the interchain coupling parameters ε and J ab . It is known that increasing J ab favors a TAF order owing to the spin-flop mechanism [17]. For fixed J ab , the enhancement of longitudinal correlations and hence the stabilization of LSDW only take place when ε is less than a critical value. As illustrated in Fig.S6 [30] the LSDW is absent for ε 0.5 at J ab /J c = 0.2. When ε → 0, however, the transverse correlations are only within each chain. Because the TAF order breaks a continuous U (1) symmetry, it can not be stabilized in this limit. In this case the QPT is from the LSDW to the fully polarized phase and belongs to the (3 + 1)D universality. But for a finite ε we always find a TAF phase before the magnetization is fully saturated (see Fig.S6 [30]). Hence the LSDW is irrelevant to the QPT, which is always controlled by the (3 + 2)D XY universality.
In what follows we discuss implications of our results to YbAlO 3 and other related Q1D quantum magnets. It is found that the intrachain exchange coupling of YbAlO 3 is almost isotropic but the interchain one is dominant by dipole-dipole interaction containing strong Ising anisotropy [26,38]. This is fully captured by the Hamiltonian in Eq. (1), where the interchain Ising anisotropy is ensured by the finite ε < 1. Taking the measured values J c ∼ 0.22 meV and g ∼ 7.6 [26], we can compare our results with experimental ones for YbAlO 3 . As shown in Fig. 1(d), the phase boundary of the model agrees qualitatively with the adapted experimental one. In particular, the LSDW state in our model naturally explains the observed unusual incommensurate AFM order. Note that INS measurement suggests a ferromagnetic interchain coupling [26]. Though for demonstration we take J ab > 0, the stabilization of the LSDW phase and the agreement on the phase boundary indicate that our model has already captured the essential physics of the system, and the results for J ab < 0 are qualitative the same. For the quantum criticality, the linear scaling of T cr and the quantum critical TLL behavior obtained in our theory are also observed in YbAlO 3 [26], while the predicted genuine (3 + 2)D XY universality, which also well applied to other coupled XXZ spin chains such as (Ba,Sr)V 2 Co 2 O 8 , can be tested by future experiments.
Our theory also predicts a TAF order for h > h 2 . In real materials, the LSDW and TAF orders may coexist or are phase separated [22,23]. This possibility and the strong anisotropic gyromagnetic tensor in YbAlO 3 complicates the detection of transverse spin correlations and the related TAF order by INS measurements [26]. In light of the theoretical results, we hereby propose to probe the dominant low-energy spin fluctuations and associated magnetic ordering by measuring NMR 1/T 1 . Even when the system is highly anisotropic such that only 1/T zz 1 or 1/T xy 1 is detectable, the peak or kink feature in the temperature dependent 1/T 1 near the ordering temperature can still tell the dominant fluctuations and the associated magnetic order, as shown in Figs. 4(a) and (b). Moreover, it would be interesting to examine the possible η inversion by a careful study on the 1/T 1 data above the ordering temperature. Such a study can also provide important information on the dominant fluctuations and the underlying magnetic ground state. Given the universal property of TLL, similar analysis can be applied on other related Q1D systems, such as (Ba,Sr)V 2 Co 2 O 8 , where the dominant fluctuations and associated long-range magnetic orders are still under debate [22,24].
In conclusion, we study the field-induced phase diagram and quantum criticality of S = 1/2 AFM Heisenberg chains with Ising anisotropic interchain couplings. We find the interchain Ising anisotropy enhances incommensurate AFM correlations, stabilizing a LSDW ground state at low fields. The transverse spin correlations dominate at high fields and a TAF ground state is stabilized. This leads to a QCP controlled by a (3 + 2)D XY universality but displaying a finite-temperature 1D-3D crossover. The calculated NMR relaxation rates show enhanced critical fluctuations at magnetic ordering, and the enhancement takes place at particular channel relevant to the underlying magnetic order. Above the ordering temperature, the system exhibits universal TLL behavior and shows an η inversion with increasing field, where the dominant spin fluctuation changes from longitudinal to transverse. These features make NMR an ideal probe for the spin fluctuations and associated ground states of coupled spin chains. Our findings thus shed light on future experimental and theoretical studies on YbAlO 3 and other Q1D quantum magnets.  Fig. 1 of the main text. At each field, the transition to an AFM state is signaled as either a peak or a kink feature. ] spin structure factors, S zz (π, π, π + q) and S xy (π, π, π + q), respectively, at t = 0.05 and h = 0.9 in the TAF phase.  Fig. 3(b) of the main text, and the extracted correlation length exponent ν agrees with the value of 3D XY universality within error bar.