The complex Langevin analysis of spontaneous symmetry breaking induced by complex fermion determinant

In many interesting physical systems, the determinant which appears from integrating out fermions becomes complex, and its phase plays a crucial role in the determination of the vacuum. An example of this is QCD at low temperature and high density, where various exotic fermion condensates are conjectured to form. Another example is the Euclidean version of the type IIB matrix model for 10d superstring theory, where spontaneous breaking of the SO(10) rotational symmetry down to SO(4) is expected to occur. When one applies the complex Langevin method to these systems, one encounters the singular-drift problem associated with the appearance of nearly zero eigenvalues of the Dirac operator. Here we propose to avoid this problem by deforming the action with a fermion bilinear term. The results for the original system are obtained by extrapolations with respect to the deformation parameter. We demonstrate the power of this approach by applying it to a simple matrix model, in which spontaneous symmetry breaking from SO(4) to SO(2) is expected to occur due to the phase of the complex fermion determinant. Unlike previous work based on a reweighting-type method, we are able to determine the true vacuum by calculating the order parameters, which agree with the prediction by the Gaussian expansion method.


Introduction
The sign problem is a notorious technical problem that occurs in applying Monte Carlo methods to a system with a complex action S. The importance sampling cannot be applied as it is since the integrand exp (−S) of the partition function cannot be regarded as a Boltzmann weight. If one uses the absolute value | exp (−S) | for generating configurations and treats the phase factor as a part of the observable, huge cancellations occur among configurations, and the required statistics grows exponentially with the system size. This problem occurs in various interesting systems in particle physics such as finite density QCD, gauge theories with a theta term or a Chern-Simons term, chiral gauge theories and supersymmetric theories.
The complex Langevin method (CLM) [1,2] is a promising approach to such complexaction systems, which may be regarded as an extension of the stochastic quantization based on the Langevin equation. The dynamical variables of the original system are naturally complexified, and the observables as well as the drift term are extended holomorphically by analytic continuation. It is known that the CLM works beautifully in highly nontrivial cases [3][4][5][6], while it gives simply wrong results in the other cases [7][8][9][10].
In the past several years, significant progress has been made in theoretical understanding of the method and the conditions for justifying the CLM. First it was realized that the probability distribution of the complexified dynamical variables has to fall off fast enough in the imaginary directions of the configuration space [11,12]. In order to satisfy this condition, a new technique called gauge cooling [13] was proposed. Using the gauge cooling, the CLM has been successfully applied to finite density QCD 3 either with heavy quarks [13] or at high temperature [20]. An explicit justification of the gauge cooling has been provided recently [21] extending the argument for justification of the CLM without gauge cooling [11,12].
It was known for some time that the CLM gives wrong results also when the determinant that appears from integrating out fermions takes values close to zero during the complex Langevin simulation. This was first realized in the Random Matrix Theory for finite density QCD [22,23] and confirmed also in effective Polyakov line models [24]. In these papers, it was speculated that the problem occurs due to the ambiguity associated with the branch cut in the logarithm of the complex fermion determinant, which appears in the effective action. On the other hand, ref. [25] pointed out that the singular drift term one obtains from the fermion determinant breaks holomorphy, which plays a crucial role in justifying the method.
A theoretical understanding of this problem and a possible cure have been given recently. First it was pointed out in ref. [26] that the branch cut cannot be the cause of the problem since the CLM can be formulated solely in terms of the weight w = exp(−S) without ever having to refer to the action S. Indeed it was found that a similar problem can occur when the action has pole singularities instead of logarithmic singularities. In the same paper, it was shown that the probability distribution of the complexified variables has to fall off fast enough near the singularities of the drift term, based on the argument for justification in ref. [11,12]. It was then proposed [27,28] that the gauge cooling can be used to satisfy this condition as well with an appropriate choice of the complexified gauge transformation. A test in the Random Matrix Theory shows that the gauge cooling indeed solves the singular-drift problem unless the quark mass becomes too small.
In ref. [29], the argument for justification with or without gauge cooling was revisited. In particular, it was pointed out that the expectation values of time-evolved observables, which play a crucial role in the argument, can be ill-defined. Taking this into account, it was shown that the CLM can be justified if the probability distribution of the drift term falls off exponentially or faster at large magnitude. This condition serves as a useful criterion, which tells us clearly whether the results obtained by the CLM are trustable or not.
In this paper, we focus on the singular-drift problem that occurs in a system with a complex fermion determinant. In many such systems, the phase of the fermion determinant is expected to play a crucial role in the determination of the vacuum. An example of this is finite density QCD at low temperature and high density, where various exotic fermion condensates are conjectured to form (See ref. [30], for instance.). Another example is the Euclidean version of the type IIB matrix model [31] for 10d superstring theory, where the SO(10) rotational symmetry is conjectured to be spontaneously broken [32][33][34][35]. When one applies the CLM to these systems, the singular-drift problem occurs due to the appearance of eigenvalues of the Dirac operator close to zero. We propose to avoid this problem by deforming the action with a fermion bilinear term and extrapolating its coefficient to zero. The fermion bilinear term should be chosen in such a way that the nearly zero eigenvalues of the Dirac operator are avoided and yet the vacuum of the system is minimally affected.
We test this idea in an SO(4)-symmetric matrix model with a Gaussian action and a complex fermion determinant, in which spontaneous breaking of SO(4) symmetry is expected to occur due to the phase of the determinant [36]. This model was studied previously by the Gaussian expansion method (GEM) [37] and the spontaneous breaking of the SO(4) symmetry down to SO(2) was suggested by comparing the free energy for the SO(2)-symmetric vacuum and the SO(3)-symmetric vacuum. The same model was studied also by Monte Carlo simulation using the factorization method 4 , and the order parameters obtained by the GEM were reproduced for both the SO(2)-symmetric vacuum and the SO(3)-symmetric vacuum [38,39]. However, the comparison of free energy for the two vacua suffered from too much uncertainty to make a definite conclusion on the true vacuum by this approach.
When one applies the CLM to this system, the singular-drift problem is actually severe because the fermionic part of the model is essentially an exactly "massless" system.
Indeed, it turns out that the gauge cooling proposed in refs. [27,28] is not sufficient to solve this problem in the case at hand. Following the idea described above, we therefore add a fermion bilinear term, which breaks the SO(4) symmetry minimally, down to SO (3). The results of the CLM show that the SO(3) symmetry of the deformed model is broken spontaneously to SO (2). Extrapolating the deformation parameter to zero, we find that the SO(4) symmetry of the original matrix model is broken spontaneously to SO (2) and that the order parameters thus obtained agree well with the prediction obtained by the GEM. We also try another type of the fermion bilinear term for the deformation and show that the final results obtained after the extrapolations remain the same, which supports the validity of our analysis. Note that we are able to determine the true vacuum directly without having to compare the free energy for each vacuum preserving different amount of rotational symmetry.
In order to probe the spontaneous symmetry breaking (SSB), we need to introduce an O(ε) symmetry breaking term in the action, on top of the deformation described above, and send ε to zero after taking the large-N limit. The singular-drift problem occurs at small ε even for the deformed model. Here, the criterion for correct convergence proposed recently [29] turns out to be useful since it tells us which data are free from the singular-drift problem and hence can be trusted. Indeed, we find that the data points in the reliable region can be fitted nicely by an expected asymptotic behavior, while the data points in the unreliable region deviate from the fitting curve. We hope that our strategy to overcome the singular-drift problem enables the application of the CLM to the type IIB matrix model and to finite density QCD at low temperature and high density.
The rest of this paper is organized as follows. In section 2, we define the SO(4)symmetric matrix model and briefly review the results obtained by the previous approaches. In section 3, we explain how we apply the CLM to the SO(4)-symmetric matrix model. In particular, we deform the action with a fermion bilinear term, which enables us to investigate the SSB without suffering from the singular-drift problem. In section 4, we present the results of our analysis. In particular, we extrapolate the deformation parameter to zero, and confirm that the SSB from SO(4) to SO (2) indeed occurs in this model. The order parameters thus obtained are in good agreement with the prediction of the GEM. Section 5 is devoted to a summary and discussions. In appendix A we give the details on how we determine the region of validity of the CLM, which is useful in making the ε → 0 extrapolations. In appendix B, we present the results obtained by deforming the action with another type of the fermion bilinear term, which turn out to be consistent with the ones obtained in section 4.

Brief review of the SO(4)-symmetric matrix model
The SO(4)-symmetric matrix model investigated in this paper is defined by the partition function [36] where the bosonic part and the fermionic part of the action is given, respectively, as Here we have introduced N × N Hermitian matrices X µ (µ = 1, . . . , 4), which are bosonic, and N f copies of N -dimensional column vectors ψ using the Pauli matrices σ i (i = 1, 2, 3). The model has an SO(4) symmetry, under which X µ transforms as a vector, whereas ψ α andψ α transform as Weyl spinors. Also, the model has an SU(N ) symmetry, under which the dynamical variables transform as Integrating out the fermionic variables for each f , one obtains the determinant of the Dirac operator which is complex in general. Thus, the partition function (2.1) can be rewritten as It was speculated that the SO(4) rotational symmetry of the model is spontaneously broken in the large-N limit with fixed r = N f /N > 0 due to the effect of the phase of the determinant [36]. In the phase-quenched model, which is defined by omitting the phase of the fermion determinant, the SSB was shown not to occur by Monte Carlo simulation [39]. We may therefore say that the SSB, if it really occurs, should be induced by the phase of the fermion determinant. Throughout this paper, we consider the r = 1 case, which corresponds to N f = N . In order to see the SSB, we introduce an SO(4)-breaking mass term in the action, where and define the order parameters for the SSB by the expectation values of where no sum over µ is taken. Due to the ordering (2.8), the expectation values obey at finite ε. Taking the large-N limit and then sending ε to zero afterwards, the expectation values λ µ (µ = 1, · · · , 4) may not take the same value. In that case, we can conclude that the SSB occurs. Explicit calculations based on the GEM were carried out assuming that the SO(4) symmetry is broken down either to SO(2) or to SO(3) [37]. For r = 1, the order parameters are given by The free energy was calculated in each vacuum, and the SO(2)-symmetric vacuum was found to have a lower value. Monte Carlo simulation of this model is difficult due to the sign problem caused by the complex fermion determinant. Among various reweighting-type methods, the factorization method [35] turned out to be particularly useful in the present case. Assuming that the SO(4) symmetry is spontaneously broken down either to SO(2) or to SO(3), the results of the GEM (2.11) and (2.12) were reproduced [38,39]. However, the calculation of the free energy difference had large uncertainties, and it was not possible to determine which vacuum is actually realized using this approach.

Application of the CLM to the SO(4)-symmetric matrix model
In this section, we explain how we apply the CLM to the SO(4)-symmetric matrix model (2.1). Including the symmetry breaking term (2.7), we can write the partition function as The drift term that appears in the Langevin equation is given by as a function of the Hermitian matrices X µ . Note that the second term in (3.3) is not Hermitian in general corresponding to the fact that the fermion determinant is complex. Thus, the application of the idea of stochastic quantization naturally leads us to complexifying the dynamical variables, which amounts to regarding the Hermitian matrices X µ as general complex matrices X µ . Accordingly, the definition of the drift term (3.3) is extended to general complex matrices X µ by analytic continuation. Then we consider the fictitious-time evolution of the general complex matrices X µ described by the discretized version of the complex Langevin equation where η µ (t) is an N × N Hermitian matrix generated with the probability proportional to e − 1 where t 0 represents the time required for thermalization and T should be large enough to achieve good statistics.
In order to justify the CLM, the probability distribution of the drift term (3.3) measured during the complex Langevin simulation should fall off exponentially or faster at large magnitude [29]. In the present model, this condition can be violated for two reasons. First, the first term in (3.3) can be large when the configuration X In order to avoid the first problem, we use the gauge cooling [13]. Note that the original theory (3.1) has the symmetry X µ → g X µ g −1 with g ∈ SU(N ), under which the drift term (3.3) transforms covariantly as v µ → g v µ g −1 and the observables (2.9) are invariant. Upon complexifying the variables, the symmetry property of the drift term and the observables enhances to X µ → g X µ g −1 with g ∈ SL(N, C). Using this fact, we can implement the gauge cooling procedure [13] in the Langevin process as where the transformation matrix g ∈ SL(N, C) is chosen appropriately as a function of the configuration X which measures the deviation of X µ from a Hermitian configuration, and choose the SL(N, C) transformation g in (3.6) in such a way that the norm is minimized. In practice, this is done by using the steepest descent method as follows. Let us consider an infinitesimal SL(N, C) transformation 9) where N × N traceless Hermitian matrices t a are the generators of SU(N ) normalized as tr (t a t b ) = δ ab . Since the norm (3.8) is invariant under SU(N ), we restrict the infinitesimal parameters ǫ a to be real. Under the infinitesimal transformation, we have Therefore, the change of the Hermiticity norm (3.8) becomes from which the gradient of the norm is obtained as Using this f a , we consider a finite SL(N, C) transformation where the real positive parameter α is chosen in such a way that the Hermiticity norm (3.8) is approximately minimized. We repeat this procedure until the norm (3.8) stops decreasing within certain accuracy. and the Langevin step-size is chosen as ∆t = 2.0 × 10 −4 unless stated otherwise. We find that the gauge cooling keeps the Hermiticity norm well under control.
Next we turn to the second problem, which is associated with the eigenvalues of the Dirac operator D close to zero. In Fig. 2, we plot the eigenvalue distribution of the Dirac operator obtained during the complex Langevin simulation for ε = 0.1 (Left) and ε = 0.5 (Right) with N = 32. We find that there are many eigenvalues close to zero for ε = 0.1, but not for ε = 0.5. This suggests that there is some critical ε, below which the results of the CLM cannot be trusted because of the singular-drift problem. It turns out that the extrapolation to ε = 0 is rather difficult in this situation.
In order to avoid this problem, we add a fermion bilinear term to the action (2.3). The partition function of the deformed model is defined as Note that the extra fermion bilinear term explicitly breaks the SO(4) symmetry of the original model (2.1). Here we choose the parameters M µ in such a way that the SO(4) We can then ask whether the SO(3) symmetry of this deformed model is spontaneously broken in the large-N limit. In Fig. 3, we plot the eigenvalue distribution of the Dirac operator (3.16) obtained during the complex Langevin simulation of the deformed model for ε = 0.1 (Left) and ε = 0.5 (Right) with m f = 1.0 and N = 32. We find that the distribution is shifted in the real direction. This is understandable since, at large m f , the eigenvalue distribution of the Dirac operator would be distributed around m f . As a result, the distribution avoids the singularity even for ε = 0.1 in contrast to the undeformed (m f = 0) case. Therefore, we can extrapolate ε to zero using data obtained with smaller ε for finite m f . Eventually, we extrapolate the deformation parameter m f to zero, and compare the results with the prediction (2.11) obtained by the GEM for the original model.

Results of our analysis
In this section, we present our results obtained by the CLM as described in the previous section. Let us recall that we have introduced an O(ε) mass term (2.7) for the bosonic matrices, which breaks the SO(4) symmetry explicitly. In order to probe the SSB, we need to take the large-N limit with fixed ε, and then make an extrapolation to ε = 0.
In Fig. 4, the expectation values λ µ ε,m f (µ = 1, 2, 3, 4) obtained for N = 16, 32, 48 with ε = 0.1 and m f = 1.0 are plotted against 1/N , where the data can be fitted nicely to straight lines. Thus we can extrapolate the expectation values to N = ∞ for each ε and m f . In what follows, we assume that the large-N limit is already taken in this way. Next we would like to make an extrapolation to ε = 0. For that purpose, it is convenient to consider the ratio This is motivated from the fact that the mass term (2.7) tends to make all the expectation values λ µ ε,m f smaller than the value to be obtained in the ε → 0 limit. By taking the ratio (4.1), the finite ε effects are canceled by the denominator, and the extrapolation to ε = 0 becomes easier. Since ε is a parameter in the action (2.7), the expectation values λ µ ε,m f and hence the ratios (4.1) can be expanded in a power series with respect to ε. By taking the ratios, the coefficients of higher order terms become smaller, and the truncation of the series becomes valid for a wider range of ε.
In Fig. 5, we plot the ratio (4.1) against ε for m f = 1.0 (Top-Left), 0.8 (Top-Right), 0.6 (Bottom-Left) and 0.4 (Bottom-Right). The data obtained at small ε suffer from the singular-drift problem, and hence cannot be trusted. Here the condition for justifying the CLM proposed recently in ref. [29] turns out to be useful since it enables us to determine the range of validity as we explain in appendix A. Taking this into account, we fit the data in Fig. 5 to the quadratic form using the fitting range given in Table  1, where we also present the extrapolated values. We find for each value of m f that ρ 1 (ε, m f ) and ρ 2 (ε, m f ) approach the same value in the ε → 0 limit, while the others approach smaller values. This implies that the SSB from SO(3) to SO(2) occurs in the deformed model.
In Fig. 6, we plot the extrapolated values lim ε→0 ρ µ (ε, m f ) obtained in this way against m 2 f . We find that our results within 0.4 ≤ m f ≤ 1.0 can be nicely fitted to the quadratic behavior, which is motivated by a power series expansion of the ex- which agree well with the results (2.11) obtained by the GEM. Here we emphasize that in the GEM, the true vacuum was determined by comparing the free energy obtained for the SO(2) vacuum and the SO(3) vacuum. In contrast, the CLM enables us to determine the true vacuum directly without having to compare the free energy for different vacua. As a further consistency check, we repeat the same analysis with a different choice of the deformation parameter M µ = (0, 0, m f , 0) in (3.16) instead of (3.17). We find that the results obtained after the extrapolation m f → 0 turn out to be consistent with the ones obtained above. See appendix B for the details.   Fig. 5 for the ε → 0 extrapolations is listed with the extrapolated values obtained by the fits.

Summary and discussion
In this paper, we have shown that the CLM can be successfully applied to a matrix model, in which the SSB of SO(4) is expected to occur due to the phase of the complex fermion determinant. The SSB does not occur if the phase is quenched, which implies that it is extremely hard to investigate this phenomenon by reweighting-based Monte Carlo methods. In the factorization method, for instance, one introduces a constraint with some parameters and extremizes the free energy with respect to these parameters. While this has been done successfully in refs. [38,39], the comparison of the free energy for the SO(2) and SO(3) vacua turns out to be subtle and a definite conclusion on the true vacuum was not reached. In contrast, we have shown by the CLM that the SSB from SO(4) down to SO(2) occurs as predicted by the GEM. For the success of the CLM, it was crucial to overcome the singular-drift problem associated with the appearance of nearly zero eigenvalues of the Dirac operator. The gauge cooling was used to suppress the excursions in the imaginary directions, but the singular-drift problem in the present case was too severe to be solved by the gauge cooling. This is understandable because the fermionic variables are exactly "massless" in the present case. Our strategy to overcome the singular-drift problem was to deform the Dirac operator in such a way that the singular-drift problem is avoided while maintaining the qualitative feature of the vacuum as much as possible. On top of this, we have to introduce an O(ε) symmetry breaking term to probe the SSB, which should be removed after taking the large-N limit. In making the ε → 0 extrapolations, the criterion for correct convergence proposed in ref. [29] turns out to be useful since it tells us the range of parameters for which the CLM is free from the singular-drift problem and the results are trustable. The order parameters obtained after extrapolating the deformation parameter to zero turn out to be consistent with the prediction by the GEM.
We have actually tried two types of deformation to avoid the singular-drift problem and confirmed that the extrapolated results agree with each other within fitting errors. While this confirms the validity of the extrapolations to some extent, we cannot exclude the possibility that something dramatic happens when the deformation parameter approaches zero. Let us recall, however, that the singular-drift problem can occur at some point in the parameter space even if the system itself does not undergo any dramatic change. For instance, in QCD at finite density, the singular-drift problem is anticipated to occur at the quark chemical potential µ m π /2, where m π is the pion mass, but the first order transition to the phase of nuclear matter occurs at µ ∼ m N /3, where m N is the nucleon mass. Nothing really happens in the wide parameter range 0 µ m N /3. This example clearly shows that the singular-drift problem has more to do with the methodology rather than the physics of the system to be investigated.
The CLM with the proposed strategy can be directly applied to the type IIB matrix model, which is conjectured to be a nonperturbative formulation of type IIB superstring theory in ten dimensions [31]. While the SO(10) symmetry of the model is expected to be spontaneously broken down to SO(4) for consistency with our 4d space-time, the GEM predicts that it is spontaneously broken down to SO(3) rather than SO(4) [40]. It would be interesting to investigate this issue using the CLM extending the present work.
We consider that the same strategy would be useful also in applying the CLM to finite density QCD at low temperature and high density, where various exotic condensates are speculated to form [30] due to the complex fermion determinant. In this case, one can deform the Dirac operator by switching on the corresponding fermion bilinear term without disturbing the vacuum significantly. Now that we have a useful criterion [29] for justifying the CLM, we can try possible deformations and see whether any of them allows us to extrapolate the deformation parameter to zero within the region of validity.

A How to determine the region of validity
In this appendix, we explain how to determine the region of validity of the CLM. When the symmetry breaking parameter ε becomes small, the singular-drift problem occurs and the results obtained by the CLM can no longer be trusted. In order to make ε → 0 extrapolations, it is important to determine the value of ε, below which the results become unreliable. Here we use the criterion based on the argument for justifying the CLM [29]. For that, we calculate the magnitude of the drift term for each configuration and obtain its probability distribution. If the tail of the distribution falls off exponentially or faster, we can trust the results obtained with those simulation parameters. We find that the finite step-size effects can modify the tail of the distribution significantly without changing the expectation values λ µ ε,m f . In order to make the plots in this section, we therefore have to decrease the step-size when it turns out to be necessary. Let us define the magnitude of the drift term by where v µ is the drift term defined by (3.3). Then, we define the probability distribution p (u) with the normalization ∞ 0 du p (u) = 1. In Fig. 7, we plot p (u) against u in the log scale for various ε with m f = 1.0 and N = 48. We find that p (u) falls off exponentially or faster for all the ε. Thus, we can trust the results obtained in this region.
In Fig. 8, we show a log-log plot (Left) and a semi log plot (Right) of the distribution p (u) for various ε with m f = 0.8 and N = 48. Since the drift term can become fairly large for ε = 0.1, we decrease the Langevin step-size to ∆t = 2.0× 10 −6 in order to probe the tail of the distribution correctly. We find that the distribution falls off exponentially or faster for ε ≥ 0.2, but a power-law tail develops for ε = 0.1. Therefore, we can trust the data for ε ≥ 0.2, but not the ones at ε = 0.1.   In Fig. 9, we show a log-log plot of p (u) for various ε with m f = 0.6 and N = 48. Here the drift term tends to become even larger than in the m f = 0.8 case, and we have to investigate the tail of the distribution more carefully. We therefore present the results obtained for two Langevin step-size, ∆t = 2.0 × 10 −4 (Left) and 2.0 × 10 −6 (Right). Indeed, we find that the behavior of the tail seems to change qualitatively by decreasing the step-size. In Fig. 10, we show a semi-log plot for ∆t = 2.0 × 10 −6 , which suggests that the tail of the distribution falls off exponentially for ε ≥ 0.3, but not for ε = 0.1. The result for ε = 0.2 is marginal. We may therefore trust the results for ε ≥ 0.3.
In Fig. 11, we show a log-log plot of p (u) for various ε with m f = 0.4 and N = 48. Here we have decreased the Langevin step-size to ∆t = 2.0 × 10 −8 , but the tail of the distribution still follows a power law for all values of ε within the region. However, the comparison of the two plots in Fig. 9 suggests a possibility that the step-size ∆t should be decreased further to see the behavior of the tail correctly. Thus for the m f = 0.4 case alone, we had to determine the lower end of the fitting range empirically from the plausibility of the fit to the quadratic behavior. Even if we omit the m f = 0.4 point in

B Results for another type of the fermion bilinear term
In this appendix, we present the results obtained by choosing the deformation parameters in (3.16) as instead of (3.17). Taking into account the ordering (2.10), we can preserve only an SO(2) symmetry with this choice. In Fig. 12, we plot the eigenvalue distribution of the Dirac operator (3.16) for ε = 0.1 (Left) and ε = 0.5 (Right) with m f = 0.6 and N = 32. We find that the distribution is separated in the imaginary direction. This is understandable since, at large m f , the eigenvalue distribution of the Dirac operator would be distributed around ±i m f . As a result, the singularity at the origin can be avoided for even smaller ε than in the case of (3.17). This enables us to extrapolate ε to zero using the data obtained in the large-N limit for finite m f .
In Fig. 13, we plot the ratios (4.1) obtained after taking the large-N limit against ε for m f = 0.6 (Top-Left), 0.5 (Top-Right), 0.4 (Middle-Left), 0.3 (Middle-Right) and 0.2 (Bottom). The data obtained for small ε cannot be trusted because of the singular-drift problem. We fit the data in Fig. 13 to the quadratic form using the fitting range given in Table 2, where we also present the extrapolated values. We find for each value of m f that ρ 1 (ε, m f ) and ρ 2 (ε, m f ) approach the same value in the ε → 0 limit, while the others approach smaller values.