Nonequilibrium-relaxation approach to quantum phase transitions: Nontrivial critical relaxation in cluster-update quantum Monte Carlo

TL;DR

This paper demonstrates that the nonequilibrium relaxation (NER) method reveals nontrivial critical relaxation behavior in quantum phase transitions analyzed with cluster-update quantum Monte Carlo, challenging previous assumptions about relaxation speed.

Contribution

The study shows that NER in quantum systems with cluster-update QMC exhibits stretched-exponential relaxation, enabling new analysis of quantum critical dynamics.

Findings

NER exhibits stretched-exponential relaxation in quantum phase transitions.

Cluster-update QMC allows effective NER analysis despite rapid relaxation.

Application to the 2D Heisenberg model confirms the method's utility.

Abstract

Although the nonequilibrium relaxation (NER) method has been widely used in Monte Carlo studies on phase transitions in classical spin systems, such studies have been quite limited in quantum phase transitions. The reason is that relaxation process based on cluster-update quantum Monte Carlo (QMC) algorithms, which are now standards in Monte Carlo studies on quantum systems, has been considered "too fast" for such analyses. Recently the present authors revealed that the NER process in classical spin systems based on cluster-update algorithms is characterized by the stretched-exponential critical relaxation, rather than the conventional power-law one in local-update algorithms. In the present article we show that this is also the case in quantum phase transitions analyzed with the cluster-update QMC, and that advantages of NER analyses are available. As the simplest example of isotropic quantum spin models which exhibit quantum phase transitions, we investigate the N\'eel-dimer quantum phase transition in the two-dimensional $S=1/2$ columnar-dimerized antiferromagnetic Heisenberg model with the continuous-time loop algorithm.

Equations

$$\mathcal{H}$$

$$J_{ij}$$

$$m_{\mathrm{N}}\equiv \frac{1}{N}\sum _i(-1{)}^i{S}_i^z,$$

$$\mathcal{H}=\sum _{⟨ij⟩\in n.n.}{J}_{ij}\left(-S_i^x{S}_j^x-S_i^y{S}_j^y+S_i^z{S}_j^z\right),S=1/2,$$

$$⟨|m_{\mathrm{N}}(t,L)|⟩\sim L^{-d/2}\exp \left(+c_m{t}^{\sigma }\right),$$

$$⟨|m_{\mathrm{N}}(t,L)|⟩L^{\beta /\nu }\sim f_m(c_m{t}^{\sigma }-\ln L^{d/2-\beta /\nu }),$$

$${\delta }_{\mathrm{c}}=0.90947(5),$$

$$\beta /\nu =0.514(2),\sigma =0.502(9),c_m=0.426(17).$$

Full text

Introduction. Dynamical Monte Carlo (MC) simulations have been widely utilized in statistical-mechanical studies on classical spin systems, where the Boltzmann weight for a MC flip is determined locally. In quantum spin systems, the Boltzmann weight cannot be determined locally unless diagonalizing the Hamiltonian matrix in principle, and treatable system sizes are strictly reduced. This difficulty was overcome by introducing the Suzuki-Trotter decomposition Trotter59 ; Suzuki76 , where noncommutable exponential operators are approximately divided into $n$ Trotter layers, and the original system is reproduced in the $n\to\infty$ limit. This procedure can be regarded as introduction of an extra imaginary-time axis.

In the original formulation of quantum Monte Carlo (QMC) Suzuki76 , the $n\to\infty$ limit is taken by numerical extrapolation with several finite-$n$ systems, and nontrivial global-update processes should be introduced additionally by hands in order to keep ergodicity. Such a complicated procedure was simplified by the continuous-time QMC algorithm Beard96 , where the $n\to\infty$ limit is taken as a part of formulation without numerical extrapolation. This QMC algorithm should be coupled with cluster-update QMC algorithms such as the loop algorithm Evertz93 , worm algorithm Prokovev98 , or directed-loop algorithm Syljuasen02 , in which ergodicity is ensured within the formulation.

The nonequilibrium relaxation (NER) method NERrev is one of the MC approaches to investigate phase transitions against the critical slowing down. In contrast to other approaches such as the cluster algorithms SW ; Wolff and extended ensemble methods Berg ; Hukushima ; WL , the critical slowing down is not avoided but rather utilized in this scheme for the evaluation of critical phenomena. The critical point is estimated from the power-law behavior of physical quantities expected there, and the critical exponents are evaluated from the exponents of such relaxation behaviors. This relaxation process is usually terminated much earlier than arrival at equilibrium. In order to avoid artifacts originating from biased initial states, completely-ordered or completely-disordered states are usually chosen as initial ones for the relaxation process.

In equilibrium MC simulations, thermal averaging is taken during long-time measurement after discarding equilibration data. In NER calculations, such averaging is replaced by that for independent random-number sequences (RNS), and all the numerical data are utilized.

Several attempts to generalize the NER scheme to above-mentioned modern QMC algorithms had not been successful so far, and cluster-update QMC algorithms had been considered “too fast” for NER. Succeeded examples of the application of NER to QMC were based on a world-line local-update QMC algorithm with finite Trotter layers Nonomura98 or a continuous-time QMC algorithm only along the imaginary-time axis Nakamura03 , both of which simply slowed down relaxation to fit for power-law NER.

Recently the present authors numerically found that the early-time nonequilibrium critical relaxation in cluster algorithms is described by the stretched-exponential simulation-time dependence, not by the power-law one in various classical spin systems Nonomura14 ; Nonomura15 ; Nonomura16 . NER framework for cluster algorithms can be constructed on the basis of this relaxation formula. Quite recently the present authors derived this relaxation formula phenomenologically in the Ising models in the Swendsen-Wang algorithm Tomita18 .

In the present article, we propose an NER analysis of quantum phase transitions based on the continuous-time loop algorithm. Although this algorithm is based on cluster update, it is not trivial if the stretched-exponential critical relaxation is observed in the present one-dimensional loop clusters, which is different from bulky ones in cluster algorithms in classical spin systems.

Model and method. In the present article we investigate the two-dimensional $S=1/2$ columnar-dimerized antiferromagnetic Heisenberg model on a square lattice,

[TABLE]

with the nearest-neighbor interaction and the columnar dimer pattern shown in Fig. 1. For $\delta=0$ (the uniform case), this model has the Néel order (reduced due to quantum fluctuations) in the ground state. For large enough $\delta$, singlet pairs on the dimerized bonds wipe out the Néel order. There exists a critical point $\delta_{\rm c}$ between the two cases, and this quantum phase transition is of conventional second order, because the dimerization does not accompany spontaneous symmetry breaking.

This is the simplest isotropic quantum spin model with quantum phase transition without frustration. This model has been intensively studied with QMC Matsumoto01 ; Wenzel09 ; Yasuda13 , and its universality class has been considered to be the same as that of the three-dimensional (3D) classical Heisenberg model. Here we analyze the early-time relaxation behavior of this model with the continuous-time loop algorithm. In the NER analysis, choice of the initial state is crucial. In the 3D classical Heisenberg model, we found that the ordering process from the perfectly-disordered state gives much smaller deviation from the stretched-exponential critical relaxation than the decaying process from the perfectly-ordered state Nonomura16 .

In classical spin systems, the perfectly-disordered state can be one of the states for $T\to\infty$. However, when NER calculations are started from such states in quantum phase transitions, they quickly converge to a classical ordered state corresponding to the basis used in the MC algorithms, and then converge to equilibrium at the quantum critical point. In the present study, we utilize the continuous-time loop algorithm formulated with the Ising basis, and measure physical quantities based on the classical Néel order. MC time evolution is based on the standard loop algorithm; loop clusters are generated in the whole system and each cluster is flipped with $50$% probability, similarly to the Swendsen-Wang algorithm in classical spin systems SW . Then, we start from the isolated dimer state, in which only singlet pairs are on the dimerized bonds, and it becomes the ground state in the $\delta\to\infty$ limit. Since the parameter $\delta$ plays a role of temperature in this quantum phase transition, this process corresponds to NER from a perfectly-disordered state.

In the present study, we analyze the absolute value of the Néel order measured on the initial Trotter layer,

[TABLE]

with abbreviations $i\equiv(i_{x},i_{y})$ and $(-1)^{i}\equiv(-1)^{i_{x}+i_{y}}$. In QMC simulations of antiferromagnets on bipartite lattices, the original Hamiltonian (1) is transformed to

[TABLE]

with the spin rotation $S_{i}^{x}\to-S_{i}^{x}$, $S_{i}^{y}\to-S_{i}^{y}$, $S_{i}^{z}\to S_{i}^{z}$ on one of the sublattices in order to refrain from the negative-sign problem. Then, the singlet state on each dimerized bond is transformed from $(|\uparrow\downarrow\rangle-|\downarrow\uparrow\rangle)/\sqrt{2}$ to $(|\uparrow\downarrow\rangle+|\downarrow\uparrow\rangle)/\sqrt{2}$ in actual simulations, and typical states consisting of the isolated dimer state in the transformed system with the Ising basis are constructed as follows:

Consider the system with $J_{ij}=(1+\delta)J$ on dimerized bonds and $J_{ij}=0$ elsewhere. 2. 2.

Assign the basis $|\uparrow\downarrow\rangle$ or $|\downarrow\uparrow\rangle$ on each dimerized bond randomly but to keep the Néel order (5) vanishing on the initial Trotter layer, on which physical quantities are measured. (When the RNS average is taken, this initial configuration is also changed.) 3. 3.

Insert gates along the imaginary-time direction QMCrev on each dimerized bond with the probability corresponding to $J_{ij}=(1+\delta)J$. The number of gates should be even on each dimerized bond in order to satisfy the periodic boundary condition along the imaginary-time direction. 4. 4.

Flip the basis $|\uparrow\downarrow\rangle$ to $|\downarrow\uparrow\rangle$ and vice versa at each gate. (When the steps 3 and 4 are skipped, we have one of the classical perfectly-disordered states.)

When the stretched-exponential critical relaxation holds, early-time behavior of the absolute value of the Néel order at the quantum critical point $\delta_{\rm c}$ is given by

[TABLE]

with the RNS average $\langle\cdots\rangle$, simulation time $t$, linear size $L_{x}=L_{y}=L$, spatial dimension $d=2$ (this size dependence originates from the normalized random-walk growth of spin clusters, and here we do not take the summation along the imaginary-time direction), quantity-dependent coefficient $c_{m}$ and relaxation exponent $\sigma$ ($0<\sigma<1$) common in all the physical quantities. Combining this formula with the equilibrium finite-size scaling relation, namely $\langle|m_{\rm N}(t=\infty,L)|\rangle\sim L^{-\beta/\nu}$, we arrive at the following nonequilibrium-to-equilibrium scaling relation,

[TABLE]

with a quantity-dependent scaling function $f_{m}$ Nonomura14 ; Nonomura16 .

Numerical results. We analyze the absolute value of the Néel order of the model (1) and evaluate the quantum critical point $\delta_{\rm c}$, critical exponent $\beta/\nu$ and relaxation exponent $\sigma$. As an example, this quantity multiplied with $L^{d/2}$ at $\delta=0.90947$ (as will be shown below, the most probable value of $\delta_{\rm c}$) is plotted versus simulation time in Fig. 2 for $L=56$ ($2.56\times 10^{6}$ RNS are averaged), $80$ ($2.56\times 10^{6}$ RNS), $112$ ($1.28\times 10^{6}$ RNS) and $160$ ($0.64\times 10^{6}$ RNS), where the system size along the imaginary-time axis is taken the same as those along the real axes, namely $\beta\equiv 1/(k_{\rm B}T)=L$. As expected, this quantity is scaled well with Eq. (7) initially, and tends to be away from this equation and to saturate as $t$ increases.

The data in Fig. 2 are scaled with Eq. (8) in Fig. 3 with $\beta/\nu=0.514(1)$, $\sigma=0.502(8)$ and $c_{m}=0.423(14)$. The exponent $\beta/\nu$ is evaluated from the scaling behavior in the vicinity of equilibrium, which little depends on the values of $\sigma$ and $c_{m}$. Details of this evaluation will be explained later together with the estimation of $\delta_{\rm c}$. The values of $\sigma$ and $c_{m}$ are evaluated so as to minimize the mutual residue of the data with the estimate of $\beta/\nu$. In Ref. Nonomura16 , we also used the criterion that the early-time behavior is described by the stretched-exponential relaxation formula (7), and the slope of the scaled data should be unity in the region $c_{m}t^{\sigma}-\ln L^{d/2-\beta/\nu}\lesssim-0.5$ in Fig. 3. Here we do not use this criterion, because validity of this formula cannot be fully justified comment .

In order to evaluate the quantum critical point $\delta_{\rm c}$, the data in the vicinity of equilibrium for $L=56$ (green symbols), $80$ (pink symbols), $112$ (blue symbols) and $160$ (red symbols) at $\delta=0.90952$ (open squares with $\beta/\nu=0.516$), $0.90951$ (full squares with $\beta/\nu=0.515$), $0.90947$ (open circles with $\beta/\nu=0.514$), $0.90943$ (full triangles with $\beta/\nu=0.513$) and $0.90942$ (open triangles with $\beta/\nu=0.512$) are plotted in Fig. 4 with $\sigma=0.502$ and $c_{m}=0.423$ estimated at $\delta=0.90947$ in Fig. 3. The data except for $\delta=0.90947$ are slightly shifted upwards or downwards for clear visualization. Although the data at $\delta=0.90951$ and $0.90943$ are still scaled on a single curve with $\beta=0.515(1)$, $\sigma=0.501(8)$, $c_{m}=0.427(16)$ and $\beta=0.513(1)$, $\sigma=0.504(7)$, $c_{m}=0.421(13)$, respectively (see Fig. 3; these data are shifted upwards or downwards for clear visualization), such a behavior is not observed at $\delta=0.90952$ and $0.90942$. Here the data for $L=80$ and $112$ are tuned to be scaled on a single curve. Then, deviation of the data for $L=56$ and $80$ and that for $L=112$ and $160$ are in opposite directions, and all the data cannot be scaled anymore. In order to include all the results, our final estimates are given by

[TABLE]

These estimates are consistent with the previous ones for the same model based on the equilibrium QMC simulations Yasuda13 , $\delta_{\rm c}=0.90947(3)$ and $\beta/\nu=0.513(9)$, and with ours for the 3D classical Heisenberg model based on the cluster NER Nonomura16 , $\beta/\nu=0.515(5)$ and $\sigma\approx 1/2$. Although the evaluation of $\delta_{\rm c}$ based on the deviation of data in Fig. 4 is rather subtle, that of $\beta/\nu$ based on the scaling formula (8) is promising, where wide scaling region results in high precision in comparison with a simple power-law fitting of equilibrium data Yasuda13 , though accuracy of the estimation is still not comparable to that of a detailed study on critical phenonena in the 3D classical Heisenberg model Campostrini02 . We made a similar analysis on the staggered susceptibility, and obtained the critical exponent $\gamma/\nu$ consistent with the hyperscaling relation suscexp .

Summary and discussion. In the present article we generalized the cluster nonequilibrium relaxation (NER) scheme to quantum phase transitions. Since modern quantum Monte Carlo algorithms such as the loop algorithm are based on cluster updates, the present scheme is indispensable for the NER analysis of quantum systems. As an example, we considered the Néel-dimer quantum phase transition in the two-dimensional $S=1/2$ columnar-dimerized antiferromagnetic Heisenberg model on a square lattice. This model is the simplest isotropic quantum spin system to exhibit a quantum phase transition with respect to the strength of dimerization $\delta$, and belongs to the universality class of the 3D classical Heisenberg model.

In the present study, numerical calculations were based on the continuous-time loop algorithm with the Ising basis and started from the isolated dimer state. Although we have numerically and theoretically clarified that physical quantities at the critical point show the stretched-exponential relaxation behavior in the early-time relaxation in cluster algorithms in classical spin systems, this behavior is not trivial in the present case, because one-dimensional loop clusters are geometrically different from the bulky ones in the Swendsen-Wang and Wolff algorithms. Then, the critical point $\delta_{\rm c}$, critical exponent $\beta/\nu$ and relaxation exponent $\sigma$ were estimated from the nonequilibrium-to-equilibrium scaling plot of the absolute value of the Néel order as in our previous study on the 3D classical Heisenberg model. The present estimates $\delta_{\rm c}=0.90947(5)$, $\beta/\nu=0.514(2)$ and $\sigma=0.502(9)$ are comparable to previous studies (even more precise than previous equilibrium QMC simulations in $\beta/\nu$). Consistency with the 3D classical Heisenberg model holds not only for $\beta/\nu$ but also for $\sigma$. These results reveal that the cluster NER scheme can be generalized to quantum phase transitions on the basis of the continuous-time loop algorithm.

As shown in the present article, large numbers of random-number sequences should be averaged for accurate data in the NER analysis instead of long-time measurements in the equilibrium Monte Carlo analysis, and such averaging can be replaced by sample averaging in random systems. This fact indicates that numerical efforts in random systems may not be so different from pure systems, which is the essential merit of the present scheme. Studies along this direction is now in progress.

Acknowledgments. Y. N. thanks K. Harada for helpful comments. The present study was supported by JSPS KAKENHI Grant No. 16K05493. The random-number generator MT19937 MT was used for numerical calculations. Most calculations were performed on the Numerical Materials Simulator at National Institute for Materials Science.