Cosmology meets functional QCD: First-order cosmic QCD transition induced by large lepton asymmetries

The lepton flavour asymmetries of the Universe are observationally almost unconstrained before the onset of neutrino oscillations. We calculate the cosmic trajectory during the cosmic QCD epoch in the presence of large lepton flavour asymmetries. By including QCD thermodynamic quantities derived from functional QCD methods in our calculation our work reveals for the first time the possibility of a first-order cosmic QCD transition. We specify the required values of the lepton flavour asymmetries for which a first-order transition occurs for a number of different benchmark scenarios.

Introduction According to the hot big bang model, the early Universe has experienced (at least) two epochs where phase transitions could occur: the electroweak transition at temperatures around T ew ∼ 100 GeV when fermions and gauge bosons became massive particles; and the transition of quantum chromo dynamics (QCD) at T QCD ∼ 130 MeV when quarks confined into hadrons.
The Standard Model (SM) of particle physics predicts both the electroweak as well as the QCD transition to be crossovers. However, various observations from cosmology and particle physics unavoidably hint towards the existence of physics beyond the SM. Many extensions of the SM would render the electroweak transition to be first-order (for the most simple extension see [1]), while there are much fewer mechanisms known that could provide a first-order QCD transition ( [2,3]). A first-order QCD transition would have an impact on the spectrum of gravitational waves and the production of exotic relics [4,5]. The details of the QCD transition furthermore have an impact on the production of primordial black holes [6][7][8][9][10] which are one of the prime candidates for dark matter.
In this work, we show for the first time that large lepton flavour asymmetries can render the QCD transition to be first-order. The fact that the value of the lepton asymmetry has an impact on the cosmic QCD epoch has already been known for a while and been discussed in a series of papers [11][12][13]. However, so far it has not been possible to provide a definite answer to the question whether a sufficiently large lepton asymmetry can indeed induce a first-order transition and what sufficiently large means quantitatively.
QCD diagram QCD matter experiences a transition from quarks and gluons to hadrons as temperature and chemical potential changes. To chart the phase diagram of QCD for non-zero temperature, charge-and baryon chemical potential has attracted a lot of interest. Lattice QCD simulations have confirmed a crossover at vanishing chemical potential [14][15][16][17]. However, due to the sign problem, it is difficult for lattice QCD to explore the region of large chemical potentials, leaving the possibility of the existence of a critical end point (CEP) and a region of a first-order transition unsettled. On the observational front, searching for signals of the QCD phase transition, especially the CEP, are the main goals of the current and the future experimental program at the relativistic heavy ion collider [18][19][20][21]. On the theoretical front, large chemical potentials are accessible by continuum QCD methods, i.e. functional QCD methods including Dyson-Schwinger equations (DSEs) [22][23][24][25] and the functional renormalization group (fRG) method [26,27]. The functional QCD methods are non-perturbative continuum methods which are capable of describing both the dynamical chiral symmetry breaking (DCSB) and the confinement simultaneously. Though limited by truncations, the functional QCD approach has made fruitful progress in studying the QCD phase structure and thermal properties [25,[28][29][30] in a large range of chemical potential.
Concerning the cosmic QCD epoch, the question about the nature of the QCD transition is actually two-fold: The first part of the question is if QCD matter under any circumstances could ever experience a first-order phase transition. Even though different functional QCD methods all predict the existence of a CEP (with varying locations) [25,27] its final confirmation is still awaiting experimental evidence. The second part is under which conditions such a first-order transition would have been realized in the early Universe. This work deals with the second part of the question by applying results from functional QCD.
Large Lepton Asymmetry When the Universe expands and cools down it follows a certain path in the arXiv:2106.11991v2 [hep-ph] 18 May 2022 QCD phase diagram, the so called cosmic trajectory. The standard cosmic trajectory is based on the assumption that matter and antimatter are (almost) equally abundant in the Universe. Indeed, measurements of primordial element abundances and the anisotropy spectrum of the Cosmic Microwave Background (CMB) show that the Universe only has a tiny asymmetry in the baryonic sector, characterized by the baryon asymmetry b = (8.7 ± 0.06) × 10 −11 (inferred from [31] and where b is going to be properly defined in eq. (2)). Explaining the creation of this tiny number necessitates theories of baryogenesis and leptogenesis. The idea of leptogenesis is to create an asymmetry in the leptonic sector and transfer this lepton asymmetry into a baryon asymmetry by sphaleron processes. Therefore, the lepton asymmetry would be of the same order of magnitude as the baryon asymmetry, or more explicitly l = − 51 28 b [32]. In this case the cosmic trajectory passes the QCD diagram at (extremely) small chemical potential where lattice calculations reveal a crossover. Observational constraints from big bang nucleosynthesis (BBN) [33] and the CMB [34] on the value of the lepton asymmetry however allow values of the lepton asymmetry as large as |l| < 1.2 × 10 −2 [34], i.e. around 8 − 9 orders of magnitude larger than the baryon asymmetry. From a theoretical point of view, there are also alternative models predicting a large lepton asymmetry [35][36][37][38][39][40]. In any case, lepton asymmetry is a key parameter to understand the origin of the matterantimatter asymmetry of our Universe.
Observationally even less constrained than the total lepton asymmetry l ≡ l e + l µ + l τ are the individual lepton flavour asymmetries l α=e,µ,τ : It has been shown [41,42] that initially different lepton flavour asymmetries are equilibrated at around T osc ∼ 10 MeV when neutrino oscillations become efficient such that finally l α ≈ l 3 [43]. This implies that CMB and BBN constraints are only sensitive to l but have no constraining power on the initial values of the individual lepton flavour asymmetries l α . From an agnostic point of view, the individual lepton flavour asymmetries should therefore be treated as free parameters (while their sum still needs to fulfill the observational constraint |l| < 1.2 × 10 −2 [34]). As recently pointed out in [44], large lepton flavour asymmetries could also explain the baryon asymmetry of the Universe.
Cosmic trajectory Let us in the following explain how we calculate the cosmic trajectory. The basic idea of the method is similar to the one outlined in [11][12][13].
Within the framework of the SM, before the onset of neutrino oscillations at T osc ≈ 10 MeV and after the electroweak transition at T ew ≈ 100 GeV electric charge, baryon number and lepton flavours are conserved in a comoving volume. Assuming as well entropy conservation in comoving volume, the baryon b, lepton flavour number l α and electric charge q asymmetry can be written as and remain constant throughout the evolution of the Universe, for temperatures T ew > T > T osc . The sum in eq.
(2) goes over all baryons and in eq. (3) over all charged particles. n i denotes the net number density (i.e. particle minus anti-particle number density) of particle species i and s denotes the total entropy density. As mentioned above, the baryon asymmetry is a known quantity and we fix it to the central value b = 8.7×10 −11 from [31]. There is furthermore good reason to believe that the Universe is charge neutral [45] such that we may set q = 0 from now on. While the standard assumption is |l α | = O(b), we treat the three lepton flavour asymmetries l α as free input parameters within observational constraints. This may be realized by some new physics violating lepton flavour around some temperature T BSM well before the QCD epoch, such that the l α in eq. (1) are conserved at T BSM > T > T osc . We assume electroweak processes to maintain kinetic equilibrium, i.e. all particles sharing the same temperature T i = T , and chemical equilibrium which allows us to find relations between the chemical potentials of different particle species, e.g. µ u = µ c etc.
Note that we assume no other but the particle content of the Standard Model throughout this work. Each conserved charge is associated with a chemical potential (i. e. µ B , µ Q and µ Lα ), which can be related to the chemical potentials of the different particle species by (see also e.g. [11] or [46]) The cosmic trajectory denotes the solution of eqs.
and given a choice for the lepton asymmetries l α . As first noted in [11], large lepton asymmetries l α do not only shift the cosmic trajectory towards large leptochemical potentials µ Lα but also towards large charge-and baryon chemical potentials.
In practice, for leptons and photons we assume Fermi-Dirac or respectively Bose-Einstein distributions for the phase-space densities, in accordance with the assumption of thermal and chemical equilibrium. Concerning the DCSB and confinement phenomena, modelling the QCD sector at temperatures close to the QCD transition is however much more complicated.
Let us briefly summarize the approach of [12,13] before we explain how we treat QCD matter in this work. In [12,13] the QCD sector was treated separately in three different (but overlapping) temperature regimes: i) at high T quarks and gluons were treated as free particles, where [13] also included contributions from perturbative QCD, ii) at T ∼ T QCD the QCD pressure was Taylor expanded and susceptibilities from lattice QCD were applied, iii) at low T QCD matter was described by the hadron resonance gas (HRG) model. As pointed out in [13], the method has two limitations. The appearance of a Bose-Einstein condensate of pions for large lepton asymmetries (see also [9]) would require a separate treatment of the low-energy modes of pions. The more stringent restriction however comes from the application of the Taylor expansion of the QCD pressure which is only justified for relatively small chemical potentials.
In this work, we therefore follow a different approach and consult results from functional QCD. In particular, we apply the thermodynamic quantities derived from DSEs in the rainbow-ladder(RL) truncation as it is currently the only QCD based method that has delivered a complete computation of the phase diagram and the related thermal properties [25,47,48]. The obtained thermal quantities are qualitatively and partly quantitatively in accordance with those from Lattice QCD simulations [25] and can be consistently extended to large chemical potentials. The existence of a CEP is apparent as a sudden drop in the number density n u/d of the u/d quarks below a certain threshold for the temperature T CEP = 125 MeV and above a certain threshold for the chemical potential µ CEP = 111 MeV (figure provided in the appendix ). Let us however note here that rather than being able to determine the exact location of the CEP the strength of the RL truncation is to capture the main properties of the QCD phase structure. Therefore, our work should also be understood as a proof-ofprinciple of the possibility of a first-order cosmic QCD induced by large lepton flavour asymmetries -and not as a calculation enabling the exact determination of cosmic trajectories during the QCD epoch.
As noted before, case (i) is constrained by CMB observations to values |l| < 1.2 × 10 −2 [34]. Cases (ii)-(vii) were chosen such that the total lepton asymmetry l always fulfills the CMB constraints. Note that there is of course an infinite number of realizations for scenarios with l e = l µ = l τ satisfying the CMB and BBN constraints. Scenarios (ii)-(vii) simply serve as benchmark models.
(1)-(3). The first result is the fact that there exist solutions to the set of equations, i.e. for large enough lepton asymmetries the Universe experienced a first-order QCD transition. Table I summarizes for which values of the lepton flavour asymmetries either the u or the d quark experiences a first-order phase transition. For the equilibrated case (i) we see that l has to be larger than allowed by observations of the CMB [34] and primordial element abundances [33]. We therefore conclude that in case of equal lepton flavour asymmetries a first-order QCD transition is excluded. On the other hand, scenarios (ii)-(vii) are by construction always compatible with CMB and BBN constraints and all allow a first-order transition. It is also interesting to compare tab. I to the findings of [13] where scenarios (i) and (ii) were studied and to [9] where a more general version of scenario (v) was investigated (|l e +l µ | = −|l τ |). The values for l α rendering a first-order transition in tab. I are remarkably close to the maximal values ensuring the reliability of the Taylor expansion in [13], namely |l µ | ≤ 4×10 −2 (case (ii)) and |l| ≤ 7.5×10 −2 (case (i)). Furthermore, for scenarios (i) and (ii) the values in tab. I exceed the values for pion condensation, i.e. |l| ≥ 9 × 10 −2 (case (i)) and |l µ | ≥ 6 × 10 −2 (case (ii)) [13]. For scenario (v) they are only slightly smaller than what is needed to potentially form a pion condensate, i.e. |l e +l µ | > 1×10 −1 [9]. We therefore conclude that in case of a first-order cosmic QCD transition induced by large lepton asymmetries the formation of a pion condensate is also likely. It is furthermore interesting to see how the cosmic trajectory looks like for such a first-order transition. Exemplary we thereby focus on scenario (ii) and a firstorder transition of the d quark. Fig. 1 shows the cosmic trajectory projected onto the (µ B , T )-plane for different values of the lepton flavour asymmetry l µ . Observe the appearance of a kink in the trajectory for increasing values of |l µ |. This discontinuous feature reflects a firstorder transition and arises for l µ values between the green (l µ = −6 × 10 −2 ) and the red curve (l µ = −7 × 10 −2 ), confirming thereby the solution l µ ≤ −6.85 × 10 −2 given in tab. I. However, on top of the jump at T CEP = 125 MeV the red and orange (first-order) trajectories also show some wiggly features at temperatures above T CEP . This (unphysical) behaviour is caused by cumulative errors appearing at large chemical potential in the computation of the entropy density with the RL truncation (see appendix). When studying the number density of the d quark, the first-order signal is however much better pronounced (lower plot in fig. 1), appearing as a sudden jump at T = 125 MeV. Finally, let us compare the method introduced in this work to the method of [12,13] for the case of equal lepton flavour asymmetries (i). In fig. 2 we show the cosmic trajectory derived by the method introduced in this work (solid lines) and by the method described in [13]. We see that at high temperatures our new method produces trajectories that are relatively similar to the ones obtained by [12,13]. In the intermediate temperature regime the solid lines show the same trend as the dashed lines (i.e. the curves bending in the same direction) but the bending of the solid lines is much stronger. At low temperatures the standard trajectory converges to the HRG line relatively smoothly. However, for all larger l values instead of converging to the dashed HRG lines the solid lines rather overshoot. This discrepancy is supposedly due to the fact that the RL truncation only allows to capture the main properties of the QCD phase structure but is not expected to deliver exact results, as stated above. In the same manner also the solutions in tab. I are not expected to be exact.
The truncation scheme can be systematically improved and in particular in [29,51,52] the CEP was estimated as (T CEP , µ CEP ) u/d (110, 200) MeV. The derivation of the thermal quantities with the method of [51] is however computationally more costly than the determination of the CEP and currently still in progress. We plan to incorporate the results into our method in a future study and expect the agreement to the method of [12,13] in fig. 2 to be improved.
Conclusions In this work, we included the thermodynamic quantities of QCD matter derived from DSEs into calculations of the cosmic trajectory during the QCD epoch. Bearing in mind the limitations of the applied RL truncation, our method offers an up-to-now unique possibility to study the cosmic trajectory in a consistent way over a wide range of temperature and chemical potential and covering the main properties of QCD matter. In particular, for the first time our work revealed that large lepton flavour asymmetries can induce a first-order cosmic QCD transition. For equal lepton flavour asymmetries l α this would require a total lepton asymmetry l that is already ruled out by observations of the CMB and BBN. For unequal lepton flavour asymmetries a firstorder transition is however feasible. Examples thereof are given by the cases (ii)-(vii) studied in this work and the corresponding values in tab. I.

APPENDIX
Computation of thermal quantities The required thermal quantities, number density and entropy density, can be computed after solving the quark gap equation of DSEs [47][48][49]. We here apply the same settings for the gap equation as described in Ref. [47], and then solve it through Broyden's iteration method in Fortran.
The number density for the u/d, s and c quarks can be directly computed from the solution of the quark gap equation after subtraction, i.e. through the quark propagator G, The entropy density can then be expressed as the integral along the chemical potential as Due to the lack of gluon pressure in the approach of Ref. [47], we here use the lattice QCD computation at zero chemical potential with N f = 2 + 1 + 1 as parametrized in Ref. [50]. The full entropy is given as: s QCD = s latt (T, µ = 0) + δs(T, µ). The data of the number density and the entropy density are generated on a grid of values for temperature and chemical potential. The grid spans over the temperature range T = 121 − 400 MeV, where we chose a grid size of 1 MeV at T < 140 MeV, a grid size of 5 MeV a 140 < T < 300 MeV and a sparse grid size of 20 MeV at 300 < T < 400 MeV. Even though the data for T < 121 MeV are available, the truncation has been much less reliable in this regime. The number density below T < 120 MeV gives an odd rising behaviour as T decreases, which can possibly be solved by improving the truncation. Nevertheless, the temperature range we chose still covers the CEP and first-order phase transition region.
The pion condensation effect on the quark propagator is not included in the current truncation. Such an effect is in the hadron resonance channel of the quark gluon vertex which is the subleading contribution [51], and hence will barely have an impact on the QCD transition. In some certain cases, the pion condensation phase transition and QCD transition can take place successively. In µ-direction the grid spans over the range 0 < µ < 1000 MeV, with a uniform grid size of 1 MeV. Both number densities and entropy densities of the different particle species were stored in data files.
The data for the number density are generally smooth, except for the characteristic jump when the CEP is crossed, see fig. 3 (top). However, almost negligible (unphysical) changes in the curvature of the number density get amplified significantly once the derivative in T is taken in order to compute the entropy density according to eq. (9). This issue becomes increasingly severe for large chemical potentials since the integral in eq. (9) cumulates small errors in the µ-dependence of the number density. Therefore -on top of the distinctive plateau that is expected for a first-order transition -the entropy density shows a wiggly behaviour for large chemical po-tentials, see fig. 3 (bottom). This unphysical feature of the entropy density causes the cosmic trajectories to look rather wiggly too for large lepton asymmetries.
Computation of the cosmic trajectory In order to calculate the cosmic trajectory we modified the c-code developed in [11][12][13] such that it reads in those data files and interpolates them. The code uses Broyden's method in order to find solutions of eqs. (1)-(3) in the main text for (µ Le , µ Lµ , µ Lτ , µ B , µ Q ), given the above specified temperature values and with l α as input parameters. In order to interpolate the data for number and entropy density in µ-direction we applied cubic spline interpolation.