Multi-Pion States in Lattice QCD and the Charged-Pion Condensate

The ground-state energies of systems containing up to twelve $\pi^+$'s in a spatial volume V ~ (2.5 fm)^3 are computed in dynamical, mixed-action lattice QCD at a lattice spacing of ~ 0.125 fm for four different values of the light quark masses. Clean signals are seen for each ground state, allowing for a precise extraction of both the $\pi^+\pi^+$ scattering length and $\pi^+\pi^+\pi^+$-interaction from a correlated analysis of systems containing different numbers of $\pi^+$'s. This extraction of the $\pi^+\pi^+$ scattering length is consistent with than that from the $\pi^+\pi^+$-system alone. The large number of systems studied here significantly strengthens the arguments presented in our earlier work and unambiguously demonstrates the presence of a low energy $\pi^+\pi^+\pi^+$-interaction. The equation of state of a $\pi^+$ gas is investigated using our numerical results and the density dependence of the isospin chemical potential for these systems agrees well with the theoretical expectations of leading order chiral perturbation theory. The chemical potential is found to receive a substantial contribution from the $\pi^+\pi^+\pi^+$-interaction at the lighter pion masses. An important technical aspect of this work is the demonstration of the necessity of performing propagator contractions in greater than double precision to extract the correct results.


I. INTRODUCTION
Multi-hadron systems, from the deuteron to heavy nuclei to neutron stars, represent a significant fraction of the universe that we observe, and for decades the phenomenological study of these systems defined the field of nuclear physics. Understanding how nuclei and nuclear interactions emerge from Quantum Chromodynamics (QCD), the underlying theory of the strong interaction, is now a central goal of modern sub-atomic physics. Since hadrons are bound-states of quarks and gluons, they do not arise at any finite order in perturbation theory and a description from QCD has proved elusive. The only known non-perturbative method that systematically implements QCD from first principles is its formulation on a discretized space-time, lattice QCD. While it is still not possible to directly calculate the properties of even the simplest nucleus from QCD, the tools and technology are gradually being put in place to make such calculations possible with lattice QCD in the near future. The impact of the successful realization of this goal cannot be overstated. For the first time, it would allow reliable calculations of strongly-interacting many-body processes that are not (or are only poorly) accessible experimentally. Important examples are hyperon-nucleon interactions that many play significant roles in the interior of neutron stars. Further, it would enable the exploration of how the properties of such systems depend upon the fundamental constants of nature, exposing the (possible) fine-tunings between the light-quark masses that give rise to the multiple fine-tunings observed in nuclear physics.
Our current understanding of nuclei requires a small but non-zero three-nucleon interaction. In the study of the structure of nuclei, we now have refined many-body techniques, such as Green's function Monte-Carlo (GFMC) [1] with which to calculate the ground states and excited states of light nuclei, with atomic number A < ∼ 14. Using modern nucleon-nucleon potentials that reproduce all scattering data below inelastic thresholds with χ 2 /dof ∼ 1, such as AV 18 [2], one fails, quite dramatically, to recover the structure of light nuclei. The inclusion of three-nucleon interactions greatly improves the predicted structure of nuclei, but at present, such interactions are difficult to constrain. At some point in the future, lattice QCD will be able to predict the interactions between multiple neutrons (and proton and mixed proton-neutron systems), bound or unbound in the same way it will be used to determine the two-body scattering parameters. A calculation of the three-neutron interaction, for instance, will be possible.
As a first step toward the study of nuclei and their interactions using lattice QCD, in this work we study systems composed of up to twelve π + 's. These multi-pion systems are conceptually and computationally the simplest multi-hadron systems that can be constructed. In addition to the two-body interactions, there are expected contributions from multi-body interactions. Such multi-body interactions are not forbidden by the symmetries of QCD, and are expected to be present with a magnitude that can be estimated with naive dimensional analysis (NDA) [3]. However, they are qualitatively (and obviously quantitatively) different from systems involving nucleons. The lowest-lying continuum state of multiple pions in a large volume is perturbatively close to each pion carrying zero-momentum, whereas the wavefunctions of systems of multiple nucleons are subject to the Pauli exclusion principle.
In a recent letter [4], we have reported results of the first many-pion calculation 1 in lattice QCD. We explored systems containing up to five pions, extracted the π + π + scattering length, and, for the first time, found indications of a non-zero renormalization group invariant (RGI) π + π + π + -interaction. Here we continue our investigations, presenting more detailed results of the previous studies and extending our work to systems containing up to twelve charged pions. We also use the recently derived expression for the ground-state energy of n-identical bosons in a finite volume at O(L −7 ) [5] in our analysis, one order beyond that at which our previous calculations [4] were analyzed.
Multi-pion systems are of interest in their own right as a strongly interacting boson gas at finite density and temperature. Such systems may be important for the late-time evolution of heavy ion collisions, such as those at RHIC and also in the interior of neutron stars [6,7]. During the last several years there have been a number of theoretical explorations of pionic systems at finite isospin chemical potential, with or without a finite baryon number chemical potential. Leading order (LO) chiral perturbation theory (χPT) has been used to study the vacuum realignment that takes place in the presence of an isospin chemical potential that exceeds the mass of the pion, leading to a charged pion condensate, and excitations about this realigned vacuum. One of the important results of this present work is the calculation of the isospin chemical potential as a function of the isospin density. Our results are in good agreement with the LO χPT result, and further, demonstrate the sizable contribution from multi-pion interactions even at moderate densities.
The structure of this paper is as follows. In Section II we review the theoretical expectations for the ground state energy of multi-pion systems at finite volume and discuss methods for extracting their interactions. In Section III we provide details of our lattice QCD measurements and analysis and in Sections IV and V we present the main results of our calculations. In Section VI we discuss the implications of our results for the equation of state of the pionic gas, its isospin chemical-potential and its pressure. Finally in Section VII we discuss the results in a global context and conclude. Certain technical details of contractions and numerical implementation are relegated to the Appendices.

II. MULTI-MESON ENERGIES : ISOLATING THE TWO-AND THREE-BODY INTERACTIONS
In recent works [5,8,9], the analytic volume dependence of the energy of n identical bosons in a periodic volume has been computed to O(L −7 ), extending the classic results of Bogoliubov [10] and Lee, Huang and Yang [11]. The resulting shift in energy of n particles of mass M due to their interactions is [5] where the parameter a is related to the scattering length 2 , a, and the effective range, r, by The geometric constants that enter into eq. (1) are and n C m = n!/m!/(n − m)!. The three-body contribution to the energy-shift given in eq. (1) is represented by the parameter η L 3 , which is a combination of the volume-dependent, renormalization group invariant quantity, η L 3 , and contributions from the two-body scattering length and effective range, where The quantity η 3 (µ) is the coefficient of the three-π + interaction that appears in the effective Hamiltonian density describing the system [5]. It is renormalization scale, µ, dependent. The quantity S is renormalization scheme dependent and we give its value in the minimal subtraction (MS) scheme, S MS = −185.12506.
2 In this work we use the Nuclear Physics sign convention for the scattering length, which is opposite to that of the Particle Physics sign convention. In this convention, the π + π + scattering length is positive.
For n = 2, the last two terms in eq. (1) vanish and the remaining terms constitute the small a/L expansion of the exact eigenvalue equation derived by Lüsher [12,13]: which is valid below the inelastic threshold. Below this threshold p cot δ(p) is the real part of the inverse scattering amplitude. The regulated three-dimensional sum is [14] S where the summation is over all triplets of integers j such that |j| < Λ and the limit Λ → ∞ is implicit. For n = 3, eq. (1) reproduces the shift of the ground state energy of three identical bosons that was recently calculated by Tan [9]. Naively, one might have expected to be able to determine the two-body effective range parameter, r, by calculating the energies of systems with different numbers of pions. However, given the scattering parameter redefinitions of eq. (2), this is clearly not possible. 3 Instead, the effective range will be extracted by calculating the energies of the excited states of, for example, the n = 2 pion system at finite volume, the lowest lying of which is perturbatively close to the state in which the pions both carry (back-to-back) one unit of lattice momentum, 2π/L.
It is useful to form various combinations of the many-body energies in order to isolate or eliminate the various contributions from two-body or three-body interactions. Further, important checks can be made regarding the convergence of the large volume expansion, as more particles are added into the volume. In particular, the combinations involving systems with n, m and two bodies vanishes at order O(L −7 ) for any n, m. 3 Writing −1/p cot δ = a = a + a 2 rp 2 /2 + ..., and evaluating it at the shifted energy of two particles in the volume at LO in the volume expansion, gives a = a + 2πa 3 r/L 3 + ....
The three-body interaction, η L 3 , can be eliminated by forming combinations of the manybody energies, allowing for various determinations of a. One such combination is which allows for the scattering length to be extracted at N 3 LO (omitting the last set of square brackets) and N 4 LO in the large volume expansion. Similarly, the three-body interaction can be isolated from combinations of the many-body energies. One such combination formed from the n-body and the two-body energies is which allows for η L 3 to extracted at two orders in the expansion. This, of course, can be straightforwardly generalized to other combinations of energies that may or may not include ∆E 2 .
The ground-state energies at finite volume, given in eq. (1), have been computed in nonrelativistic quantum mechanics, with relativistic effects added perturbatively [5]. They are given in terms of the scattering parameters, and the three-body interaction. For pionic systems in particular, it is natural to ask about the role of chiral perturbation theory in such a calculation. In χPT, the expansion parameters for this system (in the p-regime) are p/Λ χ and m π /Λ χ , where Λ χ ∼ 4πf π is the chiral symmetry breaking scale. The χPT expression for the ground state energy of n-π's in a finite volume (which remain to be determined) will be a dual expansion in 2π/(Λ χ L) and m π /Λ χ . As such, in order to extract the three-body interaction that first enters at O(L −6 ), the calculation in χPT would need to be performed at N 3 LO. While providing the complete quark-mass dependence at this order, such a calculation would involve the evaluation of three-loop diagrams, and contributions from the relevant N 3 LO counterterms. Organizing the perturbative expansion using nonrelativistic effective field theory, EFT(π /) (see e.g. Refs. [15,16]), greatly simplifies the result, but does so at the expense of ignorance of the quark mass-dependence of the EFT(π /) parameters without further calculation (such as the quark mass dependence of the scattering length).
It is important to notice that for the systems containing a large enough number of pions, the energy-shift of the ground state exceeds that required to pair-produce pions. The leading relativistic effects that are included in the analytic expression for the energy-shift include only the single particle relativistic kinematics and corrections to the two-pion scattering amplitude, they do not include contributions due to inelasticities, including pair-production. However, we conclude that such effects make a small contribution to the ground-state energy as we find no evidence for deviations from eq. (1). Nonetheless, this aspect of these systems must be explored further.

A. Lattice configurations and quark propagators
The results of the numerical computations presented in this paper were obtained using the mixed-action lattice QCD scheme developed by LHPC [17,18] and are based on the coarse MILC lattice configurations [19]. These lattices have a lattice spacing of b ∼ 0.125 fm, and a spatial extent of L ∼ 2.5 fm. They were generated using the asqtad-improved [20,21] staggered formulation of lattice fermions, taking the fourth root of the fermion determinant, and the one-loop, tadpole-improved Symanzik gauge action [22]. Herein, we assume that the "fourth-root trick" 4 recovers the correct continuum limit of QCD. These ensembles of configurations have a fixed (almost physical) strange quark mass while the degenerate light quarks were varied over a range of masses; see Table I and Refs. [38,39,40,41,42] for details.
Based on these configurations, valence quark propagators using the domain-wall (DW) formulation of the lattice fermion action [43,44,45,46,47] were computed from smeared sources on each gauge-field configuration. Hyper-cubic (HYP) smearing [24,48,49,50] was applied to the gauge links used in the domain-wall fermion action to improve chiral symmetry, and in calculating the quark propagators, Dirichlet boundary conditions were imposed to reduce the original temporal extent of 64 down to 32. This procedure is optimized for nucleon physics and indeed leads to minimal degradation of a nucleon signal, however it does limit the number of time slices available for fitting meson properties in which the ratio of signal to noise remains constant in time. Further details about the mixed-action scheme can be found in Refs. [38,51]. A summary of the lattice parameters and resources used in this work is given in Table I. In order to generate large statistics on the existing MILC configurations, multiple propagators from sources displaced both temporally and spatially on the lattice were computed.
In the continuum chiral limit the n f = 2 staggered action has an SU (8) L ⊗SU (8) R ⊗U (1) V chiral symmetry due to the four-fold taste degeneracy of each flavor, and each pion has 15  The mass and decay constant of the π + , and the energy-shift and scattering parameters in the π + π + system calculated previously [54]. The first uncertainties are statistical, the second uncertainties are systematic uncertainties due to fitting and the third uncertainty, when present, is a comprehensive systematic uncertainty [54]. l (I=2) ππ is the one-loop counterterm in χPT contributing to π + π + scattering, and δ is the phase-shift. Recall that we are using the nuclear physics convention for the sign of the scattering length. degenerate partners. At finite lattice spacing, this symmetry is broken and the taste multiplets are no longer degenerate, but have splittings that are O(α 2 b 2 ) for the asqtad staggered action. When determining the mass of the DW valence quarks there is an ambiguity due to the non-degeneracy of the 16 staggered bosons associated with each pion. One could choose to match to the taste-singlet meson or to any of the mesons that become degenerate in the continuum limit. The choice of tuning to the lightest taste of staggered meson mass, as opposed to one of the other tastes, provides for the "most chiral" domain-wall mesons and therefore reduces the uncertainty in extrapolating to the physical point. The mass splitting between the domain-wall mesons and the staggered taste-identity mesons, which characterizes the unitarity violations present in the calculation, is then given by [52,53] Simple properties of the π + have been computed to high precision on the ensembles of coarse MILC lattices that are used in this work, both with staggered valence quarks and domain-wall valence quarks. Further, the energies of the lowest-lying 2-π + states, which lead directly to the scattering lengths using Lüscher's method, have been computed relatively precisely [54]. The results of the previous mixed-action calculations on these lattice ensembles [54] are shown in Table II 5 .
The results of the present calculation are presented in lattice units (l.u.), or in terms of dimensionless quantities such as m π /f π which eliminates the requirement of scale setting. They are performed only at one lattice spacing, due to limited computer time, and as a result the continuum limit cannot be determined. Unlike the two meson system, for which mixed-action chiral perturbation theory (MAχPT) [55,56,57] has been used to include the leading order effects of the finite lattice spacing, MAχPT calculations have not yet been performed for the multi-π + systems, and therefore the leading lattice spacing artifacts in these calculations cannot be removed at present. The lattice spacing artifacts are assumed to be small, occurring at O(b 2 ), but a systematic study must be performed in the future.

B. Correlation functions
In this work we determine the π + π + and π + π + π + interactions from the ground-state energy of n < 13 π + 's (isospin stretched states). By working in the m u = m d limit and restricting the calculation to states of maximal isospin, only the simplest sets of propagator contractions are required to be performed (i.e. no disconnected diagrams) in order to form the correlation functions from which the ground-state energies are extracted.
Naively, there are (n!) 2 contractions (for large n this behaves as ∼ (2n + 1 3 )πe 2n(log n−1) ) contributing to the correlation function of n-π + 's, where π + (x, t) = u(x, t)γ 5 d(x, t). However, this correlation function can be written as 6 where and S(x, t; 0, 0) is a light-quark propagator. The object (block) Π is a 12 × 12 (4-spin and 3 color) bosonic time-dependent matrix, and η α is a twelve component Grassmann variable. Using leads to correlation functions 5 Until this point the two-body scattering length for a generic system has been denoted by a. For the π + π + system, we denote the scattering length by a (I=2) ππ . 6 We thank David Kaplan and Michael Endres for discussions on this topic. For a general approach to evaluating contractions involving a large number of fermions, see Ref. [58]. While correct, further simplifications are possible. Let us recall that for an arbitrary 12 × 12 matrix, A, where in the last line we identify the matrix A with Π. Further, Therefore, by equating terms of the same order in the expansion parameter λ in eq. (18) and eq. (19), one can recover the n-π + correlation functions in eq. (17). As an example, the contractions for the 3-π + system are Rewriting the contractions in terms of traces over the Π-blocks greatly reduces the required number of calculations, with the number of independent contributions to the correlation function equal to the partition, P (n), of n objects. An estimate of the number of operations that must be performed to generate the correlator for n-mesons is ∼ n(12 × 13 − 1) + n j=1 P (j), which for large n scales as ∼ 2n/3 using a classic result of Hardy and Ramanujan [59]. While for n = 12 there are ∼ 2.3×10 17 independent contractions that must be performed, this can be accomplished with ∼ 2 × 10 3 calculations to produce the ∼ 80 terms contributing to the contraction. Since, in this work, each contraction is performed with only a single quark propagator on each configuration, the Pauli-exclusion principle requires that the n ≥ 13 identical meson contractions vanish identically, e.g. C 13 (t) = 0 ∀ t, implying the λ 13 and higher terms in the expansion of eq. (19) vanish. Written in terms of contractions of propagators in flavor and color space, the n = 13 case of eq. (19), represents a generalized Cayley-Hamilton identity satisfied by all matrices of size less than 13 × 13. To perform calculations on systems containing more than twelve pions, additional propagators will be required. In order to calculate the n-π + correlation functions, particularly for n ≥ 8, it is necessary to use a numerical representation with precision greater than that of standard 64-bit machine precision. 7 This need arises because of the large products of propagators that must be computed and is not particular to the contractions studied here. In particular, the numerical issues impacting calculations of multi-pion systems that we have found in this work will also impact calculations of multi-nucleon systems.
Our implementation of these contractions uses arbitrary precision arithmetic based on the ARPREC library [60] which was extended for the particular operations needed here, matrix multiplications and traces. For the correlators studied here 64 decimal digit precision (approximately octupule precision) in internal operations is sufficient to give results accurate to sixteen digits. The additional overhead of using this numerical representation causes the high precision contraction code to run ∼10-50 times slower than a double-precision version but is only marginally dependent on the precision used 8 . For the n = 13 correlation function it is instructive to look at the dependence of the resulting correlator on the precision used in the computations. Since the correlation function must vanish identically for any input propagator, it is a very stringent test of the codes used herein. In fig. 2, the logarithm of the sum over time-slices of the absolute value of the correlation function as a function of the digits of precision used to perform the contractions on a representative configuration is shown. From extrapolating the results shown here, we conclude that the correlator is indeed identically zero.

D. Analysis
The correlation functions from which we extract the ground-state energy of the n-π + system are given in eq. (13), and, on a lattice with infinite extent in the time direction, behave as at large times. It is the difference between this energy, E n and n times the π + rest mass that is equated to the energy difference given in eq. (1), and which is extracted from the ratio of correlation functions where ∆E n is that of eq.(1). While there are a number of ways to extract the energy difference from the correlation function, perhaps the most visually pleasing one is to construct the effective energy difference function, defined to be ∆E eff.
In the limit of an infinite number of measurements, this function would tend to a constant equal to the ground-state energy splitting. Of course, for any real calculation, both the number of gauge fields and the number of propagators per gauge field are finite, and as such the object ∆E eff. n (t) consists of a central value and an associated uncertainty at each timeslice, t. Further, the temporal extent of the lattice is finite, giving rise to both forward and backward propagating hadrons. As such, ∆E eff. n (t) is constant (up to statistical fluctuations) only over a finite number of time-slices, in a region between where the forward propagating excited states have died-out sufficiently, and where the backward propagating states have not yet become significant. In the current situation where Dirichlet boundary conditions are imposed, the behavior of the correlator is modified by the "reflection" of forward and backward propagating states from the boundaries. As these reflections are poorly understood, the region close to the boundary is omitted in our analysis.
In our lattice calculations, multiple propagators and correlation functions are computed on each gauge configuration. These propagators and correlation functions are not statistically independent unless the sources are separated by many correlation lengths. To account for this, all the correlation functions for a fixed n computed on a given configuration are averaged (blocked) into one correlation function. The configurations are found to be statistically independent 9 , and the blocked correlators on each configuration form the basis of our statistical analysis.
In order to extract the energy difference ∆E n from ∆E eff. n (t), a fitting interval must be selected. This interval is chosen to be entirely contained in the region where the ∆E eff.
n (t) is consistent with a constant. Once the fitting interval has been selected a correlated χ 2 minimization is performed to extract the parameter ∆E n , defined in eq. (22). The covariance matrix that determines the correlated weightings of each of the values of ∆E eff. n (t) on any given time-slice is generated using a single-subtraction Jackknife procedure. 10 The central value of ∆E n is the value that minimizes the correlated χ 2 , and the standard statistical uncertainty is determined by the values of ∆E n for which χ 2 → χ 2 +1. The fitting systematic uncertainty associated with the fitting procedure is determined by varying each end of the fitting range by −2 ≤ ∆t ≤ +2, and refitting the energy-splitting.
To study the scattering length and three-body parameter, η L 3 , using eqs. (10) and (11), the appropriate ratios of the C n (t) correlators are used to define effective scattering length functions for each n (LO, NLO, N 2 LO) or each pair {n,m} (N 3 LO and N 4 LO) and effective η L 3 functions for each n. We then analyze these in the same manner as the energy differences above. This leads to multiple determinations of a (I=2) ππ and η L 3 . The effective functions defined by ζ 6,7 , eqs. (8) and (9), are studied similarly.
To make use of the full data set, we also perform a simultaneous, correlated fit of a (I=2) ππ and η L 3 to the effective masses for n = 2, 3, . . . , N max for N max = 3, . . . , 12.
In order to do this, fitting ranges, t (1), and where the covariance matrix that determines the correlated weightings of each contribution is generated using the Jackknife procedure. The standard (statistical) uncertainties on a (I=2) ππ and η L 3 are determined from the maximum and minimum values of each parameter dictated by the uncertainty-ellipse corresponding to their values for which χ 2 → χ 2 + 1. The systematic uncertainty associated with the fitting procedure is determined by repeatedly and randomly varying each end of the fitting range 9 This is tested by averaging over sets of neighboring configurations and performing analysis on the resulting blocked ensemble. For block sizes of 1, 4 and 12, no noticeable difference is seen. 10 As a check, we also performed a separate analysis using bootstrap re-sampling. The resulting energies and parameters were consistent, and, for simplicity, we focus on a single analysis in the main discussion.
for each correlator by −2 ≤ ∆t ≤ +2, refitting the parameters a (I=2) ππ and η L 3 , determining the complete range of values for each parameter associated with χ 2 → χ 2 + 1. 11 The systematic and statistical uncertainties are combined in quadrature in this work.

IV. LATTICE QCD RESULTS
Using the techniques discussed in the previous section, we now turn to analysis of the results of the lattice calculations. Our main aim is to extract the parameters a (I=2) ππ and η L 3 which we can do in a number of ways, either by forming particular combinations of energies, eqs. (10) and (11), or by a coupled analysis accounting for correlations among different n. The different methods give consistent results but we find that the most precise extraction is achieved using the latter method and consequently our final results are generated from this technique. Before we present these results, we first detail the simpler analysis using combinations of energies. As there are a large number of correlation functions that we study in this work, in some intermediate stages, we only display results for a single quark mass corresponding to m π = 291 MeV (in terms of uncertainties, this ensemble is neither the best nor the worst).
A. Multi-pion energies and energy differences A priori, it may seem surprising that the correlation function of twelve π + 's can be calculated at all. On the ensemble associated with the lightest pion mass, m π = 291 MeV, the 12-π + state has an energy of E 12 ∼ 3.5 GeV, while on the ensemble with m π = 591 MeV, the 12-π + state has an energy of E 12 ∼ 7.1 GeV. It is not immediately obvious that such a rapidly diminishing exponential can be cleanly measured but in these systems we find that it is possible. 12 As an example, for m π = 291 MeV, the effective energies (in lattice units) for systems with n = 1, 2, .., 12 π + 's are shown in fig. 3. Well-defined plateaus in the effective energy plots are seen for all systems, with the relative statistical uncertainty in the data almost constant as a function of n. As can be seen from the fits to the energies (statistical and systematic uncertainties are shown in quadrature), the precision with which the energy can be extracted is high, typically < 2%. The total relative uncertainties on the energy of the n-π + systems are shown for all data sets in fig. 4, with only a slight dependence of n apparent. An additional point is that it was not obviously the case that the Gaussian-smeared source for the light-quark, that is suitable for the single pion ground-state, would have sufficient overlap onto the multi-pion ground-state to produce useful correlation functions. However, it is clear that it did.
The energy differences which enter into eq. (1) and subsequent results can also be extracted cleanly, although with somewhat less precision than the individual energies. The effective energy difference plots, along with our fits to the energy differences are shown in fig. 5 (again for the m π = 291 MeV ensemble) while the relative uncertainties in our extrac- 11 An alternative systematic procedure of repeatedly and randomly choosing triplets of time-slices in each fit range ±1 time-slice and refitting is also used, giving qualitatively similar results. 12 In purely pionic systems there is no exponential degradation of the signal-to-noise ratio, unlike in most other hadronic systems [61]. tions are given in fig. 4. All effective energy splitting plots show behavior that is consistent with a single exponential (within statistical uncertainties) for a number of time slices. As discussed above, the region above t/b ∼ 16 is contaminated by reflections from the Dirichlet boundary at t/b = 22 and is discarded in our analysis.
B. π + π + scattering length Since the π + π + system is free from three-and higher-body hadronic interactions, it is the ideal place to extract the two body parameter, a (I=2) ππ . As is well known, this can be done without resorting to an expansion in a (I=2) For the n = 2 data, it is clear that NLO and higher extractions, fig. 7, yields the same a (I=2) ππ as the exact eigenvalue method of Lüscher. However, for the multi-pion systems, an n-dependent systematic deviation from the exact eigenvalue method of Lüscher is found at LO, NLO and N 2 LO. This is particularly clear for the lighter mass ensembles. In contrast, the extractions of a (I=2) ππ in which the π + π + π + -interaction (or at least, a term that behaves as n C 3 ) is eliminated (N 3 LO and N 4 LO) are in close agreement with the n = 2 exact eigenvalue result for all n. From this alone we conclude that the calculation suggests the presence of a significant π + π + π + -interaction.
We can further demonstrate the need for a π + π + π + -interaction by using eq. (1) to compute the energy shifts at O(L −7 ) using the value of a (I=2) ππ extracted from the π + π + -system using the exact eigenvalue method and setting η L 3 = 0. These can then be compared to the calculated effective energy differences as shown in fig. 8 for the m π = 291 MeV ensemble. The deviations between the predictions and the effective energy splittings are significant and grow with increasing n. Over the same sets of time-slices as used in the fits to the two-particle energy difference, the χ 2 /d.o.f. of such an ansatz is 8.62 and therefore η L 3 = 0 very poorly describes the results of the calculation. 13 Given the accuracy with which m π and f π have been extracted on these ensembles, the uncertainties in

C. Three body interaction
The π + π + π + -interaction can be explicitly constructed using eq. (11) for n = 3, . . . , 12. As an example, the effective η L 3 plots for the m π = 291 MeV ensemble are shown in fig. 9. A clear plateau inconsistent with zero is seen in most cases. Fig. 10 shows the n dependence of the extracted value of η L 3 for each quark mass. In general, the combined systematic and statistical uncertainty of the extractions decreases with increasing number of π + 's. This is not surprising given the combinatoric factors that appear in the expression for the energy shift.
D. Convergence: ζ 6,7 Before presenting the n-correlated analysis we briefly turn to the quantities ζ 6,7 defined in eqs. (8) and (9). A plateau in the corresponding effective-ζ 6,7 plots for a particular pair {n, m} at a value inconsistent with zero would signal the breakdown of the large-volume   Table III and their uncertainty ellipses are shown in fig. 15. For comparison with the simpler analysis above, the shaded regions in fig. 10 correspond to the values of η L 3 extracted from this correlated analysis and are seen to be consistent with the extractions made using eq. (11). Similarly, the extracted a (I=2) ππ agrees with that obtained from the exact eigenvalue method for n = 2 (and hence with all the N 3 LO extractions above).
In fig. 16, η L 3 is plotted versus m π /f π , in units of the estimate based upon NDA, η L,N DA 3 = 1/(m π f 4 π ), as discussed in Ref. [5,8]. While the three-body interaction is consistent with its NDA estimate for the lightest three pion masses, it is found to be consistent with zero at the  heaviest mass, m π = 591 MeV. It is desirable to reduce the uncertainties in this calculation to see, if in fact, the three-body interaction is decreasing with increasing pion mass. If this is found to be the case then a more detailed study in this high pion mass region is warranted. It is interesting to study how the larger n energy differences influence the extraction of the two parameters. To do so, we have performed a series of fits including only the energy differences up to a given n max . The resulting confidence regions for the parameters extracted from the m π = 291 MeV ensemble are shown in fig. 17 for n max = 3, 5, 7, 9, 11, 12. Clearly from the π + π + energysplitting (using the exact eigenvalue method) and setting η In the preceding Section the π + π + π + -interaction, η L 3 , has been extracted successfully from the lattice QCD calculations and is non-zero at the lighter quark masses used in this study. The renormalization group invariant, but volume dependent three-body interaction that arises at O(L −6 ), η L 3 , that receives perturbative corrections in the volume expansion to become η L 3 at O(L −7 ), is somewhat more problematic to extract from these particular lattice calculations due to the size of a (I=2) ππ in relation to the size of the lattice. There are two terms that must be calculated in order to determine η L 3 , one is additive and one is multiplicative, as given in eq. (4). Assuming that the effective range for π + π + scattering is not orders of magnitude larger than a (I=2) ππ , the additive term makes a very small contribution to η L 3 . To make this explicit, it is useful to rewrite eq. (4) as where The numerical values of β η 3 are shown in Table IV, and are all seen to be very small for r ∼ |a|. The multiplicative factor, α η 3 , is also shown in Table IV, and its values make clear that the perturbative expansion for the three-body interaction is converging extremely slowly at the lightest pion mass, with a ∼ 75% correction at O(L −7 ) to the O(L −6 ) result, and fails completely at both m π = 491 MeV and 591 MeV for the lattice volumes used in this work. Clearly larger lattice volumes will be required in order to determine the three- body interaction η L 3 with precision. 14 This is in stark contrast to the two-body scattering parameters which can be extracted to high precision from these same lattice volumes and the perturbative expansion in eq. (1) appears to be converging. So while it is possible to determine the "dressed" three-pion interaction, the uncertainty in the determination of the bare three-pion interaction is large, not because of the uncertainty in the lattice calculation, but due to the relatively large higher-order terms in the volume expansion when evaluated at this present lattice volume (a theoretical systematic uncertainty). In Table IV, the correction factors are also given for a lattice of spatial volume (3.5 fm) 3 for comparison. However, even in these large volumes the correction factors at one higher order in the large volume 14 We note that it is η L 3 that was extracted in Ref. [4]. expansion of ∼ 50% and clearly even larger volumes, L > ∼ 3.5 fm, are required in order to have the volume expansion of the three-pion interaction under perturbative control.
The quantity η 3 (µ) is a renormalization scheme dependent quantity that is independent of the volume, and as such is the quantity that most directly enters into the calculation of other many-body processes. It is easily extracted from η L 3 via eq. (5). However, given the present theoretical systematic uncertainties in η L 3 , we do not attempt a determination of η 3 (µ).

VI. THE EQUATION OF STATE AND THE ISOSPIN CHEMICAL POTENTIAL
The energy of the n-π + system as a function of volume and of the number of π + 's is given explicitly in eq. (1) in the large-volume expansion. From the equation of state, the isospin chemical potential is defined as which can be constructed analytical from eq. (1) or numerically from the results of the lattice calculation by using a simple finite difference approximation. The resulting ratio of this isospin-chemical potential to the pion mass for each of our ensembles is shown in fig. 18 as a function of n, and for convenience, the isospin density, ρ I (which in this case is the number density). We note that for the 12-π + system, the number density is ρ (12) I = 12/L 3 = 0.77 fm −3 . In fig. 18, the solid curves corresponds to the prediction at O(L −7 ) from eq. (1), and this prediction differs insignificantly from that at O(L −6 ). There have been recent works that perform lattice QCD calculations at finite isospin chemical potential, e.g. Ref. [62,63], where the starting point is a chemical potential term for the quarks added to the QCD action. Compelling results for the formation of the charged-pion condensate at µ I ∼ m π have been produced.
Pionic systems at finite isospin chemical potential (and temperature) have been explored theoretically in χPT [64,65,66,67], primarily as a step toward understanding finite density nuclear systems. The inclusion of the isospin chemical potential can be accomplished . For µ I < m π the vacuum state is the same as it is for µ I = 0, Σ = 1, but for µ I > m π the vacuum alignment changes and becomes, using the notation of Refs. [66,67]. Minimizing the potential energy at LO in χPT gives cos α = m 2 π /µ 2 I and a relation between the isospin density, ρ I , and chemical potential 15 , By construction, our lattice QCD calculation is in the regime of µ I > m π , and we expect that these LO χPT results should come perturbatively close to describing the properties of the n-π + systems. Implicit in this LO χPT analysis are not only the π + π + -interactions, but also multi-pion interactions, including the π + π + π + -interaction. In fact, the NDA estimate of the π + π + π + -interaction arises from the LO terms in χPT. In fig. 19, we show the isospin chemical potential versus the number of π + 's in the lattice volume, which can be directly translated into isospin density, ρ I = n/L 3 . To determine the importance of the π + π + π +interaction, we show the µ I calculated directly from the lattice calculation. We also show the curve resulting from eq. (1) with the values of a (I=2) ππ and η L 3 extracted from the lattice calculation and their associated correlated uncertainties (red curves and shaded regions), and we show the curve resulting from setting η L 3 = 0 in eq. (1) (blue curve and shaded region). It is clear from fig. 19 that the π + π + π + -interaction plays an important role in the relation between the isospin chemical potential and the isospin density. Further, in fig. 19, the dashed curve is the result of LO χPT, as given in eq. (28), which is seen to describe the result of the lattice calculation well at all the pion masses we have explored. There does appear to be a slight m π -dependent systematic difference, but the magnitude of this effect is consistent with terms higher order in χPT.
Another quantity of interest that one wishes to determine for such systems is the pressure as a function of the ρ I , Unfortunately, we cannot recover the pressure directly from our lattice calculations using eq. (1) because both a (I=2) ππ and η L 3 depend implicitly upon the volume. Further, without lattice calculations of different volumes, the pressure cannot be determined directly from the lattice calculations either. Therefore, with the present lattice calculations we have performed, while the isospin chemical potential can be determined as a function of density, the pressure cannot.
An important issue to consider is the impact of the finite lattice volume upon the relation between the isospin chemical potential and the isospin density. The periodic boundary conditions imposed in the spatial directions leave the zero-mode untouched, but discretize the non-zero modes. As such, the tree-level matrix element of the interaction Hamiltonian between initial and final states in which each pion carries zero three-momentum is unaffected by the finite volume. However, the boundary conditions will modify contributions at the one-loop level and beyond. Therefore, in the limit where the contributions from looplevel diagrams are small, as is the case for a small scattering length, the modifications of the relation between isospin chemical potential and isospin density due to the boundary conditions is expected to be perturbatively small 16 . In systems with large scattering lengths, 16 In infinite volume, the leading contributions to the ground state energy of an imperfect dilute Bose gas e.g. nucleons, the finite-volume modifications may be substantial, and this requires further study.
are calculated easily by performing a Bogoliubov transformation on the Hamiltonian, leaving the leading order mean-field contribution and a term resulting from a summation over the non-zero momentum states. The leading terms in the density expansion of the energy per particle on the ground-state are When placed in a finite cubic volume with periodic boundary conditions, the integral over non-zero momentum states that gives rise to the term that is non-analytic in the number-density, ρ, is replaced by a summation over the allowed momentum states in the volume. When the unperturbed level spacing in the volume is large compared with the energy-shift at leading order in the density expansion, the contents of the square brackets in eq. (30) becomes 1 + 128 15 recovering the leading terms in the finite-volume expansion in eq. (1) as the number of particles becomes large. See also Ref. [9].

VII. CONCLUSIONS
One of the major challenges facing nuclear physics is the solution of the hadronic many-body problem, including the effects of the multi-hadron interactions induced at the scale of chiral symmetry breaking. While the two-nucleon interaction is overwhelmingly the dominant interaction in multi-nucleon systems, a three-nucleon interaction contributes significantly to the structure and interactions of nuclei. Calculating the properties and interactions of nuclei is a major goal of lattice QCD and as a small step in this direction we have performed the first lattice QCD calculations of the properties of systems comprised of multiple-π + 's. The present paper is a detailed follow-up to the letter we recently published [4], and extends that work from studies of systems comprised of up to five pions, to systems containing up to twelve pions. As a technical detail, this required calculations of very high (64 decimal digit) precision.
We have convincingly determined a non-zero value for the dressed three-body interaction, η L 3 , defined in eq. (1), at the lightest three pion masses, m π = 291, 352 and 491 MeV, and find a value consistent with zero at the heaviest mass, m π = 591 MeV. These central results are summarized in Table III and fig.16 above. Given the lattice volumes in which the calculations have been performed, and the value of the π + π + scattering length, the connection between η (5)) is slowly converging, making a meaningful extraction of η 3 (µ) currently impractical.
Clearly, in addition to higher order calculations (O(L −8 ) and beyond) of the energy of multi-pion systems in a finite volume, lattice calculations in larger volumes, and of excited states in these volumes, are required to disentangle the bare from the dressed interactions. An important outcome of our calculations is a determination of the relation between the isospin chemical potential and the isospin density. For m π < ∼ 400 MeV the π + π + π + interaction is seen to make a sizable contribution to this relation. Such contributions are implicit in the LO χPT theoretical result, but our calculation has provided the first explicit QCD calculation of the importance of multi-pion interactions.
This work and our previous Letter represent the first study of multi-hadrons systems directly from QCD. Whilst encouraging, we conclude by reiterating that a significant amount of theoretical work remains to be performed in order to explore meaningfully more complicated systems such as those involving large numbers of baryons. The isospin chemical potential as a function of the number of pions (equivalent to the isospin density ρ I ) at a fixed volume with and without the contribution from the π + π + π +interaction, η double precision (squares) and using 64 decimal digit internal precision (stars, slightly offset to the right for clarity) for n = 9, 10, 11 and 12 for the m π = 352 MeV ensemble. The breakdown of double-precision calculations with increasing n and t, where non-zero values of the correlation functions are lost to noise, is clear. This breakdown has its origin in the large powers to which propagators are raised in the contractions required for systems with a large number of mesons. Viewing the zero three-momentum propagator as a timedependent 12 × 12 matrix, the range of numbers that can be represented in double precision is too limited and small elements of high powers of these propagators are rounded away in tracing or other matrix manipulations. Such terms are crucial for maintaining gauge invariance, and their removal leads to a serious degradation of the correlation function. For the numbers of mesons considered in this work, n ≤ 13, it is sufficient to use 64 decimal digit precision, but for n > ∼ 13, such precision will rapidly become insufficient.

Propagator Precision
Given a particular gauge field, the accuracy with which the quark propagator is computed (e.g. the tolerance in conjugate gradient (CG) algorithm) will influence the correlation functions computed from that propagator. In extreme cases, where the solution is quite inaccurate in comparison to the statistical precision of the importance sampling procedure, this residual error can persist through the ensemble averaging and lead to errors in the correlation functions. In the case of the multi-pion correlations studied here, the situation is particularly acute as the high powers to which the propagators are taken can enhance small effects. At this point in time we have not performed a detailed study of this issue. It would require generating full sets of propagators at different inversion precisions on the same configurations. We have simply studied the effect on a few representative configurations, and it is clear that for large n, significant differences in the correlation functions can arise from loose tolerance of the CG-solver. In future studies of even larger n (and very high precision studies of any observable), this issue must be investigated in more detail.