Exploring Lee-Yang and Fisher Zeros in the 2D Ising Model through Multi-Point Pad´e Approximants

,


INTRODUCTION
The knowledge of phase transitions is essential in determining the equation of state of a physical system.Phase transitions are characterised by non-analyticities that develop in thermodynamic functions in the limit of infinite degrees of freedom in a system.In finite volume systems, thermodynamic functions like the free energy are analytic functions of system parameters and the study of the onset of divergences becomes challenging.However, there are remnants of these divergences (if existing) even in finite volumes and occur as peaks in the susceptibility of the order parameter, whose height and width scale with volume.The conventional way to look for phase transitions in finite systems is to study the finite size scaling of these susceptibilities [1][2][3][4][5][6].Although these techniques are still used extensively in the study of phase transitions in the fields of condensed matter and lattice gauge theories, one of our goals in this paper is to show the power of an alternative scheme in extracting critical exponents for certain theories from their numerical simulations at finite volumes, namely, the Lee-Yang (LY) zero analysis [7][8][9][10][11].While results on the numerical extraction of LY zeros from some lattice models already exist [12][13][14][15][16][17][18][19][20][21][22][23][24], the novelty of our work is to extract them using only the leading order cumulants of the partition function, evaluated at multiple points in parameter space, using a rational function re-summation of the cumulants called the "multi-point Padé" method.
Another motivation for this work comes from our recent studies of the QCD phase diagram in the complex chemical potential plane.We have shown [25] that it is possible to gain more information from the generated Taylor coefficients, if we re-sum them into a rational function.
In the present paper, we aim at providing confidence in the claim that the multi-point Padé method can extract genuine LY zeros, illustrating it in the 2D Ising model.The Ising model is a simple physical system displaying a phase transition in dimension D ≥ 2 .Although an exact solution of the 2D Ising model exists due to Onsager in the thermodynamic limit, numerical simulations usually have to be performed for any finite size results.These have to be a fortiori accurate if we are interested in a numerical study of divergences arising in the thermodynamic limit.These simulations are based on Monte Carlo methods, and hence are plagued by statistical errors.The presence of such errors is one of the reasons we want to apply the multi-point Padé analysis to study the LY zeros; we want to check whether the tool is robust enough to live with finite accuracy of input data.A preliminary analysis leading up to this work can be found in [26,27].
We will begin our analysis by giving an overview of the complex zeros of the partition function as applied to the Ising model, putting emphasis on the well known properties of the LY zeros in Section II.In Section III we will briefly describe the simulation procedure used for the 2D Ising model and outline the phase transitions we are looking for.In Section IV, we will give a brief overview of the multi-point Padé method and the error analysis procedure used.In Section V we will describe results for the scaling analyses used.In Section VI we will conclude with some results and outlook.

II. COMPLEX ZEROES OF THE PARTITION FUNCTION
We begin by explaining the origin and nature of complex zeros of the partition function, in the context of the Ising model in the presence of an external magnetic field H.The zeros of the partition function in the complex H plane are known as the Lee-Yang (or Yang-Lee) zeros [7,8].We can also discuss the zeros of the partition function in the absence of an external magnetic field, called the Fisher zeros [10], which appear in the complex inverse temperature β plane.However, some of the properties shown by these zeros, like the Circle theorem [7,8] only hold for zeros in the complex H plane, while other properties are common to both kinds of zeros.
We will focus our attention on LY zeros in the following.
We further restrict our study to nearest neighbour ferromagnetic interaction.The Hamiltonian describing this theory is given by (with J > 0) where J is the coupling and σ i ∈ {±1} are the spins at the site i.The canonical partition function is given by: From the form of the Hamiltonian, it can be seen that up to an overall functional dependence on z = e βH , the partition function is a polynomial in e βH of order 2N , with N being the number of sites on the lattice [28].Therefore, the order of the polynomial grows linearly with the number of sites.This remark will become important in the discussion to follow on the relation of LY and Fisher zeros to phase transitions.In order to see why the partition function is a polynomial, consider the term e βH i σi .For N lattice sites, the partition function will be a sum over 2 N terms weighted by the factor e βJ ⟨ij⟩ σiσj that depends on the spin configuration.Note that this factor is invariant under the transformation of flipping all the spins on the lattice.Additionally, there will be a permutation factor associated with the number of spins pointing up or down which will be the same for configurations related by spin flips [29].
Taking these factors into account, the partition function can be written as: The symmetry of the partition function under z → z −1 is the statement that the Hamiltonian in Eq. ( 1) is invariant under the combined action of reversing the external magnetic field and flipping all the spins.A few comments can be made based on Eq. (3) above : I The partition function has the functional form of an even polynomial multiplied by a factor z −N .We could factor the partition function this way because z can never be zero.
II For real (β, H), the coefficients of the polynomial are strictly positive.This implies that complex roots always occur in complex conjugate pairs, in the complex H plane, for fixed, real β and in the complex β plane for fixed, real H.
III As the lattice volume N increases, the order of the polynomial increases, which leads to an increase in the number of complex zeros of Z N .In the limit of N → ∞, these zeros accumulate and coalesce into cuts.A phase transition does not occur when the zeros do not approach the real axis when increasing N .
IV A phase transition is said to occur at some critical temperature T crit when these zeros approach the real axis of the external field parameter, in the thermodynamic limit.The behaviour of the density of (Lee-Yang) zeros at the real axis is used to distinguish between a first and a second order phase transition.At a second order transition the complex conjugate pair gets infinitesimally close to the real axis but the density of zeros is zero on the real axis.On the other hand, a first order transition sees a non-zero density of zeros on the real axis.
V Furthermore, because of the even nature of the polynomial, if z is a root, then so is −z.
Therefore, for a finite N , the partition function is strictly positive when β and H are real.Hence, the only zeros of Z N occur in the complex plane of z or on the z < 0 axis.In the complex H plane this translates to having only complex zeros of Z N (z(H)).However, in the infinite volume limit, Z N is an infinite series which can now have non-trivial zeros on the z > 0 axis, or real zeros in H. Since real zeros of the partition function mark the onset of phase transitions, we have recovered the well-known result that phase transitions cannot occur in a finite volume.The natural question to ask now will be on how to access these complex zeros of the partition function.This is fortunately not hard to answer because thermal cumulants, like the average magnetization in the case of the Ising model, are related to the derivative of the logarithm of the partition function (see Section V).This means that the zeros of Z N (H) will appear as poles of the average magnetization as a function of H and the specific heat capacity as a function of β.The aim of the next few sections is to study in some detail these poles in the complex plane of the external magnetic field and inverse temperature.Arguments similar to those presented in this section can also be found in more detail in the existing literature [11,30], and we refer the reader to [20] for a detailed review on the subject.

III. SIMULATING 2D ISING MODEL
The Ising model has been extensively studied in the literature and can also be found in many textbooks [31][32][33] on statistical mechanics.The model has an exact solution in 1D due to Ising [34] and does not undergo any phase transition.However, its 2D version is one of the simplest systems to undergo a continuous (or 2 nd order) phase transition from a symmetry broken (ferromagnetic) phase to a symmetric (paramagnetic) phase at 2) ∼ 2.269 J.The model can also be seen to undergo a discontinuous (or 1 st order) phase transition when considering the average magnetization (Eq.5), as a function of H at T ≤ T c , across H = 0.An exact solution for this model in the thermodynamic limit exists due to Onsager [35].Hence we know the transition temperature and the various critical exponents characterizing this transition analytically.Because of this, the 2D Ising model serves as an ideal candidate for testing new numerical methods.In general, extracting critical exponents from numerical data is a non-trivial task and requires formidable amounts of statistics.Thankfully, the 2D Ising model is easy to simulate and not expensive where getting large statistics is concerned.Hence, we choose to test our multi-point Padé method on this model.The 2D Ising model can be simulated using both single spin flip [36] and cluster spin flip algorithms [37,38].It is well known that single spin flip algorithms suffer from critical slowing down [39], and since our goal is to extract LY zeros close to and at the critical temperature T c , we will use a cluster spin flip algorithm based on [40], where a modification to the original Swendsen-Wang algorithm [37] was described, to add an external magnetic field parameter in the code.We will discuss below two types of simulations performed.Note that the ferromagnetic coupling constant J has been set to 1 for the simulations that follow.

A. Simulation : To extract Fisher zeros
For the first part of our analysis, we want to study the Fisher zeros -zeros of the partition function in the complex β plane.For this we compute the energy density ⟨E⟩ and the specific heat capacity C H of the 2D Ising model, with H set to 0. These quantities are given by with E = −J ⟨i,j⟩ σ i σ j .The model is simulated at temperature values given by T ∈ {1.76, . . .2.15} in steps of 0.03, {2.17, . . .2.40} in steps of 0.01 and {2.43, . . .3.00} in steps of 0.03 with H = 0.These are then repeated for different lattice volumes at L ∈ {10, 20, 40, 64, 80}, to perform the finite volume analysis of the Fisher zeros obtained.
We refer the reader to Fig. 1 for the results of the simulation for ⟨E⟩ and C H .As an indication for the kind of statistics used, we list here the number of configurations used per lattice size.For L ∈ {10, 20, 40, 64, 80}, the following number of configurations were used respectively: {300K, 125K, 125K, 40K, 25K}.Before analysing this data, we detail the simulations required for the LY zero analysis.

B. Simulation : To extract Lee-Yang zeros
For the second part of the analysis, we want to study the zeros of the partition function in the complex H plane.For this we need to compute the average magnetization ⟨M ⟩ and susceptibility χ H of the model, given by with M = i σ i .In order to verify the volume scaling of the LY zeros, we simulate lattice volumes of sizes L ∈ {10, 15, 20, 30}.Each lattice volume was simulated for H ∈ {−0.125, . . .0.125}, in steps of 0.005 at T = T c ∼ 2.269 J.
The resulting ⟨M ⟩ and χ H , per lattice site, are shown in Fig. 2 and the details of the number of configurations are as follows: For each lattice size, 625K configurations were used to estimate the average magnetization and the resulting susceptibility.We will now proceed to describe the multi-point Padé method for extracting the zeros of the partition function from the poles of the C H and ⟨M ⟩ data.

IV. MULTI-POINT PAD É METHOD
Padé-type rational approximations have recently (re-) emerged as a reliable tool to re-sum Taylor series coefficients in the studies of lattice QCD [25,[41][42][43][44][45].Most of the literature that exists toward the existence, convergence and uniqueness of solutions for Padé approximations is limited to only a restricted class of functions to be approximated [46][47][48][49][50][51][52][53].Instead, many of the interesting results on Padé approximants are known only due to numerical experiments like those in [43,[54][55][56].However, most of the literature referred to above is based on what is called the "single point Padé" approximation.This involves using Taylor series coefficients of the unknown function about a single point and constructing a rational approximation using these coefficients.An immediate limitation of this method is the need to have relatively high number of Taylor coefficients to build rational functions of increasing order.Numerical simulations typically do not allow for the generation of such high order Taylor coefficients with reasonable accuracy.One can instead use lower order Taylor coefficients of the function evaluated at multiple points, which forms the basis of our analyis.This trade-off allows us to construct high order rational approximants.In the following we will only focus on the construction of such approximants and briefly outline the sources of error.For a more detailed discussion on this method and results of numerical experiments, we refer the reader to [25].

A. The multi-point Padé method in a nutshell
Let us consider the rational function R m n (z) with m and n being the degrees of the polynomials at numerator and denominator respectively.Writing Qn (z) = 1 + Q n (z) ensures that the rational function depends essentially on n + m + 1 parameters.Notice that one should naturally also demand that there is no point z 0 such that P m (z 0 ) = Qn (z 0 ) = 0, i.e. we should in principle exclude any (common) zero of both numerator and denominator.If this were not the case, we would have rather essentially defined the rational function R m ′ n ′ (z) with n = n ′ + l and m = m ′ + l for some integer l > 0. Despite this, we will nevertheless not exclude the possi-bility of common zeros, for a reason that will be clear when we explain why we are interested in R m n (z).Let us consider now a function f (z) and suppose we know a few of its Taylor expansion coefficients, let's say at different points {z k | k = 1 . . .N }.The number of coefficients we know can be different at different points, but for the sake of simplicity we will assume that f (s−1) is the highest order derivative which is known at each point, together with all derivatives of degree 0 ≤ g < s − 1.Then the system of equations we have to solve becomes In what follows, we will solve this system of linear equations to determine the coefficients of the polynomials P m and Q n while restricting our analysis to diagonal ([q, q]) and near diagonal ([q, q + 1]) type Padé approximants.
Other techniques like a generalized χ 2 method can also be used to estimate the coefficients of the rational functions by minimizing the distance between the measured Taylor coefficients and the required rational function, weighted by the estimated errors on the measured coefficients.This has been compared to the linear solver method in [25].

B. The method at work for the 2D Ising model
We will now focus on the results of the approximation and the singularity structure obtained from the Padé procedure outlined in the previous sub-section, for the average magnetization.We will first show the results of the approximation in Fig. 3.That we can see the rational functions approximate the data correctly is not surprising because we have essentially done a rational interpolation through the data, since our input only consisted of the average magnetization values and not the susceptibilities.Therefore, the first real success of the rational approximation is to see how well its derivatives approximate the susceptibilities, as shown in Fig. 4 for the L = 10 , 15 data.This is a nice result because the rational function was constructed assuming only the knowledge of the zeroth order Taylor coefficients in the expansion of the average magnetization and it faithfully returns the expected first order coefficients at all the input points [57].
Having gained some confidence in our multi-point Padé approximation, we can now proceed to look at the singularity structure, i.e. we will now study the zeros and poles of the rational function constructed in the complex H plane [58].
Selecting two lattice sizes at T c , we refer the reader to Fig. 5.In the figure, we depict the zeros of the numerator (black pentagons) and of the denominator (red crosses) of R m n (H) at different values of the lattice size L, i.e.L = 15 (left panel) and L = 30 (right panel).The order the rational approximant constructed is [m, n] = [25,25].The pale blue points are shown to indicate the interval of points where the value of average magnetization was used as an input to build the rational function.We can easily make a couple of key observations.
• A few zeros of the denominator are canceled by corresponding zeros of the numerator.These are not genuine pieces of information: actually their location can vary when varying e.g. the order of the Padé approximant [m, n].On the other hand, genuine pieces of information (i.e.actual zeros and poles) stay stable to a very good precision.Notice that the alternating sequence of zeros and poles along the imaginary H axis is how a Padé approximant alludes to a branch cut eventually showing up in the thermodynamic limit [59].
• In order to analyse the genuine poles more carefully, we remove the spurious poles from the singularity structure.We further remove all remaining zeros to focus the attention on poles and refer to Fig. 6.Looking at the pole structure of the figures, we can make the following comments : (1) As the lattice size L gets larger, the closest singularity H 0 gets closer to the real H axis.We claim that we verify the circle theorem, which translates to the Lee-Yang zeros lying on the imaginary H axis.However, looking at Fig. 5 the reader may notice that some of the poles for the L = 30 lattice are shifted away from imaginary H axis, and hence seemingly violating the H → −H symmetry.This is because we have used the central values of magnetization which were simulated at each quoted value of H and hence do not respect the anti-symmetry exactly, but do so within errors.In support of our claim we provide the location of the genuine poles of magnetization at the central values along with error bars in Fig. 6, obtained by performing a bootstrap procedure, i.e. for a given point we extract new data from a Gaussian distribution with mean given by central value and standard deviation given by the estimated error.Notice that only the closest pole tends to appear stable.We further observe that for the smaller L = 15 lattice, the second pole has some uncertainty easy to inspect by eye.For the L = 30, we observe a third pole which now appears with even larger uncertainty.In the following, only the closest pole is important for the discussions of scaling.(2) Referring the reader to Section II, all poles occur with their complex conjugate pairs.(3) Although this observation is highly dependent on the statistics used, which for the purposes of this work is relatively high, we can see an increase in the number of zeros as the number of lattice sites increases.All of these points seem to indicate that we are observing genuine Lee-Yang zeros.
C. On the stability of the closest poles to the real H axis Before moving on to discuss the scaling of these zeros, which will make use of the closest poles extracted for each lattice size, we would like to briefly discuss the procedure we have used to attach error bars on the locations of the poles in the complex plane of H and β.An important point to make is that in the absence of noise in data (e.g. one can construct this by discretizing a known function), if there is a genuine singularity of the function, it will appear as a stable pole of the rational approximation constructed.However, if the data is noisy, as it always is when dealing with simulation results, even if there is a genuine singularity of the function, the resulting pole will move, commensurate with the amount of noise present.
Here, we can distinguish between two types of errors the poles can have, although they are not strictly independent.
• Statistical errors: These are the errors propagated from the estimated error on the measured Taylor coefficients to the poles of the rational function approximant.The procedure used to estimate this error was to solve the system of linear equations in Eq. ( 7) repeatedly by choosing new Taylor coefficients for each solve.These coefficients are drawn from a Gaussian distribution centered at the central values of the measured coefficients and having standard deviation given by the estimated error on the corresponding Taylor coefficient.We refer the reader to Fig. 7, where the cloud of green points are the closest poles extracted for the L = 10 lattice.The scatter is from repeating the bootstrap procedure around ∼ 700 times, keeping the order of the Padé approximant fixed at [m, n] = [25,25].
• Systematic errors: These are the errors on the closest poles resulting from varying the order of the multi-point Padé and (or) changing the selection of the input points used to construct the Padé approximant, using only the central values of the input data.The idea is to change the input points by deleting data in a systematic way to construct the rational function of varying orders.We vary the order of the Padé approximant from [m, n] = [25,25] → [5,5].We now refer to Fig. 8 where we show the singularity structure for L = 15 lattice [60].Here the scatter of poles arising from changing the order of the rational function is shown as dark blue points, notice that the scatter is over only around ∼ 50 points, as the goal is only to show the stability of the closest pole.Note that the systematic errors mentioned above are correlated with the statistical errors.All in all, from the error analysis performed above, the stability of the closest poles of average magnetization gives us confidence in the fact that we are extracting genuine poles of the function, i.e. genuine zeros of the partition function and thus LY zeros.We can now proceed to analyse the scaling of these poles with the lattice volume to extract physical information like the critical exponents.

V. SCALING ANALYSIS OF ZEROS
Until now we have mainly focused on partition function zeros arising in the complex H plane (LY zeros) when considering cumulants at fixed temperature and varying H.However, looking at Eq. ( 2) we can also consider the partition function zeros in the complex inverse tempera-  3 for L = 10 , 15 plotted against the susceptibility data.Note that this is not an interpolation and no data on the susceptibility was used in the construction of the rational function.Since this is a derivative of an [m, n] = [25,25] rational function, the order is [m, n] = [24,25] ture (β) plane.This has been done numerically in [22], where the authors have studied the Fisher and Yang-Lee zeros of the 2D and 3D Ising model by using a relatively high number of cumulants in the temperature and external magnetic field variables.As explained before, instead of using such high order of cumulants, we have made use of the multi-point Padé method to study only two different cumulants as a function of temperature and external magnetic field.We refer again to the Hamiltonian of Eq. ( 1), in which we set J to unity.To draw parallel with the analysis in [22], we expand the partition function in terms of its zeros in the β plane, c being some constant and the product is over the k zeros given by {β k }.Thermal cumulants are defined by the relation which using the expansion above can be re-expressed as, Looking at Eq. ( 9) above, it is easy to see that near criticality, the closest zero to the real axis will contribute the most to the thermal cumulant.Additionally, it is possible to study the finite volume scaling of the Fisher zero following [11,61,62], and the relations describing the approach of leading zeros to critical inverse temperature can be written as and where β 0 is the Fisher zero, resulting in the closest singularity of cumulants to the real axis [63], β c is the critical inverse temperature and ν is the relevant critical exponent, which describes the divergence of the correlation length with respect to temperature, near criticality.The proportionality constants in Eq. ( 10

Im[H]
Poles, L=30 FIG. 6. Genuine poles extracted from Fig. 5 along with their error bars, shown for two lattice sizes (Top : L = 15, Bottom : L = 30).We want to highlight that these poles follow the circle theorem.They lie on the imaginary H axis (within errors) and the closest pole moves closer to the real H axis as the lattice size increases.Moreover, the number of poles increases.
for the energy density [22].

A. Extracting ν and βc
In order to determine these critical quantities, we will follow the steps (which we have previously also described in [26]): (1) we compute the n = 2 thermal cumulant (i.e. the specific heat) at various inverse temperatures β and lattice sizes L; (2) for each L we compute the rational approximant R m n (β) by our multi-point Padé method; (3) at each L we find the Fisher zero β 0 , which is obtained as the the closest singularity of the cumulant to the real The orders of the rational approximation used to obtain the blue cloud has been varied between [5,5] to [25,25].
axis; (4) we study the finite size scaling of the values of β 0 .We refer the reader to the first row of Table I for the results for ν.
• Using Eq. ( 10), we will try to extract the critical Once ν has been extracted from the data, one can fit the value of the critical inverse temperature βc given by the intercept A marked by a green star, which is reconstructed to 1% accuracy.The red star marks the intercept of Imβ0 as L → ∞ which is recovered to be zero.
exponent ν using the following fit with A being the logarithm of the proportionality constant in Eq. ( 10) and B is the exponent of 1 L that we want to extract and compare its value with ν.As can be seen in Table .I, the value of the relevant critical exponent ν is obtained with decent accuracy with a value of 1.014 (60), its exact value being ν = 1 [35].Shown in the top panel of Fig. 9 is a pictorial description of the fit described in Eq. (12).We plot Im(β 0 ) as a function of 1/L.On account of B ∼ 1, the dash-dotted line (which is the result of the fit) can hardly be distinguished from a straight line.
• Using Eqs.(10) and (11) and the fact that the exponent ν = 1, it is not hard to see that one can obtain a linear relation for Re[β 0 ] as a function of 1/L and define a fit function to extract β c as follows Fit IIb : Here ν = 1 simplifies Eq. ( 11) with Re[β 0 ] → β c in the limit L → ∞.Hence, in order to determine β c , we will fit the real part of the closest Fisher zeros to the real β axis as a function of 1/L and extract the intercept A. This intercept is shown as a green star in the bottom panel of Fig. 9. Also for this our estimate seems fairly accurate at β c = 0.4404 (19), when compared with the exact result of β c ∼ 0.4407 [35].We additionally show in the same figure, that after identifying the exponent, one can also find the intercept of the line Im(β 0 ) vs 1/L and show that it goes to zero within errors as seen with the red star on the figure.

B. Extracting βδ
After obtaining the inverse critical temperature, we can now perform simulations at β c , to study the closest zero Im(H 0 ) to the real axis in the complex H plane as a function of lattice volume L. This has been the focus of most of the previous discussions in Sections II & IV.Once again following the procedure outlined in [26], our program again entails the following steps: (1) we compute the n = 1 magnetic cumulant (i.e. the magnetisation) at β = β c and various values of external magnetic field H and lattice size L; (2) for each L we compute the rational approximant R m n (H) for the magnetisation by our multi-point Padé method; (3) at each L we find the Lee Yang zero H 0 , which is the singularity of the rational approximant for the magnetisation which is the closest to the real axis; (4) we study the finite size scaling of the values of Im(H 0 ) (as we have seen in Fig. 6 , H 0 always sits at Re(H 0 ) = 0), given by [11,22]: where, the exponent β is the well known critical exponent that describes how the average magnetization goes to zero when we approach the critical point from below T c and d is the dimension.in Eq. ( 15) is related to the infinite volume scaling function for the total magnetization.In order to extract the exponent in Eq. ( 15), we will fit the following function with A being the logarithm of the proportionality constant in Eq. ( 15) and B the exponent of L that we want to extract and compare with β/ν − d.Using the known scaling relations between the standard critical exponents, we can derive the following hyperscaling relation between β, ν, d and δ, Remembering the value of ν = 1 for the 2D Ising model, we can thus use the fit result for the parameter B to estimate βδ.As can be seen from Table I, the fit parameter B = −1.881(93)has been determined, which gives to decent precision the estimate for βδ whose exact value for the 2D Ising model is 1.875.Further, without using the hyperscaling relation given in Eq. ( 17), the fit parameter B should be compared against β/ν − d, to obtain β = 0.119(93) which compared against its exact value of 0.125 for the 2D Ising model, gives an estimate for the exponent.Finally, we show the results of the fit of H 0 we obtained for each lattice size, plotted against L β/ν−d in Fig. 10.In principle one should be able to follow these steps to estimate the critical region for QCD using Taylor expansions from lattice QCD.The relevant parameters which control the critical region in QCD will be the baryon density.We have tried to make some concrete steps in this direction recently in [64].However, being a much more complicated theory, both numerically and conceptually, we may have to wait for some time to be able to do that.(12,13,16), shown for each row respectively.For Fit I, B has to be compared with the exact value stated, whereas for Fit II, the intercept gives βc, hence the exact value has to be compared with A.
For the last fit, Fit III, the fit parameter B has to be compared to the exact value of the critical exponent product, namely, βδ.

VI. CONCLUSIONS AND OUTLOOK
As a first step we simulated the 2D Ising model using a cluster spin flip algorithm in two ways.For the LY zero analysis, we simulated the model on varying lattice sizes at a set of values of the external magnetic field.These simulations were performed at T c , and the goal was to approximate the average magnetization as a rational function of the external magnetic field and study the structure of zeros and poles that arise.On the one hand we were able to verify numerically, many properties of the LY zeros including the famous circle theorem for the Ising model, observing that only the genuine poles (uncancelled and stable) of magnetization lie on the purely imaginary H axis.It was further observed that the number of genuine poles poles increases with volume and for simulations at T c , comes closer to the real H axis.In order to verify that these were indeed physical effects, volume scaling of the zeros, using the prescription in [11,22] were performed leading to a decent estimate of the combination of critical exponents βδ = 1.881(93).The fact that LY zeros can be studied at T c is in our approach fully self-consistent.In fact, Fisher zeros were also studied by approximating the specific heat with a multi-point Padé function and studying its poles in the complex β plane.Finite size scaling of these zeros following the prescription of [22] was done to obtain again, precise values of the critical exponent ν = 1.014(60) and of the critical inverse temperature β c = 0.4404 (19).These results give us some confidence in our pursuits of studying the QCD phase diagram using lattice QCD simulations combined with multi-point Padé method.The main caveat being that in the case of the Ising model it was not very computationally expensive to reach statistics of the order of ∼ 625K configurations for each lattice, for each value of H and β.These kind of statistics are not currently realistic for lattice QCD simulations.
Another direction to pursue would be to study the universal location of the Lee-Yang edge singularities as done using the Functional Renormalisation Group (FRG) approach in [65,66] or using a suitable parametrization of lattice data as done in [67] for the O(N ) models.In order to do this using our approach, we would need accurate values of the closest LY zero to the real H axis, but at temperatures T > T c .We would like to close by inviting the reader to use the multi-point Padé method for extracting Lee-Yang and Fisher zeros in their choice of models.All data from our calculations, presented in the figures of this paper, can be found in [68].

FIG. 1 .
FIG.1.Top : Average energy (scaled by 1/L instead of 1/L 2 for better visualization) as a function of β from simulations.Bottom : Specific heat capacity per lattice site calculated from the configurations generated.

FIG. 2 .
FIG.2.Top : Average magnetization per site as a function of H from simulations.Bottom : Susceptibility per site calculated from the configurations generated.

FIG. 4 .
FIG.4.Derivative of the rational function obtained in Fig.3for L = 10 , 15 plotted against the susceptibility data.Note that this is not an interpolation and no data on the susceptibility was used in the construction of the rational function.Since this is a derivative of an [m, n] =[25,25] rational function, the order is [m, n] =[24,25]

FIG. 5 .
FIG.5.Zeros of the numerator (black pentagons) and of the denominator (red crosses) of the rational approximant R m n (H) for the magnetisation on L = 15 (Top) and L = 30 (Bottom), with [m, n] =[25,25].The pale blue circles are the points used as input for the Padé.Notice that the closest singularity to the real axis gets closer to the real H axis as L gets larger, with real parts being consistent with Re(H0) = 0.

FIG. 7 .FIG. 8 .
FIG.7.Stability of closest pole for L = 10.The green cloud represents closest poles extracted from varying the Taylor coefficients with noise drawn from a Gaussian distribution.The singularity structure shown, red crosses and black pentagons, are for the central values of the data at L = 10.This is the result of drawing the coefficients ∼ 700 times using a rational function of order [m, n] =[25,25]

Fit
FIG. 9. (Top)The scaling in 1/L of Im(β0), i.e. the imaginary part of the Fisher zero, detected as the closest singularity of the cumulant to the real axis.The correct critical exponent ν = 1 is reproduced with fairly good accuracy.(Bottom) Once ν has been extracted from the data, one can fit the value of the critical inverse temperature βc given by the intercept A marked by a green star, which is reconstructed to 1% accuracy.The red star marks the intercept of Imβ0 as L → ∞ which is recovered to be zero.

TABLE I .
Results for the fits shown in Eqs.