{'fraction_non_alphanumeric': 0.02707961731798281, 'fraction_numerical': 0.01037781741527485, 'mean_word_length': 4.571213489912677, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pushing the limits of EPD zeros method

The use of partition function zeros in the study of phase transition is growing in the last decade mainly due to improved numerical methods as well as novel formulations and analysis. In this paper the impact of different parameters choice for the energy probability distribution (EPD) zeros recently introduced by Costa et al is explored in search for optimal values. Our results indicate that the EPD method is very robust against parameter variations and only small deviations on estimated critical temperatures are found even for large variation of parameters, allowing to obtain accurate results with low computational cost. A proposal to circumvent potential convergence issues of the original algorithm is presented and validated for the case where it occurs.


Introduction
Phase transition is one of the core subjects in statistical physics. Several interesting behaviors such as superconductivity, ferromagnetism, and many others, emerge in systems after a phase transition. The concept of phase itself is somewhat difficult to broadly define but one can easily identify it in the presence of phase boundaries or phase coexistence. At those boundaries, where a phase transition occurs, it is well known that the system's free energy has a nonanalytic behavior which leads to divergences or discontinuities in its derivatives, i.e., in a phase transition it is expected to see quantities like magnetization and specific heat presenting discontinuities or divergences. To study these behaviors, several approximated methods were developed, Mean field approximations [1], the renormalization group [2], Monte Carlo [3] method and finite size scaling theory [4] are some of the most used approaches. In addition to those, Yang-Lee [5] introduced a different way to study phase transitions. They rewrote the grand canonical partition function as a polynomial in the complex fugacity plane and showed that the polynomial's roots have all the systems' thermodynamic information, such as the partition function itself. Furthermore, they showed that at the thermodynamic limit the density of roots completely determines the system critical behavior accumulating near phase transition points; those roots are called Yang-Lee zeros. The theory of zeros developed by Yang-Lee is among the few rigorous theories regarding phase transitions and is considered a cornerstone of statistical physics.
The concept of Yang-Lee zeros was extended by Fisher [6] to the canonical ensemble, where the canonical partition function was rewritten as a polynomial in the complex temperature plane. However, solving this polynomial showed to be a challenge. Even for simple systems the polynomial has a high degree easily surpassing four or more orders of magnitude. Its coefficients are given by the density of states, a quantity difficult to find and known exactly only for a few systems. Moreover, for most systems, the density of states has huge values making the polynomial even more difficult to deal with. Because of those characteristics, even for the state of the art root finder algorithms, one may expect several numerical instabilities.
The caveat of the Fisher zeros polynomial raised above was solved by Costa, Mól and Rocha in [7], where the polynomial was simplified after a transformation. The main advantage of this method lies in the fact that the transformed polynomial has as coefficients the energy probability density (EPD), being possible to judiciously discard coefficients regarding states with low probability to occur without losing relevant thermodynamic information, drastically reducing the polynomial degree.
Although the EPD zeros were successfully applied to the study of phase transitions, such as discontinuous, continuous and topological ones [7,8,9], some points still need to be clarified. Since the EPD method is highly dependent on the zero's positions, changes in the polynomial's coefficients such as coefficient discards or statistical fluctuations when estimating the EPD, may have a significant impact in quantities found by the method such as the critical temperature or critical exponent. To exemplify the problem, figure 1 shows EPDs and their respective map of zeros for different coefficient discard thresholds and different Monte Carlo steps (MCS) used to build the EPD. Although the map of zeros shown is different for each choice of coefficients discard threshold and number of Monte Carlos Steps, they all represent the same system and should give results for critical temperature and critical exponents statistically equivalents. With this in mind, in this contribution we address the effects of choosing different coefficient discard thresholds and MCS when using the EPD zeros approach in systems displaying continuous and discontinuous phase transition. In what follows we briefly describe the Fisher and EPD zeros methods and the models we considered. Results for the Ising model and the six-state Potts model are shown in sections 5.1.1 and 5.2, respectively. A convergence issue found for the Ising model is discussed in section 5.  For each threshold only states between the marked vertical lines are considered as coefficients of the polynomial. b) Map of zeros made using the EPD for 10 9 and 10 6 MCS without discarding any coefficient. c) Map of zeros made using the EPD with 10 9 MCS and with discard thresholds of 10 −4 and 10 −1 . Although the map of zeros in b) and c) are different, they represent the same system, and the EPD method should give results statistically equivalent to each other.

Fisher zeros
Following the ideas of Yang-Lee [5], Fisher [6] considered an analytic continuation of the temperature to the complex plane and rewrote the canonical partition function Z as a polynomial, where z = e −βε , β = 1/k B T is the inverse temperature, k B is the Boltzmann constant, g(E) is the density of states, the energy E is discretized by setting it to E n = ε o + nε with n ∈ N, ε 0 and ε being constants and g n = g(E n ) [10].
Since the polynomial's coefficients are positive, g(E) ≥ 0, all roots appear as conjugate pairs, and none of them are in the positive real axis. Only in the thermodynamic limit, at points where a phase transition occurs, it is expected to see a real positive root z c ; these roots are called dominant zeros or leading zeros.

EPD zeros
The EPD zeros are obtained by multiplying Eq. 1 by e −βoE e βoE = 1. Then, e −βoE is the non-normalized energy probability distribution function at β o and x = e −∆βε . We notice this corresponds to a rescaling of the Fisher zeros since x = e −βoε z. The main improvement is that since states with low probability to occur have little relevance to the system thermodynamics at β o they can be discarded. By doing that, the polynomial's degree is drastically reduced. Furthermore, as for the Fisher zeros in the thermodynamic limit, the dominant EPD zero approaches the real axis, touching it signalizing that a phase transition occurred. The EPD zeros have another important characteristic, in the thermodynamic limit, if we choose β o = β c , where β c is the critical inverse temperature, the dominant zero x c will be always at the point (1, 0). Moreover, for a β o not so close to β c , the distance from the dominant zero to the point (1, 0) can be used to estimate how far β o is from β c . Therefore, the following algorithm was proposed to iteratively approximate β o to β c even for finite systems.
2 Normalize h n (β j o ), so its maximum value is 1, and discard coefficients smaller than a given threshold.
3 Find the zeros of the polynomial with coefficients given by h n (β j o ).
The main idea behind it is that the point (1, 0) may behave as an attraction point for systems with a phase transition when coefficients are discarded and finite systems are considered. In this way, as close to the transition point the EPD is drawn, more precise the estimate of the transition temperature is, since for reasonable choices of coefficients discard thresholds only negligible information for the phase transition are disregarded.

Models
In this work, we consider two very well known spin models that serve as benchmarks in the study of continuous and discontinuous phase transitions, the Ising and Potts models, that we briefly present in what follows.
The Ising model Hamiltonian is given by, where J > 0 is a coupling constant, σ i is a spin at site i with value ±1 and < i, j > means a sum over nearest neighbors. The 2D version of this model in a square lattice has a continuous phase transition at T c = 2/ln(1 + √ 2)J/k b and a critical exponent ν = 1 [11].
The q-states Potts model Hamiltonian is given by, where σ i is a spin on the i-th site that assumes discrete values in range [1, ..., q], J > 0 is a coupling constant, δ is the Kronecker delta, and the sum is taken over the nearest neighbors. This model presents a continuous and a discontinuous phase transition for q ≤ 4 and q > 4 respectively in two dimensions. Its critical temperature is known exactly to be T c = 1/ln(1 + √ q)J/k b [11]. In special, for q = 2 this model can be mapped into the Ising model by a simple energy rescaling.

Results
The simulations were carried out as follows. For the Ising model, a single spin flip process was used and three temperatures close to the critical one were simulated using the Metropolis algorithm for lattices with lateral size ranging from L = 90 to L = 180. For the Potts model, the Wolff [12] algorithm was used in addition to the single spin flip and fifteen temperatures very close to the critical one were simulated, for lattices with lateral size ranging from L = 90 to L = 180. Moreover, multiple histograms reweighting technique [13] was applied to create the EPDs used in the algorithm described in section 3.
As will be described in more detail in the next section, to estimate the critical temperature and critical exponent, it is necessary to find the mean value of T c (L) and m(x c (L)) for each lattice size L. Therefore, for each lattice size, 5 different simulations were made at each temperature. Then, the multiple histograms reweighting technique together with the algorithm described in section 3 was used to find the values of T c (L) and m(x c (L)). With these values, it is straightforward to find the mean and the errors of T c and m(x c ) for each lattice size. This process was used for the Ising model and the six-state Potts model using 10 6 , 10 7 and 10 8 MCS with discard threshold ranging from 10 −100 to 10 −1 .
In sections 5.1.1 and 5.2, we present tables summarizing the results obtained for T c and ν with different discard thresholds and different MCS. Furthermore, images shown in the next sections, unless otherwise stated, were obtained using 10 6 MCS for each simulated temperature and with lattice size of L = 150.

Critical temperature and critical exponent
Estimates of the critical temperature were found using the following finite size scaling equation: and since the imaginary part of the dominant zero goes to zero in the thermodynamic limit, the following scaling equation is used to estimate the critical exponent, m(x c (L)) ∼ bL −1/ν .
The data used in these estimates were obtained by an automated application of the original algorithm, corrected for a convergence issue we found and discuss in the following sections. The results obtained for T c and ν were summarized in tables 1 and 2 for several combinations of thresholds and MCS. Moreover, a typical regression for the critical temperature and critical exponent is shown in figure 2, for the sake of clarity, the figure shows only the regression with 10 6 MCS. As can be seen in tables 1 and 2, all the results are close to the expected ones, ν = 1 and T c = 2.2691J/k b , independent of the threshold and the MCS chosen.

Convergence issue and the dominant zero
Considering all the partition function zeros, the dominant zero is the most important one for identifying a phase transition. However, find it among other zeros may not be an easy task for automated applications of the algorithm. Indeed, automated applications rely simply on the identification of the dominant zero as the one closer to the point (1,0) regardless of other factors. Here we found that there are some pathological situations where zeros that are not determinant to the phase transition may become closer to the point (1, 0) than the dominant zero. As a consequence, the algorithm convergence to the transition temperature is affected, since non-dominant zeros may be confused with the dominant one. As an example, in figure 3a it is shown the algorithm convergence for several coefficient discard thresholds and in figure 3b an oscillating pattern caused by the algorithm incorrectly selecting the dominant zero. It is important to highlight that systems with lattice of size L = [90, 120, 150, 180] were used, but only for L ≥ 150 the convergence issue was observed.

Solving the convergence issue
To solve this convergence issue some proprieties of the dominant zeros were explored. At this point it is worthy to note that we defined as dominant zero the zero closest to the point (1, 0). With that said, in a previous study [14], it was observed that the dominant zero is stable to random perturbations in the polynomial coefficients, i.e., applying a perturbation in the EPD following the equation, where h(β o ) is the EPD, h * (β o ) is the perturbed EPD, ε ∈ [−1, 1] is a uniform random variable and a the perturbation, does not change the dominant zero position appreciably, as can be seen in figure 4. Even for different thresholds, the dominant zero is stable, although it is not the only stable zero, as can be seen in figure 5. Therefore, only the stability it is not enough to identify the dominant zero. In fact, the zeros which causes the convergence issue for the threshold of 10 −2 are highly stable, as can be seen in figure 5. To correctly identify the dominant zero, it was added a bias to the algorithm to select among the most stable zeros, the zero closest to the line x = 1 as the dominant one. Then, the following algorithm was proposed to avoid the selection of undesired zeros in automated application of the EPD zeros approach.

Six-states Potts model
The six-states Potts model has a more complex EPD due to its discontinuous phase transition, and more care is necessary at the discard step. As can be seen in figure 7, the EPD has two large peaks, and as described in [7] states between those peaks that has low probability to occur should not be discarded. Therefore, even for higher threshold such as 10 −1 , since intermediate states have to be considered, almost no state is discarded. A reflection of such protection of states is seen in the critical temperature obtained for different threshold, table 3, where no difference in T c is observed at any threshold. Despite this, all the values found for T c are close to the expected one of T c = 0.807606J/k b independent of the coefficients discard threshold and MCS that have been chosen.
In contrast to what we observed for the Ising model, no glimpse of convergence issues was observed in this system. Indeed, as is clear from figure 7, there are no adjacent zeros to the dominant one in this case. Therefore, automated implementations are not expected to select undesired zeros.  One final remark about the six-state Potts model is that although in Ref. [7] a simple linear adjustment did not work to find T c , i.e., a correction to scaling had to be used to consistently estimate T c , in this work, using larger lattice sizes, a correction to scaling it is not necessary. Indeed, while in Ref. [7] the lattice sizes Thus, as can be seen in figure 8 the finite size scaling for the critical temperature and imaginary part of the dominant zero for larger sizes and the results obtained by the linear fits (see Table 3) shows no need of corrections to scaling.

Conclusion
In this paper, we studied how the choice of the coefficients discard threshold and MCS used to estimate coefficients impact automated implementations of the EPD method and its capability to find the critical temperature and critical exponent. To do this, and thinking of expanding our results for continuous and discontinuous phase transitions, we used 2D Ising and six-state Potts models.
The EPD method showed to be robust against the coefficients discard threshold and MCS used, given results for T c and ν close to the expected and that agree very well with each other even for considerable differences in the threshold and MCS. This allow considerable speedups in the estimation of transition temperatures since estimates are almost insensitive to reasonable choice of the parameters. Indeed, results obtained using 10 7 or 10 8 MCS agree withing error bars with minimal accuracy difference. This means that about 10 times less computational effort is necessary to estimate critical temperature and exponent with almost the same accuracy. As expected, we also find that the correction to scaling used in Ref. [7] for the first order transition of the six-states Potts model is not necessary for larger lattice sizes.
Despite the robustness of the EPD method, the Ising model showed a convergence issue for a specific choice of the coefficients discard threshold, 10 −2 , due to a misleading automated choice of the dominant zero. This behavior led to an oscillating pattern in the iteration of the EPD algorithm instead of a convergent one. To solve this issue, a bias was added to the algorithm to correctly select the dominant zero and ignore other adjacent zeros nearby that may be mistakenly chosen. Although this corrects the issue, we must point out that this oscillatory behavior for the automated convergence of the Ising model was found for this specific choice of coefficients discard threshold. The solution present here is not unique and it may even be much more advantageous to reduce the number of zeros by using higher coefficients discard thresholds without significant impact in the final estimates of the critical temperature and exponent than employing the algorithm to deal with this oscillatory behavior. Indeed, it is remarkable how the estimates are almost insensitive to variations on the coefficients discard threshold and on the number of MCS considered in the building of the EPD and map of zeros. Of course, for more complicated models, the algorithm presented here may be very useful in preventing this oscillatory behavior.