Lee-Yang-zero ratio method in three-dimensional Ising model

TL;DR

This paper applies the Lee-Yang-zero ratio method to the 3D Ising model using Monte Carlo simulations, demonstrating its effectiveness in locating the critical point and providing universal critical values.

Contribution

It introduces and validates the Lee-Yang-zero ratio method as an effective tool for identifying critical points, with advantages over traditional methods.

Findings

LYZR method is as effective as Binder-cumulant method

Universal Lee-Yang-zero ratios at criticality determined

Single Lee-Yang zero approach also useful

Abstract

By performing Monte Carlo simulations of the three-dimensional Ising model, we apply the recently proposed Lee-Yang-zero ratio (LYZR) method to determine the location of the critical point in this model. We demonstrate that the LYZR method is as powerful as the conventional Binder-cumulant method in studying the critical point, while the LYZR method has the advantage of suppressing the violation of the finite-size scaling and non-linearity near the critical point. We also achieve a precise determination of the values of the LYZRs at the critical point, which are universal numbers. In addition, we propose an alternative method that uses only a single Lee-Yang zero and show that it is also useful for the search for the critical point.

Tables

Table 1. Table 1: Results of the fit parameters T c T_{\rm c} , y t y_{t} , and r n ​ 1 r_{n1} or b 4 b_{4} obtained by the intersection analysis of the LYZR and Binder-cumulant methods for various smallest system size L min L_{\textrm{min}} .

	$L_{\min }$	$T_{\text{c}}$	$y_t$	$r_{n1}$ or $b_4$
$R_{21}$	192	4.5115240(148)	1.5338(410)	2.4179(82)
	128	4.5115185(59)	1.5874(115)	2.4195(26)
	96	4.5115161(43)	1.5947(72)	2.4206(15)
	64	4.5115166(36)	1.5891(56)	2.4208(11)
	48	4.5115159(30)	1.5890(32)	2.4213(08)
$R_{31}$	192	4.5115342(120)	1.4899(375)	3.6781(138)
	128	4.5115304(51)	1.5662(121)	3.6817(45)
	96	4.5115248(38)	1.5880(72)	3.6870(25)
	64	4.5115202(33)	1.5885(53)	3.6912(20)
	48	4.5115185(30)	1.5893(31)	3.6924(16)
$R_{41}$	192	4.5115228(21)	1.4834(475)	4.8830(332)
	128	4.5115375(104)	1.5311(185)	4.8603(130)
	96	4.5115305(78)	1.5680(113)	4.8700(82)
	64	4.5115202(63)	1.5844(84)	4.8832(56)
	48	4.5115217(56)	1.5841(51)	4.8812(44)
$B_4$	192	4.5115280(113)	1.5313(267)	1.60230(638)
	128	4.5115211(52)	1.5691(86)	1.59892(224)
	96	4.5115146(37)	1.5829(52)	1.59604(130)
	64	4.5115099(28)	1.5888(33)	1.59397(72)
	48	4.5115089(26)	1.5880(21)	1.59368(56)

Table 2. Table 2: Values of LYZRs at the CP r n ​ 1 r_{n1} . The upper row shows the results obtained by the intersection analysis with L min = 128 L_{\rm min}=128 . Shown in the middle row are the values obtained by the extrapolation to L → ∞ L\to\infty via Eq. ( 34 ) at T = T c T=T_{\rm c} in Eq. ( 2 ). Systematic uncertainties of the extrapolation are denoted by the second error. The values obtained for the CP in the three-dimensional three-state Potts model  [ 17 ] are shown in the lower panel.

	$r_{21}$	$r_{31}$	$r_{41}$
intersection analysis ($L_{\min }=128$)	$2.4195\left(26\right)$	$3.6817\left(45\right)$	$4.8603\left(130\right)$
extrapolation at $T=T_c$ in Eq. (2)	$2.4158\left(37\right)\left(5\right)$	$3.6906\left(68\right)\left(89\right)$	$4.8896\left(156\right)\left(23\right)$
3-state Potts model [17]	$2.408\left(12\right)$	$3.669\left(24\right)$	$4.861\left(36\right)$

Table 3. Table 3: Normalized curvatures C f C_{f} for the LYZR R n ​ m ​ ( T , L ) R_{nm}(T,L) and the Binder cumulant B 4 ​ ( T , L ) B_{4}(T,L) .

$f$	$R_{21}$	$R_{31}$	$R_{41}$	$B_4$
$C_f$	0.0093(17)	0.0053(17)	-0.0010(34)	0.0364(37)

Table 4. Table 4: Normalized curvatures C f C_{f} for the single-LYZ and cumulant methods, where f = h LY ( n ) f=h_{\rm LY}^{(n)} for n = 1 , 2 n=1,2 and ⟨ M n ⟩ c \langle M^{n}\rangle_{\rm c} for n = 2 , 4 n=2,4 .

$f$	$h_{\mathrm{LY}}^{\left(1\right)}$	$h_{\mathrm{LY}}^{\left(2\right)}$	${⟨M^2⟩}_c$	${⟨M^4⟩}_c$
$C_f$	0.07149(59)	0.04230(101)	-0.07209(44)	-0.17974(104)

Equations

$$H(h,L;\{\sigma \})=-J\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j-h\sum _i{\sigma }_i,$$

$$T_{\text{c}}=4.51152321(4).$$

$$t=\frac{T-T_{\text{c}}}{T_{\text{c}}}.$$

$$Z(t,h,L)=\sum _{\{{\sigma }_i\}}{e}^{-H(h,L;\{\sigma \})/T}=\sum _{\{{\sigma }_i\}}\exp (\frac{1}{T}\sum _{⟨i,j⟩}{\sigma }_i{\sigma }_j+\frac{h}{T}\sum _i{\sigma }_i).$$

$$M=\sum _i{\sigma }_i,$$

$$h=h_{\text{LY}}^{(n)}(t,L),$$

$$F_{\text{sing}}(t,h,L)=\tilde{F}(L^{y_t}t,L^{y_h}h).$$

$$y_t=1.58737472(29),y_h=2.481851194(24).$$

$$Z_{\text{sing}}(t,h,L)=\tilde{Z}(L^{y_t}t,L^{y_h}h),$$

$$R_{nm}(t,L)=\frac{\text{Im}{h}_{\text{LY}}^{(n)}(t,L)}{\text{Im}{h}_{\text{LY}}^{(m)}(t,L)}.$$

$$R_{nm}(t,L)=\frac{h_{\text{LY}}^{(n)}(t,L)}{h_{\text{LY}}^{(m)}(t,L)}.$$

$$h_{\text{LY}}^{(n)}(t,L)\underset{L\to \infty }{\to }a(t)\frac{2n-1}{L^3}(n\ge 1),$$

$$R_{nm}(t,L)\underset{L\to \infty }{\to }\frac{2n-1}{2m-1}(t<0).$$

$$R_{nm}(t,L)\underset{L\to \infty }{\to }1(t>0),$$

$$Z_{\mathrm{sing}}\left(t,h_{\text{LY}}^{(n)}(t,L),L\right)=\tilde{Z}\left(L^{y_t}t,L^{y_h}{h}_{\text{LY}}^{(n)}(t,L)\right)=0.$$

$$L^{y_h}{h}_{\mathrm{LY}}^{(n)}(t,L)={\tilde{h}}_{\mathrm{LY}}^{(n)}(L^{y_t}t).$$

$${\tilde{h}}_{\text{LY}}^{(n)}(\tilde{t})=i(X_n+Y_n\tilde{t}+Z_n\tilde{t}{}^2+\mathcal{O}(\tilde{t}{}^3)),$$

$$R_{nm}(t,L)=r_{nm}+c_{nm}{L}^{y_t}t+d_{nm}(L^{y_t}t{)}^2+\mathcal{O}(t^3),$$

$$⟨M^n(t,h,L){⟩}_{\mathrm{c}}=-T^{n-1}\frac{{\partial }^{n}F(t,h,L)}{\partial h^n}{|}_T.$$

$$⟨M^n(t,h,L){⟩}_{\mathrm{c}}=-L^{n{y}_h}{\tilde{F}}^{(n)}(L^{y_t}t,L^{y_h}h),$$

$$⟨M^n(t,0,L){⟩}_{\text{c}}=-L^{n{y}_h}{\tilde{F}}^{(n)}(L^{y_t}t,0).$$

$$B_4(t,L)=\frac{⟨M^4(t,0,L){⟩}_{\mathrm{c}}}{⟨M^2(t,0,L){⟩}_{\mathrm{c}}^2}+3=\frac{{\tilde{F}}^{(4)}(L^{y_t}t,0)}{{\left({\tilde{F}}^{(2)}(L^{y_t}t,0)\right)}^2}+3,$$

$$B_4(t,L)\underset{L\to \infty }{\to }\{\begin{array}{ll}3& (t>0)\\ 1& (t<0)\end{array}(\text{mbox}finiten,m).$$

$${\tilde{F}}^{(n)}(\tilde{t},0)={\tilde{F}}_0^{(n)}+{\tilde{F}}_1^{(n)}\tilde{t}+\frac{1}{2}{\tilde{F}}_2^{(n)}{\tilde{t}}^2+\cdots ,$$

$$B_4(t,L)=b_4+c_4{L}^{y_t}t+d_4(L^{y_t}t{)}^2+\mathcal{O}(t^3),$$

$$b_4=1.60356(15).$$

$$L^{y_h}{h}_{\text{LY}}^{(n)}(t,L)={\tilde{h}}_{\text{LY}}^{(n)}(L^{y_t}t)=X_n+Y_n{L}^{y_t}t+Z_n(L^{y_t}t{)}^2+\mathcal{O}(t^3),$$

$$\frac{⟨M^n(t,0,L){⟩}_{\mathrm{c}}}{L^{n{y}_h}}={\tilde{F}}^{(n)}(L^{y_t}t,0)={\tilde{F}}_n^{(0)}+{\tilde{F}}_n^{(1)}{L}^{y_t}t+{\tilde{F}}_n^{(2)}(L^{y_t}t{)}^2+\mathcal{O}(t^3).$$

$$N(t,h,L)=\frac{Z(t,h,L)}{Z(t,0,L)},$$

$$N(t,h,L)=\frac{⟨e^{-H(h,L)/T+H(0,L)/T_0}{⟩}_0}{⟨e^{-H(0,L)/T+H(0,L)/T_0}{⟩}_0},$$

Full text

YITP-25-128, J-PARC-TH-0324, UTHEP-810

1 Introduction

Critical phenomena are fundamental concepts across various fields of physics, including statistical mechanics, condensed-matter and high-energy physics, and cosmology. Although the physical scales involved in these fields differ by more than ten orders of magnitude, when a critical point (CP) exists, distinct physical systems exhibit identical critical behavior governed by the scaling laws and universality classes, in which symmetry and dimensionality of the system play crucial roles [1]. On the other hand, the existence and precise location of a critical point are not constrained by symmetries. Instead, the microscopic properties of the system, such as the system’s dynamics and interactions, play essential roles.

Various theoretical approaches have been developed to investigate critical phenomena, including the mean-field approximation [2, 3] and the renormalization-group methods [4, 5, 6]. However, they give exact solutions only in limited systems. Therefore, in most cases, numerical methods, such as the Monte Carlo [7, 8, 9], tensor network [10, 11], and conformal-bootstrap [12, 13] methods, serve as primary tools for locating the CP and determining its properties.

In numerical analyses of critical phenomena, it is crucial to properly account for finite-size effects, as simulations are always performed on finite volumes while the correlation length diverges at the CP. When the system size is sufficiently large, it is known that thermodynamic quantities in the vicinity of the CP obey an extended scaling law known as the finite-size scaling (FSS), where the system size acts as a scaling variable [14, 15, 9, 16]. The analysis based on the FSS is crucial for revealing properties of the CP, such as its location and the critical exponents. A useful method among them is the so-called Binder-cumulant method [15], in which the location of the CP is determined from the intersection point of the ratios of cumulants obtained at various system sizes.

In our recent study [17], we proposed a new systematic way to locate the CP based on the FSS. This method makes use of the Lee-Yang zeros (LYZs) [18, 19], i.e. zeros of the partition function on the complex-parameter space. Characteristic behavior of LYZs near the CP has been investigated in the literature [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. In particular, the edge of the distribution of the LYZ, called the Lee-Yang edge singularity [34], is known to follow a specific scaling law. This scaling property has been recently applied for the numerical search for the CP in quantum chromodynamics (QCD) at non-zero chemical potential in lattice QCD Monte Carlo simulations [35, 36, 37, 38, 39, 40]. However, the finite-size effects are not taken into account in these analyses.

In order to properly control finite-size effects exploiting the FSS of the LYZs, in Ref. \citenWada:2024qsk we proposed the use of the ratios of LYZs. We showed that the ratios of the imaginary parts of LYZs possess similar properties as the Binder cumulants, i.e. the ratios obtained on various system sizes intersect at the CP. The LYZ ratios (LYZRs) thus serve as an alternative and independent approach to locating the CP numerically through intersection analysis. In Ref. \citenWada:2024qsk, based on the FSS we also showed that this property of LYZRs applies to any CP in general systems, and the values of LYZRs at the CP are universal numbers that are specific for individual universality classes. Compared with the Binder-cumulant method, this method, which we call the LYZR method in what follows, has characteristics that the convergence of the finite-size effects is faster in general systems. We have also verified the validity of this method through a numerical analysis of the CP in the three-dimensional three-state Potts model with an external field, which belongs to the same universality class as the three-dimensional Ising model [41, 42]. It has been shown that the LYZR method can successfully determine the location and the critical exponents of the CP in this model. Preliminary results of the application of the LYZR method to QCD in the heavy-quark region are also presented in Ref. \citenWada:2024lat.

In the present study, we apply the LYZR method to Monte Carlo simulations of the CP in the three-dimensional Ising model. As the critical temperature $T_{\rm c}$ has been measured with great precision [7], this model serves as a good testing ground for the LYZR method and its comparison with other methods such as the Binder-cumulant method. We demonstrate that the LYZR method can determine the value of $T_{\rm c}$ and other parameters with nearly the same precision as the Binder-cumulant method. We also show that the method suppresses the correction to FSS and nonlinearities in the temperature dependence. Precise analyses of the values of LYZRs at the CP, which are the universal numbers, are also performed utilizing the simplicity of the model.

In addition, we propose a series of alternative methods to realize the intersection analysis of the CP with the use of a single LYZ. Whereas these methods assume that the universality class and the values of critical exponents are known, they are applicable to systems where the search for the second LYZ is difficult, such as the lattice-QCD simulations [35, 36, 37, 38, 39, 40].

This paper is organized as follows. In Sec. 2, we introduce the LYZR method and compare it with the Binder-cumulant method. The analysis with a single LYZ is also introduced. After describing the setup of our numerical simulations in Sec. 3, we present our numerical results of the LYZR and Binder-cumulant methods in Sec. 4. The numerical results for the single-LYZ method are reported in Sec. 5. In Sec. 6, we give a brief summary.

2 Theoretical framework

2.1 Ising model

Throughout this study, we investigate the ferromagnetic three-dimensional Ising (3d-Ising) model [44]

[TABLE]

with the interaction constant $J>0$ and external magnetic field $h$ on the cubic lattice of size $L^{3}$ under the periodic boundary conditions in all directions. Each lattice site $i$ has spin $\sigma_{i}$, which takes on the values $\sigma_{i}=\pm 1$, and $\langle i,j\rangle$ represents pairs of nearest-neighbor sites. In the following, we set $J=1$ to make all physical quantities dimensionless.

The 3d-Ising model (1) at nonzero temperature $T$ in the $L\to\infty$ limit exhibits a CP at $(T,h)=(T_{\textrm{c}},0)$, which is the endpoint of the ferromagnetic first-order phase-transition line at $h=0$ and $T<T_{\textrm{c}}$. The most precise value of $T_{\textrm{c}}$ known so far, to the best of the authors’ knowledge, is obtained by a Monte Carlo simulation as [7]

[TABLE]

In the following, we use the reduced temperature

[TABLE]

Conversion between $T$ and $t$ in the arguments of the following functions is understood.

2.2 Lee-Yang zeros

The partition function of the 3d-Ising model is given by

[TABLE]

While it is evident from Eq. (4) that $Z(t,h,L)$ is real and satisfies $Z(t,h,L)>0$ for real $T$ and $h$, it may vanish for complex parameters $T\in\mathbb{C}$ and/or $h\in\mathbb{C}$. Among them, the zeros on the complex-$h$ plane for real $T$ are called the LYZs [19, 18]. According to the Lee–Yang theorem, the LYZs in the 3d-Ising model are purely imaginary; that is, they always lie on the imaginary-$h$ axis. [18].

Since the total magnetization

[TABLE]

takes only even or odd integer values in the range $-L^{3}\leq M\leq L^{3}$, Eq. (4) is represented as a polynomial of the fugacity squared $e^{2h/T}$. The partition function can thus be factorized as $Z(t,h,L)\sim e^{-L^{3}h/T}\prod_{n}(e^{2h/T}-c_{n}(t))$ with $L^{3}$ factors, where $|c_{n}(t)|=1$ by the Lee-Ynag theorem [19]. This means that the number of LYZs for $-\pi<2{\rm Im}\,h/T\leq\pi$ is $L^{3}$ with possible degeneracies; hence they are discretely distributed. In the following, we denote the $n$th closest LYZs from the real-$h$ axis for ${\rm Im}\,h>0$ as

[TABLE]

which are pure imaginary functions of $t$ and $L$. Because Eq. (4) is real for $t\in\mathbb{R}$ and $h\in\mathbb{R}$, Eq. (4) has LYZs also at $h=[h_{\textrm{LY}}^{(n)}(t,L)]^{*}=-h_{\textrm{LY}}^{(n)}(t,L)$ for ${\rm Im}\,h<0$ from the Schwarz reflection principle.

Let us inspect the analyticity of Eq. (6), as this property plays a key role in later discussions. Since $Z(t,h,L)$ is an analytic function of $t$ and $h$ for finite $L$, it is suggested that $h_{\textrm{LY}}^{(n)}(t,L)$ should also be analytic in $t$; otherwise $Z(t,h,L)$ will not be analytic in $t$. This argument should hold as long as the LYZs do not have degeneracies. However, the situation is subtle when LYZs are degenerated at some $t$. For example, the partition function can contain a factor of the form $Z(t,h,L)=\big((h-h^{\prime})^{2}-k(t-t^{\prime})\big)\times\cdots$ when LYZs are degenerated at $t=t^{\prime}$. In this case, the LYZs $h=h^{\prime}\pm k^{1/2}(t-t^{\prime})^{1/2}$ are non-analytic at $t=t^{\prime}$ even though $Z(t,h,L)$ is analytic. (To satisfy the Lee-Yang theorem, $h^{\prime}$ should be pure imaginary and $k<0$. ) However, we have numerically confirmed that the LYZs near the real axis are always non-degenerate near the CP. We therefore assume that $h_{\textrm{LY}}^{(n)}(t,L)$ are analytic in $t$ in the following discussion.

2.3 Finite-size scaling

In the vicinity of the CP and for sufficiently large $L$, the free energy $F(t,h,L)=-T\ln Z(t,h,L)$ can be decomposed into the singular and regular parts as $F(t,h,L)=F_{\rm sing}(t,h,L)+F_{\rm reg}(t,h,L)$, where the singular part obeys the FSS relation [14, 9, 16]

[TABLE]

Here, $\tilde{F}$ is called the scaling function, and the critical exponents $y_{t}$ and $y_{h}$ are specific to the universality class; their values in the 3d-Ising model are known with high precision as [13]

[TABLE]

In response to the decomposition of the free energy, the partition function near the CP can be written as $Z(t,h,L)=Z_{\rm sing}(t,h,L)\,Z_{\rm reg}(t,h,L)$, where $Z_{\rm sing}(t,h,L)$ satisfies

[TABLE]

with a scaling function $\tilde{Z}$. In Appendix A, we generalize the argument to include the effects of irrelevant operators. The numerical influence of the dominant irrelevant operator on our results is studied in Sec. 4.2.

2.4 LYZ ratios

In Ref. \citenWada:2024qsk, we proposed a method to determine $T_{\textrm{c}}$ and the critical exponents using the FSS of LYZs in general systems. In this method, which we call the LYZ ratio (LYZR) method, we focus on the ratios of the imaginary parts of the LYZs

[TABLE]

In the Ising model, due to the Lee-Yang theorem, Eq. (10) can also be written as

[TABLE]

Let us examine the behavior of the LYZRs in the large volume limit $L\to\infty$. First, for $t<0$ in response to the first-order phase transition at $h=0$ it is known that the LYZs are distributed at equal distances on the imaginary-$h$ axis as

[TABLE]

where $a(t)$ is a purely imaginary function of $t$ [45]. The distribution around the real axis becomes dense for $L\to\infty$ in accordance with the discontinuity of $Z(t,h,L)$ at the first-order phase transition in this limit. From Eq. (12), one finds that the LYZRs behave as

[TABLE]

Second, for $t>0$, because $Z(t,h,L)$ is analytic at $h=0$, the LYZs are not accumulated around the real axis for $L\to\infty$. Instead, the LYZs in this limit are known to distribute continuously for $|{\rm Im}h|>h_{\rm LYES}(t)$, where the edges of the distribution $h=\pm ih_{\rm LYES}(t)$ is called the Lee-Yang edge singularity [46, 47]. As a result, the LYZRs behave as

[TABLE]

for finite $n,m$. Equations (13) and (14) show that the values of $R_{nm}(t,L)$ suddenly change at the CP at $t=0$ as depicted in Fig. 1 by the red-solid line. This property is similar to that of the Binder cumulant [15] that will be discussed in the next subsection.

Next, let us inspect the behavior of $R_{nm}(t,L)$ near the CP at finite $L$. Assuming that the $n$th LYZ is a zero of $Z_{\rm sing}(t,h,L)$, it is shown from Eq. (9) that $h_{\textrm{LY}}^{(n)}(t,L)$ satisfies

[TABLE]

Then, by defining zeros of the scaling function as $\tilde{h}=\tilde{h}_{\rm LY}^{(n)}(\tilde{t}\,)$, i.e. $\tilde{Z}(\tilde{t},\tilde{h}_{\rm LY}^{(n)}(\tilde{t}\,))=0$, the $n$th LYZ for different $L$ obey [20]

[TABLE]

As discussed in Sec. 2.2, the LYZs $h_{\textrm{LY}}^{(n)}(t,L)$, and hence $\tilde{h}_{\rm LY}^{(n)}(\tilde{t}\,)$, are analytic functions of $t$ and $\tilde{t}$, respectively. Therefore, $\tilde{h}_{\rm LY}^{(n)}(\tilde{t}\,)$ are Taylor expanded at $\tilde{t}=0$ as

[TABLE]

where $X_{n},Y_{n}$, and $Z_{n}$ are real coefficients. Combining Eqs. (17) and (16), one finds [17]

[TABLE]

with $r_{nm}=X_{n}/X_{m}$ and $c_{nm}=r_{nm}(Y_{n}/X_{n}-Y_{m}/X_{m})$. The coefficients $d_{nm}$ are also represented by the coefficients in Eq. (17).

Equation (18) shows notable features of $R_{nm}(t,L)$ around $t=0$. At the CP, the LYZRs satisfy $r_{nm}=R_{nm}(0,L)$, which is independent of $L$, while the slope scales as $L^{y_{t}}$. This means that $R_{nm}(t,L)$ obtained at various $L$’s intersect at a single point at $t=0$ as shown in Fig. 1. Therefore, the CP can be located from the intersection point of $R_{nm}(0,L)$ for various $L$.

The extension of this method to the CPs in general systems is discussed in Ref. \citenWada:2024qsk, where it has been shown that the values of the LYZRs at the CP, $r_{nm}$, are universal constants specific to individual universality classes.

2.5 Binder cumulant

Conventional quantities used to identify the CP are the cumulants of the magnetization $M$. From the free energy $F=-T\ln Z$, the cumulants are given by [48]

[TABLE]

By taking only the singular part $F_{\rm sing}(t,h,L)$ for the free energy assuming sufficiently small $t$ and large $L$, and combining Eq. (19) with Eq. (7) one obtains

[TABLE]

with $\tilde{F}^{(n)}=T^{n-1}\partial^{n}\tilde{F}/\partial h^{n}$. For $h=0$, Eq. (20) is written as

[TABLE]

To make systematic use of Eq. (21), it is convenient to consider ratios of the cumulants. For example, the ratio of the fourth-order cumulant and the square of the second-order cumulant is called the fourth-order Binder cumulant [15]. In the present study we introduce

[TABLE]

and refer to it as the Binder cumulant. These ratios share similar scaling properties to the LYZRs. First, Eq. (22) in the $L\to\infty$ limit behaves as

[TABLE]

Equation (23) for $t<0$ is obtained from the fact that the system is in a coexisting phase and the distribution of $M$ is two-peaked, while that for $t>0$ is a consequence of the Gauss distribution of $M$.

Second, near the CP at $t=0$, substituting the Taylor expansion of $\tilde{F}^{(n)}(\tilde{t},0)$

[TABLE]

into Eq. (22), one obtains

[TABLE]

with $b_{4}=\tilde{F}^{(4)}_{0}/(\tilde{F}^{(2)}_{0})^{2}+3$, $c_{4}=(F_{1}^{(4)}F_{0}^{(2)}-2F_{0}^{(4)}F_{1}^{(2)})/(F_{0}^{(2)})^{3}$, and so forth. Equations (25) shows that $B_{4}(t,L)$ shares properties similar to those of the LYZRs. It takes a unique value $B_{4}(0,L)=b_{4}$ for any $L$ at $t=0$, while the $t$ derivative scales as $L^{y_{t}}$. Therefore, the intersection point of $B_{4}(t,L)$ for various $L$ determines the location of the CP. The value at the intersection point, $b_{4}$, is a universal constant characteristic of the universality class [15]. In the 3d-Ising model, its value is known as [7]

[TABLE]

One can also readily verify that the ratios of even higher-order cumulants, such as $\langle M^{6}(t,0,L)\rangle_{\rm c}/\langle M^{2}(t,0,L)\rangle^{3}_{\rm c}$, exhibit the same property.

One can further construct quantities having the intersection property at the CP besides the LYZRs and Binder-cumulants [49, 50, 51]. An example is the one constructed from the correlation function $C(\ell)=\langle\sigma_{i}\sigma_{j}\rangle|_{|i-j|=\ell}$ [51]. It is known that its ratio at different distances, $C(2\ell)/C(\ell)$, intersects at $t=0$ and thus can be used for the determination of the CP. However, we do not discuss them further in this paper.

2.6 Intersection analysis with single LYZ and cumulant

In this subsection, we propose an alternative intersection analysis that uses a single LYZ assuming the known value of the critical exponent $y_{h}$.

For this purpose, we note that Eqs. (16) and (17) lead to

[TABLE]

around $t=0$. From Eq. (27) one easily finds that $L^{y_{h}}h_{\textrm{LY}}^{(n)}(t,L)$ for various $L$ intersect at $t=0$. Therefore, the CP can be identified through intersection analysis of $L^{y_{h}}h_{\textrm{LY}}^{(n)}(t,L)$ for different $L$. However, it should be noted that the crossing value $X_{n}$ is not universal in this case.

This idea is equally applicable to the cumulants of magnetization. A transform of Eq. (21),

[TABLE]

immediately tells us that the left-hand side, $\langle M^{n}(t,0,L)\rangle_{\rm c}/L^{ny_{h}}$, shares the same property as Eq. (27).

We refer to these methods as the single-LYZ and single-cumulant methods, respectively. These methods are tested in Sec. 5.

3 Numerical setup

We perform Monte Carlo simulations of the 3d-Ising model (1) using the Wolff cluster algorithm [52, 53]. We fix $h=0$ and generate spin configurations at three temperatures around $T_{\textrm{c}}$, $T=4.51125$, $4.5115$, and $4.51175$, for the lattice sizes $L=16,24,32,48,64,96,128,192$, and $256$. For thermalization, at least $4\times 10^{5}$ Monte Carlo updates are performed starting from random spin configurations. After thermalization, $2\times 10^{6}$ measurements are taken every ten Monte Carlo steps at each simulation point.

From the Monte Carlo history of the magnetization (5), we found that the autocorrelation time is shorter than three measurements ($30$ Monte Carlo steps) for $L\leq 64$ and about $10$ for $L=256$. We estimate statistical errors using the jackknife method with $20$ bins, corresponding to the binsize of $10^{5}$ measurements. Since the binsize is several orders of magnitude larger than the autocorrelation time, the effect of autocorrelation is well suppressed in these error analyses.

To search for LYZs on the complex-$h$ plane, we analyze the normalized partition function

[TABLE]

for $h\in\mathbb{C}$. Equation (29) is evaluated using the reweighting method [54, 55] as

[TABLE]

where $\langle\cdot\rangle_{0}$ means the expectation value on the spin configurations at a simulation temperature $T=T_{0}$. Since the denominator of Eq. (29) never diverges, the LYZs correspond to zeros of Eq. (29). To find zeros of Eq. (29), we numerically solve the simultaneous equations

[TABLE]

for fixed $t$ and $L$.

For each $L$ and $T$, we determine the LYZ by solving Eq. (31) at three simulation temperatures $T_{0}$ independently. Statistical errors are estimated by solving Eq. (31) on individual jackknife samples. The final result is obtained as the weighted average over the three $T_{0}$ results. In the analyses of the LYZRs (10), the ratios are computed within individual jackknife samples. This procedure is crucial for proper estimates of statistical errors, eliminating the correlations among LYZs. In the Binder-cumulant analysis, we apply the same procedure, computing cumulants and their ratios within individual jackknife samples.

The left panel of Fig. 2 shows the imaginary parts of the first and second LYZs, ${\rm Im}\,h_{\rm LY}^{(1)}(T,L)$ and ${\rm Im}\,h_{\rm LY}^{(2)}(T,L)$, for $L=96$, 128, 192, and 256 as functions of $T$ around $T_{\textrm{c}}$. The first (second) LYZ for each $L$ is shown by the circle (diamond) symbols. The statistical errors are depicted but much smaller than the symbols and are almost not visible. The vertical dashed line represents the location of $T_{\textrm{c}}$ in Eq. (2).

The right panel of Fig. 2 shows the same results, with both the vertical and horizontal axes rescaled according to the scaling relation Eq. (16), using Eqs. (8) and (2) for critical exponents and $T_{\textrm{c}}$, respectively. The figure shows that the first and second LYZs for various $L$ accumulate on individual universal curves, as expected from Eq. (16). We confirmed that the third and fourth LYZs exhibit the same scaling behavior.

4 Lee-Yang-zero ratios and Binder cumulant

4.1 Intersection analysis

In Fig. 3, we present the behavior of LYZRs $R_{n1}(T,L)$ for $n=2,3,4$, as well as the fourth-order Binder cumulant $B_{4}(T,L)$, as functions of $T$ for various system size $L$. The shaded bands indicate statistical errors. In this study, we consider only $R_{21}(T,L)$, $R_{31}(T,L)$, $R_{41}(T,L)$ for the LYZRs, since the statistical errors of $R_{nm}(T,L)$ typically grow with increasing $n$ and $m$ [17]. The figure shows that $R_{n1}(t,L)$ and $B_{4}(t,L)$ obtained at various $L$ intersect at almost a common point, as anticipated from Eqs. (18) and (25). However, results for smaller $L$ clearly deviate from this trend. The deviation at small $L$ is interpreted as a consequence of the correction to FSS. Comparing the four panels, the convergence toward the $L\to\infty$ limit is the slowest in $B_{4}(t,L)$, suggesting that this quantity exhibits the strongest correction to FSS. We will come back to this point in Sec. 4.2. In Fig. 3, the value of $T_{\rm c}$ in Eq. (2) is indicated by the vertical dashed lines. The temperature at the intersection is close to this value in all cases.

To determine the intersection point quantitatively, we perform the $\chi^{2}$ fits to these data. We use the following fitting function for the LYZRs motivated by Eq. (18)

[TABLE]

where $r_{nm}$, $c_{nm}$, $d_{nm}$, $y_{t}$, and $T_{\textrm{c}}$ are the fitting parameters. For $B_{4}(T,L)$, we use the same functional form

[TABLE]

with $b_{4}$, $c_{4}$, $d_{4}$, $y_{t}$, and $T_{\textrm{c}}$ the fitting parameters.

In our analysis that employs the reweighting method, the $T$ values at which the data for $R_{n1}(T,L)$ or $B_{4}(T,L)$ are referenced can be chosen arbitrarily. In this study, we use three $T$ values for each $L$ for the fits. These values should be close to $T_{\rm c}$ to suppress the nonsingular effects. However, they should be separated so that the correlations arising from the reweighting method on the same configurations are not too large. Following the analysis in Appendix B, for the three $T$ values we choose $T=4.5115$ and $4.5115\pm\Delta T$ with $\Delta T=5\times 10^{-5}$. As discussed in Appendix B, this $\Delta T$ is chosen so that the fit results are insensitive to the variation of $\Delta T$ around this value. To account for the correlations between different $T$ values, we perform correlated $\chi^{2}$ fits, which take into account the correlation between the data points by the covariance matrix.

We notice that the nonlinear term proportional to $d_{n1}$ or $d_{4}$ in Eqs. (32) and (33) is necessary for these fitting analyses. In fact, we have checked that the four-parameter fits with $d_{n1}=0$ or $d_{4}=0$ give unacceptably large $\chi^{2}/{\rm d.o.f.}$ except for $R_{41}(T,L)$. The magnitude of this non-linearity will be discussed in more detail in Sec. 4.3.

In Fig. 3, we show the fit results for the intersection point obtained by the five-parameter correlated fits with Eqs. (32) and (33) by the colored symbols. To estimate the magnitude of the correction to FSS, the fits are performed for various smallest system sizes $L_{\textrm{min}}$ used for the fits, while the largest size is fixed to $L=256$. For the largest $L_{\rm min}=192$, the number of degrees of freedom (d.o.f.) is one. The values of $\chi^{2}/{\rm d.o.f.}$ are shown in the figure, and the resulting values of $T_{\rm c}$, $y_{t}$, and $r_{n1}$ or $b_{4}$ are shown in Table 1 for $L_{\rm min}\geq 48$.

These results show that the intersection points converge well within statistics for $L_{\textrm{min}}\geq 24$, although the convergence of $B_{4}(T,L)$ may be slightly slower than $R_{n1}(T,L)$. The statistical uncertainties of $T_{\rm c}$ in these results are comparable at the same $L_{\rm min}$, whereas those of $R_{41}(T,L)$ are clearly larger than the others. The obtained values of $T_{\rm c}$, $y_{t}$, and $b_{4}$ are all consistent with the known results, Eqs. (2), (8), and (26), within about twice the standard deviation. From these results, we conclude that the LYZR method can determine $T_{\rm c}$ and $y_{t}$ in the Ising model with nearly the same precision as the Binder-cumulant method.

In the following, we use the fit results at $L_{\rm min}=128$ for the final result of the intersection analysis. The resulting values of $r_{n1}$ are shown on the top row of Table 2. In Ref. \citenWada:2024qsk, we studied the CP in the three-dimensional three-state Potts model with an external field using the LYZR method. This CP is believed to belong to the same universality class as the 3d-Ising model [42, 56]. As shown in Ref. \citenWada:2024qsk, the values of $r_{nm}$ in the same universality class should be the same. The values of $r_{nm}$ at the intersection point obtained in Ref. \citenWada:2024qsk are shown on the bottom row in the table. The table shows that they are consistent with each other. The statistical errors, however, are significantly suppressed in the 3d-Ising model.

4.2 correction to FSS

The fit analyses based on Eqs. (32) and (33) in the previous subsection assume the validity of the FSS relations (7) or (9). As Fig. 3 shows, however, our numerical results at small $L$ clearly do not intersect at the same point, indicating the correction to the scaling. In this section, we investigate this effect in more detail and perform the analysis of $r_{n1}$ after eliminating its effect.

In Fig. 4, we show the values of $R_{21}(T,L)$ and $B_{4}(T,L)$ at $T=T_{\rm c}$ determined by the intersection analysis based on Eqs. (32) and (33) at $L_{\rm min}=128$ by the blue diamonds as functions of $1/L$. For comparison, we also plot the same quantities at $T=T_{\rm c}$ in Eq. (2) by the red circles. Provided that the FSS is valid and $T_{\rm c}$ is determined exactly, these quantities should be independent of $1/L$, $R_{21}(T_{\rm c},L)=r_{21}$ and $B_{4}(T_{\rm c},L)=b_{4}$, whose values obtained by the intersection analyses in the previous subsection are shown by the horizontal solid lines in the figure. The black dashed line in the lower panel indicates the previous result of $b_{4}$ in Eq. (26).

We find that the values of $R_{21}(T_{\rm c},L)$ and $B_{4}(T_{\rm c},L)$ in the figure are consistent with each other for $L\geq 96$. However, a clear $L$ dependence is observed in both results for $L\lesssim 32$, which can be attributed to the correction to FSS. The comparison of two panels indicates that the correction to FSS is stronger in $B_{4}(T,L)$ than $R_{21}(T,L)$. This implies that the LYZRs are more effective in suppressing the effect of the correction in the intersection analysis. The red circles may also suggest a remaining weak dependence on $L$ even for $L\geq 96$.

To extract the value of $r_{21}$ while eliminating the effects of the correction to FSS, we perform an extrapolation to the $L\to\infty$ limit. In Appendix A, we investigate the form of corrections to Eq. (18) from irrelevant operators. From Eq. (49), we adopt an ansatz for the $L\to\infty$ extrapolation [7]

[TABLE]

at $T=T_{\textrm{c}}$ in Eq. (2), where $r_{21}$ and $c$ are the fit parameters. The second term represents the correction due to the dominant irrelevant operator. We use the exponent $\omega_{1}=0.8303(18)$ evaluated in Ref. \citenEl-Showk:2014dwa by the conformal bootstrap method. The same fit is also applied to $B_{4}(T_{\rm c},L)$. The red curves in Fig. 4 represent the results of the fit to $R_{21}(T_{\rm c},L)$ and $B_{4}(T_{\rm c},L)$ for four system sizes $L$ in the range $96\leq L\leq 256$. The yellow circles indicate the extrapolated values of $r_{21}$ and $b_{4}$. The value of $r_{21}$ is shown on the middle row of Table 2, together with the results of the same analyses for $R_{31}(T_{\rm c},L)$ and $R_{41}(T_{\rm c},L)$. The table shows that the extrapolated values of $r_{n1}$ are consistent with those obtained by the intersection analyses within statistics, while the extrapolated values would be more reliable, as they take account of the scaling violation.

To estimate the systematic uncertainties of the extrapolation with Eq. (34), we also perform the fits with the variation of $\omega_{1}$ in the range $0.73<\omega_{1}<0.93$ and $L_{\rm min}=64,128$. We take the largest deviation among these analyses as the systematic error of each result, shown as the second error in Table 2. From the table, one sees that the systematic errors are at most comparable to the statistical errors.

4.3 Non-linearity

In the analysis in Sec. 4.1, we found that the non-linearity of $R_{n1}(T,L)$ and $B_{4}(T,L)$ are not negligible in the range of $T$ used in the fit analyses. Such non-linearity prevents the use of linear fits with $d_{nm}=0$ or $d_{4}=0$, which are much simpler and usually employed in the literature [9, 16, 7]. To safely apply linear fits, one must choose a narrower temperature range near $T_{\rm c}$ to suppress the nonlinear effects, which, however, typically increases the statistical errors. In this subsection, we compare the magnitude of non-linearity in different methods. Throughout this subsection, we assume the validity of the FSS and neglect its correction for simplicity.

As a quantitative measure to compare the magnitude of the non-linearity in different methods, we introduce a normalized curvature

[TABLE]

with $f=R_{nm}(T,L)$ or $B_{4}(T,L)$. The meaning of this quantity is understood as follows. First, if one knows the exact value of $r_{nm}$, the critical temperature $T_{\rm c}$ is determined by solving $R_{nm}(T_{\rm c},L)=r_{nm}$. Second, if one evaluates $R_{nm}(T,L)$ at $T=T_{\rm c}\pm\delta T$, an approximate formula of $R_{nm}(T,L)$ by a linear interpolation is given by

[TABLE]

Finally, the solution of $R_{nm}^{\rm linear}(T_{\rm c},L)=r_{nm}$, i.e. an estimate of $T_{\rm c}$ based on the linear interpolation, is given by $T_{\rm c}-C_{R_{nm}}L^{y_{t}}\delta T^{2}$, as shown by the red band in Fig. 5. Therefore, $C_{R_{nm}}$ represents the deviation of the estimate of $T_{\rm c}$ in the linear interpolation with fixed $\delta T$. Using Eq. (18) or (32), Eq. (35) is given by $C_{R_{nm}}=d_{nm}/T_{\textrm{c}}c_{nm}$. The same argument also applies to $C_{B_{4}}$.

The measures $C_{f}$ are not universal quantities specific to a universality class. However, their ratios, such as $C_{R_{21}}/C_{B_{4}}$, in the $L\to\infty$ limit are invariant under the nonlinear variable transformations of $T$ and $h$, as discussed in Appendix C. In this sense, they can be used for measures to compare the magnitude of non-linearity in different methods.

In Table 3, we show the values of $C_{f}$ obtained by the nonlinear fits with $L_{\rm min}=128$ and $\Delta T=5\times 10^{-5}$. One finds that the value of $C_{B_{4}}$ is about four times larger than $C_{R_{21}}$, meaning the larger non-linearity in the Binder-cumulant method. In the LYZR method, $C_{R_{n1}}$ for larger $n$ is smaller for $n=2,3,4$. In particular, $C_{R_{41}}$ vanishes within statistics. These results, possibly due to accidental cancellations of the nonlinear terms in the LYZs, might be useful in some numerical analyses. It is also notable that different combinations of $n$ and $m$ have different characteristics in the LYZR method. For example, $R_{41}(T,L)$ greatly suppresses the nonlinear effect, while it has larger statistical errors than $R_{21}(T,L)$ and $R_{31}(T,L)$ as we have seen in Sec. 4.1.

5 Single-LYZ and single-cumulant methods

We now examine the single-LYZ and single-cumulant methods introduced in Sec. 2.6. Figure 6 peresents the rescaled LYZs, $L^{y_{h}}h_{\textrm{LY}}^{(n)}(T,L)$ for $n=1,2$, and rescaled cumulants, $\langle M^{n}(T,0,L)\rangle_{\textrm{c}}/L^{ny_{h}}$ for $n=2,4$, as functions of $T$, where we use the value of $y_{h}$ in Eq. (8). The figure shows that these quantities for various $L$ intersect around the $T_{\textrm{c}}$ of Eq. (2) indicated by the vertical dashed lines. However, the deviations of the intersection point at small $L$ are larger than those of the LYZRs and Binder cumulant in Fig. 3.

To obtain the intersection point quantitatively, we fit these data using fitting functions motivated by Eqs. (27) and (28),

[TABLE]

where $X_{n}$, $Y_{n}$, $Z_{n}$, $\tilde{F}^{(m)}$ $(m=0,1,2)$, and $T_{\textrm{c}}$ are the fitting parameters. The fitting procedure is the same as before, while in this case we fix $y_{t}$ to the known value in Eq. (8), resulting in four fitting parameters. The dependence of the fit results on $\Delta T$ is discussed in Appendix B.

In Fig. 6, the intersection points obtained from the fits are shown by the colored circles for various $L_{\textrm{min}}$. We find that, even excluding the $L_{\textrm{min}}=192$ results with large statistics errors, the dependence on $L_{\textrm{min}}$ in the single-LYZ and cumulant methods is stronger than that observed in the LYZR and Binder-cumulant methods. From the values of $\chi^{2}/\textrm{d.o.f.}$ indicated in the figure, we also find that the $\chi^{2}/\textrm{d.o.f.}$ of the intersection analysis increases more rapidly as $L_{\rm min}$ decreases. This indicates that the correction to FSS is more pronounced in the single-LYZ and single-cumulant methods. Table 4 lists the normalized curvatures $C_{f}$ for $f=h_{\rm LY}^{(n)}$ and $\langle M^{n}\rangle_{\rm c}$. Their values are significantly larger than those in Table 3, indicating that the magnitude of non-linearity is greater in these quantities.

From these results, we conclude that the single methods are less suitable for precise analyses, since they cannot sufficiently suppress systematic uncertainties arising from scaling violation and nonlinear terms. Nevertheless, they are convenient for obtaining a rough estimate of the location of a CP whose universality class is already known.

6 Summary

We performed Monte Carlo simulations of the three-dimensional Ising (3d-Ising) model to apply the recently proposed Lee-Yang-zero ratio (LYZR) method [17] to this system. The LYZR method enables an intersection analysis analogous to the Binder-cumulant method by exploiting the finite-size scaling (FSS) of the Lee-Yang zeros. We demonstrated that the LYZR method is as powerful as the Binder-cumulant method in the 3d-Ising model. Moreover, we found that both the correction to FSS and the non-linearities are suppressed in the LYZRs compared with the Binder cumulant. The values of the LYZRs at the critical point (CP), which are universal constants, were determined with high precision, as summarized in Table 2.

We also introduced alternative intersection analyses that use only a LYZ while assuming the critical exponent $y_{h}$ as an input. These methods have the advantage that the determination of the second LYZ is unnecessary, and are therefore applicable to systems in which such determination is difficult, such as QCD at nonzero chemical potential [39, 40]. We demonstrated that these methods are useful for studying the CP in the Ising model, although they exhibit stronger corrections to FSS and nonlinear effects.

The LYZR method is applicable to numerical searches for CPs in a wide range of systems, such as magnetic transitions in spin models, CPs in QCD at nonzero baryon density [39, 40] and in the heavy-quark region [58, 59, 60, 61]. It can also be combined with various numerical approaches, such as the tensor-network method, in addition to Monte Carlo simulations. Extending the applications of the LYZR method to these directions constitutes an interesting subject for future work.

{acknowledgment}

This work was supported in part by JST SPRING, Grant Number JPMJSP2110, JSPS KAKENHI (Nos. JP22K03593, JP22K03619, JP23H04507, JP24K07049), ISHIZUE 2025 of Kyoto University, and the Center for Gravitational Physics and Quantum Information (CGPQI) at Yukawa Institute for Theoretical Physics.

Appendix A Contribution of irrelevant operators to FSS

In the main text, we treated the finite-size effects of various quantities near the CP under the assumption that the singular part obeys the finite-size scaling (FSS), Eq. (7) or (9). However, these equations neglect contributions from irrelevant operators. When such contributions are not negligible, Eq. (7) is modified as [9, 16, 7]

[TABLE]

where $\{\lambda_{i}\}$ ($i=1,2,\cdots$) denotes the set of irrelevant operators whose scaling exponents $y_{i}=-\omega_{i}<0$ are universal constants within the same universality class. In this appendix, we discuss the modifications of Eqs. (18) and (25) arising from the contribution of the irrelevant operators. The obtained results are used in Sec. 4.2 to study the correction to FSS in the numerical results. Throughout this appendix, we assume that thermodynamics is dominated by the singular part and neglect the regular part for simplicity.

We first focus on the Binder cumulant [9, 16, 7]. From Eq. (39), the $n$th-order cumulant of the magnetization $M$ is expressed as

[TABLE]

where $\tilde{F}^{(n)}$ denotes the $n$th derivative of $\tilde{F}(\tilde{t},\tilde{h},\{\tilde{\lambda}_{i}\})$ with respect to $\tilde{h}$. At $h=0$, using the expansion

[TABLE]

the fourth-order Binder cumulant, Eq. (22), is expanded as

[TABLE]

around $t=0$, where the coefficients $b_{4}$ and $c_{4}$ are the same as those in Eq. (33). The values of $\lambda_{i}$ are determined by the details of the model. The expansions of $\tilde{F}^{(n)}$ with respect to $L^{-\omega_{i}}\lambda_{i}$ are justified for sufficiently large $L$ as $L^{-\omega_{i}}\lambda_{i}\to 0$ for $L\to\infty$.

Equation (42) shows that nonzero irrelevant operators lead to $L$-dependent contributions to Eq. (25) unless $\lambda_{i}=0$. In particular, the value of $B_{4}(t,L)$ at $t=0$ is $L$ dependent due to the effect of the irrelevant operators,

[TABLE]

with $\tilde{\alpha}_{4,i}=\alpha_{4,i}/b_{4}$. Equation (43) has been employed for the high-precision numerical measurements of the value of $b_{4}$ in Ref. \citenFerrenberg:2018zst. The $L$-dependent terms in Eq. (43) are suppressed for $L\to\infty$, since $\omega_{i}>0$ for irrelevant operators. In the 3d-Ising universality class, the smallest exponent, having the most dominant contribution at large $L$, is known as $\omega_{1}\simeq 0.83$ [57]. In Sec. 4.2, we use Eq. (43) for an estimation of the value of $b_{4}$ (the red line in Fig. 4).

Next, let us extend the argument to the LYZRs. Equation (39) indicates that the contributions of the irrelevant parameters modify the FSS relation for the partition function, Eq. (9), as

[TABLE]

with the modified scaling function $\tilde{Z}$. Since the LYZs, $h_{\textrm{LY}}^{(n)}(t,\{\lambda_{i}\},L^{-1})$, are given by the zeros of the partition function, they satisfy

[TABLE]

Denoting the zeros of $\tilde{Z}$ as $\tilde{h}_{\textrm{LY}}^{(n)}(\tilde{t},\{\tilde{\lambda}_{i}\})$, i.e.

[TABLE]

from Eq. (45) $h_{\rm LY}^{(n)}(t,\{\lambda_{i}\},L^{-1})$ and $\tilde{h}_{\rm LY}^{(n)}(\tilde{t},\{\tilde{\lambda}_{i}\})$ are related with each other as

[TABLE]

where on the second line we expanded $\tilde{h}_{\textrm{LY}}^{(n)}(L^{y_{t}}{t},\{L^{y_{i}}{\lambda}_{i}\})$ by $L^{y_{i}}\lambda_{i}$, and $\tilde{h}_{\textrm{LY}}^{(n)}(L^{y_{t}}{t})$ are the LYZs at $\lambda_{i}=0$ for all $i$, i.e. those introduced in the main text. Then, by adopting the expansion in Eq. (17) to $\tilde{h}_{\textrm{LY}}^{(n)}(L^{y_{t}}{t})$, we find

[TABLE]

where $r_{nm}$ and $c_{nm}$ are the same as those in Eq. (32). The coefficients from irrelevant terms $\{\alpha_{nm,i}\}$ are given by derivatives of $\tilde{h}_{\textrm{LY}}^{(n)}(L^{y_{t}}{t},\{L^{y_{i}}{\lambda}_{i}\})$ and depend on the microscopic details of the model. We note that the effect of the irrelevant operators in Eq. (48) is analogous to that in Eq. (42) for the Binder cumulant. Substituting $t=0$ into Eq. (48), one finds that

[TABLE]

is satisfied at $T=T_{\rm c}$, with $\tilde{\alpha}_{nm,i}=\alpha_{nm,i}/r_{nm}$. In Sec. 4.2, we study the $L$ dependence of $R_{nm}(0,L)$ by taking the first term of Eqs. (49) having the dominant contribution for $L\to\infty$.

We notice that the critical temperature $T=T_{\rm c}$, corresponding to $t=0$, is unique in Eq. (39). On the other hand, if one wants to estimate the critical temperature from the numerical simulation with a fixed $L$, one may, for example, use the peak position of the second-order cumulant. The value of $T_{\rm c}$ introduced in this way is $L$ dependent, which converges to the true critical temperature $T_{\rm c}$ in the limit $L\to\infty$. In the present study, we do not use such a definition of $T_{\rm c}$. Instead, we use $T_{\textrm{c}}$ determined in Ref. \citenFerrenberg:2018zst or by the intersection analysis in Sec. 4.2.

Appendix B $\Delta T$ dependence of the intersection analyses

In this appendix, we examine how the fitting results in Figs. 3 and 6 depend on $\Delta T$.

In Fig. 7, we plot the LYZRs and the Binder cumulant as functions of $T$ for $L=128,192,256$, together with the fit results of the intersection point with fixed $L_{\rm min}=128$ but for various $\Delta T$. The procedure of the fits is the same as that explained in Sec. 4.1. From the figure, one finds that $\Delta T$ dependence is well suppressed for $\Delta T<25\times 10^{-5}$, while the result for $\Delta T=50\times 10^{-5}$ is clearly deviated from the other results. The results for $\Delta T\times 10^{5}=3,5$ are almost overlapped and not distinguishable in the figure. Similar results are obtained also for other fit parameters such as $y_{t}$ and other choices of $L_{\rm min}$. In the main text, we thus employ $\Delta T=5\times 10^{-5}$.

Next, we plot the results of the same analyses for the single-LYZ and cumulant methods in Fig. 8. We find from the figure that the $\Delta T$ dependence of the fit results is nearly the same as in the previous analysis.

Appendix C Invariance of normalized curvatures

In this appendix, we discuss the invariance of the ratios of the normalized curvatures Eq. (35) under reparametrizations of thermodynamic variables. We introduce new thermodynamic variables $\tau$ and $\eta$, which are analytically related to the reduced temperature $t$ and the magnetic field $h$ near the CP $(t,h)=(0,0)$, respectively, such that they can be written as

[TABLE]

where we assume that $(\tau,\eta)=(0,0)$ corresponds to the CP and the linear terms are non-vanishing, $t_{(1)}\neq 0$ and $h_{(1)}\neq 0$. In the following, we show that the ratios of $C_{f}$ are invariant under the variable transformation $(t,h)\to(\tau,\eta)$ in the $L\to\infty$ limit. Throughout this appendix, we assume the validity of the FSS, which is justified in the vicinity of the CP in the $L\to\infty$ limit.

To prove the statement, we begin with the transformation of the cumulants $\langle M^{n}(t,0,L)\rangle_{\rm c}=\partial^{n}F/\partial h^{n}$. In terms of the new variables, the cumulants may be defined as $\partial^{n}{\cal F}(\tau,0,L)/\partial\eta^{n}$ with the free energy as a function of the new variables ${\cal F}(\tau,\eta,L)=F(t(\tau),h(\eta),L)$. They are calculated to be

[TABLE]

with $T_{n}=n(n-1)/2$, $F^{(i,j)}(t,h,L)=\partial^{i+j}F(t,h,L)/\partial t^{i}\partial h^{j}$, and $\tilde{F}^{(i,j)}(\tilde{t},\tilde{h})=\partial^{i+j}\tilde{F}(t,h,L)/\partial\tilde{t}^{i}\partial\tilde{h}^{j}$. The term on the far-right hand side $\tilde{F}^{(0,n)}(L^{y_{t}}t(\tau),0)$ is Taylor expanded as

[TABLE]

where we define $\tilde{F}^{(i,j)}$ without arguments to mean $\tilde{F}^{(i,j)}(0,0)$. Keeping only the leading contribution in each order of $\tau$, Eq. (52) becomes

[TABLE]

The fourth-order Binder cumulant in terms of the new variables is then calculated to be

[TABLE]

where $b_{4}$, $c_{4}$, $d_{4}$ are the same as Eq. (25). From Eq. (55), one finds

[TABLE]

It is not difficult to obtain the same result for the other cumulant ratios.

Next, we focus on the LYZs. We consider zeros of ${\cal Z}(\tau,\eta,L)=Z(t(\tau),h(\eta),L)$ on the complex-$\eta$ plane for $\tau\in\mathbb{R}$, and denote them as $\eta=\eta_{\rm LY}^{(n)}(\tau,L)$. Since the LYZs in terms of the original variables satisfy Eq. (16), $\eta_{\rm LY}^{(n)}(\tau,L)$ satisfy

[TABLE]

Denoting the Taylor expansion of $\eta_{\rm LY}^{(n)}(\tau,L)$ as

[TABLE]

$h(\eta_{\rm LY}^{(n)}(\tau,L))$ on the left-hand side is calculated to be

[TABLE]

Substituting Eq. (59) into the left-hand side of Eq. (57) and comparing terms in each order of $\tau$, we find

[TABLE]

The coefficients ${\cal X}_{n}$, ${\cal Y}_{n}$, ${\cal Z}_{n}$ in Eq. (58) must satisfy these equations. Equation (60) and the fact that ${\cal X}_{n}\to 0$ for $L\to\infty$ lead to

[TABLE]

Then, by taking the leading terms for $L\to\infty$ and using Eq. (63), one finds

[TABLE]

where we used $y_{t}-y_{h}<0$ to obtain the last equality.

From Eqs. (63)–(65), the LYZRs in the new coordinates are calculated to be

[TABLE]

where $r_{nm}$, $c_{nm}$, $d_{nm}$ are the same as those in Eq. (18), which means that

[TABLE]

Comparing Eqs. (56) and (67), we find

[TABLE]

that is, the ratio is invariant under the variable transformation (50) and (51).

Finally, we notice that the ratios of the normalized curvatures are not invariant under more general variable transformations that allow the mixing of two variables, i.e.

[TABLE]

with nonzero mixing elements such as $a_{12}$ and $a_{21}$.