Varying Entropy Degrees of Freedom Effects in Low-Scale Leptogenesis

We analyse in detail the effect of varying entropy degrees of freedom on low-scale leptogenesis models. As an archetypal model, we consider the Tri-Resonant Leptogensis${}$ (TRL) scenario introduced recently by the authors, where the neutrino-Yukawa coupling matrix is dictated by an approximate $\mathbb{Z}_n$ discrete symmetry (with $n=3,6$). TRL models exhibit no preferred direction in the leptonic flavour space and have the remarkable feature that leptogenesis can successfully take place even if all light neutrinos are strictly massless up to one-loop order. Most interestingly, for TRL scenarios with heavy Majorana neutrinos lighter than 100 GeV, temperature varying degrees of freedom associated with the entropy of the plasma have a dramatic impact on the predictions of the Baryon Asymmetry in the Universe (BAU), and may sensitively depend on the freeze-out sphaleron temperature $T_{\rm sph}$. We find that this is a generic feature of most freeze-out low-scale leptogenesis models discussed in the literature. In the same context, we consider heavy-neutrino scenarios realising dynamics related to critical unstable qudits in the thermal plasma and assess their significance in generating the BAU. The phenomenological implications of TRL scenarios at the intensity and high-energy frontiers are analysed.


Introduction
The existence of matter-antimatter asymmetry within the Universe, as well as the observations of the flavour mixing of neutrino species [1][2][3], provide convincing evidence of physics beyond the Standard Model (SM).In recent years, the Wilkinson Microwave Anisotropy Probe [4] and Planck observatory [5] have measured the Baryon Asymmetry of the Universe (BAU) to an unprecedented precision, η CMB B = (6.104± 0.058) × 10 −10 . (1.1) This value of η B is compatible with the one derived from Big-Bang Nucleosynthesis (BBN) bounds on the abundance of light elements [6].
A minimal and appealing solution to the origin of the observed matter-antimatter asymmetry in the Universe offers the well-studied scenario of leptogenesis [7].This is a simple framework which comfortably satisfies all three Sakharov conditions [8], and only requires a minimal extension of the SM to contain right-handed neutrino species, which are singlets under the SU(3) c × SU(2) L × U(1) Y gauge group of the SM.This minimal extension provides an additional Yukawa coupling with the left-handed neutrino species as well as a Majorana mass term, which violates the lepton number, L, by two units.This minimal extension not only provides a mechanism by which SM neutrinos may be massive but also ensures that this mass is small enough to fit with experimental observations through the famous seesaw mechanism [9][10][11][12][13].Moreover, the expansion of the FRW Universe leads to cooling, allowing particle species that are part of the thermal plasma to fall out of equilibrium, thus providing a mechanism for satisfying Sakharov's third condition.Finally, this framework introduces new sources of CP violation from the neutrino-Yukawa sector, as well as phenomena of heavy-neutrino mixing [14,15] which include coherent oscillations between neutrino flavours [16,17].These effects together permit the generation of a significant lepton asymmetry, which may be rapidly converted into a baryon asymmetry through (B + L)-violating sphaleron transitions [18].
One issue which is commonly drawn from models of leptogenesis is that for a generic model, one would typically require singlet neutrino masses of Grand Unified Theory (GUT) scales due to the lightness of the SM neutrinos.As a result, any hopes of detection at current and future experiments are lost since the interactions between light and heavy neutrino species are suppressed by the heavy neutrino masses.However, an elegant framework by which this may be evaded is the framework of Resonant Leptogenesis (RL) [15,19,20].In RL models, the CP asymmetry is enhanced through the mixing of near-degenerate heavy neutrinos, which have masses satisfying the condition Here, m Nα and Γ α are the mass and decay width of the heavy neutrino species N α , respectively.This framework can achieve successful leptogenesis at sub-TeV scales whilst maintaining agree-ment with current oscillation parameters.Consequently, RL models offer not only the possibility of generating the observed BAU but also of observing effects in near-future experiments.
In this article, we adopt a particular framework of RL with three singlet neutrino species with a mass spectrum in consecutive resonance, a structure which we refer to as tri-resonant [21], leading to Tri-Resonant Leptogenesis (TRL).In a TRL model, the CP asymmetry is maximised through the constructive interference of all three heavy neutrino species, allowing for greater amounts of baryon asymmetry to be generated with larger-scale Yukawa couplings.In particular, this is to be contrasted with the expectation from typical bi-resonant models where the CP asymmetry is generated through the mixing of two heavy neutrinos, with a third decoupled neutrino included to fit neutrino oscillation data.
A notable addition to previous research on this topic is the incorporation of the non-constant nature of relativistic degrees of freedom (dofs) which we study in detail.These enter through the energy and entropy densities of the thermal plasma in the Transport Equations (TEs).It has become a common practice in the literature that at temperatures, T , above the electroweak (EW) scale, the dofs of the plasma are taken to be constant.As shown in a previous study [21], however, even small deviations which are present at this scale can have a significant effect on the generated BAU, in particular for models where the heavy neutrino mass scale is below 100 GeV.In this work, we utilise a set of TEs which take into account the CP-violating effects from both mixing and heavy-neutrino flavour oscillations and show that the impact on the generated BAU from the inclusion of the EW-scale temperature variations of the dofs is the same as shown in ref. [21].In addition, we show that the inclusion of the dofs can induce large dependencies of BAU on the critical temperature of the sphaleron transitions, T sph , for low-scale models, in contrast with near TeV scale models where variations in T sph have minimal effect.
Finally, we consider whether the critical phenomena discussed in [22] may also be present in a thermal plasma.Our analysis shows that similar phenomena may be present in the approach to equilibrium.However, in a cosmological setting where the Universe cools down, these phenomena may be destroyed, with the out-of-equilibrium effects pushing the bath towards a fully mixed state.Conversely, in thermo-static scenarios (i.e. of a thermal bath with constant temperature), it may be possible to observe oscillations in the neutrino number density as the plasma approaches its equilibrium distribution.
The layout of this article is as follows.In Section 2, we introduce the TRL model which is based on the symmetry-motivated Z 6 Yukawa structure.We discuss some of the CP properties which arise from this structure.In Section 3, we introduce a set of coupled TEs which describe the evolution of an appreciable BAU, crucially preserving the temperature dependence of the relativistic dofs and including both the mixing of heavy neutrino species as well as coherent flavour oscillations.Section 4 investigates how the inclusion of the temperature dependence of the dofs impacts the generated BAU and the dependence on the critical temperature of the (B + L)-violating sphaleron transitions is analysed.Moreover, we examine the stability of the generated asymmetry under changes to the initial conditions.In Section 5, we consider the prevalence of critical scenarios within cosmological settings and discuss the TRL model as one which can produce critical phenomena.In this regard, we present in Appendix A the mathematical structure of an N -level unstable quantum system of the non-diagonalizable Jordan form, which may also be called a critical unstable qudit.In Section 6, we summarise our results and show numerical solutions to the TEs. 1 The accessible parameter space generated by the TRL model is delineated, and comparisons to current and future experiments are given.Finally, in Section 7, we concisely present our conclusions and discuss possible extensions and future directions.

The Tri-Resonant Leptogenesis Model
We utilise a minimal extension of the SM through the inclusion of right-handed neutrino fields, ν R , which are singlets of the SM gauge group and have lepton number L ν R = 1.The inclusion of these additional degrees of freedom leads to the extended Lagrangian where L i = (ν i L , e i L ) T with i = e, µ, τ are the left-handed SU(2) L lepton doublets, and Φ = iσ 2 Φ * is the weak-isospin-conjugate Higgs doublet of Φ.The two flavour space matrices h ν and m M are the neutrino Yukawa couplings and the Majorana mass matrix, respectively.We have reserved the boldface notation to indicate matrices with flavour structure.Without loss in generality, we may employ covariant flavour-space transformations to bring the Majorana mass matrix into the diagonal form: In the broken phase of the theory, this Lagrangian gives rise to a light neutrino mass spectrum described by the famous seesaw relation [9][10][11][12][13], where m D = v h ν / √ 2 is the Dirac-neutrino mass matrix, with v being the vacuum-expectationvalue (VEV) of the Higgs field, and m M is the Majorana mass matrix describing the spectrum of heavy neutrino species, N α , with α = 1, 2, . . ., n.In a generic model, the seesaw relation suggests the existence of GUT-scale heavy neutrinos since the mass scale of the Majorana mass matrix is determined by Here we made use of the Frobenius norm, e.g., The Lagrangian of the phenomenologically relevant W ℓN charged current interaction may be written as where ) is the left-handed chirality projection operator.Consequently, the phenomena arising from such heavy neutrinos would be out of reach in current and foreseeable experiments, since, in the broken phase, these interactions are severely suppressed by the lightto-heavy-neutrino-mixing matrix, B ≃ m D m −1 M , where the conventions of [25] were adopted.Therefore, there is great phenomenological interest in finding models where observable physics lies closer to current experimental bounds, such as at sub-TeV mass scales.
As was detailed in [21], by selecting a symmetry-motivated flavour structure for the Yukawa, the light neutrino spectrum can be made to vanish, with the observed masses and mixing observables generated through small perturbations about the symmetric Yukawa couplings.We may parameterise the Majorana mass matrix with small mass splittings as In this expression, we use the average mass for N heavy neutrino species m N .At the leading order, the expression for the light neutrino mass matrix may be given by Hence, for singlet neutrinos with a near-degenerate mass spectrum, the light neutrino mass matrix may be made to approximately vanish by demanding that This motivates the symmetric Yukawa structure where a, b, c ∈ C, and ω are the generators of the discrete group Z 6 .We note that due to the isomorphism Z 3 ≃ Z 6 /Z 2 , the generators of Z 3 also satisfy the zero mass condition, and so they may be used to produce identical results up to a change in the sign in the CP phase.For definiteness, however, we explicitly take the generators of Z 6 .
It is worth highlighting the following observation.Since the above symmetry-motivated Z 6 form appears in the flavour structure of the Dirac mass matrix, flavour covariance entails the vanishing of the tree-level mass matrix which in turn leads to the vanishing of the sum of the tree plus the one-loop contribution as well [25][26][27]: respective masses of the W , Z, and Higgs bosons.As a consequence, the symmetric structure we present gives full control over the light neutrino mass spectrum.As was described in [21], by judiciously rescaling the Yukawa couplings through the Majorana mass terms, it is possible to reinforce that the Z 6 symmetry gives exactly vanishing neutrino masses to all loop orders in perturbation theory.
The Z 6 structure not only offers naturally light neutrino masses but also has large CPphases available for generating significant lepton asymmetry.When one calculates the CP-odd invariant quantity [15,19,28,29] it can be shown that the resulting expression for the Z 6 model depends only on the imaginary part of the generator, ω, viz.
where m Nα is the mass of the heavy neutrino N α .
In order to maximise the available CP asymmetry, we make use of this symmetric structure within the RL framework [19].This framework generates enhanced CP-violating effects through the mixing of heavy neutrino species by taking a near-degenerate heavy neutrino mass spectrum, in particular, through the contributions from the absorptive part of the wavefunction [30][31][32][33].

Making use of the wavefunction and vertex coefficients
a fully consistent resummation of these two wavefunction and vertex contributions is obtained through the effective Yukawa couplings [19][20][21]34] ( hν (2.16) In the above expression, x is the Fukugita-Yanagida one-loop function [7], and (2.17) The CP-conjugate effective Yukawa couplings, associated with the N L c Φ * interaction, and represented by hν − , may be found by using the CP-conjugate tree-level Yukawa couplings, (h ν ) * , in (2.16) in place of h ν .The (rest frame) decay matrices of heavy neutrino species into LΦ and as the CP-conjugate decay to L c Φ * are then given by It has been well established in numerous works [15,19,35] that the CP asymmetry generated through the mixing of two heavy neutrino species is maximised when the mass difference between two heavy neutrino species is of the same scale as the width with such mass splitting being described as resonant.It has additionally been shown in [21] that the CP-asymmetry generated in models with three singlet neutrinos may be maximised when the heavy neutrino mass spectrum is in consecutive resonance with resonant mass splitting between not only N 1 and N 2 , but also between N 2 and N 3 .Consequently, we adopt this triresonant structure along with the Z 6 symmetric Yukawa couplings in an effort to bound from above the size of the allowed light-to-heavy-neutrino mixings whilst maintaining agreement with the observed BAU.

Transport Equations for Leptogenesis
To determine the amount of BAU generated through TRL, we require a set of TEs to track the evolution of the number densities of the heavy neutrino species and the SM leptons. 2 In this section, we lay out these TEs.
In [36], a general framework for the construction of quantum TEs was presented, where a matrix of densities was proposed as a way to track the evolution of number densities.Such a matrix of densities was defined as In the above, V = (2π)3 δ(0) is the infinite 3-volume of the coordinate space, ρ is the density matrix of the thermal ensemble, and a i (p) are annihilation operators for the relevant particle species.In contrast with flavour diagonal TEs, this approach offers the ability to account for both the number densities and quantum correlations between them.From this matrix-ofdensities approach, a set of coupled TEs may be derived, which capture CP-violating effects from both the oscillation between singlet neutrino flavours and mixing.

Heavy Neutrino Transport Equations
We follow a similar prescription to that laid out in [37] 3 , with the TEs for the distribution functions of the singlet neutrinos in the mass basis given by Where C N and C N are collision terms of the TEs.We consider the thermal bath in the relativistic limit, with the bath in kinetic equilibrium due to the rapid decay of heavy neutrinos Γ ≫ H(z).
We also neglect quantum effects such as Bose enhancement and Fermi blocking.In addition, we do not consider the contribution of back-reaction terms to the heavy neutrino TEs.With these assumptions, the pertinent collision terms are written as With these definitions, we may derive a set of TEs for the CP-even and CP-odd linear combinations of distribution functions.However, as was extensively discussed in the literature [19,37,40], we note that the subtraction of the so-called Real Intermediate States (RISs) in the 2 → 2 scattering processes contributes a term at the same scale as the decay processes, i.e.
This term changes the sign of the CP-odd inverse decay terms, allowing the TEs to reach thermal equilibrium.In detail, these TEs read4 : In the above expressions we have defined the CP-even and CP-odd components of the distribution functions as well as the CP-even and CP-odd decay terms: ) We identify the equilibrium distribution function as the distribution which makes the RHS of both of these TEs vanish, and we then see that the equilibrium distribution of the heavy neutrino species is given by To integrate these TEs, we may utilise the kinetic equilibrium assumption, which has been discussed in the literature [40,[42][43][44], where it was shown that only small deviations are to be expected.Making use of the equilibrium distribution given in equation (3.8), we find the corresponding equations for the number density of the neutrino species, ) (3.10b) In these TEs, we have defined the following thermally averaged quantities: where the integral labels indicate over which Lorentz-invariant-phase-space the integral takes place, the rest-frame decay matrices, Γ ± , are calculated in equation (2.18), and K n (z) are modified Bessel functions of the second kind.
The terms which enter the TEs (3.10) are then written as ) We would like to write the TEs given in equation (3.10) in a way which explicitly accounts for the expansion of the Universe and consequential phenomena.Following the conventions given in [20], we introduce the dimensionless run parameter z = m N /T and normalise the particle number density to the photon number density where ζ(3) = 1.202 is Apéry's constant.
From the TEs given in equation (3.10), it can be seen that it may be convenient to express the TEs in terms of CP-even and CP-odd departure-from-equilibrium matrices, which we define as respectively.Here, we use the approximate expression for the equilibrium value With all these definitions in place, we may find a set of TEs for the departure-fromequilibrium matrices, which crucially preserve the effects from changing dofs, which appear through the energy and entropy density of the thermal plasma In the TEs, the variations in the dofs are accounted for through the quantity The TEs are then written as ) In contrast with the majority of the TEs that appeared in the literature, the respective equations displayed here explicitly show the distinct phenomena that may impact the evolution of the heavy neutrino number densities, particularly the relativistic dofs and the changes in the equilibrium, η B eq .Finally, in our evaluations of the two differential equations in (3.21), we make use of the well-known expression for the Hubble parameter where M Pl ≃ 1.221 × 10 19 GeV is the Planck mass.

Lepton Asymmetry Transport Equation
As in the case of the heavy neutrino TEs, we identify the flavour-covariant TE for the lepton asymmetry in the manner described in [37].However, we additionally include the extra contributions appearing from the variations in the dofs.The TEs we use contain not only decay and inverse decay terms but also relevant 2 → 2 scatterings.Moreover, in order to avoid double counting of the decay processes, a proper subtraction of the RIS contributions is performed.The TE for the lepton asymmetry has been written in a manner such that collision decay terms do not appear, as outlined in [34,37].Consequently, the TE does not contain any negative contributions to the collision terms arising from the RIS subtraction of the ∆L = 2 interactions [45].
Since we are primarily interested in the evolution at temperatures higher than that of the critical temperature of the EW phase transition, the SM leptons are massless apart from thermally generated masses, and as a result, we neglect the oscillation of SM lepton species.We assume that the sum of lepton number density is close to thermal equilibrium, and so η L + η L ≃ 2η L eq 1.The flavour-covariant TE for the lepton asymmetry is then written as On the RHS of (3.23), the first two terms contain the decay contribution to the TE, with the rank-4 decay terms defined as The remaining terms contribute to the 2 → 2 scattering.The rank-4 scattering terms corresponding to LΦ → LΦ and LΦ → L c Φ † are defined as with the scattering matrices defined through a contraction over the final two indices Finally, the interactions with SM particles are accounted for in the decoherence matrices γ dec and δγ back dec .The former of these is a CP-even contribution, with γ dec = Γ L n L eq , and Γ L = diag(Γ L ℓ ), where as calculated in [46].In this expression, h L ℓ are the Yukawa couplings for the SM leptons, h Q t is the Yukawa coupling for the top quark, x = z M Φ (z)/m N , and M Φ (z) is the thermal Higgs mass.The CP-odd contribution may be calculated from the CP-even term as where the repeated index l is not summed over.
The total lepton asymmetry is easily extracted from δη L through its trace.Part of this lepton asymmetry is then re-processed into a baryon asymmetry through (B + L)-violating sphaleron transitions, which become exponentially suppressed at the sphaleron critical temperature T sph ≃ 132 GeV [47].In [48], the conversion factor for leptons-to-baryons is calculated to be Finally, in order to compare the generated BAU with the observed value, we should take into account the dilution of η B due to the expansion of the Universe.At temperatures below T sph , we assume that the total number of baryons is fixed, but the number density may be washed out due to the expansion of the Universe.Assuming that there are no considerable entropyreleasing processes between the sphaleron temperature and the recombination temperature, T rec , we may use entropy conservation to estimate with the overall conversion factor found to be around 1/27 [19,49].Therefore, the BAU, as measured at the recombination epoch, may be extracted from the TEs by

Impact of the Relativistic Degrees of Freedom
In order to estimate the reliability of our results in different regions, in this section we investigate how the results for the generated BAU change once the temperature dependence of the relativistic dofs is included.As shown in the literature [50,51], the contributions to the effective dofs are calculated through the energy and entropy densities of the Universe The inclusion of lattice [52] and perturbative [53] QCD effects gives additional contributions to the equation of state for the plasma.This results in an equation of state that deviates from the ideal gas assumption in a relevant manner [50] and, consequently, small variations in the dofs are present at temperatures above 100 GeV.
In the TEs, the variations in the dofs are accounted for in three places: (i) the Hubble rate, H, (ii) as a global factor of δ h on the commutator and collision terms, and (iii) as an additional term proportional to (δ h − 1).It is worth noting that these three factors do not affect our TEs in the same way.When the variations to the dofs enter into the Hubble rate, the T 2 scaling of H dominates over the small variations in g eff (T ).Similarly, the global factor δ h is overshadowed by the temperature dependence of the thermally averaged energy and decay matrices.However, the final term proportional to (δ h − 1) modifies the TEs in a significant way by introducing additional dynamics.In particular, this term provides an extra way in which the number density may be drawn away from equilibrium.With this in mind, we study the dynamics of the system in two initial-condition scenarios: (i) the case where the bath of heavy neutrinos begins in equilibrium, (ii) the case with a vanishing number density of heavy-neutrinos as the initial condition of the bath.

Equilibrium Initial Conditions
In the early evolution, we assume that thermal out-of-equilibrium effects generate the initial heavy neutrino asymmetry.Therefore, we take the initial condition ∆(z 0 ) = δ(z 0 ) = 0 when z 0 ≪ 1. Neglecting neutrino flavour oscillations and the collision terms, we may study the out-of-equilibrium dynamics alone.We consider the equation The first term in (4.2) gives the dynamics from the change in the equilibrium number density as the Universe is cooling down, while the second term describes the inclusion of changes in the entropy dofs.Taking the equilibrium initial condition, it may be shown that this equation has the solution In the above, we have defined the dimensionless parameter which will allow us to showcase better the relevant information as z > 1.In fact, the value of R provides important information on the dominant source of out-of-equilibrium dynamics.From (4.3), it can be seen that the sign of ∆ is equal to that of R.Moreover, from the definition of R, we can determine that when R < 0, the change in the entropy dofs becomes more significant than the reduction in the equilibrium number density, whereas the converse  is true for R > 0. Finally, since both h eff and η N eq decrease monotonically with z, the value of R must lie in the range [−1, 1].The extreme values represent scenarios with either constant dofs and a significant reduction in the equilibrium distribution (R → 1), or a scenario with minimal losses in the equilibrium number density but with a large reduction in the relativistic dofs (R → −1).In order to acquire a better understanding of the early behaviour of the number density, thanks to the linearity of the TEs, we may utilise R as a measure of the interplay between the two sources of out-of-equilibrium dynamics.
In Figure 1, we display the effect of the dofs on the TEs.Specifically, this figure shows the value of δ h − 1 for temperatures T > 100 GeV.If the dofs were constant, we would expect this term to vanish.As can be seen in these figures, the value of δ h − 1, albeit small, is nonvanishing, even as the temperature increases by multiple orders of magnitude.We consider two independent parameterisations of the dofs using data from Hindmarsh and Philipsen [50] (HP, 2005) and from Laine and Meyer [54] (LM, 2015).Both of these data sets are generated using lattice techniques, with the former data matched to two-loop perturbative computations and the latter using a dimensionally reduced effective field theory approach to compute corrections to thermodynamic potentials.In this context, one may notice that Figures 1(a  show where the washout of the number density due to changes in the relativistic dofs is the dominant contribution (R < 0).Solid lines indicate where the reduction in the equilibrium number density is the dominant effect (R > 0).Vertical lines show the critical temperature of the sphaleron transitions.The left panel is produced using data from Hindmarsh and Philipsen [50] (HP, 2005), and the right panel utilises the data given by Laine and Meyer [54] (LM, 2015).when 10 2 GeV ≲ T ≲ 10 4 GeV.To explore this difference, however, a dedicated study, which goes beyond the scope of this work, may be necessary.In Figure 1(a), the temperature dependence of the dofs is extracted from a discrete set (EOS C from [50]).Consequently, the value of δ h at each point was necessarily numerically estimated.The functional form is then generated from a cubic spline of these points, and we employ the derivative of this spline to find δ h .As such, some accuracy may be lost, particularly in regions where the data changes rapidly.In contrast, figure 1(b) is extracted from the data points provided in [54] using the analytic relationship where c is the heat capacity of the thermal bath, and s is the entropy density.A cubic spline is then taken of these data points to produce the continuous behaviour shown in Figure 1(b).
In Figure 2, we show the evaluation of R over the range of z pertinent to leptogenesis.We show four Majorana mass scales to represent scenarios at a very low scale and others above the electroweak scale.In the same figure, solid lines represent positive values of R, and dashed lines represent negative values.Also included are vertical lines in representative colours to show when the sphalerons become suppressed.As can be seen in Figure 2, for scenarios with Majorana masses well above the electroweak scale, we expect the variation in the dofs to have minimal influence on the dynamics.The simple reason is that at the sphaleron critical temperature, the value of R lies close to positive unity, i.e.R ≃ +1.This should be contrasted with low-scale scenarios, where R takes values closer to zero and may also be negative, i.e.R ≲ 0. The latter implies a greater interplay between the two sources of out-of-equilibrium dynamics.For 10 GeV Majorana neutrinos, we see that the freeze-out of entropy dofs is the primary out-ofequilibrium effect, which is expected to be reflected in the evolution of the number densities.If the Majorana mass scale is around 50 GeV, the variations in the entropy dofs have become the subdominant phenomena, but they still remain non-negligible.Therefore, we expect that the evolution of the neutrino number density has complex dynamics at this scale.Importantly, the R parameter we consider provides predictions for regions where the inclusion of dofs will be significant, in a way which does not rely upon computationally evaluated values of δ h .
The importance of variations within the relativistic degrees of freedom is not entirely unexpected.In the absence of interactions, the number density dilutes inversely proportional to the co-moving volume, i.e. n N ∝ a −3 .However, it is also known that the temperature dependence of the number density scales as n N ∝ T 3 .In the absence of entropy-releasing transitions, conservation of the total entropy may be used to find the condition In the standard paradigm, with constant degrees of freedom, this condition implies an inverse proportionality between the scale factor a and the temperature of the bath, T , reconciling the two scenarios outlined above.However, when variations in the relativistic degrees of freedom are included, the scale factor is related to the temperature through As a consequence, extra contributions are present, which provide another mechanism by which the number density may be moved out of chemical equilibrium.

Vanishing Initial Heavy-Neutrino Number Density
We now perform a similar analysis with the initial condition that η N (z 0 ) ≃ η N (z 0 ) ≃ 0 for z = z 0 ≪ 1. Hence ∆ ≃ −1 and δ ≃ 0. In such scenarios, it can be seen that the dominant contribution to the TEs comes from the inverse decay of leptons into heavy neutrinos, which drive the system into thermal equilibrium.Since the effects of the variations in the dofs are of significance when the bath has an appreciable abundance, the deviation in the BAU derived from variations in the dofs is directly related to how quickly the system can generate a considerable number density for the heavy neutrino species.
To make this last effect explicit, we first consider the initial evolution of the heavy neutrino number density.In addition, we may make some pertinent simplifications.Firstly, since the bath begins with η N ≃ 0, we may neglect the contributions to the TEs from variations in η N eq (z), and h eff (z).Similarly, the fact that the heavy neutrino number density is so far away from equilibrium means that we only need to consider inverse decays in the TEs since these will occur significantly more rapidly than the decay process.Finally, since the commutator term of the TEs does not increase or decrease the particle number, we may also neglect these terms.With these simplifications in mind, the TEs may then be written in the compact form Integrating the former of these expressions, we see that for z ≪ 1, For low values of z, the value of η N eq ≃ 1/ζ(3) is close to constant.We define the system as close to equilibrium by the condition Tr η N (z eq ) ≃ Tr η N eq (z eq ) .
Using this last expression and assuming a democratic flavour structure a = b = c = |h ν | and constant dofs, i.e. g eff (z = 1) ≃ 105, we may estimate the equilibrium value, z eq , to be From (4.11), we observe that the system will reach thermal equilibrium early in the evolution when the Yukawa couplings are large and the heavy neutrino mass scale is low.This result is to be expected since these two conditions give the most preferable conditions for the inverse decay of SM leptons and Higgs bosons into heavy neutrinos.In the TRL model, the neutrino-Yukawa couplings have been estimated to be of order 10 −4 [21], and therefore it is expected that 1 TeV heavy neutrinos will be close to equilibrium for z eq ∼ 10 −1 .Since the approach to equilibrium occurs faster for lower-scale heavy neutrinos, we may take this value of z eq as an upper bound.Consequently, we may conclude that in strong washout regimes such as the one we present, the variations in the dofs can still be relevant, even when the neutrino number density begins far from equilibrium.However, their effect may only be significant once an appreciable number density of heavy neutrinos has been generated through inverse decays.

Attractor Trajectories of the Transport Equations
From the discussion in Section 4.2, the rapid approach to equilibrium of the heavy neutrino number density poses the question as to whether the TEs we utilise exhibit strong attractor properties.To address this question, we first consider the unusual initial condition, η N (z 0 ) = 2η N eq 1.As opposed to the scenario given in Section 4.2, the decay of heavy neutrinos to SM leptons and Higgs bosons becomes now the dominant channel.Moreover, this over-abundance of heavy neutrinos also enhances the dynamics from the dofs.Taking these two sources into account for z ≪ 1, the heavy-neutrino transport equation may be approximated as follows: In order to assess whether or not the variations in the dofs are significant, we must compare the relative size of the two contributing terms that occur on the RHS of (4.12).In this way, we deduce the condition Taking a trace and making some appropriate approximations [49], we see that the decay of heavy neutrinos will dominate the evolution when As may be inferred from Figure 1(a), the value of the RHS of this last inequality typically lies in the range 10 −4 − 10 −3 for low values of z.Again assuming a Yukawa coupling scale of |h ν | ∼ 10 −4 and a heavy neutrino mass scale m N ≃ 1 TeV, we see that the LHS is dominant for z > 10 −2 .Consequently, we may conclude that the effect of the dofs is drastically enhanced when the neutrino number density is close to equilibrium.Otherwise, the asymmetry between the rate of decays and inverse decays is the dominant effect, pushing the bath toward the equilibrium.
By a similar argument to that given in Section 4.2, we may estimate again that z eq ∼ 10 −1 , but approaching from above rather than below.This indicates that the TEs we use may be highly attractive.Since we are now convinced that the variations in the dofs are only relevant close to chemical equilibrium, it is enough to verify that the system with no variations in the dofs exhibits such an attractor behaviour.
In all of our benchmark scenarios, we assume the Yukawa couplings to have a democratic flavour structure such that no SM lepton flavour is preferred, and flavour effects are minimised (i.e. a = b = c in (2.9)).In practice, this means that the model we adopt is similar to that of a single lepton generation.With this in mind, we study the attractor trajectories of a model with three singlet neutrino species and a single lepton flavour.In Figure 3, we show the evolution of such a model, assuming a tri-resonant mass spectrum with m N 1 = 1 TeV and a = b = c = 8.85 × 10 −4 .In addition, we do not stop the generation of baryon asymmetries at the sphaleron temperature and simply allow the freeze-out present in the figure to be achieved thermally.The colouring of the figure is determined through the quotient normalised to lie in [−1, 1].The values indicate not only the strength of the attraction but also its direction, with bluer hues identifying a suppression of η B and redder hues indicating an enhancement.From Figure 3, we see that the baryon asymmetry is driven upwards toward the attractor solution early in the evolution, z ≲ 10 −2 .Once the baryon asymmetry reaches the attractor solution, 10 −2 ≲ z ≲ 10, we see the most noticeable feature of this figure in the sharp transition from enhancement to suppression as we move across the attractor, highlighting the level to which this system of equations pushes the solution toward the attractive solution.
Hence, we find that this set of TEs is highly attractive and exhibits independence of the initial conditions.Notice that Figure 3 also identifies when thermal freeze-out occurs at z ∼ 10.Displayed through yellow shades, this region implies that |η ′ B | ≪ |η B |, and hence there is little change in the value of η B past this point in the evolution.
Finally, we identify that due to the inverse proportionality between z eq and the Yukawa coupling scale, the rapid approach to equilibrium of the heavy-neutrino number density only applies to strong washout models of leptogenesis.In the case of weak washout models, such as the freeze-in models studied in [55,56], the time taken to reach equilibrium is far longer.In such models, the out-of-equilibrium Sakharov condition relies on some unaccounted for preexisting dynamics that set the initial heavy-neutrino number density to vanish in an adhoc manner at a convenient moment of the cosmic time during the evolution of the early Universe.Therefore, due to the slow approach to equilibrium, the contributions from the relativistic dofs may be overshadowed, with the generated BAU showing little deviations away from a model with constant dofs.

Impact on the Baryon Asymmetry
In previous sections, we have shown how the variations in the dofs may be significant in many leptogenesis scenarios.In this subsection, we will now discuss how these effects can impact the predictions for the BAU, including the interplay of sphaleron interactions that enter through the sphaleron freeze-out temperature T sph .As previously noted, the largest contribution to the TEs from the variations in the dofs comes from an addition out-of-equilibrium term.This term washes out the number density through the freeze-out of the dofs.As a consequence, in regions where the variations in the dofs are dominant, we expect the neutrino number density to be under-abundant when compared to the equilibrium.Furthermore, since the transport equation for the lepton asymmetry is directly dependent on ∆, the lepton asymmetry follows a similar evolution.This leads to the less appealing feature that the generated baryon asymmetry turns suddenly negative, i.e. η B < 0, thus prompting us to resort to a contrived re-tuning of the CP phases of the model, often depending on the exact values of T sph and m N considered.
As was shown in Section 4.1, models with a Majorana mass scale below 100 GeV may be highly sensitive to variations in the relativistic dofs.Figure 4 displays how the generated BAU is altered by including the variations in the relativistic dofs.For this figure, the Yukawa couplings at each mass scale are taken such that the baryon asymmetry is reproduced by a system with no variations in the relativistic dofs, with democratic flavour structure assumed (a = b = c = |h ν ij |).These Yuakwa couplings are then utilised in a system with the variations in the dofs included to find the generated BAU. Figure 4 then shows the ratio between the two values of η B .For clarity, grey dotted lines are shown on the figure at zero and unity.Finally, in order to better gauge the significance of the sphaleron effects, we also vary the sphaleron temperature widely around the typical value T sph = 132 GeV.
As can be seen in Figure 4, when the Majorana mass scale is above 100 GeV, the generated BAU does not vary significantly from the constant dofs assumption.This occurs for two reasons.Firstly, as was discussed in Section 4.1, the lack of deviation occurs at these mass scales because the changes to η N eq become the dominant out-of-equilibrium effect, and we expect independence from the variations in the relativistic dofs.Secondly, once the heavy neutrinos are out of equilibrium, the inverse decay of SM leptons and Higgs bosons will begin to drive the neutrino number density back towards equilibrium.This latter process weakens the influence of the initial dynamics from the dofs.Therefore, the inclusion of the variations in the relativistic dofs in models with m N > 100 GeV effectively generates new initial conditions for the TEs.For comparison, we observe that for m N < 100 GeV, large deviations are apparent.In a model with such low Majorana mass scales, the washout effect from the relativistic dofs is more significant than the change in the equilibrium number density, producing a number density below the equilibrium value.Moreover, early in the evolution, the interactions with the SM leptons and Higgs bosons do not have the necessary time to restore the number density to its equilibrium value.Consequently, we may find that the baryon asymmetry at this scale drops below zero.
Figure 4 also indicates that regardless of the sphaleron critical temperature, the observed BAU drops below zero and takes on greater negative values as we approach lower mass scales.However, we point out that due to the complex interplay of the two sources of out-of-equilibrium dynamics, the deviations observed depend non-trivially on the sphaleron temperature.Due to the sign change in the generated BAU in low-scale mass regions, it is necessary to invert the CP phase in the Yukawa sector.In the Z 6 model, this is easily achieved by the replacement ω → ω * in the tree-level Yukawa couplings, which exchanges As already alluded to, depending on the parameterization for the dofs, there may be a complex interplay between the competing effects from the changes in the relativistic dofs, h eff , and changes in the equilibrium η N eq .Consequently, there may be a complicated dependence of the generated BAU on the sphaleron temperature.Figure 5 demonstrates this relationship when we use the parameterisation given in [50].To produce this figure, the Yukawa coupling necessary to reproduce the observed BAU at T sph = 132 GeV was found for each mass scale Figure 5: The dependence of the generated BAU on variations to the sphaleron temperature, T sph when the dofs are generated from [50].These figures show the ratio between the generated baryon asymmetry at a given sphaleron temperature, η B , and the generated baryon asymmetry when the sphaleron critical temperature is taken to be T sph = T * = 132 GeV, η T =T *

B
. The left panel shows this ratio for heavy neutrino masses above the electroweak scale, and the right panel shows the ratio for low-scale models below the electroweak scale.considered.This Yukawa coupling was then used to estimate the BAU generated in models with different values of T sph ∈ [125, 140] GeV, and the ratio between the two values of η B plotted.
Figure 5(a) shows the dependence of η B on T sph for high-scale Majorana masses.As expected, models at this scale show little dependence on the sphaleron critical temperature, with most scenarios showing less than ∼ 10% deviation away from the value of η B with T sph = 132 GeV.However, at low mass scales, shown in Figure 5(b), the generated BAU is significantly more sensitive to changes in the sphaleron critical temperature.Again, at these scales, the relevant contribution from the variations in the dofs produces a non-trivial relationship between the observed BAU and the sphaleron critical temperature.Here, we see that the observed deviations are significantly larger, up to a factor of 2.3 for m N = 50 GeV.In line with previous results, we also see that for lower-scale heavy neutrino masses, larger deviations are present, and in some cases, a negative value for the baryon asymmetry may be generated (indicated through the dashed line), which again may be rectified through the inversion of the CP phases.
Another source of potential deviation arises from the dilution factor, which comes from the conversion of the baryon asymmetry value at the freeze-out value at T sph to the observed value at T rec .As is shown in equation (3.31), the conversion comes from a ratio between the number of degrees of freedom at T sph and T rec , therefore a change in T sph can alter the predictions for the  The deviation in the washout factor, which converts the value of the BAU at T sph to the observed value at T rec , as a function of T sph .The value of this factor at T sph is taken as the reference value.The blue line is generated using the data from EOS C by Hindmarsh and Philipsen [50] (HP, 2005), and the red line is made using the latest data from Laine and Meyer [54] (LM, 2015).baryon asymmetry.In Figure 6, we show how this value of the conversion factor changes from the value at T sph = 132 GeV.As can be seen in the figure, when the parameterisation given in [50] is used, the variations in the predictions for the BAU vary by less than 1%.Conversely, for the data extracted from [54], these variations may increase to up to 6% at the more extreme values of the sphaleron critical temperature.
As a final point, it is worth re-iterating that the parameterisation of the variations in the relativistic dofs may be of significance.As was shown in [21], using a different data set for the dofs may affect the predictions for the BAU.In particular, when using [51], it was shown that the deviations had a lesser effect on the generated BAU, although notable features such as sudden parametric changes in the sign of η B are still present.In this work, we have made additional comparisons of the predicted Baryon asymmetry when using [50], which offers an improved data set through better modelling of the QCD quark-gluon plasma, or [54] which carefully considers distinct thermal regions, as well as a phenomenologically accurate Higgs mass.In general, these two data sets exhibit similar behaviour; however, the more recent data given in [54] seems to produce more robust predictions for the BAU with less dependence on the assumed value of the sphaleron critical temperature, as is highlighted in figure 4, demonstrating the importance of accuracy in the determination of cosmological parameters.

Critical Scenarios
In [22], the behaviour of unstable two-level quantum systems, famously called qubits, was discussed in detail.The dynamics of an unstable multi-level quantum system, also known as qudit in quantum information theory, is typically described by employing the Weisskopf-Wigner approximation (WW) [57].In the WW approximation, the states of a given quantum system evolve under a Schrödinger-like equation, with effective Hamiltonian Here, E and Γ are Hermitian matrices corresponding to the self-energy corrected energy matrix and the absorptive part of the self-energy correction to a multi-level particle system, respectively.Using this effective Hamiltonian, it may be shown that the density matrix of the multilevel system follows the evolution equation However, by taking the trace of this expression, it can be seen that this evolution equation loses probability as the states of the multi-level system decay to states not contained within the Hilbert space under consideration.As a result, the density matrix, which only contains the information pertinent to the states which have not yet decayed, is modified from the typical definition.This adjustment is done through the trace of the density matrix which we employ to define the co-decaying density matrix ρ ≡ ρ Tr ρ . (5. 3) The latter obeys the evolution equation for which it may be seen that the trace of the evolution equation vanishes provided Tr ρ = 1, implying conservation of the total probability of the co-decaying density matrix.
In the case of a two-level system, the matrices we consider may be expanded on the Pauli basis through four real-valued coefficients where , and Γ µ = (Γ 0 , Γ).Note that in this section, boldface identifies vector-like quantities rather than objects with flavour structure.
From these definitions, an equation which describes the motion of the Bloch vector, b, may be derived In this expression, two dimensionless parameters are introduced as well as the unit vectors In [22], it was shown explicitly that in the majority of cases, the system will evolve into a pure state with a stationary co-decaying Bloch vector, b ⋆ = b(τ → ∞), corresponding to the state with the longest lifetime.However, under special circumstances for which e • γ = 0 and r < 1, it may be seen that there is no preferred state of the qubit system since the lifetime for both states is identical.As a result, the co-decaying Bloch vector never reaches an asymptotic limit and oscillates between states indefinitely.Classes of models which satisfy the two conditions e • γ = 0 and r < 1 are referred to as critical scenarios and exhibit a number of interesting properties.In particular, one finds oscillations between coherent and decoherent states as well as the anomalous oscillations of the Bloch vector sweeping out unequal areas in equal time, in contrast with the typical Rabi oscillations for stable particle systems.It may be the case that the energy and decay vectors are close to orthogonal and r < 1, and so it may still be possible to observe some critical phenomena such as coherence-decoherence oscillations.However, in such cases, the system is expected to eventually relax into its longest-lived state.

Critical Scenarios in the Thermal Plasma
We now investigate whether similar phenomena can be present within cosmological settings.Comparing the TEs given in equation (3.21) with the evolution equation given in equation (5.2), we notice a significant overlap between the terms present.However, in the TEs, the matrix we are interested in would be that of the departure-from-equilibrium matrices, (∆, δ) [c.f.(3.17)], rather than the density matrix.Therefore, any phenomena we observe occur about the equilibrium number density.
In order for us to get a better insight of the relevant dynamics, it would be preferable to work with a single TE rather than two interdependent equations.Since we expect that interesting dynamics will occur about the equilibrium, we would like to retain the use of the matrices ∆ and δ.Consequently, we give a definition close to that of the distribution functions, f N and f N , which are linear combinations of the CP-even and CP-odd parts.To this end, we define in terms of which the following TE can be derived: For definiteness, we study the behaviour of the D + matrix.However, from the CP properties of ∆ and δ, we know that D + is related to D − through its transpose.In order to study the approach to equilibrium, we normalise the D + matrix by its trace Computing the TE for D+ , we find the expression where N is the number of heavy neutrino generations.From this TE, we can see that when the heavy neutrino number density is close to its equilibrium value, only the final term contributes in a significant way since the trace in the denominator enhances its contribution.As a result, the normalised matrix D+ must be close to a full mixed state Therefore, it would seem that the inclusion of terms which pull the number density away from equilibrium also drive the system toward a mixed state with no flavour preferences.
However, we notice that the remaining terms are similar to the evolution equation for a co-decaying Bloch vector.Therefore, we expect that in thermo-static backgrounds, it may be possible to observe some critical phenomena.Since we are now considering scenarios of constant temperature, the dimensionless parameter z is no longer fit for use.Consequently, we parameterise the evolution through the cosmic time.Starting from equation (3.10), we can find the evolution equation for the departure from equilibrium matrices as functions of time This expression can clearly be mapped onto Equation (5.2), with the substitutions After these replacements, we can expect the Bloch vector for the normalised matrix D+ to act analogously to the co-decaying Bloch vector of an unstable qubit.By following a similar procedure to that of the qubit system and taking definitions in line with those given in 5.6), we find that the r parameter for this system is given by In addition, the τ parameter satisfies the transformation The Z 2N models under consideration, with the assumed resonant mass splitting, naturally give energy and decay vectors which are close to orthogonal.We consider a Yukawa matrix with entries given by where ℓ i ∈ C are complex parameters and ω = e iψ is simply a phase factor, with ψ ∈ R. In this case, the tree-level decay matrix is given by (5.21) It is not difficult to see from (5.21) that the diagonal elements Γ αα of the tree-level decay matrix are all identical.Therefore, with a structure like this, the only way to induce differences in the diagonal elements is through higher-order corrections, such as through the mixing of heavyneutrino species.Consequently, we expect the coefficients of all SU(N ) diagonal generators of the decay matrix to vanish, and only the coefficients corresponding to off-diagonal SU(N ) generators will take non-zero values.For N = 2, we usually have that the component Γ 3 of the four-vector Γ µ is zero, but the other components Γ 1 and Γ 2 are non-zero [cf.(5.6)].Instead, the converse will hold true for the mass-matrix or energy components E 1,2,3 , where E 1,2 are zero, but E 3 is not.The latter is a simple consequence of the fact that the Majorana mass matrix is taken to be diagonal, and so all contributions to the energy vector E appear as coefficients of the diagonal generators of the SU(N ) group.Thus, for a model with N = 2 heavy Majorana neutrinos, this coefficient is proportional to the third Pauli matrix σ 3 .
It is important to note here that for models akin to the TRL Z 2N neutrino-Yukawa structure, the orthogonality of the vectors E and Γ is automatically satisfied by construction, i.e.E•Γ = 0.Moreover, computations of r in these models yield values that are very close to 1. Therefore, the TRL model can serve as a prototypal example of a non-diagonalisable Jordan-form model, details of which are discussed in [14,22] and Appendix A.
In Figure 7, we display the dynamics of the D + matrix.We write the D+ matrix in its Bloch decomposition and display the elements of the co-equilibrium Bloch vector, d. Figure 7(a) shows the evolution of the co-equilibrium Bloch vector in a standard cosmological setting where the Universe cools through its expansion.We use a model with two heavy neutrino species which have resonant mass splitting and Yukawa scales, |h ν | = 4.5 × 10 −4 .Since the use of Z 4 generators would lead to no lepton asymmetry, we use the Z 6 generator as the phase in the Yukawa coupling.In Figure 7(a), we observe that at the early phases of the evolution, the elements of the coequilibrium Bloch vector lie close to zero, implying the number density matrix is fully mixed with no preference in the states.But later in the evolution, one may notice some dynamics arising from the oscillation and decay terms of the TEs.The latter results from the out-ofequilibrium dynamics which becomes less significant for the evolution since the decay of heavy neutrinos will occur quickly when compared with the expansion of the Universe.In this late stage, we see that the co-equilibrium Bloch vector no longer describes a fully mixed state but instead asymptotically approaches a non-zero state.the energy and decay vectors, E and Γ, to be orthogonal to one another with r = 0.5.The evolution is shown in terms of the dimensionless run parameter, τ .Here, the dynamics of critical scenarios are clearly present, in which oscillations between states take place indefinitely and the co-equilibrium Bloch vector never relaxes into an asymptotic state.Moreover, we see that the oscillations occur only along the directions d 1 and d 2 , with the value in the final direction remaining constant, highlighting a defining property of critical scenarios in that the oscillations occur on a plane.
In the next section, we will give numerical estimates of the BAU predicted in the TRL model, including potential constraints from high-and low-energy observables of charged lepton-flavour and lepton-number violation.

Numerical Estimates
We now present numerical solutions to the TEs for the TRL model with democratic flavour structure a = b = c, and we make use of the dofs parameterisation given in [50].Additionally, we assume the central value of the sphaleron critical temperature, T sph = 132 GeV.Our analysis is limited to heavy neutrino masses above 40 GeV since, below this scale, the inclusion of thermal corrections to the mass of SM particles leads to phase-space suppression of the decay processes.As a result, the generation of lepton asymmetries must be treated with additional care in regions of the parameter space with very low-scale masses for the heavy neutrinos.Nevertheless, the results presented here will still be valid in order to describe the associated theoretical uncertainties in these regions of the parameter space.In Figure 8, we give representative numerical evaluations of the heavy neutrino TEs, with the departure of η N from the equilibrium value η N eq displayed.To carry out these evaluations, we take the masses to be in consecutive resonance with Z 6 symmetric Yukawa structure.In addition, we assume that the heavy neutrino number density is initially in equilibrium, ∆(z 0 ) = δ(z 0 ) = 0, and δη L (z 0 ) = 0, where z 0 = 10 −2 .In these two figures, the vertical dotted line indicates the point in the evolution where the sphaleron transitions become exponentially suppressed.More explicitly, Figure 8(a) shows the evolution of the singlet neutrinos when the lightest singlet neutrino is of mass 1 TeV and the Yukawa coupling scales are |h ν ij | = 2.95 × 10 −4 .As was discussed in Section 4.1, at early phases in the evolution, z < 10 −1 , we observe that variations in the dofs are the dominant phenomena.As a result, the entries of the neutrino matrix take negative values.However, once z > 10 −1 , we see that the decay and inverse decay of the heavy neutrinos push the solution towards the attractor value, and hence the effect on the generated BAU from the dofs is minimal.
In Figure 8(b), we now illustrate how the above picture changes for light singlet neutrino masses.For this figure, the lightest singlet neutrino has mass m N 1 = 50 GeV, and the neutrino-Yukawa scale is |h ν ij | = 3.1×10 −4 .In this case, the variations in the dofs have a large impact on the evolution of the different entries of the matrix ∆, producing sudden changes that arise from competing effects between the changes in the relativistic dofs and the changes in the heavyneutrino number density.Consequently, this non-trivial evolution of ∆ can source significant uncertainties when accurate determinations of the theoretical parameter space of such low-scale heavy-neutrino models are attempted based on successful leptogenesis.
In Figure 9, we demonstrate how the above effects feed through into the generated BAU.
Again, for these figures, we assume a TRL mass splitting and a Z 6 symmetric Yukawa matrix with a democratic flavour structure.The initial conditions are again chosen assuming that the heavy-neutrino number density is initially in equilibrium, so ∆(z 0 ) = δ(z 0 ) = 0, and δη L (z 0 ) = 0 with z 0 = 10 −2 .The vertical dotted line indicates the sphaleron critical temperature and the horizontal line gives the CMB central value of the baryon asymmetry η B = 6.104 × 10 −10 .For Figure 9(a), we assume the lightest singlet neutrino has mass m N 1 = 1 TeV and |h ν ij | = 2.95 × 10 −4 to align with the inputs of Figure 8(a).We see that in this high-mass regime, the generated BAU is independent of the variations in the dofs at the start of the evolution, and therefore, the attained values are protected from such phenomena.Moreover, this protection means that any variations within the sphaleron temperature are minimal, as was discussed in Section 4.4.
Let us now turn our attention to Figure 9(b), where we display the evolution of the BAU when m N 1 = 50 GeV and |h ν ij | = 3.1 × 10 −4 .As before, this figure should be considered in line with Figure 8(b).In drastic contrast with heavier neutrino models, the generated BAU at this scale is heavily dependent upon the variations in the dofs.As can be seen from Figure 9(b), there are dramatic variations in the BAU whose value changes sign on multiple occasions.In addition to this less controllable situation, the fact that the sphaleron transitions become suppressed at a significantly lower value of z means that the decays and inverse decays no longer have the time to dominate within the evolution.For such low-scale TRL scenarios, it would require a drastic re-tuning of the mass parameters and CP phases to obtain predictions that lie close to the observed BAU.Hence, we consider that such TRL scenarios are not predictive in standard cosmological settings.
We now delineate the parameter space for the TRL model and make some pertinent comparisons between this model and special cases of interest.Figure 10 shows the regions of the parameter space which may lead to successful leptogenesis.The solid green line displays the predicted values of the light-to-heavy mixing, B ≃ m D m −1 M [cf.(2.5)], which should be used to match the observed value for the baryon asymmetry η CMB B given in (1.1).The green-shaded area under the solid line represents regions of the parameter space where the generated BAU exceeds the observed value.In this area of parameter space, the observed baryon asymmetry may still be matched through some additional breaking of the Z 6 symmetry or by a softening of the tri-resonant condition, leading to sub-maximal levels of CP-asymmetry.Alongside the full TRL model, we show the adjusted parameter space when we take the dofs to be constant and another for a model with no oscillation effects.The former of these is shown in Figure 10 as brown dashed lines and exhibits a clear departure from the TRL model with fully implemented dofs when in the sub-100 GeV regime.For mass scales in excess of 100 GeV, it is apparent that the two lines overlap, in line with the discussion given in Section 4.4.Once the TEs are fully implemented, the variations in the dofs prevent the generation of baryon asymmetries in the low mass regime.In order to compensate for this, we must reduce the washout of any generated asymmetries, necessitating smaller neutrino-Yukawa couplings.This is then seen in  the parameter space as a drop-off in the available mixing.
In Figure 10 we have also presented the parameter space for the TRL model for which the number density evolution was taken to follow a set of diagonal Boltzmann equations.This is indicated with brown dotted lines in Figure 10.From this figure, we see that this line follows that of the full model, aside from a global factor, which may be numerically estimated to be around two [37,63].It would then appear that the omission of oscillation effects from the evolution seemingly permits additional CP asymmetries whilst maintaining agreement with the observed baryon asymmetry.It also implies that the CP asymmetry associated with oscillation effects interferes destructively with the CP asymmetry associated with singlet neutrino mixing in the decay [64].Such a result is not entirely unexpected.As has already been discussed, the Z 6 structure of the Yukawa couplings, when combined with a tri-resonant mass hierarchy, provides large CP phases and saturates the available mixing CP-asymmetry.Therefore, we again must reduce the scale of the Yukawa couplings to ease some of the washout effects and generate the necessary baryon asymmetry.
Figure 10(a) shows detection limits on the parameter space from charged Lepton Flavour Violating (cLFV) transitions involving muons.Evidently, the experimental bounds for many experiments still lie far above the parameter space in which we would expect to see cLFV phenomena.However, a small partion of the parameter space may be probed by the PRISM experiment in their future studies of coherent transitions within Titanium.In Figure 10(b), the bounds from collider experiments are displayed.The blue dashed line gives the projected sensitivity bound for the LHC running at √ s = 14 TeV with integrated luminosity L I = 300 fb −1 , for the process p p → N ℓ ± j j [62,65].The orange line gives the estimated sensitivity of the process Z → N ν to a 95% confidence level from the DELPHI collaboration [66].Moreover, through the decay of Z bosons to charged leptons, it may still be possible to observe some LFV processes [67], although the deviations are expected to be small in TRL models with Z 6 Yukawa structure [21].As may be seen in Figure 10(b), many of these experiments do not appear to have the required sensitivity to be able to probe the region of the parameter space which can achieve successful leptogenesis.However, there may be some future experiments, such as the Future Circular Collider (FCC), which can probe the parameter space for masses below 50 GeV [68].

Conclusions
We have studied how the inclusion of relativistic dofs influences the predictions for the BAU in the context of a low-scale Tri-Resonant Leptogenesis model whose neutrino-Yukawa couplings are dictated by a Z 6 symmetry.The TRL model, which we presented in [21] and reiterated here, is distinctive in that it can produce a vanishing SM neutrino mass spectrum whilst still maintaining large CP phases for the successful generation of the BAU.By considering resummed neutrino-Yukawa couplings, the available CP violation is maximised when the heavy neutrino mass spectrum is taken to be in consecutive resonance, as was detailed in a previous work by the authors in [21].
The effect of the relativistic dofs was analysed by introducing a set of modified Transport Equations, which include additional contributions from variations within the relativistic dofs.These additional contributions tend to wash out the number density of the heavy neutrinos and drop their number density below its equilibrium value.When this enters the TE for the lepton asymmetry, the generated lepton asymmetry will have the opposite sign to that expected in standard thermal leptogenesis models.At higher scales for which m N > 100 GeV, the variations in the dofs become a sub-dominant effect, and so thanks to their attractor properties, the TEs can return the generated lepton asymmetry to the expected value in the standard paradigm.
At low scales, however, where m N < 100 GeV, the (B + L)-violating sphaleron transitions, which produce the baryon asymmetry, become exponentially suppressed at a point where the generated lepton asymmetry carries the incorrect sign.This gives rise to the unpleasant feature that a negative BAU gets predicted unless the CP phases in the neutrino Yukawa sector are appropriately re-tuned, according to the sphaleron temperature T sph assumed, the parameterisation of dofs adopted, and the specific value of m N considered.For regions of parameter space where these re-tunings become excessive, such low-scale TRL models, as well as similar low-scale leptogenesis scenarios, lose their inherent predictive power.
We highlighted that in strong washout regimes, there are always strong attractor trajectories towards the equilibrium.This means that the generated BAU exhibits independence from the initial conditions of the pertinent number densities.Hence, variations in the dofs play a prominent role regardless of the initial conditions of the heavy-neutrino number density.Since the dofs phenomena we have studied here contribute an out-of-equilibrium effect, the observed departures from the standard paradigm may be present in many low-scale freeze-out mechanisms, including those to low-scale baryogenesis scenarios and to the generation of a light Dark Matter relic density.However, for weak washout models that rely on the freezein framework [69], the thermal bath of heavy neutrinos approaches equilibrium slowly since the initial number density vanishes; presumably through some mechanism to be specified that provides low-reheat temperature of a rather convenient scale, not too far above T sph .Apart from these strong assumptions on the initial cosmological setting, any deviations in the BAU from the variations in the dofs are expected to be small.Along the lines of our work in [22], we have extended this previous study to critical-like scenarios within a quasi-thermal bath.In particular, we have demonstrated how, in systems in which the bath is drawn out of equilibrium, the critical phenomena of interest to us may be destroyed.However, in thermo-static scenarios, there still exists the possibility that anomalous behaviours, such as coherence-decoherence oscillations, can be present.
We have delineated in Figure 10 the parameter space for light-to-heavy neutrino mixings in the TRL model which can successfully reproduce the observed BAU.We see that the inclusion of variations in the relativistic dofs results in a suppression of the light-to-heavy neutrino mixings for heavy-neutrino masses below 100 GeV.Moreover, we find that the heavy-neutrino oscillation effects contribute destructively to the generation of the BAU.Specifically, the diagonal TEs give rise to a bigger parameter space, which is increased by a factor close to two.The parameter space shown in Figure 10 indicates that many current and future experiments may still be far away from detectable phenomena.However, there still exists some possibility of finding coherent µ → e transitions within Titanium at the PRISM experiment.Simple extensions to the TRL model, such as supersymmetry, may offer additional sources of CP violation or additional contributions to the light-to-heavy neutrino mixing.These can bring the predictions for most new-physics observables close to their current experimental limits while achieving successful leptogenesis.

Figure 1 :
Figure 1: The deviation away from unity of the δ h parameter.Observe that at temperatures well above the electroweak mass scale, the relativistic dofs are almost constant, δ h → 1.As the temperature approaches this scale, SM particles (e.g. the top quark) start to decouple, and δ h decreases.In this figure, we consider two datasets.The left panel uses data extracted from EOS C by Hindmarsh and Philipsen [50] (HP, 2005), and the right panel shows the latest data from Laine and Meyer [54] (LM, 2015).Solid lines express positive values of δ h − 1, whereas dashed lines indicate negative values.

Figure 2 :
Figure 2: This figure indicates the dominant out-of-equilibrium effect in the TEs.Dashed linesshow where the washout of the number density due to changes in the relativistic dofs is the dominant contribution (R < 0).Solid lines indicate where the reduction in the equilibrium number density is the dominant effect (R > 0).Vertical lines show the critical temperature of the sphaleron transitions.The left panel is produced using data from Hindmarsh and Philipsen[50] (HP, 2005), and the right panel utilises the data given by Laine and Meyer[54] (LM, 2015).

Figure 3 :
Figure 3: The attractor properties of the numerical solution to the TEs.Red shades indicate regions where the baryon asymmetry is increasing, and blue shades indicate regions where the baryon asymmetry is decreasing.Yellow shades indicate regions where the changes to η B are small.

TFigure 4 :
Figure 4: The ratio between the baryon asymmetry generated in a model with variations in the dofs, η δ h B and the baryon asymmetry generated in a model where the dofs are constant (h eff = 105), η δ h =1 B .Different colours represent models where the sphaleron temperature takes different values.Two horizontal dotted lines are included to show unity and zero.The left panel is produced using data from Hindmarsh and Philipsen [50] (HP, 2005), whereas the right panel utilises the latest data from Laine and Meyer [54] (LM, 2015).

Figure 6 :
Figure6: The deviation in the washout factor, which converts the value of the BAU at T sph to the observed value at T rec , as a function of T sph .The value of this factor at T sph is taken as the reference value.The blue line is generated using the data from EOS C by Hindmarsh and Philipsen[50] (HP, 2005), and the red line is made using the latest data from Laine and Meyer[54] (LM, 2015).

Figure 7 :
Figure 7: The prevalence of critical phenomena in two different scenarios.The left panel shows the components of the co-equilibrium Bloch vector, D+ , in a cosmological setting where outof-equilibrium phenomena are present.The evolution is given in terms of the cosmological parameter z, since the cosmic time is not fixed through the expansion of the Universe.The right panel shows the components of the co-equilibrium Bloch vector in a thermo-staticUniverse where the temperature is kept constant.

Figure 7 (Figure 8 :
Figure 7(b) shows the dynamics in a thermo-static background.For this figure, we take

Figure 9 :
Figure 9: The evolution of the BAU at two points in the parameter space.Dotted lines show where the values are negative, and solid lines indicate positive values.In order to make pertinent comparisons, the left and right panels align with Figures 8(a) and 8(b), respectively, with identical initial conditions and model parameters.

Figure 10 :
Figure 10: The parameter space for the TRL model with limits from current experiments (solid lines) and projected sensitivity limits (dashed lines) on future experiments.The left panel exhibits the bounds on the parameter space for coherent Lepton Flavour-violating transitions.The red dashed line gives projected bounds on µ → eγ [58], the orange dashed line shows µ → eee [59], the blue solid line indicates the current bound on cLFV transitions within Gold [60], and the blue dashed line gives the projected bound on cLFV transitions within Titanium [61].The right panel presents similar bounds for collider experiments.The blue dashed line identifies the projected sensitivity for the LHC running at √ s = 14 TeV [62].The orange line specifies the current bound on Z-boson decays to singlet neutrinos at DELPHI.Finally, we include brown dashed lines to indicate the parameter space when δ h = 1, and a brown dotted line which indicates the parameter space in a flavour diagonal TRL model.