Pion condensation in the early Universe at nonvanishing lepton flavor asymmetry and its gravitational wave signatures

We investigate the possible formation of a Bose-Einstein condensed phase of pions in the early Universe at nonvanishing values of lepton flavor asymmetries. A hadron resonance gas model with pion interactions, based on first-principle lattice QCD simulations at nonzero isospin density, is used to evaluate cosmic trajectories at various values of electron, muon, and tau lepton asymmetries that satisfy the available constraints on the total lepton asymmetry. The cosmic trajectory can pass through the pion condensed phase if the combined electron and muon asymmetry is sufficiently large: $|l_e + l_{\mu}| \gtrsim 0.1$, with little sensitivity to the difference $l_e - l_\mu$ between the individual flavor asymmetries. Future constraints on the values of the individual lepton flavor asymmetries will thus be able to either confirm or rule out the condensation of pions during the cosmic QCD epoch. We demonstrate that the pion condensed phase leaves an imprint both on the spectrum of primordial gravitational waves and on the mass distribution of primordial black holes at the QCD scale e.g. the black hole binary of recent LIGO event GW190521 can be formed in that phase.

Introduction. The origin of matter-antimatter asymmetry in the Universe is unknown as yet. There are several theoretical attempts to explain this fact which has to originate from the evolution of the very early Universe [1,2]. The asymmetry can be expressed in terms of the values of charges that are conserved in the standard model: baryon number B, electric charge Q and lepton number L. These numbers are conserved during the cosmic evolution following baryo-and lepto-genesis [1][2][3][4]. Neutrino oscillations start to occur in the early Universe at T ∼ 10 MeV, therefore, at higher temperatures not only the total lepton asymmetry is conserved, but also the individual electron, muon, and tau lepton asymmetries. Conservation of these numbers leads to the evolution of chemical potentials of different particles that were present in the thermal bath and contributed to the equation of the Universe at early eras.
Recently, the LIGO experiment detected several gravitational wave (GW) events from the merger of black holes predicted by general relativity [5,6]. GWs may also have a cosmic origin due to the inflation or possible cosmic (phase) transitions [7]. Primordial gravitational waves (PGWs) can be produced from the perturbation of spacetime [8,9] by the inflationary phase in the early Universe [10]. Passing through the different stages of cosmic history like the QCD and electroweak transitions, and the matter dominated epoch will leave imprints on PGWs due to the variation of the Hubble expansion rate [11][12][13][14].
Black holes (BHs) can either form by the collapse of matter in stars or in the early Universe due to primordial density perturbations generated by inflation [15,16]. The latter ones are known as primordial black holes (PBHs) -possible dark matter candidates [17]. The formation of PBHs is caused by the collapse of inhomogeneous high density regions during the time the modes cross the horizon [18][19][20]. These processes depend on the inflationary scenario and the scales adopted, as well as on the thermal history of the early Universe, making them sensitive to the matter-antimatter asymmetry.
For an isentropic expansion of the Universe it is common to express the asymmetries in terms of the conserved charge per entropy ratios: b = n B /s, q = n Q /s, and l α = n Lα /s with α = e, µ, τ . One can associate a chemical potential to each of the conserved charges B, Q, and {L α }. The cosmic trajectory is a line in the sixdimensional space of T , µ B , µ Q , µ e , µ µ , and µ τ defined by five conservation equations: n Lα (T, µ Q , {µ α }) s(T, µ B , µ Q , {µ α }) = l α , α ∈ e, µ, τ.
The conserved charge and entropy densities entering the above equations are given as functions of the temperature and chemical potentials through the equation of state of cosmic matter. For the cosmic QCD epoch, the equation of state is mainly determined by strongly interacting matter, but also contains the contributions of leptons and photons. Naturally, non-trivial dynamics is mainly contained in the QCD part.

arXiv:2009.02309v1 [hep-ph] 4 Sep 2020
Tight constraints on the baryon asymmetry and electric charge are available: b = (8.60 ± 0.06) × 10 −11 and q = 0. The total lepton asymmetry in the standard scenario arises through sphaleron processes, giving l = −(51/28)b, equally distributed among the three lepton flavors [21]. This yields the standard cosmic trajectory where all chemical potentials are close to vanishing for the majority of the cosmic trajectory. Values of the total lepton asymmetry considerably larger than the baryon one are also possible: Ref. [22] gives the constraint of |l| < 0.012. Here l = l e + l µ + l τ . A recent analysis of Ref. [23] shows that pion condensation is unlikely to occur under this constraint if the lepton asymmetry is equally distributed among the three flavors. However, due to the absence of neutrino oscillations at T 10 MeV, the individual lepton flavor asymmetries are not strongly constrained. It has been pointed out in [23,24] that sufficient conditions for pion condensation to occur can be achieved for unequally distributed lepton asymmetries. Complementary to [24], in the present letter we determine these conditions specifically. Moreover, we point out, for the first time, signatures of a pion-condensed phase in the early Universe, namely its impacts on the spectrum of PGWs and on PBH formation.
Equation of state. Pion condensation is expected to occur if the electric charge chemical potential µ Q exceeds the pion mass. First-principle lattice QCD studies at finite isospin density do suggest pion condensation to take place at T 160 MeV and µ I m π [25,26], with µ I being the isospin chemical potential (we get back to the distinction between µ I and µ Q in Appendix B). Here we analyze the cosmic trajectories determined by Eqs. (1)- (3) at different values of l e , l µ , and l τ to determine the conditions for pion condensation to occur. Notice that the weak decays of pions are blocked in the present setting of weak equilibrium, since all outgoing neutrino states are filled due to the high lepton chemical potentials, stabilizing the pion condensate [27].
Neglecting QED interactions, the standard model equation of state is partitioned into contributions from QCD, leptons and photons: The leptonic pressure is modeled by an ideal gas of charged leptons and neutrinos, including all three lepton flavors. The photonic pressure is given by a massless ideal gas of photons.
As we focus our study on temperatures T < 160 MeV that are relevant for hadronic matter, the QCD pressure is approximated by a variant of the hadron resonance gas (HRG) model. In the standard HRG model one includes all known hadrons and resonances as free particles. The HRG model provides a reasonable description of the QCD equation of state in this temperature range when compared to the results of first-principle lattice QCD calculations [28,29]. To incorporate the pion-condensed  phase we modify the HRG model by replacing the free pion gas by an interacting pion gas, modeled by a quasiparticle (effective mass) approach [30] matched to chiral perturbation theory [31] and lattice QCD results at zero temperature. We describe the details of this model in Appendix A. Furthermore, we carry out a systematic check of the model by comparing it to lattice results at µ I > 0, identifying a validity range where the model is deemed reliable, see Appendix B for details.
The QCD pressure thus consists of the pressure of three pions species, each described by an effective mass model, and by contributions of the rest of the hadrons and resonances that are modeled as free particles: Here µ j = B j µ B + Q j µ Q with B j and Q j being the baryon and electric charge of hadron species j, respectively. The index i sums over the three pion species and the index j sums over all hadrons excluding pions. We include all established light flavored and strange hadrons listed in Particle Data Tables [32]. All the conserved charge densities and the entropy density entering Eqs. for the chemical potentials at each temperature. The numerical solution is achieved using Broyden's method [33]. The procedure is implemented within an extended version of the open source Thermal-FIST package [34]. We tested this procedure by reproducing the cosmic trajectories reported in Ref. [23] using the HRG model.
Cosmic trajectories. We fix b = 8.6 · 10 −11 and perform a parametric scan in l e and l µ . As the restriction |l| < 0.012 on the total lepton asymmetry is rather strong we shall set l τ = −(l e + l µ ), meaning that we have a vanishing total lepton asymmetry (l = 0) in all our calculations. For each value of l e and l µ , we start calculations at T = 10 MeV, where all cosmic trajectories are very similar, and gradually increase the temperature. If the cosmic trajectory enters the phase with a Bose-Einstein condensate of pions, we register the temperature T cond where the trajectory crosses the pion condensation boundary.
Our calculations reveal that T cond depends mainly on the sum l e + l µ of the electron and muon lepton asymmetries, whereas the dependence on the difference l e − l µ is mild. This is shown in Fig. 2, where we depict the dependence of the temperature T cond on the sum l e + l µ . The difference l e − l µ is varied in a range |l e − l µ | < 0.5, giving the narrow black uncertainty band in Fig. 2. At temperatures between T cond and T pc ≈ 160 MeV the cosmic matter is in a pion-condensed phase. We find that pion condensation occurs in the early Universe at T < 160 MeV if the following condition is met: Pion condensation is not observed at smaller absolute values of l e + l µ . The relation (6) can therefore be regarded as a universal criterion for pion condensation in the early Universe. Positive values of l e + l µ correspond to π + condensation, while negative l e + l µ imply π − condensation.
The temperature dependence of µ Q is shown in Fig. 1 for several different values of lepton flavor asymmetries, l e +l µ = 0, 0.1, 0.2, and 0.4. For l e +l µ = 0 one essentially recovers the standard cosmological trajectory where µ Q is very close to zero throughout and far away from the pion condensed phase. For sufficiently large absolute values of l e + l µ [see Eq. (6)], the cosmic trajectory crosses the pion condensation boundary. The kink-like structure in the cosmic trajectory, visible for the l e + l µ = 0.4 case at T ≈ 95 MeV, is associated with a rapid growth of the lepton chemical potentials.
The equation of state exhibits an interesting behavior for trajectories that enter the pion condensed phase. Of particular interest is the interaction measure, (ε−3p)/T 4 . The interaction measure is negative deep in the pioncondensed phase at moderate temperatures (see Fig. 1) -a distinctive feature of the pion condensed phase also seen in lattice QCD calculations. Figure 3 depicts the temperature dependence of (ε − 3p)/T 4 along the cosmic trajectory for the three different cases of negative l e + l µ  values discussed above. The behavior of these two quantities is significantly affected at large lepton asymmetries. For |l e + l µ | 0.3 the cosmic trajectory passes through a region with negative (ε − 3p)/T 4 , as illustrated by the magenta curve in Fig. 3 for l e + l µ = 0.4. Negative interaction measure correlates with large sound velocities that go above the conformal limit of c 2 s = 1/3. The interaction measure grows to large values (ε − 3p)/T 4 10 at larger temperatures. This drastic rise is a consequence of large lepton chemical potentials at these temperatures, which emerge from lepton flavor number conservation.
Effects on the spectrum of PGWs. Due to the presence of a nonvanishing lepton asymmetry and the possible formation of the pion-condensed phase, the equation of state before big bang nucleosynthesis (BBN) can change, which will leave an imprint on the PGW spectrum [13,14,35,36].
The evolution of each polarization λ of tensor pertur-bation h for a mode k in cosmology is given by [10,37] where the derivative is with respect to conformal time d/dη and a is the scale factor (adη = dt, t is the cosmic time). The primordial tensor perturbation can be written in terms of the transfer function X, tensor perturbation amplitude h prim k,λ and tensor power spectrum parameterized with respect to a characteristic scalek = 0.05 Mpc −1 where A T = r A S and A S , n T are scalar and tensor perturbation amplitudes, and the tensor spectral index, respectively. The tensor to scalar ratio denoted by r has an upper limit from measurements by PLANCK of r 0.07 [38,39].
To compute the temporal evolution of the scale factor one needs to solve the Friedmann equation (H 2 = (ȧ/a) 2 = (8π/3M 2 P l )ε, M P l = 1.22 × 10 19 GeV). We solve Eq. (7) for a mode k using (8) until horizon crossing 1 , i.e. when k = |k| = a(η h )H(η h ), then we use the WKB (Wentzel, Kramers, Brillouin) approximation for the PGW afterwards until today [11,12]. Using Eqs. (7) and (8) the relic density of PGWs for different frequencies ν = k/2π at today (a 0 ) can be computed from [11,12] Using the equations of state computed for different lepton asymmetry values, for which the cosmic trajectory can enter the pion condensed regime, one can estimate the PGW spectrum by using Eqs. (7)-(9). We consider entropy conservation (s a 3 = const.) and use the number of degrees of freedom after neutrino decoupling [40] to find the relation between the scale factor and the temperature. The PGW relic spectra are shown in Fig. 4. As the lepton asymmetry increases, so does the amplitude of the spectrum because the entropy, energy and pressure densities become larger. Moreover, the formation of pion condensation can enhance the PGW due to the change of equation of state. Pulsar timing arrays, as the future SKA [41,42], can measure the predicted PGW spectrum especially around the QCD phase transition if it is scale invariant (n T = 0) or blue-tilted (n T > 0). The LISA experiment [43] can also measure such effect at higher frequencies. The lepton asymmetry at BBN time and afterwards is constrained by cosmic microwave background measurements. Since nonvanishing lepton asymmetry and pion condensation before BBN can modify the PGW spectrum, GW observatories with high sensitivity are able to measure these effects in the early Universe. 1 The initial conditions that we consider at superhorizon scale (k aH) are: X(k, 0) = 1 and X (k, 0) = 0.  Impact on the formation of PBHs. The population of primordial black holes that formed in the early Universe depends on the Hubble rate and the total mass within the Hubble horizon [44][45][46][47][48][49][50]. As mentioned earlier, a nonvanishing lepton asymmetry and a pion condensed phase modify the Hubble rate thereby modifying the production of PBHs in specific range of masses. The horizon mass, defined as M h = 4π 3 H −3 ε [15,16], relates a given temperature in the early Universe to the horizon mass and later on to a typical black hole mass M BH . Figure 5 shows the fraction f PBH of PBHs with respect to total cold dark matter (CDM) abundance for different lepton asymmetry cases (the details of the calculation can be found in Appendix C). The presence of pion condensation is signalled by a modification of f PBH at masses larger than one solar mass. What we have shown is mainly the effect of uncertainty from thermal bath due to nonvanishing lepton asymmetry and pion condensation on the formation of PBH in the early Universe.
The parameter f PBH can be indirectly measured by different experiments. The fraction of PBHs with masses 10 −6 M M BH 10 3 M from some experimental constraints (OGLE, HSC, Caustic, EROS, MACHO) should be f PBH 0.05 2 [53][54][55][56][57]. The SKA [41,42] and LISA [43] can also indirectly constrain the fraction of PBHs by putting limits on the induced PGWs from curvature perturbation or using GWs produced by coalescing events [13,58,59]. Summary. The present analysis of cosmic trajectories at non-vanishing lepton flavor asymmetries reveals a simple criterion for the onset of pion condensation in the early Universe -it occurs when the total electron and muon asymmetry parameter is sufficiently large, |l e +l µ | 0.1. This results does not exhibit large sensitivity to the modeling of pion interactions. For even larger asymmetries, |l e + l µ | 0.4, the cosmic path enters a region with negative interaction measure, (ε − 3p)/T 4 < 0, located deep inside the pion condensed phase. The possible presence of such a Bose-Einstein condensed phase of pions would have significant cosmological implications such as the strong enhancement of the spectrum of PGWs and the change of the fraction of PBHs with mass larger than one solar mass. The experimental signatures of pion condensation from the early Universe can be probed by pulsar timinng and GW detectors. The recent BHs merger event of LIGO GW190521 can be from PBHs produced during the pion condensation epoch [60,61].
Pion condensation could also affect big bang nucleosynthesis. If the pion condensed phase is present, spheres of pions and leptons, pion stars, can form which are stabilized by the high density of neutrinos due to the high lepton chemical potentials [26,62,63]. Typical pion star masses will be in the range of a few solar masses when the early Universe leaves the pion condensed phase. The neutrinos will diffuse out of the pion stars on the timescale of weak interactions. The situation is similar to the one for proto-neutron stars where neutrinos leave on the timescale of several seconds. Hence, pion stars would decay around the time of BBN. The produced high energy leptons would influence the abundances of primordially produced nuclei, which could be addressed by a modified BBN simulation.

A. Effective mass model for pion condensation
We use a quasiparticle (effective mass) approach to describe interacting pions with a pion-condensed phase. Outside of the pion condensed phase, the pressure of a single pion species in the effective mass model reads [30] p EM π (T, µ π ; m * ) = p id π (T, µ π ; m * ) + p f (m * ).
Here π ∈ π + , π − , π 0 . The rearrangement term p f (m * ) is a consequence of interactions and it preserves the thermodynamic consistency in the quasiparticle model. The specific form of p f (m * ) defines the quasiparticle model.
Here we take p f (m * ) in the form chosen to match the model to chiral perturbation theory and lattice QCD results in the pion-condensed phase at T = 0 (see below). The pressure at a given T and µ π has to be maximized with respect to m * , resulting in a gap equation (∂p EM π /∂m * ) T,µ = 0: Here n id σ (T, µ π ; m * ) ≡ −∂p id π /∂m * is the scalar density of an ideal gas of pions with mass m * . A numerical solution to the gap equation determines m * at given T and µ π , allowing to calculate all other thermodynamic quantities through Eq. (11). The transition to the pion-condensed phase takes place when the effective pion mass becomes equal to the chemical potential, m * = µ π . The equation determining the transition line in the µ π -T plane reads [30] p f (µ π ) = n id σ (T, µ π ; m * = µ π ) .
The effective mass equals the chemical potential in the phase diagram region with a pion condensate, m * = µ π for µ π ≥ µ cond , as a consequence of interactions between thermal and condensed pions [64]. Here µ cond is the pion chemical potential at pion condensation boundary. Therefore, the pressure in this phase reads reads p EM π (T, µ π ) = p id (T, µ π ; m * = µ π ) + p f (µ π ).

B. Lattice simulations
Here we describe the details of our first-principles lattice QCD simulations at nonzero isospin density. On the one hand, the lattice results at (approximately) zero temperature are used to guide the construction of the effective mass model described above. Here we use our data at a single lattice spacing from Ref. [26]. On the other hand, the finite-temperature results serve to test the validity range of the model at nonzero isospin and zero baryon density. To this end we employ our data from Refs. [25,65] on four lattice spacings.
To simulate the path integral Z we take the tree-level Symanzik-improved gauge action and 2 + 1 flavors of rooted staggered quarks with physical masses [66]. The isospin chemical potential µ I enters the Dirac operator 3 via the quark chemical potentials µ u = −µ d = µ I /2, while µ s = 0. Comparing to the standard basis with baryon and charge chemical potentials, one can read off µ Q = µ I , µ B = −µ I /2. The simulations therefore correspond to a situation with a specific linear combination of baryon and charge chemical potentials, which only couples to hadron species containing an unequal number of up and down quarks (predominantly charged pions). In addition, we need to introduce an auxiliary pionic source λ > 0 that is extrapolated to zero at the end of the analysis. The role of the λ parameter is twofold. First, it triggers the spontaneous symmetry breaking corresponding to pion condensation in a finite volume. Second, it serves to stabilize the theory in the infrared by making the Goldstone boson of the pion condensed phase slightly massive [25].
To calculate the equation of state, our primary observable is the isospin density The details of the λ → 0 extrapolation of this observable are explained in Ref. [65] and in the following we work with the so extrapolated quantity. From n I , we can calculate ∆O(T, µ I ) ≡ O(T, µ I ) − O(T, 0) for any observable O. In particular, the pressure difference and the trace anomaly difference can be constructed as ∆I(T, µ I ) = µ I n I (T, µ I )+ The zero-temperature results for n I near µ I = m π are well-described by the chiral perturbation theory formula (15) with f π = 133(4) MeV [26]. This is smoothly matched by a spline interpolation for n I (µ I ) at higher values of the chemical potential. The interaction measure is determined via Eq. (18) -note that at zero temperature ∆I = I and, moreover, the first contribution to the integral in ∆I of Eq. (18) vanishes, simplifying this expression considerably. The so obtained curve is plotted in Fig. 6 as the yellow band.
For testing the effective mass model at T > 0 we concentrate on ∆I because compared to other observables it is found to contain the least amount of lattice discretization errors 4 . The integrals and the derivatives in Eq. (18) need to be evaluated numerically. To this end we fit n I (T, µ I ) via a two-dimensional spline surface. The spline nodepoints are drawn from a Monte-Carlo procedure with the goodness of the fit playing the role of the action, providing a direct estimate of systematic errors (see Ref. [67] for more details). The µ I -dependence of ∆I is plotted for two representative values of the temperature in Fig. 7. Here we include the results for our two finest lattice spacings, N t = 10 and N t = 12. (The continuum limit at constant T corresponds to N t → ∞, but we do not carry out this extrapolation here.) The model is found to capture the notable features of the lattice data qualitatively. A quantitative description is obtained if neither T nor µ I are too large.
To make the comparison between the N t = 12 lattice results and the model more systematic, in Fig. 8 we show the deviation between the two in the form of a heat plot. Here we normalize by the error σ of the lattice resultstherefore a value of n indicates a difference by n standard deviations. The plot shows substantial differences for µ I > m π at high temperatures as well as slight deviations near the boundary of the pion condensed phase. We take the contour line at 3 standard deviations as a marker and consider the model reliable in the parameter range where This range is indicated by the solid line sections of the cosmic trajectories in Fig. 1. The above comparisons were performed at nonzero isospin chemical potential µ I , where lattice results are available. For the analysis of the cosmic trajectory, the model is employed instead at nonzero charge chemical potential µ Q (as well as low baryon chemical potential µ B ). At zero temperature, µ I and µ Q can be identified  as long as the only charged states that contribute to the equation of state have zero strangeness and zero baryon number. This is the case for µ I < m K (even in this case, kaon condensation is not expected to occur if a pion condensate is already present [68]) and sufficiently low µ B as is the case for the parameters considered in this paper. Contrary to the identification µ I = µ Q at zero temperature, for T > 0 the different couplings of the two chemical potentials to hadronic states becomes relevant and the equation of state differs in the two cases. Nevertheless, in the effective mass model the pion condensation boundary expressed in µ I or in µ Q remains the same, because interactions between pions and other hadrons are neglected in the model. The difference between the critical lines, µ crit Q (T ) and µ crit I (T ) can be estimated using lattice results for the estimators of the convergence radii of the corresponding Taylor series around µ Q = µ I = 0. In particular, we consider the expansions of the pressure, and the estimators for the convergence radius for the susceptibilities χ I,Q = ∂ 2 p/∂µ 2 I,Q . We use the Taylor coefficients determined in Ref. [28] for our action and lattice spacings. The leading estimator for the isospin direction was found to give a remarkably good approximation to the true critical line, µ crit I [65]. We assume this is also the case for the expansion in µ Q . Thus we approximate the critical line in the electric charge direction by For the second factor we use the lattice results [28] above and ideal HRG with quantum statistics below a matching temperature of 105 MeV. The N t = 12 results for this approximation, together with the corresponding directly determined isospin critical line µ crit I [25,65], are plotted in Fig. 9.
Comparing to the effective mass model, considerable differences are only observed for T 160 MeV. This, and the fact that the lattice results for µ crit I and µ crit Q do not differ by more than a few percents again confirms that our model represents a reasonable approximation to the phase diagram at nonzero isospin (charge) densities.

C. Primordial Black Holes Formation
At the time of PBH formation a region of the Universe within the Hubble horizon starts to collapse due to local inhomogeneities amounting to The relation between the amplitude of the density perturbation δ, the PBH mass M BH and the horizon mass M h can be defined as [49,50] Figure 10.
Threshold of primordial black hole formation versus horizon crossing temperature for different values of the lepton asymmetry.
The fraction of PBHs Ω PBH with respect to the total cold dark matter (CDM) abundance Ω CDM reads [44] f PBH (M BH ) = 1 Ω CDM The mass or scale dependence of the density perturbation width can be assumed to be [44] The density spectral index n M can be related to the scalar spectral index n S − 1 −2n M , where n S 0.96 [38,39]. We choose the benchmark value of n M = 0 to compute the fraction of PBH from Eq. (25). The parameter f PBH for masses smaller than M increases (decreases) when n M is negative (positive). However, f PBH increases (decreases) for larger masses, respectively. For a fixed n M as lepton asymmetry increases the value of f PBH will change depending on the behavior of ω or δ c and the energy and pressure density. When δ c increases (decreases) the fraction of PBH decreases (increases).