{'fraction_non_alphanumeric': 0.04147955566060033, 'fraction_numerical': 0.016662727487591587, 'mean_word_length': 4.290715848702719, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo simulation with aspect ratio optimization: Anomalous anisotropic scaling in dimerized antiferromagnet

We present a method that optimizes the aspect ratio of a spatially anisotropic quantum lattice model during the quantum Monte Carlo simulation, and realizes the virtually isotropic lattice automatically. The anisotropy is removed by using the Robbins-Monro algorithm based on the correlation length in each direction. The method allows for comparing directly the value of critical amplitude among different anisotropic models, and identifying the universality more precisely. We apply our method to the staggered dimer antiferromagnetic Heisenberg model and demonstrate that the apparent non-universal behavior is attributed mainly to the strong size correction of the effective aspect ratio due to the existence of the cubic interaction.

We present a method that optimizes the aspect ratio of a spatially anisotropic quantum lattice model during the quantum Monte Carlo simulation, and realizes the virtually isotropic lattice automatically. The anisotropy is removed by using the Robbins-Monro algorithm based on the correlation length in each direction. The method allows for comparing directly the value of critical amplitude among different anisotropic models, and identifying the universality more precisely. We apply our method to the staggered dimer antiferromagnetic Heisenberg model and demonstrate that the apparent non-universal behavior is attributed mainly to the strong size correction of the effective aspect ratio due to the existence of the cubic interaction. Quantum phase transitions [1] are the transitions between ground states with different symmetries. They are triggered at absolute zero temperature by the change of a parameter that controls the strength of quantum fluctuations. A quantum phase transition in d dimensions, if it is of second order, is widely considered to belong to the same universality class as the finite-temperature phase transition of the (d + 1)-dimensional classical system with the same symmetry. As a concrete example, let us consider the columnar dimer model, a spin-1/2 dimerized Heisenberg antiferromagnet. The Hamiltonian of this system is written as where S i is the spin-1/2 operator at site i, J and J ′ are the positive (antiferromagnetic) coupling constants, and A (B) denotes the set of the pairs of sites connected by the thin (thick) bond shown in Fig. 1(a). It undergoes a second-order quantum phase transition at some critical value of J ′ /J. By increasing J ′ /J from 1, the ground state changes from the Néel ordered state to the dimer state [2,3], and the critical exponents of this transition are known to coincide with those of the three-dimensional (3D) classical Heisenberg [O(3)] universality [4][5][6][7].
In the meantime, however, there are a number of intensive researches that aim to find novel critical phenomena that have no classical counterparts. As for the lattice spin models, the staggered dimer model has been examined as a candidate that might exhibit such phenomena. Its Hamiltonian is the same as Eq. (1), but the only difference between the columnar and staggered dimer models is the configuration of the dimerization pattern [ Fig. 1(b)]. In Ref. 8, based on the results of the quantum Monte Carlo simulation the authors claim that the critical exponents of the staggered dimer model are different from those of the other dimerized models, such as the columnar dimer model. In the more recent study [9], on the other hand, it is pointed out that the simulation results become consistent with the conventional O(3) universality by carefully choosing the aspect ratio of the lattice. It is further discussed that the models with specific dimerized patterns, including the staggered dimer model, could exhibit apparent unconventional critical phenomena due to the presence of the weakly irrelevant "cubic term" [10,11], though the relation between such cubic term and the numerically observed large corrections to scaling is yet unclear.
Usually, quantum Monte Carlo simulations of quantum critical phenomena are carried out with a cubic geometry, e.g., L x : L y : L τ = 1 : 1 : 1 in (2+1) dimensions. Here, we denote the linear length of the system in α-direction as L α (α = x, y, or τ ). The length in τ -direction means the inverse temperature, 1/T . One should be noticed that in the case where the system has spatially anisotropic interactions, the correlation lengths, ξ α 's, generally depend on their directions. In such a case, it is natural to introduce the virtual aspect ratio, R −1 x : R −1 y : R −1 τ , by using the effective system linear length defined as the inverse of relative correlation length, R α ≡ ξ α /L α .
In principle, results of the finite-size scaling analysis do not depend on the aspect ratio chosen for a series of simulations as long as sufficiently large lattices are simulated. In practice, however, one can simulate effectively larger systems with minimal computational cost by tuning the aspect ratio so that the system becomes virtually isotropic [2], i.e., R −1 x : R −1 y : R −1 τ ≈ 1 : 1 : 1 instead of L x : L y : L τ = 1 : 1 : 1. By adopting such a geometry, one can examine the universality of quantum critical phenomena even more closely, because not only the critical exponents but also the scaling function of quantities with vanishing scaling dimension becomes universal under the virtually cubic geometry. For example, let us consider the Binder ratio that is defined as where m z is the summation of the z-component of spins.
The y-intercept of the scaling function of such a quantity, called the critical amplitude, is a useful index for identifying the universality class, because such an amplitude can usually be calculated with higher accuracy than the critical exponents [12][13][14][15]. Last of all, for the case where the virtual aspect ratio of the system changes gradually as increasing the system size, additional care must be taken, since it might cause the strong corrections to scaling. As we will see below, the staggered dimer model is the very case, and the result reported in Ref. 8 is an artifact due to the strong influence of the non-trivial system size dependence of the virtual aspect ratio.
Tuning the aspect ratio by hand is generally a difficult and complicated task. In the present paper, we propose an algorithm that optimizes the aspect ratio during the Monte Carlo simulation automatically. This method enables one to make R α the same value in all directions in order to simulate the virtually isotropic system. It is also possible to search for the quantum critical point J ′ = J ′ c (we set J = 1 without loss of generality). For example, in a (2+1)-dimensional system, we solve the equation R x = R y = R τ = R 0 for each fixed L x , where R 0 > 0 is an arbitrarily chosen constant. In this case, we have three parameters, L y , L τ , and J ′ , to be determined and three equations, which imply that we can determine all parameters. Note that the reason why R 0 is arbitrarily chosen is that in the thermodynamic limit J ′ smaller (larger) than the critical value J ′ c gives the limit R α → ∞ (0). Thus any positive finite constant R 0 leads In the present simulation, we use R 0 = 0.5664, which is an estimate of the critical amplitude of the 3D classical isotropic Heisenberg model at the critical point (1/T c = 0.693035(37) [16]) by the Wolff algorithm [17]. This is because if the model considered here belongs to the O(3) universality class, this choice of R 0 can reduce the corrections to scaling.
In the present simulation, we adopt the loop algorithm based on the continuous-time path integral representation [18,19]. The correlation length in each direction is evaluated by the second-moment method [20,21] as where C( q) is the imaginary-time dynamical structure factor of the z-component of the magnetization at wavevector q, q 0 = (π, π, 0), and δ q α = (2π/L x , 0, 0), (0, 2π/L y , 0), and (0, 0, 2π/L τ ) for α = x, y, and τ , respectively. As the estimates fluctuate statistically, the naive Newton method becomes unstable and does not work well for the present purpose. Instead, we employ a more robust method, the Robbins-Monro algorithm [22,23], from the field of machine learning. This algorithm enables us to estimate the zero of the regression function with probability unity for an observable with a finite variance. Let z(θ) be a random variable parameterized by θ with mean f (θ) and a finite variance. We assume that the regression function f (θ) increases monotonically as increasing θ and has a zero, θ = θ * . The zero θ * can be obtained by repeating the Robbins-Monro procedure: where n = 0, 1, 2, · · · is the iteration step, α some positive constant, and θ (n) the estimate of θ * at step n. It is proved that θ (n) converges to θ * with probability one [22]. In Eq. (4), the feedback coefficient, α/n, is chosen so as to satisfy (i) the summation about n diverges, while (ii) the sum of squares converges to a finite value. Condition (i) ensures that θ (n) can reach θ * irrespective of the initial value θ (0) , and condition (ii) keeps the accumulated variance to be finite. Although the choice of α affects the convergence rate and the fluctuation of θ (n) around θ * , the convergence is guaranteed as long as α is positive and finite. Extension of this algorithm to higher dimensions is straightforward [24]. At each Robbins-Monro step (RMS), we evaluate the correlation lengths by 500 Monte Carlo steps with fixed parameters, J ′ , L y , and L τ . Then, the parameters are updated according to Eq. (4). We choose the regression functions as R x + R y + R τ − 3R 0 , R y − R x , and R τ − R x , and the feedback parameter α as 2, 500, and 500 for the parameters, J ′ , L y , and L τ , respectively. Here, it should be noted that the quantum Monte Carlo simulation can be carried out only for integral values of L y , since L y is the number of lattice sites in y-direction. Furthermore, L y should be even in order to avoid negative signs. For a non-even value of L y we adopt an arithmetic mean of two independent Monte Carlo estimates for L (1) y = W and L (2) y = W +2, where W is the maximum even integer not in excess of L y . Fig. 2 shows the convergence of the parameters for staggered dimer model with L x = 64. One can see all R α 's oscillate coherently with decreasing amplitude from 20 to 50 RMS, reflecting the convergence of L y and L τ . In that region, J ′ still shows oscillatory behavior. This means the gain for J ′ in the Robbins-Monro procedure is too large. Although one might be able to optimize the gain to increase the convergence rate, such tuning only affects the number of steps before convergence. For the present example (shown in Fig. 2), the average is taken over 100 RMS with discarding first 300 RMS as the thermalization. We have to discard more steps when the system size becomes larger. For example, we discard 650 RMS in the case of L x = 384.
After estimating J ′ c (L x ) for each L x , then we extrapolate it in the thermodynamic limit L x → ∞. We found is a good fitting function, where a, b, and J ′ c (∞) are fitting parameters. We used the data with the system size ranging from L x = 16 to L x = 384 for the staggered dimer model and to L x = 192 for the columnar dimer model. From this fitting, we conclude that the critical points J ′ c = 2.51941(2) and 1.90947(3) for the staggered and columnar dimer model, respectively. They are consistent with those in the past literature, J ′ c = 2.5196 (2) for the staggered dimer model [8] and 1.9096 (2) for the columnar dimer model [2,3], but our estimates are much more precise.
Let us move on to the calculation of the Binder ratio at the critical point. For classical ferromagnetic Heisenberg models, the Binder ratio is defined as Eq. (2), where m z = i S z i . The quantum counterpart is defined in the same way except that m z s is not the simple staggered magnetization (the Néel order parameter) but the integrated staggered magnetization, which is written as as it reduces to i S z i when one maps quantum systems into classical systems. Fig. 3 shows the size dependence of the Binder ratio at the critical point calculated under the dynamic controlling of anisotropy. Fitting the data with quadratic functions of 1/L x gives the critical amplitudes We conclude that these three models share the same critical Binder ratio and thus they belong to the same 3D O(3) universality class. Note that these values are clearly different from that of other universalities, e.g., Q c ≃ 0.62 for the 3D Ising model.
We estimated the critical exponents as well. The absolute value of the Néel order parameter |m z s | and the staggered susceptibility χ behaves as |m z s | ∼ Λ −β/ν and χ ∼ Λ γ/ν at the critical point, respectively. Here, Λ is the characteristic length of the system. In the present simulation, L y and L τ are optimally tuned so that the system should be virtually isotropic, thus we define Λ as Λ ≡ (L x L y L τ ) 1/3 . For the magnetization, fitting for the largest four data gives β/ν = 0.522(3) and 0.513(9) for the staggered and columnar dimer model, respectively. Both coincide with the standard O(3) value, 0.518(1) [16,25] and exclude the possibility of β/ν = 0.545(4) from Ref. 8. For the staggered susceptibility, our estimate is γ/ν = 1.970(5) and 1.983 (13) for the staggered and columnar dimer model, respectively. They are consistent with γ/ν = 1.9750(35) in Ref. 16, but is slightly bigger than 1.963 (2) in Ref. 25. In any case, we can see no evidence that the columnar and staggered dimer model belong to different universality classes.
The apparent unconventional critical behavior observed in Ref. 8 is attributed to the strong size dependence of the virtual aspect ratio in the staggered dimer model. Fig. 5 shows the the optimized aspect ratio at the critical point J ′ = J ′ c as a function of the system size. As one can see, L y /L x of the staggered dimer model exhibits a strong and non-monotonic size dependence, though the other ratios converge to finite values quite rapidly. This fact means that in the conventional simulation with a fixed aspect ratio L x : L y : L τ , the virtual aspect ratio gradually changes as the system size increases that introduces strong corrections to scaling in the staggered dimer model.
Let us analyze this unconventional behavior in a different way. From the quantum field theory, the low-energy effective action for the present models is described by the standard where the constants which determine the length scale are left explicitly as c α . Along this action, only for the staggered dimer model, it is discussed that the action has the cubic term S 3 ∼ γ φ · (∂ x φ × ∂ τ φ) [10]. It is straightforward to show that rotating in the x-τ plane effectively pushes S 3 into the kinetic term of the action S 0 , resulting that the coefficients are renormalized as c , where [ · ] denotes the scaling dimension of ·, whereas c y remains unchanged. This simple dimensional analysis explains the behavior of the optimized aspect ratio of the staggered dimer model observed in Fig. 5, where c y /c x suffers from large corrections but c τ /c x does not. We assume the form of finite size correction as L y /L x ≈ a + bL −ω x . The exponent of correction, ω = 0.57(4), obtained by least-squares fitting is consistent with the above estimate assuming that S 3 is weakly irrelevant, i.e., [γ] is negative and has small absolute value, and [φ] ≈ 0.5. Our preliminary simulation for the herringbone dimer model (see Fig. 1(c) in Ref. 10) that is considered to have the cubic term shows a similar unconventional behavior of the aspect ratio as the staggered dimer model [26].
In this paper, we presented the finite size scaling method with controlling anisotropy of the system dynamically. In virtually isotropic systems, the corrections to scaling peculiar to the anisotropic systems is reduced, and we can compare the critical amplitudes among the classical and quantum systems. This method can give the optimal system size including L τ , which has been chosen sufficiently large value because there has been no index. We applied this method to the spatially anisotropic Heisenberg models, which was considered to be hard to judge its universality class because of the extremely large corrections to scaling. We concluded they belong to the same standard O(3) universality class based on the critical amplitudes and the critical exponents. We also revealed the optimized aspect ratio shows the nonmonotonic behavior from the cubic term S 3 of the effective action.
The simulation code has been developed based on the ALPS/looper library [21,27,28]. I acknowledge support by the Grand Challenge to Next-Generation Integrated Nanoscience, Development and Application of Advanced High-Performance Supercomputer Project from MEXT,