Reliability of Taylor expansions in QCD

We investigate the reliability of the Taylor expansion method in QCD with isospin chemical potentials using lattice simulations. By comparing the expansion of the number density to direct results, the range of validity of the leading- and next-to-leading order expansions is determined. We also elaborate on the convergence properties of the Taylor series by comparing the leading estimate for the radius of convergence to the position of the nearest singularity, i.e. the onset of pion condensation. Our results provide a handle for quantifying the uncertainties of Taylor expansions in baryon chemical potentials.


I. INTRODUCTION
The thermodynamic properties of QCD at finite temperature and density are in the focus of current research in theoretical and experimental physics and are of fundamental relevance for the structure of compact stars and the evolution of the early universe. On the theoretical side, simulations of lattice QCD are the preferred nonperturbative tool to investigate the properties of QCD at strong coupling. Monte-Carlo simulations of lattice QCD, however, are hindered for nonzero baryon density n B due to the well known complex action problem (for recent reviews see [1,2]). Despite a number of proposals for methods to potentially overcome this problem and to enable direct simulations in this regime (for reviews see Ref. [3][4][5], for instance) there is currently no method which can provide reliable results in the interesting region with T T pc at physical quark masses. Here T pc is the crossover, or pseudo-critical, temperature associated with effective chiral symmetry restoration.
Nonzero-density QCD is studied most conveniently in the grand canonical ensemble, where the baryon density is traded for the associated chemical potential µ B . One approach to circumvent the complex action problem is based on a Taylor expansion of observables in powers µ B at zero chemical potential. Following the pioneering works [6][7][8], today the Taylor expansion method is the most established approach [9,10] to investigate regions of the finite-density phase diagram relevant for heavyion phenomenology. However, since the expansion can only be carried out to a finite order n (currently typically n ≤ 8), the region of reliability of the series is a priori unknown, leaving systematic uncertainties due to higher orders difficult to estimate.
Another piece of information encoded in the Taylor expansion coefficients is the potential existence of a singularity in the complex µ B -plane -for example a phase transition at real critical chemical potential µ B,c . Due to the non-analyticity at µ B,c , the phase transition cannot be described by a series expansion in one of the adjacent phases. In turn, this shows up as a finite radius of convergence for the series expansion of observables [11]. This method has been applied extensively in QCD to probe the presence of a possible second order critical endpoint in the µ B − T plane, see, e.g., Refs. [12][13][14][15].
A similar expansion can also be applied in the case of QCD at finite isospin chemical potential µ I [12,16], which is also realized in the aforementioned physical systems. One advantage of QCD with a pure isospin chemical potential (i.e. µ I = 0 but µ B = 0) is that the complex action problem is absent and the theory can be simulated with standard Monte-Carlo methods [17]. Consequently, QCD at pure isospin chemical potential can serve as a realistic test system to investigate the range of applicability of the Taylor expansion method.
After the initial lattice studies of QCD at µ I > 0 and µ B = 0 at finite lattice spacings and unphysical pion masses [17][18][19][20][21][22][23], we have recently determined its continuum phase diagram [24]. It has a rich structure: besides the chirally broken and restored regions at low chemical potential, it exhibits a Bose-Einstein condensed (BEC) phase of charged pions beyond a critical chemical potential µ I,c (T ). According to our findings, the boundary of the BEC phase is at µ I,c ≈ m π /2 for temperatures up to about 150 MeV. This is followed by a pronounced turn and a saturation at around T ≈ 160 MeV for chemical potentials µ I ≤ 120 MeV. The appearance of the BEC phase is accompanied by the spontaneous breaking of the residual U τ3 (1) symmetry, remaining from the chiral SU V (2) symmetry group at finite µ I . Consequently, the phase transition to the BEC phase is expected to be of second order in the O(2) universality class [25], which is consistent with the finite volume-dependence and the critical scaling of the lattice results [24].
In this letter, we extend the simulations of [24] to test the performance of Taylor expansion in µ I . In particular, we investigate the applicability of the Taylor expansion method for a broad range of temperatures and study its capability to determine the BEC phase boundary via estimates of the radius of convergence.

II. LATTICE SETUP
We employ the tree-level Symanzik improved gluon action and 2+1 flavors of rooted staggered quarks with two-levels of stout smearing at physical quark masses, following the line of constant physics from [26]. The continuum limit is approached using lattice ensembles with N t = 6, 8, 10 and 12, corresponding to lattice spacings of arXiv:1810.11045v2 [hep-lat] 1 Nov 2018 a = 0.20, 0.15, 0.12 and 0.10 fm around the zero-density crossover temperature T pc (µ I = 0). To enable the observation of the spontaneous breaking of the U τ3 (1) symmetry in finite volumes and to regulate the theory in the infrared, we introduce a pionic source λ in the fermion matrix for the light quark masses. This source term leads to an unphysical explicit breaking of the U τ3 (1) symmetry and physical results are obtained in the limit λ → 0. For a more detailed discussion see [24].
Our main observable is the isospin density which is free of ultraviolet divergences and, thus, does not require renormalization. Its computation in terms of lattice operators is discussed in detail in App. A. The most difficult task for a reliable computation of n I is the extrapolation in λ. Similarly to our experience with other observables [24], the λ-dependence of n I is very pronounced, so that a naive extrapolation cannot be performed in a controlled manner. In Ref. [24] we introduced an improvement program for the λ-extrapolations, using the singular values of the massive Dirac operator. Similar improvements can be applied to n I as well, and we discuss the details in App. A. The dependence on λ is reduced substantially, allowing for fits of the data to a constant or a linear function in λ 2 (note that n I is an even function of λ due to U τ3 (1) symmetry). The Taylor expansion for the isospin density with respect to µ I /T is given by where c 2 and c 4 are the associated Taylor coefficients. In the following we will consider the leading order n I LO (including c 2 ) and the next-to-leading order n I NLO (including c 2 and c 4 ) series. For our action and temporal extents N t , the Taylor expansion coefficients have been computed in Ref. [16], albeit at different temperatures and, in some cases, on slightly different volumes. To arrive at the temperatures used in our study, we have performed a cubic spline interpolation of the associated results. In addition, we found the volume dependence of the Taylor coefficients to be sufficiently small, so that the effects due to the slight differences in volume are negligible. The details of the interpolations and the study of volume effects are provided in App. B.

III. TESTING THE TAYLOR EXPANSION AGAINST DIRECT RESULTS
To perform a detailed comparison between the direct results for n I and the Taylor expansion in a wide range of µ I , we have extended our existing results [24] with new data up to µ I 325 MeV. This value is still sufficiently far away from the saturation region for N t ≥ 8, ensuring that lattice artifacts remain under control (for all the results shown below, the chemical potential in lattice units satisfies µ I a < 0.3). The comparison between the direct data for n I and the results from the Taylor expansion is shown in Fig. 1 for N t = 6. For the lower temperature, T = 124 MeV, the data reaches the BEC phase boundary at µ I,c ≈ m π /2. Up to this point the data shows remarkable agreement with both the LO and the NLO Taylor expansion. Slight differences between the data and the LO expansion become apparent close to µ I,c (see also Fig. 7 in App. B at T = 113 MeV). The lattice data starts to deviate from the Taylor expansion curve for µ I > µ I,c . This is certainly expected, since the Taylor expansion cannot capture the change of dynamics at the transition. In contrast, the results for T = 176 MeV are above the BEC phase boundary. In this regime the agreement with Taylor expansion persists up to larger values of µ I and starting at around µ I /m π ≈ 0.6 one can clearly see that the data favors NLO Taylor expansion over the LO expansion. Furthermore, Taylor expansion at NLO fails to describe n I within the current uncertainties at around µ I /m π ≈ 1.6, where higher orders become important for this temperature. Comparing the behavior of the data points between the two temperatures reveals another characteristic of the Taylor expansion. For low temperatures, the NLO expansion tends to underestimate the results for n I , while it overestimates them for higher temperatures. Due to continuity there will be a region, where the agreement between the NLO expansion and the full result is almost perfect. Outside the BEC phase, where the expansion converges to the correct result, this region is related to the suppression of higher order terms, most dominantly of c 6 (T ), which is indeed expected to cross zero as T increases. Inside the pion condensation phase this agreement is merely accidental and does not reveal any information about the interior of the BEC region.
To quantify the regions in parameter space where Taylor expansion at a given order (LO or NLO) starts to become unreliable, we look at curves in parameter space with constant difference between the full results and the Taylor expansion. In the following we focus on the high temperature region, T 150 MeV, to be able to draw conclusions about the applicability range of the Taylor method in the absence of the BEC phase transition.
The contour lines are determined using a twodimensional spline fit to ∆, where the nodepoints have been generated via a Monte-Carlo analysis (for a description of our fit strategy, see Ref. [27]). In the spline fit we include the constraint that ∆ = 0 for µ I = 0. Note that we expect a rapid change of the data for ∆ at the BEC phase boundary in the thermodynamic limit.
Continuum extrapolation for the contour line ∆ NLO /m 3 π = 0.61. The yellow curve corresponds to the continuum extrapolation and the points show (every third of) the results from the individual lattices that were included in the fit. We have slightly shifted the data horizontally to enhance visibility.
Our N t = 8 results for the contour lines are shown in Fig. 2 for various values of ∆. The figure also includes the BEC phase boundary, which we extended to higher values of µ I compared to Ref. [24], see App. C for details. Most of the contour lines have positive slopes, indicating the general tendency that the expansion performs better and better as the temperature increases. This is partly due to the fact that the actual dimensionless expansion parameter is µ I /T , cf. Eq. (2) -however, the contour lines differ from the simple µ/T = const. lines considerably (see below). The exception is the contour line with ∆ LO /m 3 π = 0.12, which is roughly insensitive to the temperature. In addition, the results clearly reflect that the NLO expansion has a broader reliability range than the LO one, with contour lines shifted to considerably higher values of µ I .
Eventually we aim at investigating the range of applicability of the Taylor expansion in the continuum. To this end we perform a continuum extrapolation of the contour lines of ∆, using a parameterization in terms of polynomials in (T −T 0 ) with lattice spacing dependent coefficients, setting T 0 = 140 MeV. In the continuum extrapolation we focus on ∆ NLO , for which the contours are well described by a second order polynomial for T ≥ 161 MeV. In all of the cases the N t = 6 results were found to be outside of the scaling region and have thus been excluded from the fit. One of these extrapolations is visualized in Fig. 3 for ∆ NLO /m 3 π = 0.61. Finally, the contour lines of ∆ NLO in the continuum limit are plotted in Fig. 4, this time versus µ I /T . Once more, the gray area indicates the BEC phase in the continuum with the updated phase boundary from App. C. As discussed above, the naive expectation for the contours outside of the BEC phase are lines with µ/T = const., i.e. vertical lines in this plot. While this is approximately the case for large ∆ NLO , the contours with small values of ∆ NLO show clear deviations from this expectation with the tendency to shift to larger values of µ I /T with increasing temperature.

IV. TESTING THE RADIUS OF CONVERGENCE
As mentioned in the introduction, the radius of convergence of the Taylor expansion has been used extensively in the literature to extract information on the possible phase transitions of the theory for µ B > 0. The current setup with µ I > 0 is ideal to test the performance of this method in QCD, since the phase diagram features a second order phase transition [24,25] comparably close to the µ I = 0 axis. We have already seen the breakdown of the expansion close to the phase boundary (cf. Fig. 1). We will now test whether the leading estimator for the radius of convergence also indicates the presence of this phase boundary.
A possible definition for the radius of convergence r for the Taylor series of n I from Eq. (2) is given by but note that for a general singularity in the complex µ Iplane this limit is not guaranteed to exist (see Ref. [28] for a counter-example). Here c n are the coefficients of the expansion of the pressure defined in Sec. II. Note that while the same radius of convergence r is encoded in the Taylor series of other observables, the estimators at finite n can be quite different. In particular, comparing the series for the pressure p, the density n I and the susceptibility χ I = ∂ n I /∂µ I gives r n (χ I ) = n − 1 n + 1 r n (n I ) = n(n − 1) (n + 2)(n + 1) r n (p) . The leading-order estimators for the radius of convergence using the series for different observables on our Nt = 8 ensembles. The results are compared to the boundary of the BEC phase (gray area) and to the contours of ∆ NLO (colored bands). These indeed agree in the limit n → ∞, but differ at finite n.
With two coefficients at hand, only a single estimator can be constructed for r and we cannot investigate the n → ∞ limit systematically. The estimators for r obtained from the different observables are shown in Fig. 5 for N t = 8. For comparison we also included the boundary of the BEC phase and the contour lines of ∆ NLO in the figure. Both the estimators r n and the n-th order contour lines are expected to fall on top of the phase boundary in the limit n → ∞ -as long as the BEC onset is the singularity closest to µ I = 0.
The results for r 2 (χ I ) are observed to lie surprisingly close to the phase boundary for low temperatures, while r 2 (n I ) and r 2 (p) significantly overestimate the radius of convergence. While such a perfect agreement for r 2 (χ I ) is likely accidental, similar tendencies were also found in an NJL-type model [29], and in toy models of QCD with imaginary chemical potentials [14], suggesting that the series of the susceptibility gives estimators with the fastest convergence rate. The estimators also reveal a considerable change of slope around 150−160 MeV, close to the upper boundary of the BEC phase and in agreement with the qualitative trend that the two curves will agree in the limit n → ∞. Finally, a qualitative agreement is also observed between the behavior of r 2 and the contours lines of ∆ NLO .
Up to now we have ignored a subtle issue regarding the estimators of the radius of convergence in finite volumes. In a finite volume V , phase transitions are smoothed out, but the partition function has Lee-Yang zeroes at complex µ I that approach the real axis as V → ∞ [30], according to criticality [31]. How this is connected to the V -dependence of the estimators r n is highly non-trivial. We find that the leading order Taylor coefficients depend only mildly on the volume. This is discussed in App. B together with the finite size effects of the full data.

V. CONCLUSIONS
In this letter we have presented a detailed comparison between Taylor expansion and full simulations of QCD at nonzero isospin chemical potential µ I . This is a theory with a second-order phase transition between the normal phase and a phase with Bose-Einstein condensation of charged pions, enabling us to observe the breakdown of the Taylor expansion at the critical chemical potential. Up to the boundary of the BEC phase the full data for the isospin density n I is well described by Taylor expansion -both by the leading-(LO) and the next-to-leading order (NLO) series.
To test the reliability of the expansion outside the BEC phase, we extended our lattice ensembles generated in [24] with simulations at higher chemical potentials. The update for the BEC phase boundary up to µ I ≈ 325 MeV is presented in App. C. To quantify the performance of the LO and NLO expansions in this region, we introduced the deviation ∆ between the full and the Taylor-expanded results, see Eq. (3). The contour lines of ∆ LO and ∆ NLO are shown in Fig. 2 for our N t = 8 ensembles. Taking the ∆ = 0.2 m 3 π contour (which corresponds to deviations of 3−6% in the isospin density) as an indicator, the LO expansion is reliable up to µ I /m π ≈ 0.5 to 0.6, while the NLO series performs reasonably well up to µ I /m π ≈ 1.0 to 1.5. The continuum extrapolation of the contour lines is visualized in Fig. 4. Contrary to the naive expectation that the contours lie along lines of constant µ/T outside of the BEC phase, we find considerable deviations from this behavior, with a tendency towards larger values of µ/T with increasing temperature.
We have also compared the estimator r 2 for the radius of convergence, obtained from the first two coefficients (c 2 and c 4 ) of the Taylor expansion, to the critical chemical potential µ I,c known from the full simulations. We find that both r 2 and µ I,c change similarly with temperature, signaling the expected agreement when higher order estimates for the radius of convergence are taken into account. Concerning possible different definitions of r, we find that the estimate obtained from the susceptibility χ I is closest to the phase boundary. Whether this remains true for higher-order estimators remains to be seen. Our findings demonstrate that already the leading estimators for the radius of convergence are sensitive to the phase transition. We emphasize, however, that a more detailed study including higher order estimates for r is clearly desired and mandatory to be able to draw definite conclusions. Our study can be used to quantify the uncertainties of Taylor expansions in baryon chemical potentials and to guide the interpretation of results obtained for the radius of convergence of such series. This will also be of relevance for comparison to low energy models of QCD.
where n N I (λ) is the operator from Eq. (A2) with a singular value sum truncated at n = N . This truncated difference δ N does not contribute in the λ → 0 limit, allowing us to write As indicated, the λ > 0 value of the operator is determined using stochastic estimators, while the correction term δ N is calculated in the spectral representation (A3). In addition to the improvement of the operator, we also employ the leading order reweighting discussed in section III.4 of Ref. [24]. This approximates the λ = 0 distribution of the lattice ensembles and brings the expectation value n I closer to its λ = 0 limit. As discussed in Ref. [24], the optimal (or minimal) value of N necessary to achieve sufficient improvement to obtain a controlled λ-extrapolation will in general depend on the operator. We find that the behavior of n I with N is similar to the one of the chiral condensate, so that N ≈ 100 is usually sufficient to obtain reasonably flat extrapolations. A particular example for the λ-extrapolations is provided in Fig. 6. The plot indicates the tremendous increase in reliability due to our improvement scheme, enabling a well controlled λ-extrapolation and precision results.

Appendix B: Interpolation of Taylor expansion coefficients and finite size effects
The Taylor coefficients are combinations of derivatives of the pressure, with respect to the isospin chemical potential at µ I = 0. These can be rewritten using the quark chemical potentials µ u and µ d -in particular, c 2 and c 4 are given by and where ∂ f stands for the derivative with respect to µ f /T . The results for the coefficients c 2 and c 4 from Ref. [16] for the 24 3 × 8 lattice are shown in the top panel of Fig. 7 together with a cubic spline interpolation.
To check finite size effects in the temperature range of interest we also include the results from the 32 3 × 8 ensemble from Ref. [16] in the top panel of Fig. 7. Apart from some visible but not significant effects for the c 4 coefficient at T 150 MeV finite size effects are absent. To lend further support to the statement that finite size effects are negligible, we show the results for n I obtained on 16 3 ×6 and 24 3 ×6 lattices in comparison to the Taylor expansion with coefficients from a 18 3 × 6 lattice in the bottom panel of Fig. 7. For µ I < µ I,c , we see that the lattice results agree within uncertainties and are in mutual agreement with the NLO expansion. Also evident are the expected strong finite size effects at and just above µ I,c , outside the applicability region of the expansion.
The main part of our study has been done on 24 3 × 6, 24 3 × 8, 28 3 × 10 and 36 3 × 12 lattices. In contrast the Taylor coefficients have been computed on 18 3 ×6, 24 3 ×8, 32 3 × 8, 32 3 × 10 and 32 3 × 12 lattices in Ref. [16]. Given the magnitude of finite size effects visible in Fig. 7, we can thus conclude that those effects are irrelevant within the present accuracy. . The gray band is the part of the continuum extrapolation from Ref. [24] which enters the fit for the purpose of matching to the phase boundary for µI < 120 MeV.
Appendix C: The pion condensation phase boundary for large µI For the present study we extended the range of chemical potentials compared to Ref. [24], enabling a determination of the BEC phase boundary for higher values of µ I . We have performed new temperature scans in the range 120 MeV < µ I < 325 MeV, allowing to locate the critical temperature T c (µ I ), where the pion condensate vanishes. The results for three lattice spacings are shown in Fig. 8. Compared to our previous results [24] we observe a slight increase in the critical temperature with all data points approximately lying along a constant line. To capture this behavior and to approach the continuum limit, we fit all points by the function d 1 +d 2 /µ 2 I with a 2dependent coefficients d 1 and d 2 . To smoothly connect to our previous result, the data of [24] for T < 161 MeV and 90 MeV ≤ µ I ≤ 120 MeV for each value of N t and in the continuum are also included in the fit. The resulting continuum extrapolation is also shown in Fig. 8.