Griffiths Singularities in the Random Quantum Ising Antiferromagnet: A Tree Tensor Network Renormalization Group Study

TL;DR

This study uses a tree tensor network renormalization group approach to analyze the quantum phase transitions and Griffiths singularities in a disordered antiferromagnetic quantum Ising chain, revealing an infinite-randomness critical point and the effects of longitudinal fields.

Contribution

The paper introduces a new matrix product operator representation for high-order moments, enabling efficient detection of quantum phase transitions and Griffiths effects in disordered quantum chains.

Findings

Identifies an infinite-randomness quantum critical point with large dynamical exponent z>1.

Demonstrates pronounced Griffiths singularities due to rare regions with slow fluctuations.

Shows that longitudinal fields suppress Griffiths effects and alter the phase transition dynamics.

Abstract

The antiferromagnetic Ising chain in both transverse and longitudinal magnetic fields is one of the paradigmatic models of a quantum phase transition. The antiferromagnetic system exhibits a zero-temperature critical line separating an antiferromagnetic phase and a paramagnetic phase; the critical line connects an integrable quantum critical point at zero longitudinal field and a classical first-order transition point at zero transverse field. Using a strong-disorder renormalization group method formulated as a tree tensor network, we study the zero-temperature phase of the quantum Ising chain with bond randomness. We introduce a new matrix product operator representation of high-order moments, which provides an efficient and accurate tool for determining quantum phase transitions via the Binder cumulant of the order parameter. Our results demonstrate an infinite-randomness quantum critical point in zero longitudinal field accompanied by pronounced quantum Griffiths singularities, arising from rare ordered regions with anomalously slow fluctuations inside the paramagnetic phase. The strong Griffiths effects are signaled by a large dynamical exponent $z>1$, which characterizes a power-law density of low-energy states of the localized rare regions and becomes infinite at the quantum critical point. Upon application of a longitudinal field, the quantum phase transition between the paramagnetic phase and the antiferromagnetic phase is completely destroyed. Furthermore, quantum Griffiths effects are suppressed, showing $z<1$, when the dynamics of the rare regions is hampered by the longitudinal field.

Equations

$$\rho (\epsilon )\sim (\ln \lambda {)}^{-1}{\epsilon }^{\omega },$$

$$\epsilon \sim L^{-1/(\omega +1)}\sim L^{-z}$$

$${\chi }^{\text{loc}}(T)\sim T^{-1+1/z},$$

$${\epsilon }_{\text{typ}}\sim \exp (-c{L}^{1/2}),c=\text{const.}$$

$$\overline{\epsilon }\sim \exp (-c^{\prime }{L}^{1/3}),c^{\prime }=\text{const.}$$

$$H=\sum _i{J}_i{\sigma }_i^z{\sigma }_{i+1}^z-\sum _i{\Gamma }_i{\sigma }_i^x-\sum _i{h}_i{\sigma }_i^z,$$

$$H_{\text{cl}}=\sum _{i,n}{K}_i{s}_{i,n}{s}_{i+1,n}-\sum _{i,n}{K}_i^{\prime }{s}_{i,n}{s}_{i,n+1}-\sum _{i,n}{h}_i^{\prime }{s}_{i,n},$$

$$K_i$$

$$e^{-2K_i^{\prime }}$$

$$h_i^{\prime }$$

$$\epsilon \sim \left|1-\Gamma \right|,$$

$$\xi \sim {\left|1-\Gamma \right|}^{-1},$$

$${\Gamma }_c\sim 1-a{h}^2$$

$${\Gamma }_c\sim b(2-h)$$

$$\overline{\ln J}=\overline{\ln \Gamma }.$$

$$-\ln C_{\text{typ}}(r)\sim \sqrt{r}.$$

$$\overline{C}(r)\sim \frac{1}{r^{2-\varphi }},$$

$${\overline{C}}_{1,L}\sim \frac{1}{L}$$

$$\delta =\frac{\overline{\ln \Gamma }-\overline{\ln J}}{\text{var}(\ln h)+\text{var}(\ln J)},$$

$$\overline{{\left(\frac{J}{\Gamma }\right)}^{1/z}}=1,$$

$$z\approx \frac{1}{2\left|\delta \right|},\left|\delta \right|\ll 1.$$

$$\begin{array}{rl}{W}_{H_p}^{[i]}& =\left(\begin{array}{cccc}\text{mathds}1& 0& 0& 0\\ 0& \text{mathds}1& 0& 0\\ {\sigma }_i^z& 0& 0& 0\\ -{\Gamma }_i{\sigma }_i^x-h_i{\sigma }_i^z& 0& J_i{\sigma }_i^z& \text{mathds}1\end{array}\right),\\ & 2\le i\le L-1,\end{array}$$

$$W_{H_p}^{[1]}=\left(\begin{array}{cccc}-{\Gamma }_1{\sigma }_1^x-h_1{\sigma }_1^z& J_L{\sigma }_1^z& J_1{\sigma }_1^z& \text{mathds}1\end{array}\right),$$

$$W_{H_p}^{[L]}=\left(\begin{array}{c}\text{mathds}1\\ {\sigma }_L^z\\ {\sigma }_L^z\\ -{\Gamma }_L{\sigma }_L^x-h_L{\sigma }_L^z\end{array}\right).$$

$$\begin{array}{rl}{W}_{H_o}^{[i]}& =\left(\begin{array}{ccc}\text{mathds}1& 0& 0\\ {\sigma }_i^z& 0& 0\\ -{\Gamma }_i{\sigma }_i^x-h_i{\sigma }_i^z& J_i{\sigma }_i^z& \text{mathds}1\end{array}\right),\\ & 2\le i\le L-1,\end{array}$$

$$W_{H_o}^{[1]}=\left(\begin{array}{ccc}-{\Gamma }_1{\sigma }_1^x-h_1{\sigma }_1^z& J_1{\sigma }_1^z& \text{mathds}1\end{array}\right),$$

$$W_{H_o}^{[L]}=\left(\begin{array}{c}\text{mathds}1\\ {\sigma }_L^z\\ -{\Gamma }_L{\sigma }_L^x-h_L{\sigma }_L^z\end{array}\right).$$

$$O^n={\left(\frac{1}{L}\sum _{i=1}^L{O}_i\right)}^n,$$

$$O^n=\frac{1}{L^n}{W}_{O^n}^{[1]}{W}_{O^n}^{[2]}\cdots W_{O^n}^{[L]},$$

$$W_{O^n}^{[1]}=\left(\begin{array}{ccccc}{O}_1^n& \left(\genfrac{}{}{0ex}{}{n}{1}\right)O_1^{n-1}& \left(\genfrac{}{}{0ex}{}{n}{2}\right)O_1^{n-2}& \cdots & \text{mathds}1\end{array}\right),$$

$$W_{O^n}^{[L]}=\left(\begin{array}{c}\text{mathds}1\\ O_L\\ O_L^2\\ ⋮\\ O_L^n\end{array}\right),$$

Full text

I Introduction

The interplay between quantum fluctuations and randomness often leads to drastic disorder effects, resulting in some exotic phenomena at and near zero-temperature phase transitions. For example, there exists one class of disordered systems where strong Griffiths singularities, Griffiths ; McCoy ; McCoy_Wu which arise from the formation of rare strongly coupled region, lead to the divergence of the dynamical exponent at quantum critical points. The random transverse-field Ising chain is the most prominent example for such spectacular phenomena and its quantum critical point is, in the renormalization-group language, described as an infinite-randomness fixed point, where the system appears more and more disordered under coarse graining. Fisher_PRL ; Fisher_Ising ; Motrunich There are many other one-dimensional random quantum systems that also exhibit infinite-randomness phases, such as the random singlet phases of SU(2)k anyonic chains for all $k\geq 2$ (including the case $k\to\infty$ which corresponds to the random spin-1/2 Heisenberg chain) Fisher_AF ; Bonesteel and the random J-Q spin-1/2 chains with multi-spin couplings. SDRG_JQ It is noteworthy that the infinite-randomness fixed points can control more than just one dimension; in particular, the quantum phase transitions in higher dimensions with Ising (Z2) symmetry have been numerically demonstrated to be of infinite-randomness type. QMC_2D ; Lin_2D ; Istvan_2D Experimentally, Griffiths singularities have been observed in various quantum phases with quenched disorder. Neto98 ; Aeppli07 The signature of an infinite-randomness phase with the diverging dynamical critical exponent has also been detected in recent experimental studies on random Heisenberg chains RS_exp ; RS_exp1 and the superconductor-metal quantum phase transition of a two-dimensional superconducting system. SC_exp

The Griffiths singularity was first discussed nearly fifty years ago in the context of the randomly dilute Ising ferromagnet; Griffiths it arises because there are in general rare but arbitrarily large spatial regions which contain no vacancies and therefore show local magnetic order even when the system is globally in the paramagnetic phase. Such locally ordered regions produce non-analyticity in the free energy at zero field, known as the Griffiths singularity, in the paramagnetic phase below the critical temperature of the clean model. The same rare region effect can be generalized to non-dilute random-bond systems. The Griffiths singularity of the free energy in the context of a thermal transition is in general an essential singularity, Wortis ; Harris_G ; Shankar and is thus too weak to be observed experimentally. By contrast, effects of rare regions on dynamical properties are more promising to detect in experiments or numerical simulations; Palmer in particular, in a quantum system, where both dynamics and statics are inextricably coupled, the rare regions are strongly correlated along the time-like direction, which enhances the Griffiths phenomena. Thill ; Fisher_PRL ; Fisher_Ising ; Young_QCP1996 ; Young_GS1996 ; Vojta_Rev

Below we use a simple quantum spin chain with binary disorder distribution to explain the origin of the quantum Griffiths singularity. Generalizations to other disorder distributions and higher dimensions are straightforward. Thill ; Vojta_Rev We consider here a transverse-field Ising model with a constant transverse field $\Gamma=1$ perpendicular to the Ising axis and with random bonds ($J$) drawn from the binary distribution: $\pi(J)=p\delta(J-\lambda)+(1-p)\delta(J-\lambda^{-1})$, where $\lambda\gg 1$ and $p<1$, corresponding to a quantum paramagnet. RW ; Frechet The probability of finding a strongly coupled region of $n$ strong bonds ($J=\lambda$) decreases exponentially with the volume as $P(n)\approx p^{n}$. The energy gap of such a rare region is also exponentially small in its volume $\varepsilon\approx 1/\lambda^{n}$, corresponding to an exponentially large relaxation time. Thus, by change of variables we obtain a power-law density of states for the small energy gaps:

[TABLE]

with $\omega=-\ln(p)/\ln(\lambda)-1$. In a finite chain of length $L$, from the integrated density of states $N_{\varepsilon}=\int_{0}^{\varepsilon}d\varepsilon^{\prime}\rho(\varepsilon^{\prime})\sim L^{-1}$ we arrive at the relation between the length scale and the energy scale:

[TABLE]

where we have introduced the dynamical exponent $z=(\omega+1)^{-1}$, which varies continuously with the disorder parameter $p$ and becomes larger and larger as the quantum critical point ($p=1$) is approached. The power-law density of low energy excitations is responsible for singular low-energy behavior of various observables away from the critical point; for example the average local susceptibility varies with the temperature as

[TABLE]

which diverges at zero temperature $T\to 0$ if $z>1$, even though the system is in the paramagnetic phase.

Many quantitative and exact asymptotic results for the scaling behavior in the Griffiths phase and in particular at the random quantum critical point can be obtained by a powerful renormalization group method, called the strong-disorder renormalization group (SDRG), which was first introduced by Ma et al. Ma ; Dasgupta and extended by Fisher Fisher_Ising and others. SDRG_rev In the SDRG picture, the rare regions that give rise to quantum Griffiths singularities are strongly correlated sites connected with effective interactions generated during the action of the renormalization. Analytical SDRG results for the transverse-field Ising chain show that with any bond randomness the length dependence of a typical energy gap at the quantum critical point has, in contrast to the power law in Eq. (2), an exponential form: Fisher_PRL ; Fisher_Ising

[TABLE]

corresponding to activated dynamical scaling, and $z\to\infty$; on the other hand, the average gap scales as Fisher_Young

[TABLE]

where the overbar denotes an average over disorder realizations. This unconventional dynamical scaling along with the distinction between average and typical (i.e. most probable) values characterize an infinite-randomness fixed point, where energy gaps become broadly distributed even on the logarithmic scale. Away from the infinite-randomness critical point in the Griffiths phase, the distribution of low energy excitations and the dynamical exponent $z$ can also be obtained within the SDRG framework. Igloi In the paramagnetic phase, where the rare regions are localized, the distribution of their smallest energy gaps is identified as a Fréchet-type distribution. Fisher_Young ; Frechet Interestingly, many analytical SDRG results for one-dimensional systems, including the relations given in Eqs. (4) and (5), can also be obtained from a variety of statistical models. Monthus ; McKenzie ; Monthus_RW ; Z

Here we investigate strong disorder effects on the antiferromagnetic Ising spin chain in transverse and longitudinal magnetic fields using a generalization of SDRG represented as a tree tensor network. TSDRG We introduce a novel matrix product operator representation of high-order moments of the staggered magnetization and evaluate the Binder cumulant to identify the quantum critical points. In the absence of disorder, the quantum spin chain is antiferromagnetic in weak fields at zero temperature, and paramagnetic in strong fields; the two phases are separated by a critical line ending at a classical multicritical point where the transverse field vanishes. In the presence of bond randomness, in addition to the infinite-randomness critical point and Griffiths regions in zero longitudinal field, we find that the phase transition is destroyed and quantum Griffiths effects are weakened upon a longitudinal field is applied.

The paper is organized as follows. In Sec. II we define the model and summarize some known properties of the zero-temperature phases. In Sec. III we describe the matrix product operator representations used for our calculations; in particular, we propose a new representation of high-order moments. In Sec. IV we extend the technique used in Ref. TSDRG, and introduce the tensor network SDRG for the quantum Ising chain with bond randomness plus on-site disorder, and provide the results. We conclude in Sec. V with a summary and discussion.

II The model

We study the antiferromagnetic Ising chain in the presence of both longitudinal and transverse magnetic fields, described by the Hamiltonian:

[TABLE]

where $\sigma_{i}^{z}\,,\sigma_{i}^{x}$ are Pauli spin operators, $h_{i}$ ($\Gamma_{i}$) is the longitudinal (transverse) field applied to site $i$, and $J_{i}>0$ is a nearest-neighbor antiferromagnetic interaction which favors staggered magnetic ordering along the $z$ axis.

The Hamiltonian of the quantum spin chain defined in Eq. (6) can be mapped to a classical two-dimensional Ising model using the path integral formalism, where the extra dimension is regarded as imaginary time Sachdev and has extent $1/T$. In this mapping the quantum model effectively corresponds to the classical model described by

[TABLE]

where $s_{i,n}=\pm 1$ and $n$ is the index in the imaginary time direction, running from 0 to $1/(T\Delta\tau)$ with $\Delta\tau$ being the width of a time slice; the corresponding couplings and fields are given by

[TABLE]

To suppress the errors from the Suzuki-Trotter expansion Suzuki ; Trotter used in this formalism, one requires $\Delta\tau\ll 1$, implying strong ferromagnetic couplings $K^{\prime}_{i}\gg 1$ along the imaginary time direction. As a result, the classical model is strongly anisotropic, equivalent to ferromagnetic Ising chains weakly antiferromagnetically coupled with each other.

II.1 Clean system

We consider first the model with a homogeneous (i.e. site-independent) coupling, $J_{i}=1,\,\forall i$, and homogeneous applied fields $\Gamma_{i}=\Gamma,\,h_{i}=h,\,\forall i$. With $h=0$, the model is exactly solvable and is equivalent to the ferromagnetic quantum Ising chain through a gauge transformation. At absolute zero temperature, the model with zero $h$ exhibits two phases: an antiferromagnetic ordered phase for $\Gamma<1$ and a paramagnetic phase for $\Gamma>1$. The phase transition at $\Gamma=\Gamma_{c}=1$ between these two phases is characterized by a vanishing energy gap and a diverging correlation length. As $\Gamma$ approaches the critical value, the lowest energy gap vanishes as Sachdev

[TABLE]

and the correlation length diverges as

[TABLE]

implying the dynamical exponent $z=1$ and the correlation length critical exponent $\nu=1$. This zero-temperature phase transition belongs to the two-dimensional Ising universality class.

With finite value of $h$, the model is distinct from the ferromagnetic system. In the ferromagnetic case the quantum phase transition occurs only in zero longitudinal field $h=0$, while in the antiferromagnetic model the phase transition remains at $h\neq 0$, yet the longitudinal field can further destabilize antiferromagnetic order. As a result, the critical value $\Gamma_{c}$ at $T=0$ decreases as the strength of the longitudinal field increases and the critical points terminate at $\Gamma=0,\,h=2$, where a classical first-order transition takes place, corresponding to a multicritical point. The quantum critical line in the $\Gamma-h$ plane except for the two end points at $(\Gamma,\,h)=(1,\,0)$ and $(0,\,2)$ is not exactly solvable; nevertheless, it falls into the two-dimensional Ising universality class. AF_DMRG

For clarity, we present here the zero-temperature phase diagram in Fig. 1, reproduced by quantum Monte Carlo (QMC) and density matrix renormalization group (DMRG) calculations. The QMC algorithm is an unbiased method working in terms of the formalism given in Eq. (7) in the limit of $\Delta\tau\to 0$. Cluster Some details of the DMRG algorithm are given in Sec. III. Perturbation theory predicts that the critical line drops quadratically with increasing $h$ in the limit of $h\to 0$, and is linear in the vicinity of the multicritical point. AF_DMRG Our data agree with this asymptotic behavior: the critical points for $h\leq 0.8$ are well described by a quadratic function

[TABLE]

with $a\approx 0.19$; near the multicritical point, the critical line drops linearly as

[TABLE]

with an estimated prefactor $b\approx 1.4$.

The classical multicritical point at $(\Gamma=0,\,h=2)$ is a first-order transition point separating an antiferromagnetic ground state for $h<2$ and a ferromagnetic ground state for $h>2$. The ground state is extensively degenerate at the critical field $h=2$, where in a staggered $\uparrow\downarrow\uparrow\downarrow\cdots$ configuration any spin pointing in the opposite direction to the field can be overturned without change of energy. Domb

II.2 Random-bond chain with $h=0$

Consider now the model with quenched disorder in which the interactions $J_{i}$ and/or the transverse fields $\Gamma_{i}$ are independent random variables. In the space-time lattice model described in Eq. (7), the randomness introduced into $J$ and $\Gamma$ produces inhomogeneities extending only in spatial coordinates but it is perfectly correlated in the imaginary time direction. The random-bond quantum model in zero longitudinal field is effectively a solvable two-dimensional anisotropic Ising model; McCoy_Wu ; Shankar it undergoes a $T=0$ phase transition when Shankar

[TABLE]

The critical behavior is governed by an exactly solvable infinite-randomness fixed point within the framework of SDRG. At this fixed point, in addition to the energy gap, spin correlations $C(r)=(-1)^{r}\langle{\sigma_{i}^{z}\sigma_{i+r}^{z}}\rangle$ are also very broadly distributed, Fisher_Ising with the typical behavior

[TABLE]

The average correlation function, which is experimentally detectable, is however dominated by rare pairs of widely separated spins that have a correlation of order unity, much larger than the typical value. Thus, different from the typical correlations, the average correlations decay as a power of $r$ at the critical point, Fisher_Ising

[TABLE]

where $\phi=(1+\sqrt{5})/2$ is the golden mean. For an open chain of length $L$, one is often interested in the end-to-end correlation function ${C}_{1,L}=(-1)^{L-1}\langle{\sigma_{1}^{z}\sigma_{L}^{z}}\rangle$, which considers correlations between two end spins; the average of ${C}_{1,L}$ is also dominated by rare samples for which the two end spins are almost perfectly correlated with each other, and it decays algebraically as

[TABLE]

at criticality. Fisher_Young ; RW

For a random chain, the deviation from criticality is often parameterized by Fisher_Ising

[TABLE]

where $\text{var}(x)$ stands for the variance of the random variable $x$. Slightly away from criticality $\delta\neq 0$, the distinction between the average and the typical correlations leads to two correlation lengths. Fisher_Ising ; Fisher_Young In the disordered phase $\delta>0$, the average correlations, both $\overline{C}(r)$ and $\overline{C}_{1,L}$, decay exponentially with the true correlation length $\xi\sim 1/\delta^{2}$, implying the correlation length exponent $\nu=2$ for the random chain. The typical correlations decay even faster and its associated correlation length, $\xi_{\text{typ}}$, varies as $\xi_{\text{typ}}\sim 1/\delta$, defining another exponent $\nu_{\text{typ}}=1$.

The dynamical exponent $z$ which characterizes the Griffiths singularity in the off-critical region varies with the distance $\delta$ from the critical point. In fact, this nonuniversal exponent is the positive root of the equation:

[TABLE]

which is an exact expression in the entire Griffiths region. Z ; Igloi In the vicinity of the critical point, the dynamical exponent to the leading order is given by Fisher_Ising ; Igloi

[TABLE]

II.3 Random-bond chain with $h\neq 0$

In addition to the clean system and the random chain in the absence of longitudinal fields, we also consider the random chain in a finite longitudinal field, which to our knowledge has not been previously studied.

We will show in Sec. IV that with finite $h$ the quantum phase transition between the paramagnetic ground state and the antiferromagnetic ground state is completely smeared out. Also the pronounced Griffiths effects seen in the zero $h$ regime are strongly reduced upon a longitudinal field is applied.

III Matrix product operator representations

We have studied the ground-state properties and the zero-temperature phase diagram of the quantum antiferromagnetic Ising chain using tensor networks. All calculations in our study involving tensor networks were performed using the Uni10 library. uni10 The starting point of our calculations is to construct the matrix product operator (MPO) representations Scholl_Rev ; McCulloch of the Hamiltonian and observables. The Hamiltonian described in Eq. (6) for a chain of $L$ spins with periodic boundary conditions (PBC) can be written as an MPO with operator-valued matrices ${W}_{H_{p}}^{[i]}$ (where $H_{p}$ stands for the Hamiltonian with PBC) at sites $i$ in the bulk:

[TABLE]

and at the first and last sites:

[TABLE]

[TABLE]

On an open chain, the Hamiltonian (denoted by $H_{o}$) is encoded by the following MPO tensors:

[TABLE]

[TABLE]

[TABLE]

The Binder cumulant, Binder which involves the second and the fourth moments of the order parameter, is a commonly used quantity in numerical studies for determining the critical point and critical exponents. The calculation of the expectation value of a high-order moment in a tensor network state is in general not straightforward. Wei Here we introduce an MPO representation for the high-order moments of an observable, which allows an efficient computation of the Binder ratio.

Consider an $n$-th moment of some observable $O$ of the form

[TABLE]

where $i$ runs over all the sites in the system and $O_{i}$ is the on-site operator at $i$. Here we introduce a generic MPO form for $O^{n}$ with arbitrary $n$,

[TABLE]

where the edge matrices are given by

[TABLE]

and

[TABLE]

with the binomial coefficient $\binom{n}{k}$, and bulk matrices $W_{O^{n}}^{[i]}$ at sites $i=2\cdots L-1$ read

[TABLE]

Expressing an $n$-th moment in terms of this MPO, the bond dimension grows linearly with $n$. One can directly use this explicit MPO construction for any power $n$ without using MPO compression to reduce the bond dimensions, which is the advantage of this MPO representation over previously proposed constructions. McCulloch ; Generic_MPO

Using the density matrix renormalization group (DMRG) technique White_DMRG in terms of matrix product states McCulloch ; Scholl_Rev and the MPO representations described above, we determine the ground state of the chain with PBC in the absence of disorder. For a fixed $h$, we find the critical value of $\Gamma$ by studying the Binder cumulant $U$ of the staggered magnetization, defined as Sandvik_Rev

[TABLE]

where the staggered magnetization $m_{s}$ is given by

[TABLE]

and $\langle{\cdots}\rangle$ denotes the expectation value in the ground state. The curves of the Binder cumulant as a function of $\Gamma$ at fixed $h$ for different system sizes $L$ cross each other, as illustrated in Fig. 2. The critical values $\Gamma_{c}$ can then be extracted from the crossings of the Binder cumulant. In the inset of Fig. 2(b), we have rescaled the Binder cumulant for $h=1$ horizontally, multiplying $(\Gamma-\Gamma_{c})$ by $L^{1/\nu}$ with $\nu=1$, to achieve data collapse; the same scaling form applies for other critical points plotted in the phase diagram in Fig. 1 except for the classical multicritical point at $h=2$, where no crossing point is found in the curves of the Binder cumulant $U(\Gamma)$ for different system sizes. We have examined the accuracy of the DMRG results by comparing the locations of the critical points in Fig. 1 with those obtained by the continuous-time quantum Monte Carlo (QMC) algorithm. Excellent agreement can be observed between the results obtained from the two methods.

We conclude this section with the note that the MPO representation for high-order moments and cumulants introduced here circumvents tedious calculations involving sums of correlators $O_{i_{1}}O_{i_{2}}\cdots O_{i_{n}}$. This new MPO requires no additional truncation to represent an $n$-th moment with a bond dimension linear in $n$, and provides an efficient way to determine the critical point with great accuracy, also for the disordered system discussed in the next section.

IV Tree tensor network strong disorder renormalization group

The DMRG method can be applied to systems with quenched disorder, too. However, a practical implementation of the DMRG for disordered systems is computationally expensive because large bond dimensions and many DMRG sweeps are often required to avoid getting stuck in local minima. Hubbard1 ; Hubbard2

The Hamiltonian of the quantum spin model in Eq. (6) with randomness has an intrinsic separation in energy scales, which allows us to find its ground state using the SDRG technique. The basic strategy of the SDRG method, introduced in Ref. Ma, and Ref. Fisher_Ising, , is to find the largest term in the Hamiltonian successively and lock the spins associated with this term into their local ground state. For example, if the largest term is a field, $\Gamma$ or $h$, on spin $i$, then that spin is frozen in the $x$- or $z$-direction; if, on the other hand, the largest term in the Hamiltonian is an interaction $J_{ij}$, the two spins associated with this term are combined together into an antiferromagnetic cluster, in the state $\left|{\uparrow_{i}\downarrow_{j}}\right\rangle$ or $\left|{\downarrow_{i}\uparrow_{j}}\right\rangle$. The largest term in the Hamiltonian is thus effectively eliminated and energy corrections in terms of effective interactions or fields are generated using perturbation theory. This RG process yields an effective Hamiltonian with gradually fewer degrees of freedom and lower energy scale.

Although the SDRG procedure is approximate, it yields asymptotically exact results (up to logarithmic corrections SDRG_JQ ; XX_Log ) in the low-energy limit for the scaling behavior of systems at quantum critical points governed by an infinite-randomness fixed point, where the widths of the distributions of the logarithmic effective couplings diverge. To justify the accuracy of this approximate RG in general cases, one should carefully examine whether the condition of validity of perturbation theory is satisfied, that is, whether the largest energy scale to be decimated is much larger than those of the neighboring terms (the perturbation). This condition is often not satisfied in the early stage of RG when the disorder distribution is not very broad. One approach to refine the perturbative approximation is to include higher energy states, rather than only the lowest multiplet, at each RG step; Hikihara this extended SDRG method was further reformulated in terms of tree tensor networks in Ref. TSDRG, for the random-bond Heisenberg chain. Below we demonstrate that the tree tensor network SDRG is applicable to the random quantum Ising chain also with site disorder.

IV.1 Numerical scheme

For clarity of presentation, here we review the formulation of SDRG as a tree tensor network as proposed in Ref. TSDRG, .

A natural way to include the effect of excited states of local Hamiltonians in the RG process is to partition the system into blocks. Hikihara In each RG iteration, we compute the energy spectrum of each nearest-neighbor two-block Hamiltonian and identify the energy gap $\Delta_{\chi}$, which is measured as the difference between the highest energy of the $\chi$-lowest energy states that would be kept and the higher multiplets that would be discarded. We then combine the pair of blocks with the largest energy gap $\Delta_{\chi}$ into one block. The RG process is repeated until the chain is described by one block.

The SDRG procedure formulated as tree tensor networks can be described as follows TSDRG (see Fig. 3):

Decompose the chain into a set of $n$-site blocks and construct the block MPO tensors $W_{B}$.

• 2.

Obtain the energy spectrum of the two-block Hamiltonian for each pair of nearest-neighbor blocks.

• 3.

Merge the pair of blocks (say $B$ and $B^{\prime}$) with the largest energy gap $\Delta^{\text{max}}\equiv\max\{\Delta_{\chi}\}$ into a new block and contract the MPO tensors to form $W_{BB^{\prime}}$.

• 4.

Build a rank-3 isometry tensor $V$ with the $\chi$ lowest energy states of new block as column vectors; use $V$ and $V^{\dagger}$ to truncate $W_{BB^{\prime}}$ to $W_{\tilde{B}}$.

• 5.

Repeat steps 2-4 until the system is represented by one single block.

Finally, we diagonalize the Hamiltonian $H_{B_{f}}$ of last one block resulting from the RG process to obtain the eigenvalues, which represent the low-energy spectrum of original system. The ground-state expectation value of some observable can be obtained by contracting the MPO with the tree tensor network composed of the ground state of $H_{B_{f}}$ and the isometric tensors $V$ created during the RG.

IV.2 Results

Using the tensor network SDRG we investigate the $T=0$ phase diagram of the random-bond antiferromagnetic Ising chain in longitudinal and transverse magnetic fields, which to our knowledge has not been previously studied. We use open boundary conditions and the following distributions of the random bonds and random transverse fields:

[TABLE]

The width, $\Gamma_{0}$, of the field distribution is therefore the single control parameter in our simulations. We have achieved convergence for system sizes up to $L=128$ using four-site blocks ($n=4$) and $\chi=50$. In addition, for each case more than 10000 independent disorder realizations (samples) were used to obtain the disorder averages with sufficiently small error bars. In the following, we discuss the results for the case with zero longitudinal field ($h=0$) and with nonzero longitudinal field ($h\neq 0$), separately.

IV.2.1 Zero longitudinal field

We begin first with the random model in the absence of a longitudinal field. The location of the critical point for this case is determined by the relation given in Eq. (15) or Eq. (19). For the distributions of Eq. (35), the critical point occurs when $\Gamma_{0}=1$, which can be verified by using the Binder cumulant MPO described in Sec. III for finite-size systems. To obtain the disorder average of the Binder cumulant we first average the expectation values of the moments over samples and then evaluate the ratio, i.e.

[TABLE]

in this way we reduce errors in evaluating the ratio. Kawashima As shown in Fig. 5, the curves of the Binder cumulant for different system sizes cross each other, although the crossings exhibit some significant drift for the finite open chain considered here; the finite-size scaling plot in the inset of this figure shows that an extrapolation of the crossings to $L\to\infty$ based on data sets for sizes $L$ and $2L$ indeed gives the critical value $\Gamma_{0}^{c}=1$. We note that computing the Binder cumulant using the conventional SDRG method does not appear to have a straightforward implementation, which shows a clear advantage of the tensor network approach.

The random quantum critical point with $h=0$ is an infinite-randomness fixed point, where the dynamical exponent diverges $z\to\infty$, characterized by a stretched exponential scaling form of energy gaps (see Eq. (4)). We have determined the energy gap from the lowest-lying excitation, denoted by $\varepsilon_{1}$, of the final normalized block. The distribution of energy gaps at the critical point and a scaling plot using the scaling variable $x=\ln(\varepsilon_{1})/\sqrt{L}$ are shown in Fig. 6(a) and Fig. 6(b), respectively, where the scaling function agrees well with the analytical expression given in Ref. Fisher_Young, , which has a limiting form as $P(x)\approx 2\pi^{-1/2}\exp(-x^{2}/4)$ for large $\left|{x}\right|$ and approaches $P(x)\approx-2\pi x^{-3}\exp(-\pi^{2}/(4x^{2}))$ for small $\left|{x}\right|$.

We have also calculated the end-to-end correlation function $C_{1,L}$. Like the distribution of energy gaps, the distribution of the typical correlations $C_{1,L}$ (shown in Fig. 6(c)) at the critical point is extremely broad even on a logarithmic scale and its width grows with increasing $L$. Fig. 6(d) shows a scaling plot with the rescaled variable $x=\ln(C_{1,L})/\sqrt{L}$; here the scaling function is also found to be in good agreement with the analytic prediction: $P(x)=-x/2\exp(-x^{2}/4)$. Fisher_Young While the typical correlations decay exponentially with the distance $L$, the average value of $C_{1,L}$, shown in Fig.6(e), falls off according to $1/L$, Fisher_Young and is dominated by rare events in which the two end spins have strong correlations of order unity.

Away from the quantum critical point in the paramagnetic phase where $\Gamma_{0}>1$, the dynamical exponent $z$ for the whole region is also known from analytical results Z ; Igloi and is given by the root of Eq. (20); for the distribution of Eq. (35), the dynamical exponent is then determined by the equation

[TABLE]

Furthermore, deep in the paramagnetic phase $\delta\gg 1$ where the low-energy excitations are rare, the distributions of $k$th low-lying excitation energies $\varepsilon_{k}$ for a finite chain of size $L$ have the following scaling form Young_GS1996 ; Igloi_3z ; Frechet

[TABLE]

where the scaling function $\tilde{P}_{k}(u_{k})$, based on extreme-value statistics, Extreme is described by the Fréchet distribution given by Frechet

[TABLE]

in terms of the variable

[TABLE]

with a nonuniversal constant $\alpha$. In Fig. 7 the distributions of the first and second excitation energy of the final RG block are scaled according to Eq. (38) with the dynamical exponent $z\approx 1.3354$ obtained from Eq. (37). The large dynamical exponent $z>1$ is a signature of the quantum Griffiths singularity. In fact, all solutions of Eq. (37) for $\Gamma_{0}>1$ give $z>1$, the value of which increases with decreasing distance to the quantum critical point $\Gamma_{0}=1$.

IV.2.2 Nonzero longitudinal field

Now we turn to the random chain in the presence of a longitudinal field and with the rectangular distributions of Eq. (35) for the bonds $J_{i}$ and transverse fields $\Gamma_{i}$. Without loss of generality, we consider the case where the strength of the longitudinal field is a site-independent constant $h>0$.

Figure 8 shows the disorder-averaged Binder cumulant, $\overline{U}$, of the staggered magnetization as a function of the cut-off of the transverse field, ranged from $\Gamma_{0}=0$ to 1, at a small longitudinal field $h=0.02$. The curves for different system sizes have no crossing, giving evidence that the staggered magnetization does not show a phase transition. Furthermore, all values of $\overline{U}$ as well as the squared antiferromagnetic order parameter $\overline{\langle{m_{s}^{2}}\rangle}$ (shown in the inset) decrease with increasing $L$, indicating vanishing antiferromagnetic order in this region where the longitudinal field is finite but only slightly away from zero. We recall that the chain in the classical limit $\Gamma_{0}=0$ is antiferromagnetic only in the state which consists of alternating spins pointing in opposite directions; in the chain with bond randomness ($J_{i}\in(0,1]$ in our case), a finite magnetic field will break up the staggered order of the spins if some bonds are weaker than the field, which explains the destruction of the antiferromagnetic order in the whole regime where $h>0$, including the ”quantum” phase with finite transverse fields.

Next, we analyze disorder effects in the finite-$h$ paramagnetic phase near the infinite-randomness fixed point at $\Gamma_{0}=1$. In Fig. 9 the distribution of first energy gaps at $\Gamma_{0}=1,\,h=0.02$ for different system sizes is shown. Here the distribution of $\ln(\varepsilon_{1})$ looks similar to the distribution in the Griffiths phase at $\Gamma_{0}>1,\,h=0$ (see Fig. 7); the low-energy tail of the distribution for different sizes is just shifted horizontally relative to each other and shows power-law behavior. According to Eq. (39), the power of the low-energy tail of $P(\ln(\varepsilon_{1}))$ is given by $1/z$ (here $k=1$). We estimate $z\simeq 0.8$ from the slope of the tails in Fig. 9, and find good agreement between the scaling function, in terms of the scaling variable $\varepsilon_{1}L^{z}$, and the Fréchet distribution described in Eq. (39). The power-law tail and the Fréchet distribution of lowest gaps show that these low-energy excitations are localized in this phase. The small dynamical exponent $z<1$ does not lead to divergence of the linear susceptibility, although it can produce divergence in the non-linear susceptibility at $T=0$ when $z>1/3$. Young_GS1996 ; Igloi_3z In comparison with the situation in the zero-$h$ paramagnetic phase, the dynamical exponent here in the vicinity of the quantum critical point is much smaller. The strongly reduced Griffiths effects can be understood as follows. A rare region in the equivalent classical system corresponds to a quasi-one-dimensional ferromagnetic Ising model along the time-like direction. An applied longitudinal field will align the spins in the rod-like rare region and suppress fluctuations of the order parameter, so that the rare region becomes essentially static, leading to a suppression of quantum Griffiths effects.

V Discussion

Using the DMRG and the SDRG algorithms based on matrix product operators, we have explored the zero-temperature phases of the antiferromagnetic Ising chain in longitudinal ($h$) and transverse ($\Gamma$) magnetic fields. We have introduced a new matrix product operator technique for calculating expectation values of high-order moments, which allows us to compute the Binder cumulant of the order parameter efficiently and locate the quantum critical points with high accuracy, both in the clean system and in the system with disorder.

Quenched disorder and rare regions of local order give rise to drastic effects on the zero-temperature phases. The random quantum critical point triggered by the transverse field in the absence of the longitudinal component is of unconventional type governed by an infinite-randomness fixed point; also pronounced Griffiths effects exist in the zero-$h$ paramagnetic phase, characterized by a large dynamical exponent ($z>1$) in the entire region.

We have unambiguously found that the aforementioned infinite-randomness quantum critical point in the disordered system are destroyed by the longitudinal component of the applied field. A schematic phase diagram of the random chain is shown in Fig. 10. The chain with bond randomness exhibits no antiferromagnetic order even in the classical limit (with zero transverse component of the field) when the applied field disturbs the perfect staggered spin configuration by flipping any spin over; this explains the destruction of antiferromagnetic order and quantum criticality in the random chain with finite $h$ that we consider. Furthermore, the longitudinal field hampers the quantum tunneling of rare ordered regions by pinning, leading to the suppression of quantum Griffiths behavior in the finite $h$ regime. Similar phenomena also occur in metallic Ising magnets. Millis_PRL ; Millis_PRB ; Vojta_rounding

The model we consider in this paper is interesting in many aspects. First, the model in the presence of a longitudinal field is nonintegrable even in the clean case, offering a playground to explore how dynamics depends on the integrability and the nonintegrability. clean_dyn ; clean_dyn_f Second, the model with bond-randomness has unconventional ground-state properties as discussed in this paper, and due to its strong disorder effects the model provides a theoretical paradigm for the problem of localization. dyn2012 ; Altman2014 ; RSRG-X ; Imbrie Furthermore, the model is experimentally realizable using ultracold spinless bosonic atoms confined in a tilted optical lattice with a magnetic field gradient along one direction, where each Ising spin is encoded in the motional degree of freedom of a single atom. cold_atom In such an ultracold atomic system, the longitudinal and transverse field are represented by the lattice tilt and tunneling of the atom, respectively; the spin-spin interaction arises from a nearest-neighbor constraint in the site occupation number. cold_atom ; cold_atom1 The paramagnetic-antiferromagnetic transition in the spin model can then be mapped to a transition between a phase with singly occupied sites at small tilt and a density wave phase at large tilt, which is detectable via occupation measurements using high-resolution imaging or noise correlation measurements. cold_atom ; noise_cr Disorder (spatial inhomogeneity) is naturally present in optical lattices and a near-homogeneous system can be realized by selecting a smaller sample size; cold_atom it is then possible to experimentally investigate effects of controlled disorder on quantum phase transitions in ultracold atomic lattice gases.

The higher-dimensional random transverse-field Ising antiferromagnets have similar zero-temperature phases to the one-dimensional case. Infinite-randomness fixed points accompanied by pronounced quantum Griffiths singularities also govern the low-temperature physics in higher dimensions with discrete Ising symmetry. Motrunich ; Lin_2D ; Istvan_2D In the presence of a longitudinal field, the quantum dynamics of rare ordered regions and quantum Griffiths effects are expected to be suppressed, too. In principle, there is no constraint on the implementation of tree tensor network SDRG in two or more space dimensions and with different types of symmetry (such as SU(2) symmetry), as long as disorder effects are relevant. Damle ; Bilayer ; Refael_2D Also the MPO technique proposed in this paper for calculating higher moments can be applied in high dimensions. A potential difficulty in implementing high-dimensional SDRG is to deal with a large number of block-block couplings resulting from fusion of blocks under the action of the RG. Nevertheless, for systems governed by strong disorder, it should be possible to find a criterion for discarding some block couplings without affecting the low-lying energy spectrum.

Acknowledgements.

We would like to thank Anders Sandvik, Matthias Vojta, and Frank Pollmann for useful discussions. Y.C.L. is grateful to Z.-R. Hsu for previous collaboration. This work was supported by the Ministry of Science and Technology (MOST) of Taiwan under Grants No. 105-2112-M-002-023-MY3, 105-2112-M-004-002, 104-2112-M-002 -022 -MY3. We also acknowledge support from the NCTS.