Criticality and scaling corrections for two-dimensional Heisenberg models in plaquette patterns with strong and weak couplings

TL;DR

This study uses quantum Monte Carlo simulations to precisely determine critical points and exponents in two-dimensional Heisenberg models with various plaquette arrangements, confirming their universal critical behavior and highlighting the importance of scaling corrections.

Contribution

The paper provides improved critical point estimates for three plaquette models and demonstrates their universal critical behavior within the 3D O(3) class, emphasizing the role of scaling corrections.

Findings

Critical points for three models are precisely determined.

All models exhibit 3D O(3) universality class behavior.

Effective scaling corrections are crucial for accurate analysis.

Abstract

We use the stochastic series expansion quantum Monte Carlo method to study the Heisenberg models on the square lattice with strong and weak couplings in the form of three different plaquette arrangements known as checkerboard models C$2\times2$, C$2\times4$ and C$4\times4$. The $a\times b$ here stands for the shape of plaquette consisting with spins connected by strong couplings. Through detailed analysis of finite-size scaling study, the critical point of C$2\times2$ model is improved as $g_{c}=0.548524(3)$ compared with previous studies with $g$ to be the ratio of weak and strong couplings in the models. For C$2\times4$ and C$4\times4$ we give $g_{c}=0.456978(2)$ and $0.314451(3)$. We also study the critical exponents $\nu$, $\eta$, and the universal property of Binder ratio to give further evidence that all quantum phase transitions in these three models are in the three-dimensional O(3) universality class. Furthermore, our fitting results show the importance of effective corrections in the scaling study of these models.

Tables

Table 1. Table 1: The finite-size scaling results of the critical point g c subscript 𝑔 𝑐 g_{c} using different quantities for the C 2 × 2 2 2 2\times 2 (top), C 2 × 4 2 4 2\times 4 (middle), and C 4 × 4 4 4 4\times 4 (bottom) models. Here we fit all the data points in Fig.  3 using scaling formula Eq. ( 16 ) with 1 / ν = 1.406 1 𝜈 1.406 1/\nu=1.406 .

	$g_c$	$\omega$	${\chi }^2/d.o.f$
$R_2$	$0.548532\left(6\right)$	1.14(2)	0.89
${\chi }_{u}L$	$0.548522\left(8\right)$	0.83(4)	0.88
${\rho }_{s}L$	$0.548521\left(5\right)$	0.68(2)	1.09
$R_2$	$0.456985\left(6\right)$	1.08(2)	0.63
${\chi }_{u}L$	$0.456972\left(8\right)$	0.86(3)	0.67
${\rho }_s^{x}L$	$0.456975\left(5\right)$	0.70(2)	0.92
${\rho }_s^{y}L$	$0.456983\left(6\right)$	0.60(3)	0.95
$R_2$	$0.31446\left(1\right)$	1.08(5)	1.11
${\chi }_{u}L$	$0.314441\left(9\right)$	0.87(5)	0.80
${\rho }_{s}L$	$0.314449\left(6\right)$	0.67(3)	0.92

Table 2. Table 2: Estimate results of the critical binder ratio R 2 ​ c subscript 𝑅 2 𝑐 R_{2c} for the C 2 × 2 2 2 2\times 2 , C 2 × 4 2 4 2\times 4 , and C 4 × 4 4 4 4\times 4 models. The fitted curves are shown in Fig.  6 .

	$R_{2c}$	$\omega$	${\chi }^2/d.o.f$
C$2\times 2$	$2.2549\left(5\right)$	1.163(9)	1.00
C$2\times 4$	$2.2549\left(7\right)$	1.17(1)	1.29
C$4\times 4$	$2.2542\left(9\right)$	1.02(2)	0.84

Equations

$$H=J_1\sum _{⟨i,j⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{⟨i,j{⟩}^{\prime }}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$R_2=\frac{⟨(m_s^z{)}^4⟩}{⟨(m_s^z{)}^2{⟩}^2},$$

$$m_s^z=\frac{1}{N}\sum _i^N{S}_i^z(-1{)}^{x_i+y_i},$$

$${\chi }_u=\chi (0,0)=\frac{\beta }{N}⟨{\left(\sum _{i=1}^N{S}_i^z\right)}^2⟩,$$

$$\delta f=\frac{1}{2}\rho (\nabla \theta {)}^2=\frac{1}{2}\rho (\Phi /L{)}^2$$

$${\rho }_s^a=\frac{3}{2\beta N}⟨(N_a^{+}-N_a^{-}{)}^2⟩,$$

$$Q(g,L)=L^{\kappa /\nu }f(\delta L^{1/\nu },{\lambda }_1{L}^{-{\omega }_1},{\lambda }_2{L}^{-{\omega }_2},\dots ),$$

$$Q(g,L)=f_Q^{(0)}(\delta L^{1/\nu })+L^{-{\omega }_1}{f}_Q^{(1)}(\delta L^{1/\nu }),$$

$$\begin{array}{r}{f}_Q^{(0)}(g^{\ast }{L}^{1/\nu })+L^{-{\omega }_1}{f}_Q^{(1)}(g^{\ast }{L}^{1/\nu })\\ =f_Q^{(0)}(g^{\ast }{L}^{{}^{\prime }1/\nu })+L^{{}^{\prime }-{\omega }_1}{f}_Q^{(1)}(g^{\ast }{L}^{{}^{\prime }1/\nu })\end{array}$$

$$g^{\ast }=\frac{f_Q^{(1)}(0)}{f_Q^{(0{)}^{\prime }}(0)}\frac{b^{-{\omega }_1}(b^{{\omega }_1}-1)}{b^{1/\nu }-1}{L}^{-{\omega }_1-1/\nu },$$

$$g_c(L)=g_c(\infty )+\frac{f_Q^1(0)}{f_Q^{(0{)}^{\prime }}(0)}\frac{b^{-{\omega }_1}(b^{{\omega }_1}-1)}{b^{1/\nu }-1}{L}^{-{\omega }_1-1/\nu }.$$

$$Q(L)=Q_{\infty }(g_c)+\frac{b^{-{\omega }_1}(1-b^{1/\nu +{\omega }_1})f_Q^{(1)(0)}}{b^{1/\nu }-1}{L}^{-{\omega }_1}.$$

$$\frac{1}{\nu (L)}=\frac{1}{ln(b)}\left(ln\frac{S(L^{{}^{\prime }})}{S(L)}\right),$$

$$S(L)=\frac{dQ(g,L)}{dg}\bigg{|}_{g=g_{c}(L)},$$

$$\frac{1}{\nu (L)}=\frac{1}{\nu }+a{L}^{-\omega }$$

$$g_c(L)$$

$$Q_c(L)$$

$$⟨m_s^2⟩\propto L^{-(1+\eta )}(1+a{L}^{-\omega }).$$

$$\eta (L)=\frac{ln[⟨m_s^2(L)⟩/⟨m_s^2(2L)⟩]}{ln(2)}-1.$$

$$\eta (L)=\eta +d{L}^{-\omega },$$

Full text

I Introduction

The $S=1/2$ Heisenberg antiferromagnetic model with different interactionsoverall ; Murg2009 has always been a very interesting topic in both theoretical and experimental fields because of its rich ground states and close relations to cuprate superconductorsHorsch1988 ; cup91 ; Keimer1992 ; Lee2006 , Bose-Einstein condensation of magnonsbec08 ; Rakhimov2012 , etc. One of the best studied two-dimensional (2D) Heisenberg models is the dimerized modelMatsumoto2001 ; andersbook ; Jiang2012 ; accurate with inter- and intradimer antiferromagnetic couplings on the square lattice, which bring in quantum fluctuations to destroy the N$\acute{e}$el ground state and make the model undergo a quantum phase transition (QPT)always from antiferromagnetic (AFM) to quantum paramagnetic (QPM)cubic ; accurate ; different ; wangl06 ; Syljuasen2006 ; gc ; yao10 . Based on field analysis mapping to a nonlinear $\sigma$ model, this QPT belongs to the three-dimensional (3D) O(3) universality class3dmap , which is also proved by several separate numerical results with high accuracydifferent ; ma18 ; suwa15 ; wangl06 .

Apart from those well studied dimerized models, a QPT from AFM to QPM can also be realized by introducing strong and weak couplings which favor the formation of singlets in quadrumerized or other patterns, which connect more spins as long as there are an even number of strong couplings in the unitKoga1999 ; wangl06 ; Mambrini2006 ; different ; troyer08 . These patterns are referred to as the checkerboard patterns, which were proposed to explain the experiments of real-space structures observed in Bi2Sr2CaCu2O8+δ and Ca2-xNaxCuO2Cl2experiment1 ; experiment2 ; Shen2005 . The quadrumerized Heisenberg model on the square lattice with $2\times 2$ spins connected by stronger couplings can also be very helpful in the study of the Shastry-Sutherland modelssmodel ; wessel02 , which explains the critical properties of SrCu2(BO3)2ssexperiment . Recently there has been a very popular discussion about the order of QPT in the Shastry-Sutherland modelssdqcp ; bowen . However, except for the quadrumerized Heisenberg model, no numerical study has ever been done on plaquette models in which larger numbers of spins have been connected by strong couplings. Even for the quadrumerized model, the previous best estimate for the critical point is $g_{c}=1.8230(2)$, with $g$ being the ratio of strong and weak couplings in the systemwenzel09 , whose accuracy is at least one order of magnitude larger than that of the dimerized model (e.g., $g_{c}=1.90951(1)$ in the columnar dimerized model (CDM)ma18 ) or the classical 3D Heisenberg modelclassical . We note that the coupling ratio $g$ in this work is weak couplings divided by strong couplings while in dimerized models it is reciprocal as mentioned above. Besides, a recent work concerning QPT from AFM to QPM shows that different local symmetries may bring in different critical corrections at QPTsma18 . It has answered a long standing issue that the quantum Monte Carlo (QMC) simulation results of critical exponents in a staggered dimerized model (SDM) are not the standard O(3) valuesdifferent ; cubic . In their work they compared the critical exponents and correction forms of the SDM with those of the CDM to show that they belong to the same universality class with different corrections. However, the previous work that had claimed to find different exponents in the SDM also compared it with the quadrumerized Heisenberg model, whose correction form has not been carefully studied yet.

In this paper we study a series of plaquette antiferromagnetic Heisenberg models on a square lattice with the Hamiltonian

[TABLE]

where ${\bf S}_{i}$ denotes an $S=1/2$ spin operator at lattice sites $i$, and $\langle i,j\rangle$ and $\langle i,j\rangle^{\prime}$ are the nearest-neighbor sites connected by corresponding coupling strengths, which are represented by strong couplings $J_{1}$ and weak couplings $J_{2}$. According to different checkerboard patterns of the arrangements of $J_{1}$ and $J_{2}$ shown in Fig. 1, we refer to these plaquette models as the C$2\times 2$, C$2\times 4$, and C$4\times 4$ models. We set $J_{1}=1$ and define the ratio of weak and strong couplings to be $g=J_{2}/J_{1}$. When $g=1$, the model becomes an isotropic Heisenberg plane which has an antiferromagnetic ground state with long-range order. When $g=0$, the ground state turns into a disordered phase with no magnetism. It is a product state of singlets that differs in different plaquette modelsXu2019 . In this case, for $0<g<1$ there is a critical point $g_{c}$ at zero temperature where a QPT from AFM to QPM would happen. This QPT is in the O(3) universality according to the nonlinear $\sigma$ mapping analysis class3dmap . We use the stochastic series expansion (SSE) QMC methodandersbook and a finite-size scaling (FSS) study to estimate the critical points and exponents in the thermodynamical limit for all three models.

There are mainly two purposes to studying these three plaquette models. The first one is that we want to obtain the critical point of the C$2\times 2$ (quadrumerized) model with higher accuracy by large scale QMC calculation and a detailed FSS study on several different variables. The final result for the C$2\times 2$ model is $g_{c}=0.548524(3)$, which is improved quite obviously compared with $g_{c}=0.54854(1)$ as the previous best estimate. Except for being helpful to the exploration of finite-T quantum criticality by minimizing the quantum regime in the C$2\times 2$ model itself, the increase of the statistical accuracy here would be useful to study the related models (i.e., Shastry-Sutherland model) toowessel02 . We also obtain $g_{c}=0.456978(2)$ and $g_{c}=0.314451(3)$ for the C$2\times 4$ and C$4\times 4$ models, respectively, for the first time. This could be very necessary in future study of a certain material with this kind of real-space structure in the experiment. These results of critical points offer a very good benchmark for further tests or development of numerical techniques of FFS methods as well. Second, from the FSS of criticality of the C$2\times 2$ model, we can determine its correction form and correction exponent to compare with those of the SDM as a complement to the comparison in Ref. [ma18, ]. The nonmonotonic scaling behavior is not found in the criticality with critical exponents in the O(3) universality class in any of these three models, which is similar to the CDM. Thus, our results give more examples that having local symmetry lacking cubic couplings which can bring in corrections not present in the standard O(3) class. The scaling correction forms and correction exponents in different models are the same as each other, which implies that the different local $Z_{4}$ (C$2\times 2$ and C$4\times 4$) or $Z_{2}$ (C$2\times 4$) symmetries do not result in any difference. These results could be helpful in understanding the QPTs with irrelevant field.

In addition to the main purposes above, the standard O(3) value $1/\nu=1.4061(7)$ is chosen for all three models from the 3D classical Heisenberg modelclassical and $\chi^{2}/d.o.f$ ($\chi^{2}$ per degree of freedom) close to 1 for all fits. Through these scalings, all critical points obtained from different quantities agree with each other in one system, which offers computational evidence to confirm that QPTs in all three models belong to the O(3) universality class. Besides $\nu$, we also compare the anomalous dimension $\eta$ and the dimensionless Binder ratio at critical points as further evidence to prove the predicted universality class.

The rest of the paper is organized as follows. In Sec. II we introduce the physical quantities calculated in this work and the finite-size scaling method that we utilize to analyze data from simulations. In Sec. III simulation and FFS results of criticalities for all three models are presented with detailed analysis. In the end we give a brief summary and discussion in Sec. IV.

II Observables and Finite-size scaling

We use the SSE QMC simulation method with an operator-loop updating algorithm to study all plaquette models in our work. This computing method is based on sampling of the diagonal elements of the Boltzmann operator $exp(-\beta H)$, with $\beta$ being the inverse temperature. In order to rule out the effect of temperature on the scaling function near the quantum critical point, $\beta$ is always chosen as $\beta\sim L^{z}$. The QPTs in the plaquette models studied here are believed to be in accordance with O(3) behavior so that $z=1$Troyer1997 in these models; therefore we consider $\beta=L$ in our study with simulated system size $L$ up to $160$. In all calculations we use $~{}10^{4}$ Monte Carlo samplings to obtain average values of the observables.

II.1 Observables

In order to study the criticality of the certain spin model, we chose to measure several important physical quantities in our work. The first one is the Binder ratio defined as

[TABLE]

where

[TABLE]

where $N=L\times L$ is the total number of spins on the square lattice and $(x_{i},y_{i})$ are the coordinates of the corresponding spin $S_{i}$. The Binder ratio is dimensionless and universal regardless of the detailed structures and couplings of the model. However it does depend on the boundary conditions and effective aspect ratios from previous studiesHasenbusch2001 ; Selke2007 ; Selke2009 ; Qian2016 . Here we use periodic boundary conditions on these three models and the effective aspect ratio of time-space is related to the critical spin-wave velocity.

Another quantity is the uniform susceptibility

[TABLE]

whose scaling form at $g_{c}$ is $\chi_{u}\sim L^{z-d}$, giving $\chi_{u}\sim L^{-1}$ and $\chi_{u}L$ to be dimensionless in our case.

The last physical observable calculated in our work is the spin stiffness $\rho_{s}$. The stiffness $\rho$ is covered in the calculation

[TABLE]

in the continuum field theory with $f$ being the density of free energy, $\Phi$ the boundary twist, and $\theta$ the order parameter field. In the Heisenberg model, $\rho_{s}$ is the spin stiffness determined by the twist $\Phi$ directly to the Hamiltonian, which in the SSE procedure can be obtained through the calculation

[TABLE]

where $N_{a}^{+}$ ($N_{a}^{-}$) represent the total number of $S^{+}_{i}S^{-}_{j}$ ($S^{-}_{i}S^{+}_{j}$) operators in the sampling along the $a$ ($x$ or $y$) direction of the square lattice. When the system is isotropic, the lattice $\rho_{s}^{x}$ is the same as $\rho_{s}^{y}$, while for the anisotropic system they are different. So in the C$2\times 2$ and C$4\times 4$ models we only calculate $\rho_{s}=(\rho_{s}^{x}+\rho_{s}^{y})/2$, and for the C$2\times 4$ model both $\rho_{s}^{x}$ and $\rho_{s}^{y}$ are recorded separately. However they all have the same scaling form at the critical point as $\rho_{s}\sim L^{2-d-z}$, with $\rho_{s}\sim L^{-1}$ in our models, which means that $\rho_{s}L$ is a size-independent dimensionless quantity.

II.2 Finite-size scaling

After all the mean observable values mentioned above are obtained from the simulations, we need to deal with all these data using the finite-size scaling study method to estimate the critical properties in the thermodynamical limitandersbook ; Jiang2012 ; accurate ; cubic ; different ; wangl06 ; ma18 ; troyer08 ; Shao2016 ; Liu2018 . From the renormalization group theory we know that the physical quantity $Q$ near its critical point obeys

[TABLE]

where $\kappa$ is the critical exponent of $Q$, $\nu$ is the correlation length exponent, and $\delta=g-g_{c}$. The set $\{\lambda_{i}\}$ refers to all irrelevant fields with their correction exponents $\{\omega_{i}\}$, which are arranged as $\omega_{i+1}>\omega_{i}$. Usually at most one irrelevant field is supposed to be considered in the FSS analysis, but there are still some special cases where more than one field is necessaryma18 . Here we start with one correction exponent to the first order of the dimensionless quantities ($\kappa=0$) so that Eq.(7) can be written as

[TABLE]

in which $L^{-\omega_{1}}$ is regarded as a deviation value of the theoretical scaling function $f_{Q}$ near the critical point. Ignoring irrelevant items, the dimensionless quantity $Q(g,L)$ does not depend on the size of the system at the critical point $g_{c}$ because $g=0$ then. Thus, $Q(g,L)$ values for different sizes cross at the critical point in this simplified situation, but here we need to take irrelevant items into consideration and $Q(g,L)$ values for different sizes would cross at $g_{c}(L)$, which is near to the real $g_{c}$ with a correction to the order $L^{-\omega_{1}}$. For two different simulated sizes $L$ and $L^{{}^{\prime}}$, using Eq.(8), we have

[TABLE]

at the cross point $g_{c}(L)$, with $g^{*}=g_{c}(L)-g_{c}$. Expanding $f_{Q}^{(0)}$ and $f_{Q}^{(1)}$ to the first order of $L^{-\omega_{1}}$ with $L^{{}^{\prime}}=bL$ we can get

[TABLE]

which is more easily understood as

[TABLE]

Insert Eq. (10) into Eq. (8) and again expand $f_{Q}^{(0)}$ and $f_{Q}^{(1)}$ to the first order of $L^{-\omega_{1}}$, then we have

[TABLE]

Besides, with the definition of

[TABLE]

where

[TABLE]

we can also obtain the scaling of critical exponent $\nu$ combining Eqs.(11) and (13) as

[TABLE]

with a free parameter $a$. For simplicity the scaling forms of the coordinates of crossing points $(g_{c}(L),Q_{c}(L))$ are written as

[TABLE]

with $b$ and $c$ to be fitted as free parameters. From Eqs.(16) and (15), we know that by using crossing points from the $g$ dependence of the dimensionless quantity for two sizes $(LandbL)$, the extrapolation value when $L\rightarrow\infty$ can give the quantum critical point $g_{c}$ and the critical exponent $\nu$ at the thermodynamical limit. In our work, we use $b=2$ in obtaining all crossing points for different models.

III Simulation results and data analysis

We performed the SSE QMC simulations on the C$2\times 2$, C$2\times 4$, and C$4\times 4$ models and obtained the average values of all the observables $R_{2}$, $\chi_{u}$, and $\rho_{s}$ (or $\rho_{s}^{x}$ and $\rho_{s}^{y}$ especially for C$2\times 4$). One example of the simulation results for C$2\times 4$ is illustrated in Fig. 2 to show the crossings of different sizes for the four dimensionless quantities $R_{2}$, $\chi_{u}L$ , $\rho_{s}^{x}L$, and $\rho_{s}^{y}L$. Similar figures can also be obtained from the SSE data for the C$2\times 2$ and C$4\times 4$ models. The obvious shift of crossings from different sizes implies that it is necessary to take the correction into account in the scaling analysis.

III.1 Critical points and corrections

After all crossing points are extracted from the raw data we use the finite-size scaling method to estimate the critical points of our models. By fitting all the data points in Fig. 3 using the function in Eq. (16) separately for each quantity we obtain all the $g_{c}$ results shown in Table 1. In all the fits we use the standard O(3) value $1/\nu=1.406$ with one correction exponent $\omega$ to the first order. For the C$2\times 2$ model, $g_{c}$ from each quantity is the same considering one error bar and agrees with the former result $g_{c}=0.54854(1)$. It is also true for the other two models with $\chi^{2}/d.o.f$ close to 1, implying the credibility of the fits. These results give further evidence to show that plaquette models with different checkerboard patterns all belong to the O(3) universality class.

In order to obtain a better estimation of the critical points, we continue to deal with the crossing points by joint fits as all size dependencies of $g_{c}(L)$ for different quantities should converge to the same value in one system. Therefore we fix $g_{c}(\infty)$ to be the same in each curve and fit all data together with other parameters being independent and $1/\nu=1.406$. The fitting results are shown in all curves in Fig.3 with $g_{c}=0.548524(3)$ in C$2\times 2$, $0.456978(2)$ in C$2\times 4$, and $0.314451(3)$ in C$4\times 4$. Our result for the C$2\times 2$ model is fully consistent with the value in Ref. [wenzel09, ] with higher precision. By comparing these critical point values we find $g_{c}$ gets smaller from the C$2\times 2$ model to the C$4\times 4$ model, indicating that our models more easily turn into the QPM state with less strong couplings in the unit cell. Therefore we deduce it is a universal rule in QPTs of any C$L_{x}\times L_{y}$ models.

From the separately fitting results in Table 1, we find that with only one correction term included the correction exponents $\omega$ are not the same for different quantities in the same model, while they are the same for same quantity in different models within at most two error bars with taking average of $\rho_{s}^{x}$ and $\rho_{s}^{y}$ in the C$2\times 4$ model. The difference shows that $\omega$ calculated here is more likely to be an “effective correction” including higher orders. However, fitting including $2\omega$ or higher order is very difficult and challenging with too many free parameters. Here we did not find any nonmonotonic behavior in the size dependence of all crossings in the plaquette models as shown in Fig. 3; therefore, one correction term can also give convincing criticality analysis, which is also confirmed by the $\chi^{2}/d.o.f$ of each fit. The joint fitting results still share the same rule as the separate ones. Thus, we can estimate the effective $\omega$ by taking the weighted average values of the joint fitting results of $\omega$ from three models. We have the effective correction exponent $\omega=1.058(7)$ for $R_{2}$, for $\chi_{u}L$ $\omega=0.834(7)$. For $\rho_{s}L$ we first get the average $\omega$ from the correlated results of $\rho_{s}^{x}$ and $\rho_{s}^{y}$ in the C$2\times 4$ model and we take the larger error of them as the error, which gives $\omega=0.65(1)$ in the end. Then taking the weighted average of all three values gives $\omega=0.66(1)$ for $\rho_{s}L$. Comparing with the standard correction exponent $\omega\approx 0.78$classical in the O(3) universality class, we can see that the system sizes included in our fits are still not large enough to rule out the affection of higher order corrections even with $L$ up to $160$. Therefore the value of the effective $\omega$ becomes very important in the FSS study to obtain the critical point and critical exponents.

III.2 Universal quantities at critical points

As discussed above, in order to study the critical point we use the fixed standard O(3) value $1/\nu=1.406$. The goodness of all fitting results also proves the theoretical prediction. In this section we consider further tests by studying some universal properties to give further evidence of the universality class of the phase transitions.

We start with two critical exponents $\nu$ and $\eta$, which are two sensitive universal quantities derived from the Binder ratio and magnetization at the critical point, respectively, and share scaling behavior similar to that of the physical quantities as we discussed in Sec. II. To begin with, the correlation length exponent is calculated using the scaling of $1/\nu(L)$, which is defined as Eq. (13), in Eq. (15). The simulation and scaling results are shown in Fig. 4 for all three models. Fitting values of $1/\nu$ are the same for all models considering error. The weighted average of all three $1/\nu$ value is $1.404(4)$. Compared with the best estimate of $1/\nu=1.4061(7)$ (reciprocal value of $\nu=0.7112(5)$ in Ref. [classical, ]) in O(3) it is proved again that the QPTs here are in the same universality class as the CDM, the SDM, and the 3D classical Heisenberg model. However, the accuracy of the estimation using scaling of $1/\nu(L)$ in our work is much less compared to the previous results. Usually $\nu$ can be obtained from the data collapse together with the critical point and corrections. Here we use the combination of two sizes together at once in order to lease the influence of the corrections. It does help as the fitting results of $\omega$ are very large in our study, which means that $1/\nu$ converges very fast with the increase of $L$. But the values of $1/\nu(L)$ obtained from simulations have much larger errors compared with other quantities studied before, which brings in larger error to the final extrapolation value. Much more computational effort is needed in order to obtain better estimation of $\nu$. We just stop here in this paper as it is not a key point of our work, but we want to point it out for other studies using this procedure.

Another critical exponent considered here is the anomalous dimension $\eta$. Once the critical point $g_{c}$ is obtained, we can study the scaling of order parameter at $g_{c}$ using

[TABLE]

Similar to $1/\nu(L)$, we can also define $\eta(L)$ from the scaling of pairs of size $L$ and $2L$ as

[TABLE]

In this way, the size dependence of $\eta(L)$ is

[TABLE]

with correction to the first order. This time our fits use the best-known estimation of $\eta=0.0375(5)$classical and leave the other parameters in Eq. (20) free. The fitting results shown in Fig. 5 again imply that it is correct to set $\eta=0.0375(5)$ here as all QPTs are in the O(3) universality class. Furthermore, all corrections are the same in the three models considering error and the weighted averaged $\omega=0.78(1)$ is the same as the first correction exponent $\omega_{1}=0.782(13)$ in the O(3) modelclassical . This shows that $\eta$ could be a good quantity in testing the correction exponent in the FSS study once a good estimation of $g_{c}$ is obtained. And these scaling results also give us more confidence in the accuracy of $g_{c}$ here for all three plaquette models.

At last we test the universal quantity Binder ratio $R_{2}$ in the thermodynamic limit. With crossing points extracted from two different sizes $(Land2L)$ of the $g$ dependence for $R_{2}$ near the critical point, we can obtain the scaling of $R_{2c}(L)$ in Fig. 6. The results of fits using Eq. (17) in Table 2 indicate that $R_{2c}$ converges to the same value within one error bar in three models with different checkerboard patterns. Except for a further proof of the same universality class, the same $R_{2c}$ in all cases implies that different plaquette models might have the same aspect ratios as well. Taking the weighted average of all $R_{2c}$ values gives $R_{2c}=2.2547(4)$ for a series of plaquette models. Thus, we predict that in all C$L_{x}\times L_{y}$ models according to our model definition, the Binder ratios $R_{2}$ would all converge to $2.2547(4)$ as long as there is an even number of spins in one unit cell with simulated $\beta=L$.

IV Summary and discussions

In this paper we carried out the FSS study on data with high-precision using the SSE QMC method. The criticality of three $S=1/2$ Heisenberg models on the square lattice with strong and weak couplings in plaquette patterns C$2\times 2$, C$2\times 4$, and C$4\times 4$ was studied using the Binder ratio, uniform susceptibility and spin stiffness. By the joint fits combining the scalings of crossing points from all three quantities we have obtained the most accurate estimates of the critical points $g_{c}$ for three plaquette models up to now. Our scaling analysis implies the importance of corrections in FSS, and with only one correction term the value of $\omega$ is more likely to be an effective one. The effective $\omega$ does not change in different models as long as it describes the scaling behavior of the same physical observable. The calculation of $1/\nu$, $\eta$ and $R_{2}$ at the critical point shows that QPTs in all three models are in the O(3) universality class as predicted. The scaling of $\eta$ using the order parameter at the critical point also gives $\omega\approx 0.78$, the same as $\omega_{1}$ determined in the 3D classical Heisenberg model, which further supports the estimate of the critical points.

The fitting results of $\omega$ using different quantities in these three models help us to understand the influence of corrections in the scalings. With system sizes up to $L=160$ the correction exponent $\omega$ is still an effective one that differs with different variables. However, fitting including higher orders of correction terms would be quite challenging and difficult. Here we find that for models with detailed difference structures in our case the effective $\omega$ does not change for the same quantity. This might be helpful in further FSS studies on other similar models. We also obtain the value of the universal quantity $R_{2}$ at the critical point, and we suggest that any C$L_{x}\times L_{y}$ models might have the same critical Binder ratio value with $L_{x}L_{y}$ being even and the same $\beta/L$. This value could be a very important referee in later study of quantum phase transitions, since the same Binder ratio value could be very convincing supporting evidence to show whether a new phase transition belongs to the O(3) universality class.

Acknowledgements.

We thank A. W. Sandvik for useful discussions and a careful reading of the manuscript. The work of X.R and D.X.Y was supported by Grants No. NKRDPC-2018YFA0306001, No. NKRDPC-2017YFA0206203, No. NSFC-11574404, and No. NSFG-2015A030313176; the National Supercomputer Center in Guangzhou; and the Leading Talent Program of Guangdong Special Projects.