{'fraction_non_alphanumeric': 0.04324086006918538, 'fraction_numerical': 0.018992064030387303, 'mean_word_length': 4.143380429094715, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Equation of State for Spin Systems with Goldstone Bosons: the 3d O(4) Case

We propose an improved parametric form for the equation of state of three-dimensional O(N) spin systems. The proposed form is a series expansion with two sets of terms, which contribute (mainly) separately to the description of the high- and low-temperature regions of the phase diagram. Our goal is a better description of the low-temperature phase at zero magnetic field (i.e. the coexistence line), characterized by singularities induced by Goldstone modes. We test our proposed form by comparison with existing Monte Carlo data for the N=4 case, which is of interest in studies of the QCD phase transition and for which the Goldstone-mode effects are quite pronounced. We find that the description of the numerical equation of state is indeed improved with respect to other fitting forms. In all cases considered we determine the coefficients nonperturbatively, from fits to the data. As a consequence, we are able to obtain a very precise characterization of the pseudo-critical line for the model.


Introduction
The O(N) (or, more specifically, the N-vector) spin models correspond to a generalization of the Ising model to the case of the continuous symmetry of rotation. The spin variables S i are taken as vectors on a sphere of unit radius in an N-dimensional space. We consider N ≥ 2. The Hamiltonian is defined in terms of the scalar product of nearest-neighbor spins on a three-dimensional square lattice as where J > 0 represents the ferromagnetic coupling and H the external magnetic field. These models are of general interest for the statistical mechanics of phase transitions [1]. The N = 2 case (also known as the XY model) describes the superfluid transition in liquid helium and the N = 3 case corresponds to the classical version of the Heisenberg model for ferromagnets. 1 Moreover, it is believed that the N = 4 case describes the chiral phase transition in finite-temperature QCD with two degenerate light quark flavors, which makes this class of models interesting to high-energy physics as well. 2 In this case, the magnetization and the magnetic field of the spin model correspond respectively to the chiral condensate and to the quark mass for the QCD analogue of the transition [2,3,4].
The O(N) symmetry is exact in the Hamiltonian for H = 0, just like the reflection symmetry for the Ising model. The main difference with respect to the Ising case is the possibility of configurations where the spins are locally aligned but for long distances this alignment is lost, yielding a null average for the magnetization. Such configurations -called spin waves -possess arbitrarily low energy and tend to destroy the order of the system even at low temperatures. In fact, as opposed to the Ising model, the O(N) models do not display a phase transition with spontaneous magnetization 3 in d = 2. In d = 3 a phase transition occurs, with the presence of spontaneous magnetization below the critical temperature. The breaking of the (continuous) rotational symmetry at low temperatures, signaled by the spontaneous magnetization, is associated with Goldstone modes, the spin waves. These modes cause the divergence of the zero-field susceptibility not only at the critical temperature, but for the entire low-temperature phase [1,5]. Note that the magnetic field defines a privileged direction in spin space and the magnetization M is the expectation value of the spin component along H. There are thus N − 1 massless Goldstone modes, corresponding to the N − 1 transverse spin components.
Spin models in the O(N) class have been extensively studied using analytic and numerical methods (see [6] for a recent review). In particular, the nonperturbative study by Monte Carlo simulations is very efficient for these models due to the Swendsen-Wang cluster algorithms [7], which can be applied to the continuous-spin case by means of the embedding technique introduced by Wolff [8]. This study is important to test the perturbative predictions and to investigate cases for which these predictions are not available, or cannot be done with great accuracy. These problems include properties of the models in the presence of magnetic field and the direct calculation of longdistance observables such as the correlation length. For example, the predicted singular behavior of the longitudinal susceptibility for vanishing H -mentioned above and induced by Goldstone modes at low temperatures -was directly observed in Monte Carlo simulations of the cases N = 2, 4, 6 respectively in Refs. [9], [10] and [11].
Here we consider the determination of the magnetic equation of state, which gives the relation between applied field, temperature and magnetization for the system. The equation of state has been determined perturbatively for general N by ǫ-expansions (see [1,Chapter 29] and references therein) and for the cases N = 2, 3, 4 by matching a high-temperature expansion (with coefficients obtained from perturbation theory) to a parametric form incorporating the leading Goldstone-mode behavior [6,12]. Of course it is interesting to compare these expressions to Monte Carlo data for the equation of state. One can also test the various forms used in the perturbative expansions (or new proposed forms) by fitting them to the Monte Carlo results and obtaining nonperturbative coefficients. This has been done (see e.g. [10]) using an interpolation of the low-temperature (Goldstone-mode) form derived in [13] with a high-temperature form determined by analyticity conditions. This method has the advantage of a clear low-temperature form, with several orders in the Goldstone-mode expansion, but has the disadvantage of needing an interpolation with the high-temperature form.
In the present paper we carry out fits using instead a variant of Josephson's parametrization [1,14], a polynomial parametric representation for the equation of state. The resulting representation is valid above and below the critical temperature and automatically satisfies the analyticity conditions mentioned above. In addition to the leading (multiplicative) Goldstone-mode contribution, we consider explicitly the higher-order terms, which are important in the low-temperature region. Our proposed form contains two sets of coefficients, which will be separately more relevant for the description of the high-or low-temperature regimes. We argue that the use of this double set of coefficients enables a better characterization of the two regimes, leading to better fits in the comparison with numerical data. This claim is verified by an application to the N = 4 case, for which the Goldstone-mode effects are fairly high, using the data reported in [10]. As mentioned above, this case is of interest for comparison with data from numerical simulations of the phase transition in two-flavor QCD. In particular, the prediction of universal behavior in the O(4) class has been confirmed for lattice-QCD data in the Wilson-fermion case [15], but not for the staggered-fermion formulation, which is believed to be the appropriate formulation for studies of the chiral region. (At the same time, some recent numerical studies suggest that the transition may be of first order [16].) We plan to extend our analysis to the N = 2 case, for which we are generating new data [17].
The paper is organized as follows. In Section 2 we describe the usual parametric representation for the equation of state, as well as our proposed form. In Sections 3 and 4 we consider the determination of important universal properties that can be obtained from the equation of state: some critical amplitude ratios and the characterization of the pseudo-critical line (respectively in Section 3 and in Section 4). Finally, in Sections 5 and 6 we present our results and conclusions.

Scaling equation of state
The magnetic scaling equation of state is given [1,Chapter 29] by where t and h are the reduced temperature t = (T − T c )/T 0 and magnetic field h = H/H 0 . We fix the normalization constants T 0 and H 0 by requiring unit critical amplitudes in the behavior of the magnetization along the coexistence line (given by t → 0 − , h = 0) and along the critical isotherm (given by h → 0, t = 0), corresponding respectively to M = (−t) β and M = h 1/δ . The equation of state can also be written as where Note that the coexistence line and the critical isotherm are given respectively by x = −1 and x = 0. The corresponding normalization conditions are thus For large values of x (i.e. in the high-temperature region of the phase diagram) the behavior of f (x) is described by Griffiths's analyticity condition [1] f (x) = ∞ n=1 a n x γ−2(n−1)β .
As said in the Introduction, at low temperatures there appear divergences in the zero-field magnetic susceptibility, due to transverse fluctuations from the massless Goldstone modes [1,5,18,19,20]. To leading order the divergence in the longitudinal susceptibility is proportional to h −1/2 and the equation of state has the leading behavior for x → −1. We note that the Goldstone-mode divergences cancel out and the equation of state is divergence free. This is seen order by order in the ǫ-expansion [21] and for fixed-dimension perturbation theory [22]. This is also observed nonperturbatively in the Monte Carlo data (see e.g. [10]). The corrections to the leading behavior are incorporated explicitly in the expression proposed by Wallace and Zia [13], which is inferred from the ǫ-expansion for the equation of state deduced in [21]. For d = 3 the expression corresponds to an expansion in powers of y 1/2 This form describes well the Monte Carlo data from the low-temperature region until around the critical temperature. The coefficient associated with the H −1/2 divergence of the susceptibility for H → 0 is c 2 . Note that the expression of c 2 derived in [13] increases with N, i.e. models with larger N should display stronger Goldstone-mode effects.
In References [9], [10] and [11] the Monte Carlo data for the equation of state have been fitted to the expression x(y) = x 1 (y) y n 0 y n 0 + y n + x 2 (y) where x 1 (y) is given in Eq. (8) above and corresponds to the first two terms of Eq. (6). (The parameters y 0 and n are chosen to ensure a good interpolation.) This interpolation of low-and high-temperature behaviors describes well the data, but of course it would be nicer to have a form valid in both regions, such as the parametric form introduced in [14]. This type of form is described in the next section. We then comment on the previous use of this parametrization and propose a new variant that is especially well suited for fits.

Parametric representation
Let us consider the polynomial parametric representation introduced in [14], in which one writes M, t and H in terms of the variables R and θ (see e.g. [1,6]) Here m(θ) and h(θ) are odd functions 4 of θ, regular at θ = 0 and θ = 1. This ensures that Griffiths's analyticity conditions are satisfied. The coexistence line is given by θ 0 , the smallest positive zero of h(θ). [From Eq. (12) it is clear that θ 0 must be greater than 1.] Without loss of generality, we may take m(θ) = θ. The equation of state then becomes f The relation between x and θ is shown schematically in Fig. 1 together with the respective ranges considered. Note that we must have θ 2 0 < 1/(1 − 2β) for the above mapping to be invertible [24]. For the O(4) case β ∼ 0.4 and we have roughly θ 2 0 < 5. With the parametrization (11)-(13), the singular part of the free energy F s can be written as where the function g(θ) is the solution of the first-order differential equation that is regular at θ = 1 [23]. (This follows from the relation H = ∂F s /∂M , where the derivative is taken at fixed t.) We note that the parametrization above was first used in perturbative studies of the equation of state for the Ising model [24]. We discuss below its application to the N-vector (Goldstone-mode) case.
In accordance with Eq. (7) the leading behavior for θ → θ 0 must be This combined with the requirement that h(θ) be an expansion in odd powers of θ suggests the general form This form is used in [6,12], in their "scheme B". They also define another scheme with a similar expression for m(θ). In both cases the differential equation becomes with coefficients a i depending on the exponent β, the root θ 0 , the coefficients c i and on the scheme considered. One can easily check that the solution of the differential equation that is regular at θ = 1 is given by Clearly, this solution is a function of the values of θ 0 and of the parameters c i , i = 1, . . . , n. In [6,12] the authors have considered the cases n = 1, 2, with parameters θ 0 , c i obtained from perturbation theory. We comment on their results in Section 5.2.
In the next subsection we introduce a more specific parametric expression for h(θ), as a combined expansion around θ = 0 and around θ = θ 0 , in order to isolate the two regions of the phase diagram.

Improved parametric form
We consider here a variant of the parametric function h(θ) above For consistency, the expression is normalized so that the contribution from the d j 's is equal to 1 at θ = 0. This normalization factor does not affect the equation of state, since h(θ) enters in f (x) only as a ratio. Let us note that this is still an odd function of θ and is equivalent (as an expansion in θ) to equation (19) with a rearrangement of terms. In particular, we may compare the series using only terms with c i coefficients (i.e. with all d j = 0) with the one using only d j coefficients. The relation between the two cases is given by Thus, a single coefficient c i corresponds to a sum of d j 's. Conversely, if we considered an expansion around θ ≈ θ 0 each d j would correspond to a sum of c i 's. Since the roles of the terms of the two series are different, considering both series may be important for truncated sums, such as the ones we use for the fits. The consideration of two types of coefficients (c i and d j ) is done for gaining better control over the description of the two distinct regions of the phase space, the lowand high-x regions. In fact, although the two sets of coefficients give rise to a similar expansion in powers of θ, the determination of the single coefficients (by fits to the numerical data) is more stable when each of the two regions is separately associated with a set of coefficients. More precisely, since we write the series as a product of two sums of terms corresponding respectively to an expansion around the high-x region (θ ≈ 0, coefficients c i ) and the low-x region (θ ≈ θ 0 , coefficients d j ), we can expect fits of the data for each of these two regions to be more sensitive to the corresponding set of coefficients, since the other set's main contribution will be a constant. An indication of this property can be seen from a "quick" expansion of the parametric form in powers of ǫ ≈ 0 for the two regions Note that we show only the leading order from each multiplicative contribution. It is interesting that the correction to the leading behavior is O(ǫ 2 ) in the first case and O(ǫ) in the second.
To be more precise, we can associate the various c i , d j coefficients with the coefficients in the separate expressions for the low-and high-x regimes used in [10], respectively equations (8) and (10) above. For the high-x behavior we expand f (x) around θ → 0 (corresponding to large x). We start by writing θ as a function of x and inverting Eq. (14) consistently in powers of θ (correspondingly in powers of x −1/β ). We then get where The equation of state then becomes This expression is of the form (6) and can be inverted and compared to Eq. (10), giving Note that the leading coefficient a contains θ 0 and sums of the coefficients c i and d j , whereas the expression for the next orders will contain isolated contributions from the c i 's (e.g. the coefficient c 1 in the expression for b) but not from the d j 's, which appear always as a sum. Analogously, for the low-x region we expand the expressions of θ, h(θ) around θ → θ 0 (corresponding to x → −1) and substitute the results into the expression for f (x). Defining we write x as a function of ǫ, invert this expression to get ǫ(x) and then obtain h(θ) in terms of x, as done above for the large-x case. The expressions are with The equation of state then becomes where This form may be inverted to give an expression of x as a series of powers of y 1/2 as in Eq. (8). We obtain the coefficients We see that in this case it is the d j coefficients that appear as single contributions, whereas the c i 's appear always as sums. Thus, the qualitative feature observed in Eqs. (24) and (25) is confirmed by a more careful expansion, i.e. the c's are more relevant for the high-x region and vice-versa. This will also be seen directly from the fits in Section 5.1. All calculations above were checked using Mathematica.

Amplitude ratios
Just like other critical properties of statistical systems (e.g. critical exponents), certain ratios of critical amplitudes are universal [25]. The amplitude ratios are taken as dimensionless combinations of critical amplitudes above and below T c for various quantities. For example, for the singular part of the specific heat one has where t ∝ (T − T c ). The ratio A + /A − is then universal. Similarly, by considering the behaviors of • the susceptibility along the critical isochore (t > 0, H = 0) • the magnetization along the critical isotherm (t = 0, H = 0) • the magnetization on the coexistence line (t < 0, H = 0) one may construct the universal ratios These and other universal ratios may be obtained directly from Monte Carlo simulations (as done e.g. in [26]) or indirectly from the equation of state. In the case of the Josephson parametrization discussed above, the universal amplitude ratios of quantities defined at zero momentum are given in terms of g(θ) by [24] A Note that in order to evaluate the function g(θ) one has to solve the differential equation (17), i.e. determine the solution (21). Our results for the above ratios are reported in Section 5.3.

The pseudo-critical line
Another important property that can be extracted from the equation of state is the characterization of the so-called pseudo-critical line, defined by the points where the susceptibility χ shows a (finite) peak for H = 0. This corresponds to the rounding of the divergence observed at the critical point, i.e. for H = 0 and T = T c . More precisely, one looks for a peak in the scaling function of the susceptibility, given by [28] where Clearly, at each fixed h the peak in χ is given by t p = z p h 1/βδ , and we have Thus, the behavior along the pseudo-critical line is determined by the universal constants z p , f M (z p ), f χ (z p ). Determining this line is important for systems where a study at H = 0 is not possible (and consequently the critical value T c is not known with accuracy), such as for the chiral transition of QCD at finite temperature. In fact, the knowledge of these universal constants allows an unambiguous normalization of QCD data (using the observed scaling along the pseudo-critical line), as done in [27]. The pseudo-critical line has been studied for O (2) and O(4) models in [28,29]. For the N = 4 case, it is found that the susceptibility peaks are given by z p = 1.33 (5). Since this value is close to the interpolating point of the equation of state in [10], it is very important to work with the smooth parametrization considered here, especially when using the derivative of f M (z) as in Eq. (59) below.
The expression for f χ (z) can be easily obtained from the equation of state, given by f M (z) or f (x). Using the original parametrization we obtain [28] f In terms of the parametric representation this gives Our results for f χ (z) and the determination of z p are shown in Section 5.3.

Results
The fits have been done using a conjugate-gradient minimization [30] of χ 2 -without considering the gradient of the function f (θ) -with a numerical inversion of Eq. (14) in order to find θ for any given value of x. For the critical exponents we used ν = 0.749(2) [31] and δ = 4.824(9) [32], implying the values β = 0.386(1), γ = 1.476(5) and the upper bound θ 2 0 ≤ 4.38 (5). We refer to these values as the first set of exponents. We note that the corresponding exponent δ from [31], 4.789(6), has slightly smaller error bars. However, we choose to use the one in [32] because it is obtained directly from (infinite-volume) simulations at nonzero magnetic field. These two exponents are not in agreement within error bars. We will also present below for comparison a few quantities obtained using the exponent δ from [31]. We refer to the resulting values as the second set of critical exponents.
The data for the magnetization are taken from [10]. In addition to the statistical errors, we have included errors due to the critical exponents, the critical temperature and the normalization constants H 0 and T 0 . These constants have been rederived using the first set of exponents above (with errors), yielding The errors reported in parentheses in all the tables below are Monte Carlo (MC) errors, obtained with 2000 MC iterations. In particular, in Section 5.1 we not only vary the y variable but also consider the uncertainties in the exponents γ and δ appearing in the fitting function [i.e. in Eqs. (14) and (15)]. The same is true for the errors reported in Section 5.2. In Section 5.3 the errors comprise the error bars in the input parameters and also the errors in the critical exponents. (7) 40.1 1.09(2) -1.18 (4) 27.3 1.07 (1) 3.3(1) 20.2 2.69 (2) 154(4) -111(2) 18.7  (7) 31.5 1.614(7) 0.58 (3) 19.

Fits
As a first step, we tried to fit the data separately in the high-and low-x regimes using only c i or only d j parameters, in order to confirm that the c i 's are more important at high x and the d j 's at low x, as suggested in Section 2.2. As one can see from Tables  1 and 2, this is indeed the case. At high x the fits using c i parameters work better that the fits using d j parameters, as can be seen in the case with one parameter plus θ 2 0 . When using two parameters plus θ 2 0 the values of χ 2 /d.o.f. obtained in the two cases coincide, but in the case with d 1 and d 2 one obtains the unphysical value θ 2 0 ≈ 7. Moreover, if one tries to do a fit using θ 2 0 , c 1 , c 2 and d 1 the fit is not better than the one reported in the third row of Table 1 and the value of d 1 is very close to 0. In the low-x region the fits using d j parameters work much better than the corresponding fits using c i parameters. Again, if one tries a fit using θ 2 0 , c 1 , d 1 and d 2 the result is not better than the one reported in the last row of Table 2. Thus, we see clearly that the 2.17(4) 0.9(1) -0.62 (7)    coefficients c i and d j are more relevant respectively at high and low x, as suggested in Section 2.2. As a second step, we checked that the fit of all the data using only the parameters c i does not work very well (see Table 3). In particular, even with four parameters c i one cannot get a large improvement in the value of χ 2 /d.o.f., compared to the case with only the parameter c 1 . The situation is slightly better when considering only d j parameters (see Table 4). Notice, however, that we have only a few data points with very large x and that the low-x expression used in [10] describes well the data up to x ≈ 2.
Finally, fits of all the data with both c i and d j parameters (see Table 5) work very well, giving a value of χ 2 /d.o.f. about a factor two smaller than the best result obtained in Table 3 (see last row). In particular, we find it interesting that for the three fits considered we get (within errors) the same value for θ 2 0 . By averaging over the three results we find θ 2 0 = 2.16 (2) .
In Fig. 2 we show a plot of the data together with the curve corresponding to the case on the second row of Table 5. We have also tried fits with 6 parameters, without significant f. may be related to remaining systematic effects in the data. This is especially true in the low-temperature regime, where the finite-size effects are very strong due to the effect of Goldstone-mode-induced singularities. It would be interesting to test our parametrization using the higher-precision data recently produced in [32]. In any case, if we use the second set of exponents above (i.e. with δ from [31]) the values of χ 2 /d.o.f. are significantly worse, as can be seen in Table 6. We note that, in order to consider this second set of exponents, we have reevaluated the normalization constants H 0 and T 0 , the values of x and y and the data errors for this case.

Comparison with other parametrizations
We now compare our results with previous expressions for the O(4) equation of state. In Reference [12], the scheme B considered by the authors corresponds to all d j = 0 and only θ 0 , c 1 nonzero. Their values for these coefficients are θ 2 0 = 2.4(2) , c 1 = 0.065 (30) .
Note that their value of θ 2 0 is consistent with ours within error bars. Using these coefficients as a "fit" of the data (considering the first set of critical exponents above), one obtains a χ 2 /d.o.f. of 268. We can also use our second set of data to evaluate χ 2 /d.o.f., but this yields the value 688.

Universal quantities
As discussed in Sections 3 and 4, we use the fits obtained above to evaluate several interesting universal quantities, such as critical amplitude ratios and the characterization of the pseudo-critical line in the phase diagram. In Table 7 we show the results obtained for the ratios A + /A − , R c , R χ using our preferred fits (reported in Table 5). The three fits give consistent results within error bars. Averaging over the three cases yields These values are in agreement with the ones reported in [12, Table 3]. (Note, however, that our values take into account the errors due to the uncertainty in the critical exponents.) We also show, in Table 8, the same quantities using the fits for our second set 1.7(2) 0.26 (2) Table 7: Results for the universal amplitude ratios using the fits reported in Table 5.
of data (from Table 6). These results show a little more fluctuation, but are essentially in agreement with the ones in Eq. (67) above. Note that the ratio R χ can also be evaluated directly [9] from the coefficient a in the interpolated form, given in Eq. (65). In this case we get R χ = a γ = 1.105 (15), in agreement with our result above and with Ref. [32].
We now turn to the numerical characterization of the pseudo-critical line (see Section 4). Using Eqs. (60) and (61), we draw the parametric plot of the scaling function for the susceptibility versus z (see Fig. 3). The peak corresponds to the pseudo-critical line and can be determined numerically from the two equations by varying θ. The peak coordinates thus obtained are reported in Table 9, where we used our preferred fits. The values are consistent within errors, yielding θ p = 0.587(2) , z p = 1.29(1) , f χ (z p ) = 0.341 (1) .
The results are in agreement with previous determinations of z p and f χ (z p ), made in Refs. [28] and [12], but our error for z p is much smaller. In Table 10 we present these quantities in the case of our second set of data. Again, the determinations are in agreement.

Conclusions
We have introduced an improved parametric form for the description of the equation of state of 3d O(N) models. This form is based on the parametrization used perturbatively  Table 5.
in [24] for the Ising model, but takes into account terms associated with the effects of Goldstone-mode fluctuations. Such effects are present in O(N) models along the coexistence line, i.e. at low temperatures and small magnetic field (or equivalently, at low values of the variable x). These new terms are included by means of the d j coefficients, associated with an expansion around the coexistence line. (The d j 's are considered in addition to the usual c i coefficients, related to the high-temperature/high-x behavior.) We show that the new parametric form indeed provides a better fit to the numerical data as compared to previous parametrizations. In particular, the consideration of the d j coefficients is essential for a good description of the Monte Carlo data in the whole range of values of x. Also, we were able to verify clearly the different roles played by c i and d j parameters in the high-and low-x regions. We note that -in the case where all d j = 0 -our parametrization is equivalent to the scheme B discussed in [6], used perturbatively by the authors for general O(N) models. We find that our value of θ 0 is consistent with their perturbative determination for the O(4) case, presented in [12]. However, we do not confirm their conjecture that fit (c i + d j ) 0.605 (2) Table 10: Results for θ p , z p and f χ (z) using fits reported in Table 6.
the c i 's get smaller with increasing i. We also stress that, in addition to providing a better fit to the numerical data, the expression considered is a continuous function, needing no interpolation between the two x regions. This is particularly useful for the determination of the pseudo-critical line, since the interpolating form introduced in [10] is unstable precisely in this region. As a result, our determination of z p is very precise in comparison to the previous estimates from the interpolated form and the perturbative equation of state.