Revisiting the Lee-Yang singularities in the four-dimensional Ising model: a tribute to the memory of Ralph Kenna

We have studied numerically the Lee-Yang singularities of the four dimensional Ising model at criticality, which is believed to be in the same universality class as the $\phi_4^4$ scalar field theory. We have focused in the numerical characterization of the logarithmic corrections to the scaling of the zeros of the partition function and its cumulative probability distribution, finding a very good agreement with the predictions of the renormalization group computation on the $\phi_4^4$ scalar field theory. To obtain these results, we have extended a previous study [R. Kenna, C. B. Lang, Nucl. Phys., 1993, B393, 461] in which there were computed numerically the first two zeros for $L\leqslant 24$ lattices, to the computation of the first four zeros for $L\leqslant 64$ lattices.


Introduction
Since their introduction, the study of complex singularities of the free energy (equivalently the zeros of the partition function) in magnetic field [1,2] and in temperature [3,4] has played a role of paramount importance in understanding criticality in statistical mechanics.
Regarding the complex singularities of the Ising model in a magnetic field, we can recall some important results obtained by Yang and Lee [1,2].They showed that the partition function of the (ferromagnetic) Ising model defined on a regular graph of  vertices (sites) in a complex magnetic field  c has its zeroes on the unit circle of the fugacity variable ( = e  c ), implying a pure imaginary magnetic field  c = i, with real .They also found that in the paramagnetic phase these imaginary zeros condense at a finite value of  in the thermodynamic limit ( → ∞), known as the Yang-Lee edge.Furthermore, this edge reaches  = 0 at the critical point (when the Ising model is above its lower critical dimension), and this is the origin of the singularity of the free energy at criticality.The use of these techniques, based on the study of complex singularities, can be cited both in numerical simulations to characterize phase transitions in statistical mechanics and in quantum field theory (e.g., see [5][6][7]) and analytically, where, for example, they have been employed to develop scaling relations in the presence of logarithmic corrections [8][9][10].
Additionally, the emergence of these logarithmic corrections, particularly in two-dimensional systems and models at their upper critical dimensions [8,9,11,12], has undergone thorough investigation.This research has been facilitated by techniques grounded in the properties of complex singularities.
An illustrative case of these investigations involved characterizing the critical behavior of the fourdimensional Ising model.This model is thought to belong to the same universality class as the  4  4 lattice field theory [13].The significance of  4  4 in high-energy physics is paramount.For example, we can cite its role in the triviality problem [14].
Nevertheless, recent studies have raised questions about this established scenario.In [15], the simulation of the Ising model in both the canonical and microcanonical ensembles suggests the possibility of the specific heat of the four-dimensional Ising model being discontinuous, akin to mean-field behavior, in contrast to its logarithmic divergence in the  4  4 -theory.Furthermore, in [16], employing tensor renormalization group techniques, a weak first-order transition scenario has been proposed1.
The primary objective of this paper is to numerically compute the Lee-Yang (LY) singularities in the four-dimensional Ising model, aiming to comprehensively understand the onset of logarithmic corrections in the scaling of zeros at criticality as a function of the lattice size.In pursuit of this goal, we intend to expand upon the initial work conducted by Kenna and Lang in [18], extending their research to larger lattices and incorporating the computation of two additional zeros (the third and fourth ones).
Furthermore, given that we have computed the first four zeros, we can analyze the cumulative probability distribution of the Lee-Yang zeros and check its associated logarithmic correction.
The overarching objective is to assess whether the numerically characterized logarithmic corrections align well with the analytical predictions derived for  4  4 using the renormalization group (RG).This paper aims to provide support for the conventional understanding that the Ising model and the  4 field theory belong to the same universality class in four dimensions.
The paper is organized as follows: in section 2, we provide an overview of the model and observables.Following this, we present relevant theoretical results in section 3. Subsequently, in section 4, we elaborate on the numerical simulations conducted for this study.In section 5, we present our findings, and the paper concludes with a section summarizing our conclusions.

The model and observables
We have considered the Ising model defined on a four-dimensional cubic lattice with periodic boundary conditions, linear size  and volume  =  4 .The Hamiltonian of the model is given by where   are Ising variables and the sum in equation (2.1) runs over all pairs of lattice nearest-neighbors.
As stated in the Introduction, in order to compute the LY zeros we add a pure imaginary magnetic field to the model, the new Hamiltonian being with  ∈ R. For further use, we define the magnetization as The partition function of the Hamiltonian H  [equation (2.2)] is that can be written as with ℎ ≡  and the average ⟨(• • • )⟩ is computed in absence of the magnetic field (i.e., using the Hamiltonian of equation (2.1)).This fact implies ⟨sin(ℎ)⟩ = 0 for all values of ℎ and that we will restrict our computation to the case ℎ > 0. Hence, the pure imaginary complex singularities in the magnetic field (YL singularities) are determined by the solutions of the equation ⟨cos(ℎ)⟩ = 0 . (2.6) We will denote the -th solution of this equation as ℎ  , ordered in such a way that ℎ +1 > ℎ  .To solve this equation, we have used the bisection method.
Once we have computed the LY zeros (ℎ  ), the cumulative distribution function of zeros can be easily computed as [19] where  is the dimensionality of the space ( = 4 in this paper).

Some theoretical results
In this section we collect some analytical results obtained on the  4 4 field theory using RG, relevant to the scaling of the LY zeros [8,9].
The scaling behavior of the -th LY zero with the reduced temperature ( = ( −   )/  ,   being the critical temperature) is in the paramagnetic phase (i.e., for  > 0).The specific heat also behaves as At the critical point the cumulative distribution function of the zeros scales (in the thermodynamic limit) as [8,9] By inverting this equation, we obtain Since  ∼  − , the logarithmic correction of ℎ versus  is the same as that of ℎ versus , and the behavior at the infinite volume critical point of the LY zeros in a finite box of size  is ℎ  () ∼  −Δ/(2− ) (log ) −  . (3.5) The values of the reported exponents in the (one-component)  4  4 theory are:  = 0, Δ = 3/2, α = 1/3, Δ = 0 [8,20].Hence, we can write the following asymptotic relations for the four-dimensional Ising model (which will be tested in the numerical part of this paper): and Inverting the previous equation [or using equation (3.4)], we obtain Since the main aim of this paper is to characterize the logarithmic corrections of the LY zeros, we assume the RG values for Δ and  exponents and write equation (3.6) as: in this way we try to evaluate the  exponent analyzing the numerical data and check if the data support the RG prediction:  = 1/4.

Numerical simulations
We have investigated the model defined in equation (2.1) through equilibrium numerical simulations.Specifically, we have brought to thermal equilibrium our systems by combining cluster and local update algorithms in order to compute the LY singularities.
In particular, we have used a combination of the Wolff's single cluster algorithm [21] with Metropolis updates [22].Our elementary Monte Carlo step on a lattice of size  is composed by  2 Wolff's singlecluster updates and subsequently by a full sweep Metropolis actualization of the lattice.
Furthermore, we have performed all our simulations at the infinite volume critical inverse temperature   = 0.14969383(6) [15].Other previous estimates were:   = 0.149694(2) [17,23] (based on numerical simulations and analysis of high-temperature series analysis) and that of Kenna and Lang   = 0.149703 (15) [18] (based on the numerical study of complex singularities)2.
To check the thermalization of our runs, we have monitored some non-local observables as the second-moment correlation length and the susceptibility as a function of the Monte Carlo time.We have simulated a large number of runs on the same lattice (different initial configurations and random numbers)3.In addition, we have computed the statistical error using the jackknife procedure [24] merging the different pseudosamples in ten groups.
Finally, details about the numerical simulations can be found in table 1.

Results
In the first two parts of this section we analyze the scaling of the four first zeros.In the third part we present the analysis of their cumulative probability distribution.

Scaling of the zeros
In table 2 we present the numerical values of the first four LY zeros for the seven lattice sizes simulated4.In figures 1 and 2 we show the dependence with  and the associated logarithmic corrections [see equation (3.6)] by plotting ℎ   3 (log ) 1/4 (that should be constant) versus 1/.2In a recent paper [16] a very large four-dimensional Ising model, 1024 4 , was numerically analyzed using the tensor renormalization group reporting a critical inverse temperature of   = 0.1503677(1) (statistically) very different from all the previous published values.
3This strategy allows us to optimize the use of clusters with a large number of processors: the obvious drawback is the need of thermalizing the different pseudosamples.However, the final outcome is positive.
4It is not easy to compare numerically our values reported in table 2 with those of Kenna and Lang [18].First, we use a slightly different value of the inverse temperature.And second, they used the spectral density method to find the zeros in near values of the temperature which could bias the final results.Overall, we found a good agreement with their computed zeros.Only for the higher-order zeros (  = 3 and  = 4) we obtain ℎ   3 (log ) 1/4 ≃ const for all the lattice sizes simulated (see figure 2) taking into account the statistical errors of the data.This initial analysis shows that the behavior of ℎ   3 (log ) 1/4 ≃ const holds, hence, strongly supporting the RG predictions.Notice that for the lower-order zeros, the points of figure 1 show a small variation of about 2-3%.
As did by Kenna and Lang [18], we can try to obtain the  exponent [defined in equation (3.9)].To do that, we plot in figure 3, log(ℎ   3 ) versus log log  to extract the  exponent via a linear fit.
The different slopes plotted in this figure are reported in table 3. Leaving free the  exponent, we have obtained good fits (see the first four columns of table 3).However, the value of the  exponent is clearly different from the RG prefiction for the lower-order zeros.In addition, fixing the exponent  to 1/4 in the fit to equation (3.9) worsens the situation for the lower-order zeros (see also figure 1).
The rationale for this behavior is as follows.We have been able to compute with high precision the lower-order LY zeroes, while the higher-order ones present, in comparison, greater relative errors.This means that the possible subleading behaviors (for example, corrections to scaling or, see below, improved scaling for the leading term) are hidden inside their error bars.However, the very small statistical errors of the lower-order zeros allow us to start to see and analyze (see below) these additional behaviors.
In the next subsection we try to improve this analysis.3. The size of the error bars is smaller than the ones of the symbols.

Two improved scalings of the zeros
In the previous subsection, to analyze the scaling of the zeros, we have used the asymptotic scaling form given by equation (3.6).In this subsection we re-analyze the data using two kinds of improvements.The first one is to write the scaling relation improving the dependence on the logarithm (as suggested by the theory), and the second one is based in the introduction of the leading correction to scaling in the fits.
The starting point in our first approach is to note that the integration of the -function5 for the coupling  in the  4  4 -model will provide a factor (1 +  log ) instead of its asymptotic form given by log  for large  ( being the scaling factor in the RG approach) [25]: 1/() ∝ 1 +  log .Using this fact, we 5 () ≡ d/d log .For the four-dimensional Ising model, the leading order of this function is  = − 2 , with  > 0 [25].
Table 3. Value of the  exponent as a function of the order of the zero using equation (3.9).The first column is the order of the zero, the second column is the minimum value of the lattice size used in the fit, the third column reports the  exponent obtained in the fit, and in the next column we report the goodness of the fit via the  2 /d.o.f.(where d.o.f.stands for number of degrees of freedom of the fit).In the last two columns, we report the minimum lattice size and the goodness of the fit  2 /d.o.f. by fitting the data to equation (3.9) with  fixed to its RG value ( = 1/4).For all the fits reported in this paper, we have computed the minimum value of , in such a way the -value of the least-squares fit is greater than 5%.Notice that for the first zero and fixed  = 1/4, we have been unable to find an acceptable fit (i.e., we have obtained a -value < 5%).

Free 𝜃
Fixed for a suitable constant  that should be independent of the order of the zero6.
In figure 4 we present the behavior of ℎ   3 versus  to check this improved scaling behavior and the fits to equation (5.1) with  fixed to its RG value of 1/4 (i.e., the improved leading prediction of the RG).In addition, in table 4 we report (see the last two columns) the statistical quality of these four fits and the small lattice size for which we have obtained acceptable statistical fits.Notice that the four zeros behave accordingly as stated by the improved RG prediction (even the lower ones).Hence, this improved scaling is capable of capturing the behavior even for the lower-order zeros despite their small statistical errors.
However, we have found that the value of  is roughly independent of the order of the zero ( = 2.0(1) for  = 2 and 3.5 (7) for  = 3, unfortunately, for  = 4 the relative error for  is greater than 100%) but not for the first zero where  = 0.64 (1).The explanation of this fact is that we would need to take into account the leading scaling correction of the model (see below).
Finally, we have tried to obtain the value of the  exponent using this improved scaling.We present the results of these fits in the 2nd to the 4th columns of table 4. The obtained  exponents are compatible with the RG prediction.However, being a three parameter fit, the statistical errors associated to the different estimates of  are large.Table 4. Value of the  exponent as a function of the order of the zero using the improved scaling relation given by equation (5.1).The first column is the order of the zero, the second column is the minimum value of the lattice size used in the fit, the third column reports the  exponent obtained in the fit, and in the next column we report the goodness of the fit via the  2 /d.o.f.In the last two columns, we report the minimum lattice size and the goodness of the fit  2 /d.o.f. by fitting the data to equation (5.1) with  fixed to its RG value ( = 1/4).

Free 𝜃
Fixed  = 1/4 6We thank an anonymous referee for pointing out to us this improved scaling relation.).The size of the error bars is smaller than the ones of the symbols.Now, we analyze the data using the second method of improvement based in considering the effect of the leading scaling correction (at the infinite volume critical point) which induces a correction term proportional to7  ( 1/2 ) ∝ 1/ √︁ log .Hence, the scaling of the -th zero can be written as In table 5 we report the results of our fits.Firstly, we have tried to perform the fits to equation (5.2) with a free value of  (see the 2nd, 3rd and 4th columns of the table 5).Notice that the best fits for the lower-order zeros provide lower values of the  exponent, as in the previous approaches.
In addition, we have fitted the data against equation (5.2) with  = 1/4 (the RG prediction with the leading correction to scaling term).The results of these fits can be seen in columns 5th and 6th of table 5. Notice that the dependence on the lattice size of all the zeros can be described by the RG prediction, including the leading correction to the scaling term.
To illustrate this point, we plot this fit for the first zero (the most problematic case) in figure 5.The scaling is very good for the first zero taking  ⩾ 16, and we remark that the data for this zero support the RG prediction.
To sum up, both improved methods produce similar (good) effects on the different fits fixing  = 1/4 (see the last two columns of tables 3, 4 and 5).Obviously, both must be present in the analysis.However, the fits using both improvements in the same equation yield fitting constants with higher errors.
7The susceptibility can be computed using the behavior of the lowest zero (see, for example, [10]):  ∼  − ℎ =1 ( ) −2 .Hence, the corrections to scaling of the LY zeros are the same as those of the susceptibility.The analysis of the leading correction to scaling for the susceptibility at the infinite volume critical point can be found, for example, in equations (4.24) and (6.28) of [25] on the four-dimensional Ising model or in equation ( 39) of [26] on the four dimensional diluted Ising model adapted to the pure case.Notice that if one considers the next leading order contribution,  ( 3 ), to the -function, an additional correction term  (log log /log ) appears, this behavior was obtained by R. Kenna in 2004 [27].5).
Table 5. Value of the  exponent as a function of the order of the zero using the improved scaling relation given by equation (5.2).The first column is the order of the zero, the second column is the minimum value of the lattice size used in the fit, the third column reports the  exponent obtained in the fit, and in the next column we report the goodness of the fit via the  2 /d.o.f.In the last two columns, we report the minimum lattice size and the goodness of the fit  2 /d.o.f. by fitting the data to equation (5.2) with  fixed to its RG value ( = 1/4).

Density of the zeros
In this subsection, we check the RG prediction (with no improvements) in a different way, by introducing a  parameter defined as8 If the prediction of RG holds, then  = 1.
In figure 6 we plot ℎ  against  3/4 /(log()) 1/4 , where we have computed  using the first four zeros and the seven lattice sizes.The scaling of all the zeros and sizes is very good, all the points falling in a common curve.
We have performed a quantitative analysis of the slope (in a double logarithmic scale) in order to compute .The results are reported in table 6.
Notice that to avoid the use of the full covariance matrix in the  2 minimization fits (the data of the zeros belonging to the same lattice size are correlated), we have performed independent fits for all the 8Notice that we test equation (3.8) instead of equation (3.7).The reason is that it is simpler to fit the data points ( , ) with the statistical in the -variables.In this case, we have computed the statistical error on the ℎ  variables, so we use them as -variables in the fit and consequently we test equation (3.8).Finally, we remark that if the numerical data behave as  (ℎ) ∼  1 ℎ  2 with  2 ≃ 1 (for small ℎ), we will be in presence of a first order phase transition with a discontinuity in the magnetization given by  1 .However, our numerical data are compatible with  2 = 4/3 (modified by logarithmic corrections with a positive power) [19] and we are unable to see a trend towards a first order transition in the analysis in the cumulative probability distribution data.(Colour online) ℎ  (for the first four zeros and the seven lattice sizes simulated) versus the analytical prediction  3/4 /(log()) 1/4 ,  is the cumulative probability distribution of the zeros.The fit has been computed (continuous line) using only the fourth-zero data in a double-log scale.The slope of the straight line is the -parameter.We report in table 6 the values of the -parameter as a function of the order of the zero.The size of the error bars is smaller than the ones of the symbols.

Conclusions
By numerically computing the first four LY zeros of the four-dimensional Ising model, we have been able to extend previous numerical results but also to unveil the onset of the appearance of the logarithmic corrections associated with these observables.
Specifically, the analysis of the scaling behavior of these zeros with the lattice size indicates that the logarithmic corrections, as foreseen by the renormalization group approach in a  4  4 continuous field theory, manifest distinctly within the simulated lattice sizes (where  ⩽ 64).In particular, the scaling is much better using the improved versions of the scaling of the LY zeros (changing log  by 1 +  log  and the introduction of the leading scaling correction) and fixing  = 1/4.
Another way to study the scaling of the LY zeros is to analyze the properties of the cumulative probability distribution.Our findings demonstrate once again a robust alignment with the predictions of the field theory.Notably, as in the ℎ  () analysis, we observed that even for relatively modest lattice sizes, the zeros exhibit the expected logarithmic power dependencies.We remark that the analyses of ℎ  () and ℎ  () are not independent.
In particular, the behavior of ℎ  () is inconsistent with a second order phase transition with a specific heat discontinuity at the transition point ( = α = 0) and also with a weak first order phase transition.Notice that we have tested the following combination of critical exponents (just in the scaling of ℎ  versus  or in the scaling of ℎ  versus ) Δ − Δ α/(2 − ) which mixes the odd (Δ, Δ) and even (, α) sectors of the theory, and this combination of critical exponents is sensitive to the behavior of the specific heat (exponents  and α).Furthermore, we would like to remark that we have been unable to detect a departure of the scaling of the zeros (as a function of the lattice size) from the RG predictions.
Finally, the results presented in this paper are fully compatible with a second order phase transition for the four-dimensional Ising model with the exponents and logarithmic exponents predicted by the RG analysis of the  4  4 continuous field theory.

Figure 1 .
Figure 1.(Colour online) We show ℎ   3 (log ) 1/4 versus 1/ in order to check the logarithmic correction [see equation (3.6)] for the -th zero.If the data are in the asymptotic regime and the analytical prediction is correct, then ℎ   3 (log ) 1/4 ≃ const.In this figure we plot only the behavior of the first two zeros.To have the data of these two zeros in the same figure, we have shifted vertically the  = 1-points by an amount of 1.55.Notice that the tiny error bars allow us to see the subleading corrections (see the text).

Figure 2 .Figure 3 .
Figure 2. (Colour online) We show ℎ   3 (log ) 1/4 versus 1/ in order to check the logarithmic correction [see equation (3.6)] for the -th zero.In this figure we plot only the behavior of the third and fourth zeros.To have the data of these two zeros in the same figure, we have shifted vertically the  = 3-points by an amount of 1.05.Notice that all the lattices follow the analytical prediction (no dependence with 1/) taking into account the error bars.

3 LFigure 4 .
Figure 4. (Colour online) We show ℎ   3 versus  in order to check the improved scaling relation [see equation (5.1)] for the first four zeros.The continuous lines are fits to equation (5.1) with  = 1/4 (see table4).The size of the error bars is smaller than the ones of the symbols.

G 3 / 4 / 4 Figure 6 .
Figure 6.(Colour online) ℎ  (for the first four zeros and the seven lattice sizes simulated) versus the analytical prediction  3/4 /(log()) 1/4 ,  is the cumulative probability distribution of the zeros.The fit has been computed (continuous line) using only the fourth-zero data in a double-log scale.The slope of the straight line is the -parameter.We report in table 6 the values of the -parameter as a function of the order of the zero.The size of the error bars is smaller than the ones of the symbols.

Table 1 .
Parameters of the numerical simulations.isthe lattice size,  therm is the number of thermalizations (composed by 2Wolff updates complemented by a Metropolis sweep),  steps is the number steps after the thermalization (in which we measure every 10 steps) and  pseudo is the number of pseudosamples (see the text).Hence, for a given lattice, we have performed  pseudo ×  steps /10 measures.

Table 2 .
Values of the first four computed zeros for all the seven simulated lattice sizes.

Table 6 .
The  parameter as a function of the order of the zero.The meaning of the different columns is the same as in table3.which are uncorrelated), obtaining four estimates of the  parameter.Overall, we find that the values of  point to  = 1, accordingly with the RG prediction.