{'fraction_non_alphanumeric': 0.04095810116506259, 'fraction_numerical': 0.019502134742021855, 'mean_word_length': 4.150764070070816, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, '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}}
Numerical study of the overlap Lee-Yang singularities in the three-dimensional Edwards-Anderson model

We have characterized numerically, using the Janus computer, the Lee-Yang complex singularities related to the overlap in the 3D Ising spin glass with binary couplings in a wide range of temperatures (both in the critical and in the spin-glass phase). Studying the behavior of the zeros at the critical point, we have obtained an accurate measurement of the anomalous dimension in very good agreement with the values quoted in the literature. In addition, by studying the density of the zeros we have been able to characterize the phase transition and to investigate the Edwards-Anderson order parameter in the spin-glass phase, finding agreement with the values obtained using more conventional techniques.


Introduction
In two seminal papers, T. D. Lee and C. N. Yang [1,2] introduced a new tool to understand the origin of a phase transition by studying the complex singularities of the free energy, or, equivalently, the zeros of the partition function. In particular, they showed that all the zeros are located on the unit circumference on the complex activity plane (taking as variable z = e −2h , where h is the external magnetic field). They also proved that if the zeros condense onto the real axis when V → ∞ a phase transition takes place. Finally, they related, in the low-temperature phase, the density of zeros with the discontinuity in the order parameter (remember that the Ising model experiences a firstorder phase transition when h changes at a fixed temperature below the critical one). This approach was subsequently extended to the temperature zeros by Fisher [3,4,5,6].
We have Lee-Yang like theorems only for a limited class of non-disordered systems (such as Ising models). However, it is possible to develop a scaling theory by assuming that asymptotically the complex singularities (wherever they lie) touch the real axis (thus generating the phase transition). Hence, despite the lack of formal theorems it is still possible to apply Lee and Yang's main results to a wide class of systems (e.g., Potts models [7]). In this class of systems the zeros do not live on a circle as stated by the Lee-Yang theorem, but they still control the critical properties of the model. We will only assume this last fact irrespectively of the form of the locus of the zeros in the complex plane [7].
In Refs. [8,9,10] the analysis of complex singularities was applied successfully to diluted systems (in particular diluted Ising models in two and four dimensions). The key point for the applicability of the standard results, well tested in non-disordered systems, is to compute the complex singularities individually for each disorder realization (called sample) and then use the mean of the individual zeros (sample zeros) in order to test the scaling properties of the zeros and to study the properties of the integrated density of zeros. In this work we will also introduce the analysis of the median.
Nowadays we are interested in gaining a deeper understanding (from the point of view of the complex singularities) of the properties of an interesting frustrated and disordered system: the three-dimensional Ising spin glass. The magnetization, while very interesting in off-equilibrium dynamics and in experiments, plays no role in the critical behavior and in the understanding of the low-temperature properties in a finitedimensional spin glass. The observable that controls this spin-glass phase is the overlap. Hence, in this work we have focused on the numerical study of the complex singularities linked with the overlap in order to study the phase transition and the properties of the spin-glass phase.
In the past, Lee-Yang and Fisher zeros were obtained in spin glasses by means of the numerical evaluation of the partition function on small lattices [11,12,13,14]. This methodology was also applied to models defined on Bethe lattices [15]. Finally, some calculations were performed with the help of replicas [16].
More recently, the complex singularities linked with the external magnetic field were studied for the two and three-dimensional Ising spin glass model in the interesting reference [17], which focuses on the Griffiths singularity and computes all the zeros for small lattices.
In particular we are interested in characterizing the scaling of the individual zeros at the critical point (which will allow us to compute the anomalous dimension exponent) and checking the scaling in the spin-glass region. In addition, we want to study the properties of the density of zeros in the critical and spin-glass region: the behavior of this observable will clearly signal the phase transition. Finally, we will show how this density of zeros can be used to compute the Edwards-Anderson order parameter. However, the spin-glass susceptibility presents strong scaling corrections (even on an L = 32 lattice and β = 1.4), which induce strong corrections on the density of zeros allowing us (from the numerical point of view) only to test our density of zeros against the values of q EA found in the literature, rather than attempting a direct numerical computation of the order parameter. We want to stress that in cases in which the spin glass susceptibility reaches the asymptotic value, the method we propose will be able to provide directly the order parameter (q EA ) giving an additional method to those used nowadays [18,19,20,21].
Let us finally mention that we have obtained the data presented in this work from the analysis of the configurations produced in parallel tempering runs [20,21] using the Janus computer [22,23,24].

Model and observables
We have studied the three-dimensional Edwards-Anderson model with dynamical variables σ i . These variables are Ising spins and are placed on the nodes of a cubic lattice of linear dimension L and volume V = L 3 . The Hamiltonian of the system is where ij indicates that the sum is over the nearest neighbouring sites. The couplings J ij are random quenched constants with bimodal probability distribution, that is, J = ±1 with 50% probability. Every realization of the couplings is called a sample. Due to the fact that we have a random Hamiltonian, we have to deal with a double average: first the thermal average, which we will denote by (· · ·) , and then the average over the samples, which we will denote by (· · ·). We have simulated several real replicas of the system, so we can define the local overlap where σ (1) x belongs to the first replica and σ (2) x belongs to the second one. The spin overlap is defined from this local overlap as where the sum runs over the whole volume (V ). In addition, we define q ≡ Q/V . These observables allow us to define some new quantities, for example the non-connected spinglass susceptibility Let us now rewrite the Hamiltonian adding a new perturbation Q and including the two replicas explicitly, This Hamiltonian looks like that of the Ising model in a magnetic field We can write the partition function, whose zeros we want to study, as Let Z 0 the partition function of the non-perturbed system, so and we have to find the zeros of the function cos(β Q) since in absence of a magnetic field sin(β Q) is zero. The algorithm to find them is quite easy: we start from the list of individual measurements of Q for each sample (see Section 4) and evaluate the average cos(β Q) , increasing in small steps ∆ . When the function changes signs from one step to the next, we have found a zero in this interval. Obviously, the smaller ∆ the better the precision of the zero that we have found, but also the slower the analysis, so we have to be careful with the error estimates. We have analyzed the first four zeros of this function.

Finite-Size Scaling
One can obtain the expected behavior of the LY zeros by means of (see for example Ref.
therefore, the finite-size dependence, at the critical point, of the Yang-Lee zeros can be expressed as: where and D is the dimensionality of the system, being D = 3 in this work. If corrections to scaling are taken into account, the above relation becomes where x 2 is the leading correction-to-scaling exponent, x 2 = ω.
In the broken symmetry phase, where the non-linear susceptibility diverges as the volume of the system, we expect the following behavior: We can take scaling corrections into account, as in the critical point, and where x 3 is the leading correction-to-scaling exponent in the broken phase. ‡. In order to discuss the density of zeros we need to describe some known properties of the Hamiltonian defined in Eq. (5). This Hamiltonian was introduced in the past [25,26]. In particular it experiences a first-order phase transition in , below the critical temperature of the uncoupled model. Hence, the overlap is discontinuous: being the discontinuity at the transition just 2q EA .
We can also introduce the density of zeros and its integrated version which takes the following value computed for a given zero: where j labels the j-th zero (j = 1, 2, . . .). In order to deal with the discontinuous behavior of G( ) at the zeros, we follow the recipe of references [27,28] and use the mean between two consecutive plateau values (j − 1 and j). Anyhow, the asymptotic value of the integrated density computed in the j-th zero is j/V . We will discuss this point again in subsection 5.3. This integrated density is very useful to characterize a phase transition. In general it behaves as ‡ Both droplet and RSB predict algebraic decays for the connected correlation functions in the spin glass phase (the spin glass phase is critical in both models). In the droplet model the exponent of the decay is y (sometimes denoted as θ), so one can show that x 3 = y. In RSB depending of the value of q we have different decays (of the q-constrained correlation functions), denoting the decay exponents as θ(q). So the leading correction exponent can be shown to be the smallest of the different θ(q). See Refs. [20,21] for a detailed discussion on θ(q). and we can extract a great amount of physical information from these three numbers (a 1 , a 2 and a 3 ): • In the symmetric phase a 3 < 0. In a broken phase a 3 > 0.
• In the onset of a first-order phase transition, varying as it is our case: a 2 = 1 and a 3 = 0. In addition we can extract the order parameter of the broken phase: • At the critical point, a 3 = 0 and a 2 is related with the anomalous dimension η by means:

Simulation details
We have run simulations for several lattice sizes on the Janus supercomputer [22,23,24] (for L = 16, 24, and 32) and on conventional computers (for L = 8 and 12). These simulations were originally reported in [20], which gives full details on the chosen parameters and the thermalization protocol. In this section we give only a brief summary.
We have used the parallel tempering algorithm [29,30], choosing the temperatures to maintain an acceptance around 20% in parallel tempering updates. Besides, since Janus needs far more time to do a parallel tempering update than a heat-bath one, we have chosen to do one parallel tempering update every 10 heat-bath ones. In table 1 one can find a summary of the simulations parameters. In order to choose the simulation length, we have assessed thermalization on a sample-by-sample basis, using the temperature random walk technique [31,20] (table 1 gives the average number of lattice updates for each L).
In general, each of the single processors (FPGAs) of Janus takes care of the simulation of one replica of the system. However, some samples have such a slow dynamics that even with this algorithm the simulation would be too long (more than six months of continuous running time), so we would need to accelerate it. For these few cases we have created a special low-level code that is in charge of the parallel tempering in the control FPGA of a board of Janus. This allows us to spread the simulation § In Lee and Yang's paper, the starting point is the Hamiltonian βH h = βH h=0 + hM , where M is the total magnetization of the system. In terms of the fugacity z = e −2h , they obtained the following result (valid below the critical temperature) for the density of zeros (in the fugacity variable that we will denote as µ z (z)): µ z (0) = m sp /(2π), where m sp is the spontaneous magnetization below the critical temperature. In order to transfer this result to our notation we remark that our "magnetic field" is β , q EA plays the role of m sp and we need to use the standard law of the transformation of the probability densities (µ z (z) = µ ( )|d /dz|, where z = exp(−2β )), obtaining: Notice that near = 0 we can identify a 1 with µ (0). over several processors running only a subset of temperatures in each FPGA, thus accelerating the simulation by increasing the parallelism. Table 1. Summary of the simulations. N T is the number of simulated temperatures (evenly spaced between T min and T max ); N mes is the number of Monte Carlo steps (updates of the whole lattice) between measurements; N med HB is the average simulation time (since we use the random-walk technique the simulation time depends on the sample); N sam is the number of simulated samples. We have simulated four real replicas for each sample.  We have saved on disk every individual measurement of the overlap. Since we have simulated four real replicas of the system, for each sample we have a total of 6N HB /N mes values of Q. Given the variable N HB , this ranges from 1.2 × 10 5 to 2 × 10 7 measurements for our largest lattices, so we have a very good precision for computing the zeros of the partition function. We have discarded the first half of the measurements for equilibration. We want to study the behavior of the system in the critical temperature and in the low-temperature phase of the system, analyzing the scaling of the zeros. Therefore, we need to compute the zeros of different system linear sizes, L, but at the same temperature. Since we have not simulated the same temperatures for every lattice size, we have interpolated, using cubic splines, in order to estimate the zero at each of the chosen scaling temperatures.

Results
In this section we will study the behavior of the zeros as a function of the lattice size, both in the critical and in the spin-glass phase. Finally, we will compute the density of zeros and extract the η exponent from the analysis at the critical temperature and the Edwards-Anderson order parameter from the scaling in the low-temperature phase. Table 2. Fits of the zeros to j (L) = A j L −x1 , for L ≥ L min . As we can see, with L min = 8 the χ 2 per degree of freedom is acceptable only for j = 1, 2, but with L min = 12 all the zeros have a reasonable fit. However, the value of x 1 grows with j, an indication that we have to consider corrections to scaling (see text).  (7) 2.49/2

Scaling at the Critical Point
We first consider the scaling at the critical point and use it to determine the anomalous dimension, as in (10). Our simulations were optimized to investigate the lowtemperature phase, for large system sizes, rather than to obtain the critical parameters. Therefore, we take the value of β c = 0.902(8) from [32], which features many more samples and small sizes to control scaling corrections but does not reach the lowtemperature phase, and will also use this reference to check our value of η. Let us first consider a fit of the individual zeros, leaving aside corrections to scaling, i.e., following (10). For the j-th zero, we fit to We report the results of these fits in table 2. We see that the first and second zeros follow (22) very well for L ≥ 8, but for j > 2 we need to restrict the fit to L ≥ 12. However, there is an inconsistency in the results: the value of x 1 should be the same for all zeros, but we see that it increases with j. Moreover, at least for the larger j, x 1 is incompatible with the expected value of x 1 = 2.688(5), (taking η = −0.375(10) from [32]) This hints that corrections to scaling should be taken into account, as in (12).
In order to do so, we consider all values of j at the same time and perform a global fit, enforcing data from different zeros to share the same x 1 and x 2 . As points coming from a given L are correlated, the full covariance matrix has to be considered. We label our set of points { j (L a )} by their L and their j: we have data for L = 5 different values of L (L 1 = 8, L 2 = 12, L 3 = 16, L 4 = 24, L 5 = 32) and for j = 1, 2, 3, 4. The If we combine the critical exponents of [32] with the Janus simulations studied herein, we obtain a compatible value of β c = 0.905(7) [33]. appropriate goodness-of-fit estimator is, therefore, where σ (ia)(jb) is the covariance matrix of the set of zeros (which is block diagonal, since data for different L are uncorrelated). Unfortunately, we do not have enough data to determine x 2 and x 1 at the same time (the resulting error in ω would be greater than 100%). Instead, we take x 2 = ω = 1.0(1) from [32] and fit only for x 1 and the amplitudes. The resulting fit for L ≥ 12, shown in figure 1, gives where the error in square brackets accounts for the uncertainty in ω. Our determination of x 1 is now compatible with the expected value of x 1 = 2.688 (5). Therefore, the scaling of the zeros is consistent at the critical point.

Scaling in the low-temperature phase
Now we consider the scaling of j (L) in the low-temperature phase. This time, we expect, from (13), with x 1 = D. Table 3. Scaling of the zeros in the low-temperature phase. For the two considered temperatures (β = 1.2, 1.4) we first show a fit without corrections to scaling for L ≥ 16, that is j (L) A j L −x1 . As explained in Section 5.1, this is a global fit for the four zeros, considering their full covariance matrix. We then consider the same fit with corrections to scaling, trying different values for ω (see the text for more details). In all cases x 1 is smaller than the expected value We have fitted the data for β = 1.2 ( Figure 3) and β = 1.4 ( Figure 2) to (25). ¶ The results, in Table 3, show a value of x 1 incompatible with x 1 = D = 3. We have also included corrections to scaling, using both ω = 1 (Goldstone-like correction) [34] and ω = 3 (Ising ordered correction) [34], ω = y = 0.255 (droplet) [35,36,37], ω = θ(0) = 0.39 (replicon and also 1/ν which controls the scaling correction of q EA (L) [21]), ω = 0.79 = 2θ(0) (twice the replicon [21]) and ω = 0.65 = θ(q EA ) [21] but in neither case is the asymptotic x 1 = D behavior recovered (see Table 3). In addition, we have forced the fits with x 1 = 3 and leaving free ω and the statistical quality of the fits was bad.
The origin of this discrepancy with the standard theory can be understood using Eq. 9. Notice from this equation that the scaling of the zeros depends strongly on the behavior of the non-connected spin glass susceptibility, so only with a divergence of this observable as the volume, we can recover x 1 = 3. However, for these two temperatures ¶ The crossover length L c which marks the change between the criticality induced by the critical point at T c (L < L c ) and that of the spin glass phase (L > L c ) has been computed for different values of β in reference [20]. In particular, we know that L c (β = 1.2) 6 and L c (β c = 1.4)) 2.5. Hence all the data presented in this section belong to temperatures which lie deep into the spin glass phase. In other words, we can only see the critical effects induced by the spin-glass phase itself, which is critical, not those of the critical point at T c .
(β = 1.2 and 1.4) this is not the case (see Fig. 4). Notice that q 2 = χ/V has not reached the plateau asymptotic value + : Hence at these temperatures the spin-glass susceptibility does not yet diverge as the volume.

Behavior of the integrated density of zeroes
We will start our analysis of the integrated density of zeros by plotting this density at the critical point in Fig.5. One can see that the largest lattices follow a pure power law as predicted by the theory. The slope, on a log-log scale, of this straight line should correspond with an exponent a 2 . Fitting only the L = 32 points we obtain a 2 = 1.16 (2) in good agreement with the theory a 2 = 1.116(2) (using Eq. 21 and η = −0.375(10)).
To obtain this figure we have discarded in the fit the first zero. * + In a spin glass phase, both the droplet as the RSB theory predict power law corrections on the lattice size, so the approach to the infinite volume values is really slow. * This phenomenon has been previously found in the literature. For example, the authors of [28] studied the anisotropic Ising model at the critical point and found a different behavior of the first zero in the study of the integrated density. This model shows a spreading distribution of the zeros in the fugacity complex plane. The authors suggest that the effect of this spreading distribution of the zeros is modifying the behavior of the first zero. We have not computed the complete distribution of the For large L we should expect a good collapse of all points in the same power law curve: the non collapsing part (small L in the figure) is due to the presence of scaling corrections (which we also found in the previous sections). Now we will check the theoretical predictions for the integrated density of zeros in the broken phase, which predict a linear behavior in the perturbing parameter . Notice that in our case the margins between the critical point and the broken phase are tight since in the infinite volume limit we will see a behavior 1.116 at the critical point which changes just below T c to (of course, this is due to the value of the η exponent).
In Figs 6 and 7 we show that the data nearly follow a linear behavior of the integrated density deep in the spin glass phase (more concretely at β = 1.2 and β = 1.4), in particular for L ≥ 24. The non-collapsing part of the curve is produced by the presence of scaling corrections as at the critical point.
However, it is easy to show that if the zeros do not follow (for the lattice sizes simulated), in the broken phase, a scaling as the inverse of the volume, then the complex zeros (only in the straight line i ), nevertheless, we know from reference [17] that the zeros spread in the magnetic field complex plane, so it is quite natural to assume that we will have a similar (spreading) spatial distribution of the zeros in . Another possible explanation is that the behavior of the integrated density of zeros as j/L 3 is only asymptotic. These anomalies affect only the lower order zeros. Notice that this phenomenon affects only to the pre-factor of the power law of the smallest zeros. We have seen in subsection 5.1 that the first zero scales with the right power law. We thank R. Kenna for interesting comments regarding this behavior.  integrated density of zeros does not follow exactly a linear behavior, since a 2 = D/x 1 .
We have discussed at the end of Sec. 5.2 that this lack of 1/V behavior is related to a susceptibility that is not yet diverging as the volume. In sec. 5.2 we have found an exponent x 1 = 2.842(11) for β = 1.2 and x 1 = 2.844(10) for β = 1.4, which implies that a 2 = 1.056(4) and a 2 = 1.055(3) for β = 1.2 and β = 1.4 respectively. In Fig. 8 we show the behavior of the integrated density of zeros computed for our largest lattice (L = 32) and lowest temperature, (β = 1.4). Notice the points are not lying on a straight line. A fit to a 1 a 2 works well, with a 2 consistent with the value computed from x 1 (a 2 = 1.068(10)). So we have obtained, numerically, a 2 = 1. 16(2), at the critical point which has changed to 1.068 (10) in the broken phase.
In this situation, we cannot compute the order parameter directly from the linear behavior of the integrated density since we are not observing a fully linear behavior. Hence, we confront our numerical data for G( ) against the theoretical prediction for really small , which is G( ) = (βq EA /π) . It is interesting to note (see Refs. [27] and [34]) that we can recover the exact slope for a given lattice size if we substitute the value We can do the same analysis with the x 1 exponent: we have obtained at the critical point x 1 = 2.67 (7), which should change in the broken phase to x 1 = 3, although we actually see with our numerical data x 1 = 2.842 (11). of the order parameter computed for this lattice size. We have followed this advice, and we show in Figs. 10 and 11 our data for G( ) at β = 1.2 and 1.4 showing L = 32 data. In addition we have plotted the asymptotic slope using the order parameter (q EA ) computed for L = 32 lattices for these two temperatures in Ref. [20]. Notice that we have a slow approach to the right slope, but also that the overall picture seems to be correct.
In order to gain a better understanding of this effect, we have computed the density of zeros not with the average of the sample zeros but with the median of the probability distribution of the zeros. † † We show in Fig. 9 the histogram of the 1000 first zeros computed on the L = 32 lattice at β = 1.4. Notice from this figure the asymmetry of the histogram and the presence of events at large values of the zeros (which induces a large and strongly fluctuating value of the mean).
One can see in Figs. 10 and 11 that the median data produce an improved scaling, compared with those obtained from the mean, when comparing the data with the † † The probability distributions one usually finds in disordered systems present long tails due the presence of rare events, hence, the study of the median of this kind of distributions is also very useful (see for example [38,33]). ǫ Mean: L=32 Median:L=32 theor. slope at origin (q EA =0.7114) Figure 10. Integrated density of the zeros, for the largest lattice L = 32 and temperature β = 1.2 using the average of zeros. We have also plotted the median values. We have marked the expected slope at the origin, using the Edwards-Anderson order parameter computed in Ref. [20] for the L = 32 lattice.
analytical prediction (slope provided by q EA ).
For the sake of completeness, we can cite that the integrated density of zeros using the medians does not behave completely linearly but with a law 1.06(2) (for β = 1.4).

Conclusions
By studying the complex singularities linked with the overlap we have obtained a clear picture of the critical region and of the low temperature phase fully compatible with that obtained by other more standard approaches.
In particular, we have studied the behavior of the individual zeros as well as the integrated density at the critical point. In both cases we have obtained good values for the η exponent and we have seen that the data are compatible with the corrections to scaling published in the literature [32].
Finally, we have checked the scaling laws in the spin-glass phase, obtaining strong scaling corrections as found previously [20]. In addition we have obtained, by monitoring the behavior of the integrated density, a compatible picture using the zeros with that obtained from the order parameter of the model (q EA ) computed in finite volumes with standard methods. We have also shown that the use of the median instead of the mean improves the overall picture. ǫ Mean: L=32 Median:L=32 theor. slope at origin (q EA =0.7114) Figure 11. Integrated density of the zeros, for the largest lattice L = 32 and lowest temperature β = 1.4 using the average of zeros. We have also plotted the median values. We have marked the expected slope at the origin, using the Edwards-Anderson order parameter computed in Ref. [20] for the L = 32 lattice.