Field-induced incommensurate ordering in the Heisenberg chains coupled by Ising interaction: Model for the ytterbium aluminum perovskite YbAlO$_3$

TL;DR

This study models YbAlO$_3$ as coupled Heisenberg chains with Ising interactions, using DMRG to explain experimental magnetic and neutron scattering data, revealing field-induced incommensurate order and quantum criticality.

Contribution

It introduces a detailed spin model for YbAlO$_3$ and demonstrates its effectiveness in reproducing experimental magnetic and neutron scattering results.

Findings

Quantitative agreement with experimental magnetization curve

Identification of phase transition at 0.35T

Description of quantum critical behavior at 1.21T

Abstract

We study isotropic antiferromagnetic Heisenberg chains coupled by antiferromagnetic Ising interaction as an effective spin model for the ytterbium aluminum perovskite YbAlO$_3$. Using the density-matrix renormalization group (DMRG) method we calculate the magnetization curve, local spin, central charge, and dynamical spin structure factors in the presence of magnetic field. From the fitting of the experimental magnetization curve, the effective intrachain and interchain couplings are estimated as $J=2.3$K and $J_{\rm ic}=0.8$K, respectively. We can quantitatively explain the experimental observations: (i) phase transition from antiferromagnetic to incommensurate order at field 0.35T, and (ii) quantum critical behaviors at the saturation field of 1.21T. Furthermore, the low-energy excitations in the experimental inelastic neutron scattering spectra can be well described by our DMRG results of the dynamical structure factors.

Figures

Equations

$$H$$

$$+H\sum _{i,j}{S}_{i,j}^z,$$

$$⟨S_{n,j}^z⟩=\pm A_L(-1{)}^n\cos (2\pi Mn+\varphi )+M,$$

$$N_{\mathrm{c}}=\frac{3\left[S_L\left(\frac{L}{2}-1\right)-S_L\left(\frac{L}{2}\right)\right]}{\ln \left[\cos \left(\frac{\pi }{L}\right)\right]},$$

$$S^{\gamma \bar{\gamma }}(\vec{q},\omega )$$

$$S_{\vec{q}}^{\gamma }=\sqrt{\frac{2}{L+1}}\sum _l{e}^{i\vec{q}\cdot \vec{r}}{S}_{\vec{r}}^{\gamma },$$

$$H=J\sum _i{\vec{S}}_i\cdot {\vec{S}}_{i+1}+(-1{)}^i{h}_z{S}_i^z.$$

$$H=J\sum _i\frac{\Delta }{2}(S_i^{+}{S}_{i+1}^{-}+S_i^{-}{S}_{i+1}^{+})+S_i^z{S}_{i+1}^z+(-1{)}^i{h}_z{S}_i^z.$$

$$S^{\gamma \bar{\gamma }}(\vec{q})=\sum _{i,j}{e}^{i\vec{q}\cdot ({\vec{R}}_j-\overrightarrow{R_i})}⟨S^{\gamma }{S}^{\bar{\gamma }}⟩.$$

Full text

I Introduction

In low dimensional strongly correlated electron systems, the effects of interactions and quantum fluctuations are maximized. This situation provides us with an ideal playground to study quantum phase transitions. In fact, quasi-one-dimensional (quasi-1D) Heisenberg magnets reveal a variety of exotic phases Vasiliev et al. (2018).

Quantum phase transitions at zero temperature are characterized by a quantum critical point (QCP) Coleman and Schofield (2005); Sachdev and Keimer (2011), where the critical fluctuations are scale-invariant and the system belongs to a universality class characterized by critical exponents, independent of the microscopic details of the system Moriya (1985); Binney et al. (1992); Timusk and Statt (1999). Physical properties over a wide range of temperatures above a QCP can be influenced by critical fluctuations and physical quantities such as magnetic susceptibility and specific heat obey simple scaling laws in the critical exponent.

In the pure 1D limit, the isotropic spin-$\frac{1}{2}$ Heisenberg model belongs to a universality class, the so-called Tomonaga-Luttinger liquid (TLL), exhibiting a gapless mode in the elementary excitations at zero temperature. This means that the 1D system is extremely fragile to external perturbations such as further interactions. Typically, when the 1D chains are coupled, they crosses over from the exotic TLL to more conventional 2D or 3D order. In some cases, weaker interchain couplings can drastically change the magnetic properties and yield the dominant contribution to physical quantities Nishimoto et al. (2011). Therefore, even in the case of small interchain coupling, an uncoupled 1D chain is not always a good approximation and dimensional crossovers are a challenging and open problem.

Very recently, TLL behavior was reported in the ytterbium aluminum perovskite YbAlO3 Wu et al. (2019). Below $T_{\rm N}=0.88$ K the system is antiferromagnetically ordered, namely, in a Néel state. With applying external magnetic field, a Néel to incommensurate (IC) phase transition occurs at $H_{\rm c}=0.35$ T, and the spins fully saturate at $H_{\rm s}=1.13$ T. Scaling behaviors associated with the TLL universality class were experimentally observed as the magnetic susceptibility $dM/dH\propto\tilde{H}^{-0.51}\phi(T/\tilde{H}^{1.04})$ and the specific heat $[C(H)-C(H_{\rm s})]/T\propto\tilde{H}^{-0.5}\psi(T/\tilde{H})$ where $M$ is magnetic field, $M_{\rm s}$ is saturation field, and $T$ is temperature. Moreover, the spinon confinement-deconfinement transition was confirmed by the inelastic neutron scattering measurements, and it was theoretically explained within a single Heisenberg chain in the presence of staggered magnetic field Wu et al. (2019).

In this paper, we propose isotropic antiferromagnetic (AFM) spin-$\frac{1}{2}$ Heisenberg chains coupled by antiferromagnetic Ising interaction as an effective spin model for YbAlO3. The density-matrix renormalization group (DMRG) method is employed to study this spin model. From the fitting of experimental magnetization curve, the effective intrachain and interchain couplings are estimated as $J=2.3$ K and $J_{\rm ic}=0.8$ K, respectively. To see the stability of AFM/IC orders and to confirm a TLL criticality at the saturation field, we calculate local spin and central charge in the presence of magnetic field, respectively. Then, we can the quantitatively explain the experimental observations: (i) phase transition from AFM to IC order at field 0.35 T, and (ii) quantum critical behavior with a saturation field of 1.21 T. Furthermore, we find that the low-energy excitations in the experimental inelastic neutron scattering spectra can be well described by our DMRG results of the dynamical spin structure factors.

The paper is organized as follows: In Sec. II our Hamiltonian for the effective spin model for YbAlO3 is explained and the applied numerical method is described. In Sec. III we present our numerical results and discuss the relevance for the experimental observations. Finally we give a conclusion in Sec. IV.

II Model and method

II.1 Heisenberg chains coupled by Ising interaction

We here propose an effective spin model to describe the fundamental magnetic properties of YbAlO3. YbAlO3 is insulating and the magnetic properties come from the Yb3+ ions. The four Kramers doublets of the Yb3+ ions are split by the crystal field and the doublet ground state $m_{J}=\pm\frac{7}{2}$ is well separated from the excited levels. Therefore, the low-energy magnetic properties can be described by a pseudo-spin-$\frac{1}{2}$ model Radhakrishna et al. (1981). Magnetically, the system consists of two 1D sublattices (see Fig. 1). The magnetic coupling between the two sublattices, i.e., interchain coupling, is similar to a dipole-dipole interaction and it leads to a Ising-type coupling $S^{z}S^{z}$. On the other hand, the intrachain coupling is dominantly created by an AFM super-exchange process. Furthermore, it is known that the exchange interaction of rare earth ions has the typical form $\vec{S}\cdot\vec{S}$ Wu et al. (2016). Thus, despite the highly anisotropic doublet state of Yb3+ ions, one may assume an isotropic exchange coupling along the chain direction. The validity of this isotropic exchange coupling is also discussed in Appendix A.

We are thus dealing with isotropic Heisenberg chains coupled by Ising interaction as an effective spin model for YbAlO3. The Hamiltonian is written as

[TABLE]

where $S^{\alpha}_{i,j}$ is the $\alpha$-component of spin-$\frac{1}{2}$ operator $\vec{S}_{i,j}$ at $i$-th site in $j$-th chain. The lattice structure is sketched in Fig. 1. Since the interchain couplings form a bipartite square lattice in the $a$-$b$ plane and no exchange processes are allowed between the chains, the extension of the cluster in the $a$-$b$ plane can be reduced to a two-chain problem. Hence, in this paper, we study two isotropic Heisenberg chains coupled by Ising interaction. We note that, in order to count the strength of the interchain coupling consistently with the material, the interchain Ising coupling should be replaced by $\tilde{J}_{\rm ic}=4J_{\rm ic}/N_{\rm c}$, where $N_{\rm c}$ is the number of neighboring chains, i.e., $N_{\rm c}=1$ for two chains.

II.2 Density-matrix renormalization group

To examine the ground state of the system (1) we employ the DMRG technique which is a powerful numerical method for various (quasi) 1D quantum systems White (1992). For the calculation of static properties, we use the standard DMRG method. Either open or periodic boundary conditions are chosen dependent on the calculated quantity. We study clusters with length up to $L\times 2=600\times 2$ and keep up to $m=4000$ density-matrix eigenstates in the renormalization procedure. In this way, the maximum truncation error, i.e., the discarded weight, is less than $1\times 10^{-13}$. This high accuracy is naively expected because no quantum fluctuations are allowed between chains coupled by Ising interaction.

For the calculation of dynamical properties, we use the dynamical DMRG method which has been developed for calculating dynamical correlation functions at zero temperature in quantum lattice models Jeckelmann (2002). Since the dynamical DMRG algorithm performs best for open boundary conditions, we study a open cluster with length $L\times 2=50\times 2$. The dynamical DMRG approach is based on a variational principle so that we have to prepare a ‘good trial function’ of the ground state with the density-matrix eigenstates. Therefore, we keep $m=1200$ to obtain the ground state in the first ten DMRG sweeps and keep $m=600$ to calculate the excitation spectrum. In this way, the maximum truncation error, i.e., the discarded weight, is about $1\times 10^{-5}$, while the maximum error in the ground-state and low-lying excited states energies is about $10^{-2}J$.

III Results and discussion

III.1 Magnetization

First, to estimate the effective coupling parameters for YbAlO3, we perform a fitting of the experimental magnetization curve observed at $T=45$ mK under applied field $H$ parallel to the $a$-axis [see Fig. 2(a)]. The experimental magnetization curve exhibits; (i) an explicit gap between singlet and triplet states, the magnitude of which corresponds to the critical field $H_{c}=0.35$ T, and (ii) a sharp increase near the saturation $H_{s}=1.21$ T, which is a typical signature of strong quantum fluctuations in 1D Heisenberg systems. (Note that, although the saturation field has been estimated to be $H_{\rm s}=1.13$ T from the peak position of $dM/dH$ at $T=0.05$ K, it may be slightly shifted to higher field in the limit $T\to 0$.) These two features can be reproduced by taking $J>0$ and $J_{\rm ic}>0$ in the system (1): Due to the AFM Ising interchain coupling, a spinon propagation costs more at low field and quantum fluctuations are allowed only within each chain along the $c$-axis. This is consistent with the spinon confinement picture suggested in Ref. Wu et al., 2019.

At present, low-temperature magnetization measurement are available only for $H\parallel a$. As seen in Fig. 2(a), the observed saturation moment $M_{\rm s}=3.8\mu_{\rm B}$/Yb leads to a large $g$-factor $g^{a}=7.6$ for $H\parallel a$. Since this value is much larger than those for the other directions $g^{b}\simeq g^{c}=0.46$, the magnetization is dominated by the $a$-component. Thus, the most reliable fitting is chosen by using the magnetization curve for $H\parallel a$. The fitting result is given in Fig. 2(a), where a periodic cluster with $L\times 2=100\times 2$ sites is used for the DMRG calculations. The best fitting is obtained by assuming that the crystallographic $a$-axis is parallel with the $z$-axis ($z\parallel a$, $g^{z}=g^{a}$) in our model and by setting $J=2.3$K and $J_{\rm ic}=0.8$K. If we assume the interchain coupling to be ferromagnetic (FM), i.e., $J_{\rm ic}<0$, only a rough fitting is possible (see Appendix B). We note that the AFM and FM Ising interchain couplings give quantitatively the same results except for the anti-phase or in-phase chains at $M=0$; and, qualitatively similar results for $M>0$. Because of this fit to the experimental magnetization, we mostly consider the case of AFM Ising interchain coupling hereafter.

For confirmation, let us now examine the validity of the Ising-type interchain coupling. To test it, additional exchange terms $\Delta_{\rm ic}J_{\rm ic}(S^{+}_{i,j}S^{-}_{i,j^{\prime}}+S^{-}_{i,j}S^{+}_{i,j^{\prime}})$ are added to the Ising interchain coupling of the Hamiltonian (1). In Fig. 2(b) we show the magnetization curve for several values of $\Delta_{\rm ic}$, where the values of $J$ and $J_{\rm ic}$ are tuned to keep the ratio between $H_{\rm c}$ and $H_{\rm s}$. We find that the magnetization is lifted up from the Ising limit ($\Delta_{\rm ic}=0$) even by small $\Delta_{\rm ic}=0.2$ at the intermediate field $0.5~{}{\rm T}\lesssim H\lesssim 1~{}{\rm T}$. This means that finite $\Delta_{\rm ic}$ only brings a further deviation from the experimental magnetization. This confirms that the simple AFM Ising interchain coupling gives the best fit.

III.2 $z$-component of local spin

Experimentally, an AFM-IC phase transition has been observed at $H=0.35$ T Wu et al. (2019). To investigate the nature of the AFM-IC transition and the stability of the IC ordering, we calculate the $z$-component of local spin, $\langle S^{z}_{n,j}\rangle$. Here, open boundary conditions are applied. This allows us to directly observe a translation-symmetry broken state due to Friedel oscillations. In Fig. 3(a-d) we show $\langle S^{z}_{n,j}\rangle$ as a function of the position $n$ in either of the chains for several strengths of the magnetic field. We use an open cluster with $600\times 2$ sites for the system (1) and the result for the central $100$ sites is plotted. Since the system (1) is unfrustrated, by considering the Fermi wavenumber after Jordan-Wigner transformation, a possible form of $\langle S^{z}_{n,j}\rangle$ is derived as

[TABLE]

where $A_{L}$ is the amplitude of oscillation and the $\pm$ sign refers to the two sublattices ($j=1$ or $2$).

At $H=0$, Eq. (2) leads to $\langle S^{z}_{n,1}\rangle=\pm A_{L}(-1)^{n}$. As shown in Fig. 3(a), this is a pure AFM ordering and remains unchanged up to the critical field $H=0.35$ T. In consistency with the experimental observation, the IC modulation appears at $H>0.35$ T. At $H=0.36$ T, which is slightly higher than the critical field, an IC modulation with long wavelength is clearly seen [Fig. 3(b)]. By fitting the modulation with Eq. (2), we estimate the wavelength $\Delta n\sim 100$ ($M\sim 0.01$) and $A_{600}=0.376$. With increasing field, $\Delta n$ becomes shorter and $A_{L}$ decreases. At $H=0.44$ T [Fig. 3(c)], we estimate $\Delta n\sim 20$ ($M\sim 0.05$) and $A_{600}=0.216$. In the low-field region, the fitting with Eq. (2) is rather straightforward because the wavelength is much longer than the lattice spacing, i.e., $\Delta n\gg 1$. However, if the wavelength is comparable with the lattice spacing, i.e., $\Delta n\sim{\cal O}(1)$, the fitting is not very simple. In Fig. 3(d) the profile of $\langle S^{z}_{n,1}\rangle$ at $H=0.89$ T is shown. It looks unclear how to fit it with Eq. (2). This problem can be solved by plotting, for example, $\langle S^{z}_{n,1}\rangle$ for odd $n$ and $\langle S^{z}_{n,2}\rangle$ for even $n$, namely, $\langle S^{z}_{n,\lambda}\rangle$ with $\lambda=\frac{3+(-1)^{n}}{2}$. Fig. 3(e) shows $\langle S^{z}_{n,\lambda}\rangle$ at $H=0.89$ T. In this way, the fitting with Eq. (2) can be handily performed and we estimate $\Delta n\sim 5$ ($M\sim 0.2$). This method works in the whole range of IC phase.

It is important to examine whether the AFM and IC modulations are long-range ordered. To determine it, we perform finite-size scaling analysis of the amplitude $A_{L}$. We look first at the stability of the AFM order as a function of the Ising interchain coupling. In Fig. 4(a) the scaling of $A_{L}$ in the AFM phase, i.e., at $M=0$, is shown for several strengths of the interchain coupling. Note that the FM and AFM Ising couplings give exactly the same results in this case. The finite-size effect seems to be very small once the system size is extended beyond the critical length. The extrapolated values of $A_{L}$ in the thermodynamic limit $L\to\infty$ are plotted as a function of $J_{\rm ic}$ in Fig. 4(b). We can confirm that the AFM long-range order is stabilized only if the chains are coupled through Ising interaction. In the AFM phase ($0<H<0.35$T) our YbAlO3 parameter $J_{\rm ic}/J=0.35$ leads to $A_{\infty}=0.393$, which is slightly reduced from the full value $0.5$ by the intrachain quantum fluctuations.

Next, we turn to the IC state, i.e., at $0<M<0.5$. In Fig. 4(c) the finite-size scaling is shown for several strengths of magnetic field. For all $H$, the best scaling function seems to be linear by choosing the horizontal axis as $1/\sqrt[3]{L}$. This might be related to the three-dimensional ordering of the IC modulation. Fig. 4(b) shows the extrapolated values of $A_{\infty}$ as a function of magnetic field. At the critical field $H=0.35$ T, $A_{\infty}$ drops down from 0.393 to 0.042, indicating a first order transition between the AFM and IC phases transition. With further increasing magnetic field, the value of $A_{\infty}$ decreases gradually and becomes zero at the saturation field $H=1.21$ T, indicating a second order or continuous transition. Thus, we find that $A_{\infty}$ is small but finite at $0.35~{}{\rm T}<H<1.21~{}{\rm T}$. This means that the IC modulation is long-range ordered in the whole range of $0<M<0.5$. To support this statement, we consider the central charge in the next subsection.

III.3 Central charge

In general, the saturation field at $T=0$ is a quantum critical point. In YbAlO3, the TLL behavior has been experimentally observed near the saturation field. To investigate the consistency with our model (1), we calculate the central charge $c$ which provides definitive information on the universality class of $(1+1)$ dimensional system Cardy (1996). A system in the TLL phase belongs to the Gaussian universality class ($c=1$) and $c<1$ is expected for the gapped phase from the renormalization in the massive region. The central charge can be numerically calculated through the von Neumann entanglement entropy $S_{L}(l)=-{\rm Tr}_{l}\rho_{l}\log\rho_{l}$, where $\rho_{l}={\rm Tr}_{L-l}\rho$ is the reduced density matrix of the subsystem with length $l$ and $\rho$ is the full density matrix of the whole system with length $L$. Using conformal field theory (CFT), the relation between $S_{L}(l)$ and $c$ has been derived Affleck and Ludwig (1991); Holzhey et al. (1994); Calabrese and Cardy (2004): $S_{L}(l)=\frac{c}{3}\ln\left[\frac{L}{\pi}\sin\left(\frac{\pi l}{L}\right)\right]+s_{1}$, where $s_{1}$ is a non-universal constant. A prime objective of using this formula is to estimate the central charge Laflorencie et al. (2006); Legeza et al. (2007). For the system (1) under periodic boundary conditions, the central charge is obtained via Nishimoto (2011)

[TABLE]

where $N_{\rm c}$ is the number of chains considered.

For reference, let us first look at the central charge for an isotropic SU(2) Heisenberg chain. It is plotted as a function of the magnetization $M(=S^{z}_{\rm tot}/(N_{\rm c}L))$ for several system lengths $L$ in Fig. 5(a). As expected, the system is always critical $c=1$ at $0\leq M<0.5$. Near the saturation $M=0.5$, it goes up steeply. In the thermodynamic limit, there would be a jump at $M=0.5$, which signifies the quantum criticality at the saturation field. In Fig. 5(b) a finite-size scaling analysis of the central charge with keeping $S^{z}_{\rm tot}=\frac{L-2}{2}$ ($M=\frac{L-2}{2L}$) as a function of $1/L^{2}$ is performed. The thermodynamic limit $L\to\infty$ corresponds to the saturation limit $M\to 0.5$. Thus, we obtain $c=1.215854$ in the limit of $M=0.5$. The origin of this $c$ value is still unknown and it should be further studied in the future.

Then, we turn to our model (1). In Fig. 5(c) the central charge of the system (1) is plotted as a function of the magnetization for several system lengths. At $M=0$, we find $c=0$ as expected from the gapped ground state. At $0<M<0.5$, interestingly unlike in the case of isotropic SU(2) Heisenberg chain, the central charge approaches the Ising universality class ($c=0.5$) with increasing $L$, as the SU(2) symmetry is reduced to U(1) symmetry due to the Ising interchain coupling. Near the saturation $M=0.5$, the convergence to $c=0.5$ seems to be very slow with $L$ since the system is nearly critical with the friable IC order. Nevertheless, the crossing point with $c=1$ for $L=64$ is already $M\sim 0.43$ close to the saturation $M=0.5$, namely, $c<1$ at $0<M<0.43$. Therefore, it is most likely that the system belongs to the Ising universality class over the whole region of $0<M<0.5$ in the thermodynamic limit. The central charge jumps from $c=0.5$ to $\approx 1.5$ at $M=0.5$. CFT Popkov and Salerno (2005); Castro-Alvaredo and Doyon (2011) predicts the value of $c=1.5$ for the ferromagnetic point. In Fig. 5(d) the central charge at fixed $M=(L-1)/(2L)$ for a system with length $L$ is plotted as a function of $1/L$. As $1/L\to 0$, i.e., $M$ approaches the saturation value $M=0.5$, the central charge approaches $1.5$. Therefore, we can confirm that the saturation field is a quantum critical point. This is consistent with experimental observations Wu et al. (2019).

III.4 Dynamical spin structure factor

For calculating dynamical properties, we use the dynamical DMRG method. To examine the low-energy excitations and their development with the magnetic field, we calculate the dynamical spin structure factor, defined as

[TABLE]

where $\gamma$ is $z$ or $-(+)$, $|\psi_{\nu}\rangle$ and $E_{\nu}$ are the $\nu$-th eingenstate and the eigenenergy of the system, respectively ($\nu=0$ corresponds to the ground state). Under open boundary conditions, we define the momentum-dependent spin operators as

[TABLE]

with (quasi-)momentum $\vec{q}=(\pi Z_{x}/(L+1),\pi Z_{y})$ for integers $1\leq Z_{x}\leq L$ and $Z_{y}=0,1$. We use an open cluster with $50\times 2$ sites. Due to the AFM correlation between the neighboring chains, we restrict ourselves to the case $\vec{q}=(q_{c},\pi)$ where $q_{c}$ is momentum along the chain ($c$-axis). In Fig. 6 DMRG results of the dynamical structure factors are shown for several strengths of magnetic field. The left, middle, and right panels corresponds to $S^{zz}(q_{c},\omega)$, $S^{+-}(q_{c},\omega)$, and $S^{-+}(q_{c},\omega)$, respectively.

At low magnetic fields ($H=0-0.35$ T), the three spectra show an explicit gap $\sim 0.22$ meV at $q_{c}=\pi$. The longitudinal structure factor $S^{zz}(q_{c},\omega)$ has a dominant $\delta$-peak at $(q_{c},\omega)=(\pi,0)$ indicating the AFM order. The weight of this peak is about 90$\%$ of the total weight $\int S^{zz}(q_{c}=\pi,\omega)d\omega$ and it is denoted by an open circle in Fig. 6. The thick (or blur) dispersion of $S^{zz}(q_{c},\omega)$ is a typical signature of Ising-type spin anisotropy, which is effectively induced by the Ising interchain coupling. The two-spinon continuum is seen but seems to be somewhat suppressed in comparison to those of the SU(2) Heisenberg chain. These features are clearly visible in the experimental spectra [see Ref. Wu et al., 2019]. In the AFM phase, $S^{+-}(q_{c},\omega)$ and $S^{-+}(q_{c},\omega)$ are equivalent: Their dispersions are very thin and the two-spinon continuum is further dilute. Since the AFM order is long ranged and stable, a magnon excitation picture may be a good approximation. In fact, the lower bound of the spectra is well described by the spin wave analysis $\omega(q_{c})=J^{\prime}+J\sin^{2}q_{c}$, which is shown with white dotted lines in Fig. 6.

As soon as the system goes into the IC phase, the spectra are drastically changed as a consequence of the first-order transition. Let us see $S^{zz}(q_{c},\omega)$. At $H>0.35$ T, the single large peak at $(q_{c},\omega)=(\pi,0)$ in the AFM phase is split into two peaks at $(q_{c},\omega)=((1\pm 2M)\pi,0)$; thus, $q_{c}$ is shifted from $\pi$ to [math], $2\pi$ with increasing the magnetic field. This $q_{c}$ value corresponds to the propagation number of the IC long-range order. The magnetization $M$ is $0.09$, $0.18$, $0.27$, and $0.36$ for $H=0.44$, $0.67$, $0.89$, and $1.03$ T, respectively, and lead to $q_{c}=(1\pm 0.18)\pi$, $(1\pm 0.36)\pi$, $(1\pm 0.54)\pi$, and $(1\pm 0.72)\pi$. They perfectly agree with our DMRG spectra. Remarkably, the total weight of $\int S^{zz}(q_{c}=\pi,\omega)d\omega$ seems to remain comparable to that of the IC peaks even at high field, although the IC order must be dominant. This weight comes from the staggered oscillation of local spin $\langle S^{z}_{i}\rangle$ induced by the Ising interchain coupling [see Fig. 3]. Therefore, the IC state could also be interpreted as the coexistence of two spin-density waves with $q_{c}=\pi$ and $(1\pm 2M)\pi$. Similarly to the AFM phase, the two-spinon continuum seems to be weakened due to the stability of the long-range order. Our DMRG spectra in the IC phase also agree with the experimental ones reported in Ref. Wu et al., 2019.

We then turn to the transverse structure factors $S^{+-}(q_{c},\omega)$ and $S^{-+}(q_{c},\omega)$. Their low-energy part can be interpreted by considering those of a single Heisenberg chain as the $x$- and $y$-component of spins are uncoupled between chains. It is known that AFM $S^{x}_{i}S^{x}_{j}$ and $S^{x}_{i}S^{x}_{j}$ correlations grow with increasing the magnetic field in the Heisenberg chain Affleck and Oshikawa (1999). Accordingly, a large peak appears at $(q_{c},\omega)\approx(\pi,0)$ in the whole region of the IC phase. The magnon dispersion (lower bound of the continuum) is explained by that of the XY chain under magnetic field along the $z$-axis. The dispersion of $S^{-+}(q_{c},\omega)$ is a sine-like function with nodes at $q_{c}=\pi$, $2\pi M$, and $2\pi(1-M)$ and width $J[1-\sin{\frac{\pi}{2}(1-2M)}]$; that of $S^{+-}(q_{c},\omega)$ is a sine-like function with nodes at $q_{c}=\pi$, $2\pi(1+M)=2\pi M$, and $-2\pi M=2\pi(1-M)$ and width $J[1+\sin{\frac{\pi}{2}(1-2M)}]$. Furthermore, it is worth noting that in the higher energy range around $\omega\sim 0.5-0.8$ meV there exist another continuum in $S^{+-}(q_{c},\omega)$. They are nearly separated from the low-energy structures by $\sim J^{\prime}$; and thus, reveal correlations between the two sublattices. If the transverse structure factors could be experimentally observed, more sophisticated study on YbAlO3 would be allowed.

III.5 Effective 1D model: Heisenberg chain with staggered field

The effect of the Ising interchain coupling $S^{z}_{i,j}S^{z}_{i,j^{\prime}}$ may be mimicked by a self-consistent staggered magnetic field $(-1)^{i}S^{z}_{i}$ Schulz (1996); Bocquet et al. (2001). If this is the case, our model (1) could be mapped onto a single Heisenberg chain with a staggered magnetic field:

[TABLE]

In fact, this model was used to analyze some experimental observations in Ref. Wu et al., 2019. A simplest mapping of the model (1) onto the model (6) can be performed by comparing the low-lying energies of a small cluster. We take a 4-site plaquette for (1) and a 2-site dimer for (6): The ground state energies in the $S^{z}=0$ sector are $-\frac{J}{2}-\sqrt{J^{2}+\frac{J^{\prime 2}}{4}}$ and $-\frac{J}{4}-\sqrt{\frac{J^{2}}{4}+h^{2}}$; and those in the $S^{z}=1$ sector are $-\frac{J}{2}$ and $-\frac{J}{4}+\sqrt{\frac{J^{2}}{4}+h^{2}}$, respectively. By comparing the energy differences between two spin sectors, i.e., $\sqrt{J^{2}+\frac{J^{\prime 2}}{4}}=2\sqrt{\frac{J^{2}}{4}+h^{2}}$, it leads to $h=\pm\frac{J^{\prime}}{4}$. Thus, we estimated the parameters for the effective single spin chain as $J=2.3K$ and $h=0.8K$. These values agrees well to those ($J=2.4$ K, $h=0.66$ K used in Ref. Wu et al., 2019.

IV Conclusion

Using the DMRG technique we studied isotropic AFM Heisenberg chains coupled by AFM Ising interaction as an effective spin model for the ytterbium aluminum perovskite YbAlO3. The effective intrachain and interchain couplings were estimated as $J=2.3$ K and $J_{\rm ic}=0.8$ K, respectively, from the fitting of the experimental magnetization curve. At $0.35~{}{\rm T}<H<1.21~{}{\rm T}$ ($0<M<0.5$), a long-ranged IC order is stabilized by the Ising interchain coupling; accordingly, the system belongs to the Ising universality class ($c=0.5$) in the IC phase. At the saturation field, the central charge jumps from $c=0.5$ to a FM value $c=1.5$ which is a signature of QCP. Our results fully agree with the experimental observations reported in Ref. Wu et al., 2019. Furthermore, we calculated the dynamical structure factors. Our calculations quantitatively reproduce the low-energy excitations in the experimental spectra. In this paper, we focused on the case of $H\parallel a$, as the low-temperature experimental data are available only for this case. For more sophisticated theoretical analysis, measurements for other field directions are required.

Acknowledgements

This work is supported by the DFG through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 39085490) and through SFB 1143 (project-id 247310070). We thank S. E. Nikitin for fruitful discussion and the experimental magnetization data and U. Nitzsche for technical assistance.

Appendix A Intrachain XXZ anisotropy

Generally, the exchange interaction of rare earth ions has a typical form $\vec{S}\cdot\vec{S}$. Thus, despite the highly anisotropic doublet state of Yb3+ ions, we assumed that the exchange coupling along the chain direction is isotropic, i.e., an XXX chain in the main text. To verify this assumption, we consider the XXZ-anisotropy dependence of the dynamical spin structure factor. We then introduce the XXZ anisotropy into the effective single chain (6) as

[TABLE]

In Fig. 7 we show the result for the dynamical spin structure factor for several $\Delta$ values with a periodic 24-site cluster. In the experimental spectra, the lowest excitations, i.e., gaps, at $q_{c}=0$ and $q_{c}=\pi$ are nearly the same for the AFM state. As seen in Fig. 7, this situation is realized at $\Delta\sim 1$. Therefore, it is reasonable to assume the intrachain interaction to be isotropic.

Appendix B Magnetization with ferromagnetic Ising interchain coupling

As shown in the main text, the best fit of the experimental magnetization was achieved by taking the Ising interchain coupling to be AFM in our effective spin model (1). If we assume it to be FM, the fitting is less accurate (see Fig. 8): In the case of FM Ising interchain coupling, (i) a sharp rise of the experimental $M$ at the AFM-IC phase transition $H=0.35$ T can not be reproduced, and (ii) the increase of $M$ around the saturation field is too steep in comparison to the experimental curve. Therefore, we decided to focus on the case of AFM Ising interchain coupling in the main text. If the low-temperature magnetization curves for the other field directions are made available, more sophisticated fitting analysis may be performed.

Appendix C Static structure factor

To see the magnetic structure in the ground state, we calculate the static spin structure factors. They are defined by a Fourier transform of the spin-spin correlation functions

[TABLE]

We plot the results of static spin structure factors in Fig. 9. For $H=0-0.35$ T, the longitudinal structure factor $S^{zz}(\vec{q})$ has a large peak at $\vec{q}=(\pi,\pi)$ indicating the AFM order. As soon as the system goes into the IC phase, this peak is abruptly lowered and the peak position is shifted to small $q_{c}$ with increasing the magnetic field. On the other hand, a FM peak at $\vec{q}=(0,0)$ in $S^{zz}(\vec{q})$ begins to grow around $H=0.67$ T. The transverse structure factor $S^{+-}(\vec{q})$ always has a peak at $\vec{q}=(\pi,0)$ and $\vec{q}=(\pi,\pi)$; the peak becomes sharper with increasing the magnetic field because of the enhancement of AFM $S^{x}_{i}S^{x}_{j}$ and $S^{x}_{i}S^{x}_{j}$ correlations.