Critical and tricritical singularities from small-scale Monte Carlo simulations: the Blume–Capel model in two dimensions

. We show that accurate insights into the critical properties of the Blume–Capel model at two dimensions can be deduced from Monte Carlo simulations, even for small system sizes, when one analyses the behaviour of the zeros of the partition function. The phase diagram of the model displays a line


Introduction
Global warming is nowadays an essential concern for everyone, and the indigence of the most powerful governments which do not take measures to meet the challenges leads more and more people, the youngest in particular, to seize the environmental issues.The message of science is unequivocal: there is no collective solution without a considerable reduction in greenhouse gas emissions.Academic circles play their part in this awareness, some in a very active manner, such as Scientists Rebellion.At universities, researchers wonder about the compatibility between their professional practices and reducing the carbon footprint.In the field of intensive computing, supercomputing centres enormously generate greenhouse gases and perhaps the use of (massively parallel) large-scale simulations and high-performance computing must be questioned.We consider here an approach which goes into the opposite direction of small-scale simulations in the field of critical phenomena and demonstrate that it is possible to obtain relatively precise results at a significantly lower computational cost.
In particular, we focus on the determination of critical exponents and other universal amplitudes of a well-known lattice model in statistical physics which features tricritical behaviour.
Thermodynamic quantities in the vicinity of a second-order phase transition are indeed characterized by power laws in the reduced temperature t and external field h.In the thermodynamic limit, at h = 0, the divergence of the correlation length is responsible for the critical singularities, following from ξ ∞ (t) ∼ |t| −ν .Additionally, one has C ∞ (t) ∼ |t| −α for the specific heat, χ ∞ (t) ∼ |t| −γ for the susceptibility, and m ∞ (t) ∼ |t| β for the order parameter or magnetisation of the system at T < T c .The critical behaviour may be identical for distinct systems, therefore the set of exponents determine a universality class.These universality classes depend only on very general properties, like the dimensionality of the system, the range of interactions, and the symmetries of the ground state of the system.
Our aim here is to show that the analysis of the scaling of the zeros of the partition function, popularized by the seminal work of Lee and Yang [1,2], is a very accurate tool for extracting the critical behaviour of spin models.This is of course not new, but we move here in the opposite direction than the usual trend.Ordinarily, critical exponents are extracted from extensive Monte Carlo simulations via finite-size scaling (FSS) methods of thermodynamic quantities like the susceptibility or the order parameter.This requires simulations at different sizes for a range of temperatures to determine first the pseudocritical parameters where the diverging quantities (e.g., the susceptibility) display their maximal values, then to analyze the size-dependent power-law behaviours of these thermodynamic observables at pseudocriticality.The larger the sizes, the better the accuracy of the critical exponents measured, especially in cases where corrections to scaling are strong and call for particular attention [3].
The zeros of the partition function on the other hand can provide an alternative to these commonly used thermodynamic quantities.In fact, their behaviour upon approaching the critical point also encodes the critical exponents needed to characterize a universality class.We show here, by comparing the two approaches outlined above, that, indeed, the FSS analysis of the partition function zeros does not require as large system sizes as those needed in the case where the usual thermodynamic quantities are involved in order to achieve a given accuracy for the estimation of the basic critical exponents.In this context, our work designates the following unexpected finding: An excellent precision can be achieved via scrutinizing the partition function zeros of very small system sizes, of the order of L ≤ 8, where L denotes the linear dimension of the lattice; thus, as few as 64 spins or less for a two-dimensional (2D) lattice considered here.To support our claim we choose as a platform the two-dimensional Blume-Capel model [4,5], one of the most well-studied fruit-fly models of statistical physics, which has a rich phase diagram featuring first-and second-order phase transition lines that merge at a tricritical point.Of course, the numerical results presented here do not consist a rigorous proof of our main conclusion, but provide a clear and strong empirical support to it, an aspect that has not been recorded previously in the relevant literature.As a side remark, the analysis of the zeros that we provide is to the best of our knowledge the most accurate ever performed for the Blume-Capel model.
The rest of the paper is laid out as follows: In section 2 we introduce the model and in section 3 we shortly outline the implemented numerical methods.Our numerical results and FSS analysis are presented in sections 4, 5, and 6.Finally, a summary of our conclusions is given in section 7.

Blume-Capel model
The Blume-Capel model is a particular case of the Blume-Emery-Griffiths model [4,5] which has greatly contributed to the study of tricriticality in condensed-matter physics.It is also a model of great experimental interest since it can describe many physical systems, including, among others, ultracold quantum gases [6], 3 He-4 He mixtures, and multi-component fluids [7,8].In fact, experimental studies have focused on systems with a tricritical point, such as FeCl 2 for example [9].The spin-1 Blume-Capel model is introduced through the Hamiltonian where J > 0 denotes the ferromagnetic exchange interaction coupling, the spin variables σ i ∈ {−1, 0, +1} live on a square lattice with periodic boundaries, and ⟨i, j⟩ indicates summation over nearest neighbors.∆ represents the crystal-field strength that controls the density of vacancies (σ i = 0).For ∆ = −∞ vacancies are suppressed and the Hamiltonian reduces to that of the Ising ferromagnet.
The Blume-Capel model has been studied extensively, mainly using numerical approaches [10][11][12][13][14][15]; for a recent review, we refer the reader to reference [16].showing the ferromagnetic (F) and paramagnetic (P) phases that are separated by a continuous transition at small ∆ (solid line) and a first-order at large ∆ (dotted line).The line segments meet at a tricritical point (∆ t , T t ) marked by a black rhombus.Numerical data shown are selected estimates from previous studies, including the critical point at zero crystal field which we denote hereafter for the needs of our work as T c (∆ = 0), or simply as T c .The numerical simulations reported in this paper were carried out at these two special points (∆, T ) = {(0, T c ), (∆ t , T t )}.
A reproduction of the phase diagram of the square-lattice Blume-Capel model in the (∆, T )-plane is shown in figure 1: For small ∆ there is a line of continuous transitions between the ferromagnetic and paramagnetic phases that crosses the vertical axis at (∆ c = 0, T c ≈ 1.6929) [12].For large ∆, on the other hand, the transition becomes discontinuous and it meets the T = 0 line at ∆ 0 = zJ/2 [4], where z = 4 is the coordination number (here we set J = 1 and k B = 1 to fix the energy and temperature scale).The two line segments meet in a tricritical point estimated to be at (∆ t ≈ 1.966, T t ≈ 0.608) [13,17].It is well established that the second-order transitions belong to the universality class of the two-dimensional Ising model [16,18].We should note here that similar phase diagrams showcasing a tricritical point can also be found in the Baxter-Wu model in transverse magnetic field [19], the spin-1 Baxter-Wu model [20][21][22] and the bond-diluted 4-state Potts model at three dimensions [23,24].Note that the two-dimensional disordered 4-state Potts model always displays a continuous phase transition [25,26].

Hybrid and Wang-Landau simulations
Most conventional Monte Carlo schemes generate a canonical distribution P (E, β) = g(E)e −βE at a fixed temperature β (to be precise β refers to the inverse temperature and should not be confused with the critical exponent of the order parameter).Then, typically as many simulations as necessary are performed in order to probe the thermodynamic observables under study within a sufficient temperature range.For the particular case of the Blume-Capel model and due to the presence of the zero-state spins, we employ here a combination of a Wolff single-cluster update [27] of the ±1 spins and a single-spin-flip Metropolis update [28][29][30] to account for the vacancies, denoted herewith as the hybrid Metropolis-Wolff (MW) approach.
We also performed an additional array of canonical Monte Carlo simulations implementing the Wang-Landau (WL) entropic sampling method [31][32][33], an efficient protocol that samples the density of states g(E).As g(E) does not depend on β (or ∆), it is possible to construct the canonical distribution at any temperature and then compute the partition function.Using the Wang-Landau method the energy and magnetisation (in zero magnetic field) are computed as follows where g(E J , E ∆ ) and g(E J , E ∆ , M ) are the appropriate densities of states and the absolute value of M breaks the Z 2 symmetry as usual (otherwise the simulations would deliver a vanishing magnetisation in zero magnetic field).We should point out at this stage that the main drawback of the Wang-Landau algorithm lies in its heavy computational time cost which limits simulations to rather moderate system sizes.Of course, several improvements of the algorithm have been presented over the last years facilitating the simulations of much larger systems [34][35][36][37][38], some of which however with unclear consequences in the numerical accuracy [39].In lieu, the hybrid approach [28][29][30] allows the simulation of considerably larger system sizes in moderate computational times but lacks direct access to the density of states.However, this problem can be surpassed with the help of histogram reweighting [40].This is a well-known technique that allows the safe extrapolation of data emerging from a Monte Carlo simulation at a fixed temperature to a nearby temperature window (within a certain level of accuracy).
For both methods (hybrid and Wang-Landau), the specific heat (C), susceptibility (χ), magnetocaloric (m T ), and "mixed" (χ 12 and χ 2 ) coefficients (see also below) are accessible through the usual fluctuation-dissipation relations and Note that in the above formulas, the notation E refers to the total energy, whereas E ∆ when only the contribution of the crystal-field term is taken into account [equations (7) and ( 8)].
As already mentioned above, we carried out several series of simulations for the square-lattice Blume-Capel model with periodic boundary conditions.In particular, for the hybrid method we run simulations at the zero crystal-field critical point [12] and the tricritical point [13], located along the phase boundary of figure 1 at respectively, establishing the targeted simulation points of our numerical endeavour.
Our results were then extrapolated to nearby temperature values benefiting from the histogram reweighting technique.In our protocol the hybrid algorithm performs 1500 × N Monte Carlo Steps per spin (MCSs) to reach a steady state and 9 × 1500 × N MCSs to collect the data; N = L × L denotes the total number of spins on the lattice.A MCS consists of 3N Metropolis steps followed by L Wolff steps, a practice suggested as optimum in references [16,41].For the Wang-Landau method now, the number of iterations increased until the modification factor f , which controls the convergence of the simulation, was less than f final = 1 + 10 −8 [32].Finally, errors were determined using the Jackknife method [42] for the hybrid algorithm, while 10 simulations with different initial configurations were carried out to compute errors of the relevant thermodynamic quantities accumulated during the Wang-Landau simulations.

Zeros of the partition function
The partition function of the Blume-Capel model on a finite D-dimensional lattice may be given in terms of the density of states g(E J , E ∆ , M ) The microstates are labelled by the energies E J (integer values {−DN, DN }), the crystal-field energies E ∆ (positive integer values {0, N }), and the magnetisation values (integer values {−N, N }).The sum is explicitly written as where h = βH and δ = β∆ -we stress that there should be no confusion with the critical exponent δ here, as the inverse temperature should not be confused with the critical exponent β.
In order to find Fisher zeros -which are the zeros of the partition function in the complex plane of the (inverse) temperature variable β = β r +iβ i -the partition function in zero magnetic field takes the form where g(E J , E ∆ ) is the associated density of states.One defines a normalised version of equation ( 12) by Z(β r , β i ) = Z(β)/Z[Re(β)], which selects only the sine and cosine parts The zeros of the partition function, labelled by the index n, are the values β rn +iβ in of the complex temperature where both real and imaginary parts of Z(β r , β i ) are equal to zero.In order to find the first zero β r0 + iβ i0 , simulations were performed in a range of complex temperatures (β r , β i ) in the vicinity of the real values of the parameters at the transition (i.e., on the transition line) to get the intersection in the complex plane between two curves: Those corresponding to the independent vanishing of the real part ⟨cos[β i (JE J − ∆E ∆ )]⟩ βr = 0 on one hand and of the imaginary part ⟨sin[β i (JE J − ∆E ∆ )]⟩ βr = 0 on the other hand.When both conditions are met a zero is found.The closest to the transition line is the first zero.At the tricritical point -see equation ( 9) -and for sizes L = 4 − 8, the zeros were found using the Wang-Landau algorithm working exactly at ∆ t = 1.966.For the larger sizes considered within the range L = 16 − 64, the simulations were performed with the hybrid algorithm enhanced by the histogram reweighting technique to explore the phase space around the simulated temperature T t with the minimal effort.An analogous strategy was also followed at the zero crystal field critical point.
A consistency check regarding the validity of the method and the numerical accuracy between the different algorithms is provided in figure 2. This plot depicts the curves in the complex temperature plane at ∆ c = 0, where either the real or the imaginary parts of Z(β r , β i ) vanish.Comparative results from hybrid Metropolis-Wolff and Wang-Landau methods are shown for a system with linear size L = 4.It is evident that there is a perfect matching between the numerical data produced by the two algorithms.The first zeros, given by the intersections of the two types of curves, are located in the vicinity of the critical temperature (shown by the vertical solid line).In order to spot an accurate location of the zeros, we follow a two-step process: During the first step the approximate position of the zeros is graphically identified by plotting the sign changes of Re(Z) and Im(Z) in the vicinity of the critical (or tricritical) temperature.In the second step and once these positions are well approximated, a polynomial fit over a few points in the vicinity of the intersection using the standard Newton algorithm determines with high precision where the real and imaginary parts of the function are simultaneously equal to zero; see figure 3 for a graphical illustration of the described approach.A single hybrid Metropolis-Wolff simulation was performed at the critical temperature T c = 1.6929.The two sets of blue and red circles show respectively changes in the sign of the real and imaginary parts of the partition function.The solid lines are fourth-order degree polynomial fits and the zero sought is enclosed within the central losange.Right panel: A zoom in of the interesting area for the benefit of the reader.

Some comments on the goodness of fits
For the fitting procedure discussed below and in order to determine an acceptable fit we employed the standard χ 2 test for goodness of fit -this χ 2 parameter here should not be confused with the observable χ 2 !Specifically, the Q value of our χ 2 test is the probability of finding a χ 2 value which is even larger than the one actually found from our data [43].Recall that this probability is computed by assuming (1) Gaussian statistics and (2) the correctness of the fit's functional form.We consider a fit as being fair only if 10% < Q < 90%.Note that if Q < 0.001 the model should be well verified to find the cause but can still be acceptable [44].At the critical point, for systems of finite size L, the singular part of the free energy density is a generalized homogeneous function [45], expressed according to the scaling with y t and y h the renormalization-group eigenvalues associated to t and h.Firstorder derivatives with respect to the temperature t and magnetic field h define the internal energy e and magnetisation, respectively, while second-order derivatives provide the specific heat C, susceptibility χ, and magnetocaloric coefficient m T [cross derivative, m T = ∂ 2 f sing (t, h)/∂t∂h; see equation ( 6)].Homogeneous forms analogous to equation ( 14) also follow for these quantities, leading at the critical point t = 0, h = 0 to the following FSS behaviour ) Note the renormalization-group eigenvalues y IM t = 1 and y IM h = 15/8 for the case of the two-dimensional Ising model [46].These exponents govern the critical behaviour along the second-order transition line of the phase boundary and the associated critical exponents fall into the Ising universality class, where In the above equation only the last two observables diverge with the system size in the vicinity of the transition and are then adapted for a numerical study, in particular for the determination of the pseudocritical point, since otherwise a possible regular contribution (e.g., in the case of the energy and of the specific heat) may hide partially the singular signal.
Ising model universality class Along the second-order transition line the scaling field g = ∆ − ∆ c associated to the crystal field (here ∆ c = 0, but it would take a non-zero value along the transition line) is a supplementary parameter which can also be used in the numerical analysis as it allows to prescribe other useful response functions, such as see also equations (7) and (8).
As discussed in Section 2 the phase diagram of the Blume-Capel is well established.Our objective here is to study first the convergence of the exponents to the Ising universality class.The next subsection will be devoted to the tricritical Ising model universality class.Simulations were therefore carried out at T c = 1.6929 (∆ c = 0) [12] where histogram reweighting was performed in the parameters T or h.
Figure 4 displays the typical shape of some thermodynamic averages.In particular, the magnetocaloric coefficient and the fourth-order magnetisation cumulant for small system sizes (L = 5 − 8) were obtained by temperature reweighting, comparing the two algorithmic approaches used.In a finite system, thermodynamic quantities do not present any singularity; fluctuations are limited by the size of the system.Thus, the maxima of these quantities lie close to the critical temperature in the thermodynamic limit (vertical line exactly at T c = 1.6929).This temperature of the maxima, T L , can be identified as the pseudocritical temperature of the corresponding quantity.Pseudocriticality may also be characterised in other directions of the parameters space, e.g., the magnetic-field direction, as the value ±H L of the magnetic field for which diverging quantities are maximised.
Further simulations were then performed up to sizes L = 64, which is still something fairly easy to achieve.Then, we proceeded to the ordinary FSS analysis, fitting the  quantities to simple power laws in order to test equation ( 16).Since we do not take corrections-to-scaling into account, the fits are linear.However, it is clear that finitesize effects are present.In order to get a qualitative imprint of their presence while approaching the asymptotic limit, we performed for each quantity of interest fits within three intervals of system sizes: (i) The small lattice-size regime (L = 5 − 8), (ii) the complete available lattice-size regime (L = 5 − 64), and (iii) the large lattice-size regime (L = 16 − 64).
As expected, the fits taking only the largest sizes into account provided the best estimates for the critical exponents -the expected values are presented in table 1 for comparison.We note that fits with only the smallest sizes systematically gave results that are not very accurate and we then proceeded to consider possible improvements with the inclusion of corrections-to-scaling terms.It should be noted however, that these effects are counterbalanced by the large sizes when all the sizes of the system are taken into account.The best results are those stemming from the magnetocaloric coefficient and susceptibility.The study by Zierenbeg et al. [16] reported a value of γ/ν = 1.750(3) for sizes within the range L = 32 − L = 128, while Malakis et al. [47], using Wang-Landau simulations and sizes L = 50 − 100 reported γ/ν = 1.748 (11) at the crystal field value ∆ = 1 (still in the second-order regime of the phase boundary where the model is governed by the Ising universality class).These results suggest that in our case scaling-correction effects appear for such sizes, even for the case L ≥ 16.One can thus try to improve the fits by considering correction-to-scaling terms via the exponent ω, which is known to take the value ω = 1.75 for the two-dimensional Ising universality class [48,49].However, this strategy appeared to only slightly improve the extrapolations, when sizes up to 64 were taken into account, but failed at improving small sizes estimates; see table 2 for more details.

Scaling hypothesis and finite-size scaling at the tricritical point
In the case of the Blume-Capel model, in order to describe the approach to the tricritical point an additional field enters the description and it is more appropriate to introduce the scaling fields t and g as [8] with the temperature T , crystal field ∆, magnetic field H, and a coefficient a.The third field h = H describes the deviations from the tricritical point outside the symmetry plane, while g and t describe the approach within the symmetry plane.The free-energy density is then expressed as It is more convenient in our case to work with the variable ∆ rather than g since it appears directly in the Hamiltonian.At the tricritical point, critical exponents have been first conjectured from the dilute Potts model [50][51][52] and determined by conformal field theory [46,53,54] as y tri t = 9/5, y tri h = 77/40, y tri g = 4/5 [55].The expected FSS forms in this case read as Since the crystal field ∆ is a linear combination of the scaling fields t and g, derivatives with respect to ∆ will lead to the following corrections-to-scaling for the two responses which involve the crystal-field fluctuations.
The same procedure as the one carried out at the critical point is also presented below, this time at the tricritical point; see also figure 1.The FSS analysis of the thermodynamic quantities under study leads to the exponents disclosed in tables 3 and 4. In table 3 the analysis is based on the maxima of the diverging quantities when the histogram reweighting is performed along the temperature direction.Therefore the quantities are measured strictly at ∆ t and at the corresponding pseudocritical temperatures T L .In table 4 we follow the exact reversed strategy, sitting at the corresponding pseudocritical values of the crystal field ∆ L .

Fisher, crystal-field, and Lee-Yang zeros
An alternative perspective to the study of critical phenomena is through the zeros of the partition function.These are fundamental quantities of special interest since their FSS properties encode universal features of underlying phase transitions [56,57].They appear when a certain parameter of the partition function (temperature, magnetic, or crystal field, or associated fugacities) attains complex values.The approach was first developed by Lee and Yang [1,2] who studied the lattice gas model via the grand partition function as a polynomial with a complex fugacity.Fisher then tackled the problem via the zeros in the complex temperature field, where, with increasing system size, the position of the zeros would get closer and closer to the real axis and the critical temperature.This line of research was proven to be quite successful for the study of criticality in the Ising model leading to high-accuracy results even on the basis of smallscale studies and has also been successfully applied to other spin models as well, such as the Potts model [58,59] and the Blume-Capel model considered here.
Fisher zeros are the zeros in the complex temperature plane, while Lee-Yang zeros are those in the complex magnetic field plane.According to equation ( 13), the partition function can be written for a complex temperature β = β r + iβ i , as   9).Note the inverse critical temperature β c = 0.5907 and that H c = 0.
In the complex plane of the crystal field ∆ = ∆ r +∆ i , the crystal-field partition function zeros will follow from the vanishing of Eventually, in the complex h-plane corresponding to the external magnetic field H, we obtain the Lee-Yang partition function It is important to note that for the present model under study the Lee-Yang theorem holds, i.e., all Lee-Yang zeros are purely imaginary, H rn = 0, and lie on a unit circle in the complex fugacity plane e −2h .This theorem holds for all dimensionalities and/or lattice sizes.Figure 5 illustrates the temperature (Fisher), crystal-field, and magnetic-field (Lee-Yang) zeros for the Blume-Capel model on a 6 × 6 square lattice in zero magnetic field at the critical point specified in equation ( 9).The zeros are located in the area where the real part (shown in blue) and the imaginary part (shown in red) of the partition function are both zero for the same complex parameter (T , ∆, or H), so that the partition function vanishes and therefore provides the position of a zero.The figure shows only the changes in the sign of Re(Z) and Im(Z) from which we can deduce the value of the zeros of the partition function.As it is well-known, the FSS of the Fisher and Lee-Yang zeros gives access to the thermal exponent, y IM t , and the magnetic exponent, y IM h , respectively Following our previous practice in the framework of FSS analysis, we carried out power-law fits over different size intervals.At this point we should note that we were  only interested in the imaginary part of these zeros; of course the real part can also be used, but as the results were less interesting we excluded them for the needs of this study.
For each type of zero we focused on the first zero β i0 or H i0 , i.e., the zero closest to the real axis.In particular, one of the best estimates that we retrieved refers to the value y IM h = 1.876 (2), which is in excellent agreement with the expected value y IM h = 1.875.Inspecting table 5 that summarises the results from all our fitting attempts, it becomes evident that both Fisher and Lee-Yang zeros give critical exponents very close to the expected values for each case.What is more, the estimated exponents are remarkably accurate even when considering the FSS analysis over the small lattice-size regime, reflecting the fact that within this approach, finite-size corrections are negligible.This is in sharp contrast to the case of other thermodynamic averages considered before.Note that we have not shown here any fitting analysis based on the zeros of the crystal field at the critical point (only at the tricritical point; see below) as the results for the fits were found to be of low quality.
Next, we contemplate the case of the tricritical point.Graphically the zeros are shown in figure 6.Compared with the critical point at the same size, the Fisher zeros are much less dense and the first zero has a much larger imaginary component.On the other hand, the crystal field and Lee-Yang zeros are very dense, even at small sizes.These first zeros almost form a vertical line.In fact, the crystal-field zeros have a similar behaviour to that of the Fisher zeros for the critical point, and vice versa.

Tricritical Ising model universality class
Zeros In table 6, we summarise the results from the FSS analysis of the imaginary parts of the zeros in the complex plane corresponding to the temperature, crystal and magnetic fields, according to We only focused on the crystal-field and Lee-Yang zeros since corrections were found to be very strong for the Fisher zeros and it was not possible to obtain a good exponent from these, unless a large proportion of the system sizes was discarded (as an example we quote the value y tri t = 1.72(1) for fits within the range L = 16 − 64 of system sizes).The opposite is not true for the crystal-field and Lee-Yang zeros.In particular, at small sizes the exponents found were already those expected.As at the critical point, these zeros indicate that finite-size corrections are negligible.Nevertheless, we also did test powerlaw fits with corrections which improved the exponent too.We would like to remind the reader here that the zeros of the two-dimensional Blume-Capel model have already been studied in the complex temperature plane.Ayat and Care [60] tackled the problem along the fugacity complex plane u = e −β∆ and the temperature complex plane t = e −β and reported for a system with linear size L = 22 the values y tri t = 1.814 and y tri t = 1.782, respectively.Kim and Kwak [61] also studied the zeros of the partition function in the complex temperature plane using Wang-Landau simulations and suggested the tricritical exponent y tri t = 1.83.Although the conventional approach to probe efficiently the critical behaviour of spin models using numerical simulations of finite sizes L dictates the need for L → ∞, the results presented in this section advocate a different scenario: the critical properties of the Blume-Capel model may be accurately established from small-scale numerical simulations after employing the FSS analysis of the partition function zeros.Other studies focusing on the zeros of the partition function for small system sizes have shown a similar trend, i.e., minor finite-size corrections.Using the cumulant method, Deger et al. [62,63] determined the critical exponents of the Ising model in two and three dimensions with good accuracy via the Lee-Yang and Fisher zeros using sizes up to L = 15 in two dimensions and L = 10 in three dimensions, respectively.The same authors [64] also studied the phase transition in a molecular zipper with relatively small lengths up to N = 30, establishing the expected exponents.The zeros were also suggested as an alternative route to probe logarithmic corrections in statistical physics.In fact, Kenna and Lang [65] have been able to numerically extract the logarithmic corrections, which were analytically predicted by the renormalization-group theory for the Φ 4 model in four dimensions, using systems of linear sizes L ≤ 24.Additionally, through the density of zeros one may gain important information and determine the order and strength of a transition, as Janke and Kenna [66] have shown for several systems, again using relatively small system sizes.Interestingly, the study by Alves et al. [67] contrasted the Janke and Kenna approach [66] with another classification scheme proposed by Borrmann et al. [68] which explores the linear behavior for the limiting density of zeros proposed, and suggested that it is indeed possible to characterize the order of the transition in the four-and five-states Potts model via a small-size scaling analysis.On the other hand, there are no real explanations or hypotheses provided in the vast literature around this topic.One suggestion that we bring forward here, in support of the above observations that the finite-size corrections to the zeros of the partition function are small, is that the zeros may not comprise of regular contributions, or at least are strongly dominated by the singular part.

Impact angle and amplitude ratios
In the complex plane of the temperature, the zeros in the vicinity of the critical point β c follow a line which should pinch the real axis at the critical temperature T c .The real axis of β and this line form an impact angle φ which is related to the exponent α and the amplitude ratio of the specific heat via the relation if one defines the impact angle φ between the positive sense of the axis and the zeros.The critical exponent α and the amplitude ratio A + /A − are universal quantities, thus the impact angle is itself universal.For the two-dimensional Ising universality class the estimate of the amplitude ratio A + /A − = 1 from Monte Carlo simulations [69] suggests that φ 2D = π/2.The values are the same also for the ϕ 6 model [70].The impact angle is usually attained in the complex plane of the temperature, and we show here that it is also possible to extract this angle, and therefore the amplitude ratio A + /A − , in the crystal-field plane at the tricritical point.Thus, at the tricritical point the study was undertaken at the complex crystal-field plane and the exponential plane u = e −β∆ .There are several ways to retrieve the impact angle.We elaborate here on three different ways for each system (i) the angle between the first zero and the real crystal field, (ii) the angle between the second zero and the real critical crystal field, and (iii) the angle between the first and second zero with respect to the real axis.Table 7 recaps the results for system sizes up to L = 64.
The angle φ j,k is the angle between the real ∆-axis and the line from ∆ (j) (L) and ∆ (k) (L).Undoubtedly, for sizes in the range L = 32−64 all the angles are asymptotically close to the expected value φ 2D = 90°, allowing us to safely extract an estimate for the average angle φ = 89.68(9), taken as an average over the u-plane estimates.From that one may extract the value A + /A − = 1.07 (6) for the amplitude ratio using equation (33) and the exact value of α.To the best of our knowledge this is the first numerical estimation of the tricritical amplitude ratio in the two-dimensional Blume-Capel model.

Partition function zeros: Energy and magnetisation cumulants
Another powerful method to extract the partition function zeros from fluctuations of thermodynamic observables, such as the energy and the magnetisation, has been developed by Deger et al. [63,64,71,72].The advantage of the so-called cumulant method is that it does not require the computation of the full density of states, only usual simulations at fixed external parameters.The cumulant method has already been successfully tested to extract the Fisher and Lee-Yang zeros of several models, including a molecular zipper [64], the Ising model at two and three dimensions [62,63], and the Curie-Weiss model [72].The goal in this section is to test the cumulant method on the basis of the Lee-Yang zeros at the critical and tricritical points of the two-dimensional Blume-Capel model and compare the outcomes with our previous results.

From cumulants to zeros
The partition function Z(q) contains fluctuations of thermodynamic quantities which are determined by the derivatives of the free energy with respect to a control parameter q of a thermodynamic observable Φ ⟨⟨Φ n (q)⟩⟩ = (−1) n ∂ n q ln Z(q)/N.
For instance, for the total energy one has Φ ≡ E and if the control parameter is the inverse temperature, then and the energy cumulants are generally expressed as non-linear combinations of central moments The first cumulant gives the mean while the second cumulant is the variance, related respectively to the average energy and specific heat.One then has, up to n = 4 where ⟨⟨E n ⟩⟩ ≡ κ (n) .Note that one can find expressions of the cumulants up to the tenth order in reference [73].The relation between the cumulants and the partition function zeros is given by where q k are the zeros of the thermodynamic observable in the complex plane of the control parameter q.High-order cumulants encapsulate precious information about zeros, especially the first zero, i.e, the one closest to the real axis, since it dominates the sum.The contributions from subleading zeros are suppressed with the cumulant order n [64].
The phase transition of the Blume-Capel model depends on three control parameters: (i) the inverse temperature q = −β, (ii) the reduced external magnetic field q = βH, and the reduced crystal field q = β∆.Thus: (1) For the case of the energy Φ ≡ E and its conjugate variuable q = −β one obtains the Fisher zeros.
(2) For the case of the magnetisation Φ ≡ M , and its conjugate variable q = βH one may extract the Lee-Yang zeros.

The case of Lee-Yang zeros
For the case H = 0, odd cumulants ⟨⟨M 2n+1 ⟩⟩ vanish due to the symmetry properties.Hence, equation (37) simplifies to Monte Carlo simulations for systems with linear sizes L = 16 − 64 were performed at the critical (∆ c = 0, T c = 1.6929) and tricritical (∆ t = 1.966,T t = 0.608) points for n = 3. Figure 7 compares the results from the cumulant method and those from the partition function analysis secured in the previous section.The values of the zeros obtained by the two methods are very close, as one can inspect from figure 7 where the difference between the two fits is barely perceptible.At the critical point the estimate y h = 1.8802(1) is computed in agreement with the expected value y IM h = 1.875, while at the tricritical point the exponent value found y h = 1.9240 (61) is in perfect agreement with the expected result y tri h = 1.925.Therefore, we may safely conclude that the cumulant method is a very efficient way to extract the exponent y h , both at the critical and tricritical points, involving at best the computation of the sixth-and eight-order cumulant of the magnetisation.
Another interesting aspect worthy of investigation relates to the crossover phenomena at the tricritical point and the sensitivity of the Lee-Yang zeros on the transition.Let us now focus only on the tricritical point and explore its critical behaviour by arbitrarily selecting two additional values of the temperature in close proximity to the tricritical value T t = 0.608.Namely one temperature below T t , T = 0.6075, and another one above T t , T = 0.6083.The crystal field is fixed at its tricritical value ∆ t = 1.966 [13].Figure 8 illustrates the fits for each temperature.Three distinct behaviours can be identified and several comments are in order at this point: (i) Below the tricritical temperature the exponents found quickly converge to the first-order characteristic exponent y h = D. We therefore present evidence of the first-order phase transition.(ii) Above the tricritical temperature we see that the exponent found is remarkably close to that of the Ising ferromagnet, falling back into the critical behaviour of the second-order transition line and the Ising universality class.(iii) These results highlight that the Lee-Yang zeros at the tricritical point are extremely sensitive quantities and it is therefore extremely important to know precisely the location of the critical line and tricritical point of the model in order to extract securely the critical exponents.(iv) Finally, scanning the phase diagram from T = 0.6075 to T = 0.6083, it is possible to characterize all the richness of the different transitions of the two-dimensional Blume-Capel model and to extract the exponent y h in a remarkably accurate manner.

Finite-size scaling of the cumulants
Once we have computed the cumulants, an ordinary FSS analysis is also welcome, and we expect the following behaviours ⟨⟨M n ⟩⟩ ∝ L −D+ny h .(43) Figure 9 shows as an example the power-law behaviour of the magnetisation cumulants at the critical point.We only focus on the even cumulants since the odd ones vanish.As one can see, the exponents are really close to the expected ones, both at T c and at T t for the large lattice-size regime L = 16 − 64.At the critical point, we extract ⟨⟨M 2 ⟩⟩ ∼ L 1.7581 , ⟨⟨M 4 ⟩⟩ ∼ L 5.5204 (3) , ⟨⟨M 6 ⟩⟩ ∼ L 9.2808 (4) , and ⟨⟨M 8 ⟩⟩ ∼ L 13.0411 (5) where one expects the values 1.75, 5.50, 9.25, and 13.00, respectively.The exponent estimates are already excellent for the small lattice-size regime, except for the case n = 2. Indeed, one finds ⟨⟨M 2 ⟩⟩ ∼ L 1.5922 , ⟨⟨M 4 ⟩⟩ ∼ L 5.500 (6) , ⟨⟨M 6 ⟩⟩ ∼ L 9.251 (9) , and ⟨⟨M 8 ⟩⟩ ∼ L 13.00 (2) .This consists strong evidence that the convergence is really robust for these higher-order cumulants and is also in contradiction to the scaling of the susceptibility (χ ∼ L 1.661 (9)  instead of ∼ L 1.75 ) and the magnetocaloric coefficient (m T ∼ L 0.956 (7) instead of ∼ L 0.875 ).

Conclusions
We have analysed the Fisher, crystal-field, and Lee-Yang zeros of the two-dimensional Blume-Capel model at the critical and tricritical points, employing Wang-Landau and hybrid simulations of Metropolis and Wolff type.The main conclusion of our work is that the zeros of the partition function allow us to find exponents that are remarkably close to their expected values even by the simulation of very small system sizes.When compared with the exponents extracted from basic thermodynamic quantities, we notice that the zeros give more accurate results, without the need to include any scaling corrections.At the tricritical point, the crystal-field zeros allow us to obtain a very good value for the thermal exponent, unlike the Fisher zeros, which require the inclusion of correction terms.At the critical point, the opposite situation occurs.Our work calls attention to the fact that the zeros are highly sensitive to the external parameters of the system, especially around the area of the tricritical point of the phase boundary where strong crossover effects are expected to obscure any attempt for finite-size scaling.An alternative technique which measures the zeros benefiting from high-order cumulants of the energy and the order parameter has also been implemented in the present work and it was shown that this method also leads to concrete results of similar numerical accuracy, in reduced computational times.Finally, we have extended the standard analysis of the impact angle of the accumulation of Fisher zeros in the complex temperature plane to the complex crystal-field plane, the first numerical computation of this specific amplitude ratio at tricriticality.The extension to crystal-field zeros is also a novelty of the present work.

Figure 1 :
Figure 1: Phase diagram of the square-lattice Blume-Capel model in the (∆, T )-planeshowing the ferromagnetic (F) and paramagnetic (P) phases that are separated by a continuous transition at small ∆ (solid line) and a first-order at large ∆ (dotted line).The line segments meet at a tricritical point (∆ t , T t ) marked by a black rhombus.Numerical data shown are selected estimates from previous studies, including the critical point at zero crystal field which we denote hereafter for the needs of our work as T c (∆ = 0), or simply as T c .The numerical simulations reported in this paper were carried out at these two special points (∆, T ) = {(0, T c ), (∆ t , T t )}.

Figure 2 :
Figure 2: Fisher zeros at the zero crystal-field critical point specified in equation (9) in the complex temperature plane for L = 4.The two sets of blue and yellow circles show respectively changes in the sign of the real and imaginary parts of the partition function obtained from hybrid Metropolis-Wolff algorithm, while the sets of black and red crosses those obtained from Wang-Landau simulations.A crossing of the lines is a zero of the partition function.The first two Fisher zeros are denoted as [β r1 , β i1 ] and [β r2 , β i2 ].

Figure 3 :
Figure 3: Left panel: Fisher zeros at ∆ c = 0 for a system with linear size L = 4.A single hybrid Metropolis-Wolff simulation was performed at the critical temperature T c = 1.6929.The two sets of blue and red circles show respectively changes in the sign of the real and imaginary parts of the partition function.The solid lines are fourth-order degree polynomial fits and the zero sought is enclosed within the central losange.Right panel: A zoom in of the interesting area for the benefit of the reader.

Figure 4 :
Figure 4: Magnetocaloric coefficient (left panel) and fourth-order magentisation cumulant (right panel) vs. temperature.Comparison of numerical data in the small lattice-size regime at ∆ = 0 between hybrid Metropolis-Wolff and Wang-Landau simulations.

Figure 5 :
Figure 5: Fisher (left panel), crystal-field (centre panel), and Lee-Yang (right panel) zeros from Wang-Landau simulations for a system with linear dimension L = 6 in the vicinity of the critical point, see equation (9).Note the inverse critical temperature β c = 0.5907 and that H c = 0.

Figure 6 :
Figure 6: Similar to figure 4 but for the tricritical point.Note the inverse tricritical temperature β t = 1.644 and that h t = 0.

Figure 7 :
Figure 7: FSS behaviour of the Lee-Yang zeros at the critical (left panel) and tricritical (right panel) points.Comparative results from the cumulant and partition function (PF) methods are shown, via the hybrid Metropolis-Wolff algorithm.We quote the expected values y IM h = −1.875 and y tri h = −1.925for the of the reader.Lines are simple linear fits in a double-logarithmic scale.

Figure 8 :
Figure 8: FSS behaviour of the Lee-Yang zeros at three different temperatures around the tricritical point, as indicated in the panel, extracted by the cumulant method.using the hybrid Metropolis-Wolff algorithm.Lines are simple linear fits in a doublelogarithmic scale.

Figure 9 :
Figure 9: FSS behaviour of the magnetisation cumulants ⟨⟨M n ⟩⟩ at the critical (left panel) and tricritical (right panel) points according to the scaling laws ⟨⟨M n ⟩⟩ ∼ L −D+ny IM h and ⟨⟨M n ⟩⟩ ∼ L −D+ny tri h , respectively.Numerical results for n = 2, 4, 6, and 8 are shown, via the hybrid Metropolis-Wolff protocol.Note the double-logarithmic scale.

Table 1 :
Summary of critical exponents for the two-dimensional Blume-Capel model approaching the second-order zero crystal field transition point along the line ∆ c = 0. Data produced via the Wang-Landau (L = 5 − 8) and hybrid Metropolis-Wolff (L > 8) schemes.
Ising model universality class (the role of corrections-to-scaling) Quantities χ 12 ∼ L y IM

Table 2 :
Similar to table 1 but this time including corrections-to-scaling (L −ω ) terms in the fits.

Table 3 :
Summary of critical exponents for the two-dimensional Blume-Capel model at the tricritical point from FSS at the pseudocritical temperatures.The numerical approaches used for the different size spectrum are analogous to those described in table1.

Table 4 :
Similar to table 3. Exponent estimates obtained from fits at the pseudocritical crystal fields.

Table 5 :
Summary of critical exponents y IM t (Fisher zeros) and y IM h (Lee-Yang zeros) for the two-dimensional Blume-Capel model at the critical point, of equation (9), for h c = 0. Data produced via the Wang-Landau (L = 5 − 8) and hybrid Metropolis-Wolff (L > 8) schemes.

Table 6 :
Similar to table 5 but this time critical exponents were obtained crystalfield and Lee-Yang zeros at the tricritical point.

Table 7 :
Estimation of φ at different system sizes and angles (first, second, and third zeros) for the two-dimensional Blume-Capel model at the tricritical point.Expected value: φ 2D = 90°.