Impact of high-scale Seesaw and Leptogenesis on inflationary tensor perturbations as detectable gravitational waves

We discuss the damping of inflationary gravitational waves (GW) that re-enter the horizon before or during an epoch, where the energy budget of the universe is dominated by an unstable right handed neutrino (RHN), whose out of equilibrium decay releases entropy. Starting from the minimal Standard Model extension, motivated by the observed neutrino mass scale, with nothing more than 3 RHN for the Seesaw mechanism, we discuss the conditions for high scale leptogenesis assuming a thermal initial population of RHN. We further address the associated production of potentially light non-thermal dark matter and a potential component of dark radiation from the same RHN decay. One of our main findings is that the frequency, above which the damping of the tensor modes is potentially observable, is completely determined by successful leptogenesis and a Davidson-Ibarra type bound to be at around $0.1\;\text{Hz}$. To quantify the detection prospects of this GW background for various proposed interferometers such as AEDGE, BBO, DECIGO, Einstein Telescope or LISA we compute the signal-to-noise ratio (SNR). This allows us to investigate the viable parameter space of our model, spanned by the mass of the decaying RHN $M_1 \gtrsim 2.4\times 10^8\;\text{GeV} \cdot \sqrt{2\times 10^{-7}\;\text{eV}/\tilde{m}_1}$ (for leptogenesis) and the effective neutrino mass parameterizing its decay width $\tilde{m}_1<2.9\times 10^{-7}\;\text{eV}$ (for RHN matter domination). Thus gravitational wave astronomy is a novel way to probe both the Seesaw and the leptogenesis scale, which are completely inaccessible to laboratory experiments in high scale scenarios.

1 Introduction The standard model (SM) of particle physics predicts that the neutrinos are massless, but due to the observation of neutrino oscillations for solar [1][2][3][4][5][6][7][8], atmospheric [9,10] and reactor [11][12][13][14] neutrinos we now know that they are massive and the flavor states mix due to the propagation of multiple mass eigenstates. Moreover the β-decay experiment KATRIN [15] has provided us with the first direct limit of the neutrino mass scale m ν < 0.8 eV.
Cosmology offers an indirect probe of this scale and demands that the sum of all neutrino masses satisfies i m ν i < 0.12 eV [16,17] in order to be consistent with the predictions for the Cosmic Microwave Background (CMB) radiation, Large scale structure (LSS) formation and Big Bang Nucleosynthesis (BBN). The accelerated expansion at the beginning of the universe provided by cosmic inflation, which was postulated in order to solve the horizon and the flatness problems and is responsible for quantum generation of the primordial fluctuations seeding the large scale structure of the universe, is thought to be driven by a scalar field known as the inflaton (see [18] for a review). In this paper, we will be concerned with the primordial Gravitational Waves (GW) background of such inflationary origin [19][20][21] (see [22] for a review on this topic). These inflationary GWs can act as a logbook of the expansion history of our universe throughout its entire evolution [23][24][25][26][27][28][29][30]. Particularly, the detailed time evolution of the Hubble rate during the expansion determines the transfer function that describes how gravitational waves at different frequencies are red-shifted to the present day. This property turns primordial GWs into a powerful tool that grants access to the thermal history of our universe prior to BBN. Primordial GWs offer, e.g. an opportunity to measure the reheating temperature after inflation [31][32][33][34][35][36][37][38]. Similarly, with help of these inference can be drwan of the equation of state during the quark-hadron phase transition in quantum chromodynamics [39,40] or constrain properties of the hidden sectors beyond the Standard Model (BSM) of particle physics [41,42].
The observed baryon asymmetry of the universe (BAU) is longstanding puzzle in particle physics and cosmology [16,43]. While the universe is expected to start in a matterantimatter symmetric phase, any primordial asymmetry set due tothe initial conditions is expected to get diluted by the exponential expansion phase during cosmic inflation. The BAU is often quoted in terms of the baryon to photon ratio measurement which, according to the latest Planck 2018 data, is given by [16] η B = n B − nB n γ = 6.1 × 10 −10 (1. 1) and agrees with the value extracted from BBN [44] as well. Similar to the BAU, there has been another question related to the presence of a mysterious, non-luminous form of matter, popularly known as dark matter (DM), giving rise to approximately 26% of the energy density in the present universe. In terms of density parameter Ω DM and h = H 0 /(100 km s −1 Mpc −1 ) with H 0 being the observed present day Hubble parameter, the current DM abundance is conventionally reported to be [16] Ω DM h 2 = 0.120 ± 0.001 (1.2) at 68% CL. Apart from cosmological evidence, the presence of DM has also been suggested by several astrophysical implications [45][46][47]. While none of the standard model particles satisfy the criteria of a particle DM candidate, the SM also does not to satisfy the criteria to dynamically generate the observed BAU, known as Sakharov's conditions [48], in adequate amounts. This has led to several BSM possibilities offering intriguing solutions to these puzzles: The Type I Seesaw mechanism [49][50][51][52][53][54], where the SM is augmented with three right handed SM gauge singlet neutrinos (RHN), may explain both the observed neutrino masses (from neutrino oscillation experiments) as well as the baryon asymmetry of the universe via first generating an asymmetry in the dark leptonic sector [55][56][57][58][59] and subsequently getting transferred to the visible baryonic sector via the electroweak sphaleron transitions [60]. Among the BSM proposals for DM, the weakly interacting massive particle (WIMP) [61] produced as a thermal relic is perhaps the most widely studied one (see [62] for a review). However due to the absence of any WIMP related signals in nuclear and electron recoil DM direct detection experiments, there has been growing interest in other (non-thermal) production modes: some examples are the well-known super-WIMP scenario [63], where frozen out WIMP decays to the actual DM, FIMPs [64] (see [65] for a review) that have such tiny couplings to the SM plasma that they never thermalize, or non-thermal production from inflaton decays [66] during the process of the formation of the radiation bath known as reheating. In leptogenesis models the RHN might also have the decay modes to other SM singlets that can be good DM candidates [67,68], which is why we will adopt this framework. Since the RHN decays out-of-thermal equilibrium the DM will be non-thermal.
We will demonstrate that the same RHN decay responsible for both the generation of the primordial baryon asymmetry via leptogenesis, as well as the production of non-thermal dark matter and a possible component of dark radiation, leaves its vestige on the primordial spectrum of inflationary GWs. In particular we consider an epoch of intermediate matter domination [61,69,70] from the lightest RHN, which decouples from the plasma while relativistic and is very long-lived compared to the characteristic time scale of the cosmic expansion. Since the decay occurs far away from thermal equilibrium it will release a large amount of entropy, which dilutes the energy density of primordial GWs that enter the horizon before the decay.
Although the Seesaw mechanism ties leptogenesis to the observed light neutrino masses, the mechanism itself is notoriously difficult to test in laboratory based experiments, as the heavy right-handed neutrino mass scale has to be above ≳ 10 9 GeV (see [71]). One should keep in mind that this bound can be evaded, see for example [72] and with some fine tuning, it is also possible to bring down the scale of the non-resonant thermal leptogenesis to as low as 10 6 GeV [73]. However indirect tests for high scale leptogenesis of course exist as well. These are primarily based on neutrino-less double beta decay scenarios [74,75], meson decay scenarios [76][77][78], and via CP violation in the neutrino oscillation [79,80], the structure of the leptonic mixing matrix [81], or via considering theoretical constraints from the demand of the SM Higgs vacuum does not become unstable in early universe [82,83]. Therefore, it is necessary, although very challenging to find newer and complementary tests of such heavy neutrino seesaw physics and consequently the leptogenesis mechanism. Recently it has been proposed to complement these indirect tests with the observations of GWs of primordial origin such as that from cosmic strings [84], domain walls [85] and other topological defects [86] or from nucleating and colliding vacuum bubbles [87,88], graviton bremmstrahlung [89] and primordial black holes [90,91]. These previous studies on GW [35,84,86,[92][93][94][95][96][97][98] focused on the stochastic GW background from the dynamics of the scalar field, whose vacuum expectation value is responsible for the RHN mass, whereas (when it comes to leptogenesis) we only extend the SM by adding nothing more than three RHNs with hard mass terms. In order to ensure a thermal population of the lightest RHN, which can not be established by the Yukawa couplings we consider, we have to assume that the RHNs are produced from inflaton decays or additional gauge interactions. In this paper we propose the imprint of the RHN decay on the inflationary first-order tensor perturbations as a novel probe of the minimal high-scale leptogenesis mechanism.
The paper is organized as follows: In the subsection 2.1 of section 2 we discuss the Seesaw model, then how the decay of the lightest right handed neutrino (RHN) leads to an intermediate era of matter domination in 2.2, and we elaborate on the generation of baryon asymmetry via leptogenesis from the decay of the lightest RHN in 2.3. We also discuss the production of non-thermal dark matter and dark radiation from such heavy RHN decays in 2.4. In section 3 we discuss the generation and propagation of inflationary tensor perturbations as Gravitational Wave signals and show how RHN decays leave their imprint on the GW spectrum. We discuss the GW detection prospects in 4.1 of section 4 and translate such experimental sensitivities into the reach for probing the parameter space and scale of leptogenesis via computing the signal-to-noise ratio (SNR) in 4.2. We end with the conclusions in section 5.
2 Decays of a long-lived RHN

Type I Seesaw mechanism
We start with a conventional Type I Seesaw [49][50][51][52][53][54] with three right handed neutrinos N where σ 2 is the second Pauli matrix and assume without loss of generality that the symmetric right handed neutrino (RHN) mass matrix is diagonal without making any assumptions about the mass spectrum yet. After Integrating out the RHN and electroweak symmetry breaking with ⟨H⟩ ≡ v = 174 GeV the active neutrino mass matrix reads at leading order in the Seesaw expansion Using the Casas-Ibarra parameterization in the basis where the charged lepton mass matrix is diagonal one finds [99] where U PMNS is the leptonic equivalent of the CKM matrix. R describes the mixing and CP-violation in the RHN sector and is expressed as a complex, orthogonal matrix that reads in terms of 2 × 2 rotation matrices R (ij) in the ij-plane with an angle z ij .

Conditions for intermediate matter domination
The lightest RHN N 1 has the tree level decay width summed over all SM lepton flavours of For T ≫ M j the decay in the plasma is suppressed by a time dilation factor of M 1 /T [100], which goes to one for T ≤ M 1 . It is customary to define the effective neutrino mass mediated by N 1m which appears when comparing the decay rate to the characteristic time scale of cosmic expansion H(T ) −1 , where H(T ) is the Hubble rate during radiation domination This effective mass only coincides with the physical mass (m j = m j ) for R ji = 0, ∀i ̸ = j. A small effective massm 1 implies that N 1 is weakly coupled to other two RHN. One can show that this effective mass is larger than the lightest active neutrino mass [101] We find that the N 1 decays after it has become non-relativistic (K 1 ≪ 1) as long as The energy density of the non-relativistic RHN redshifts slower than radiation, so it overtakes the radiation component and becomes the dominant contribution to the energy budget of the universe at [102] T dom. = 7 4 where we used that the number of relativistic degrees of freedom above the electroweak crossover is g * (T dom. ) = O(100). Once   . (2.12) The decay takes place after the onset of early matter domination for [102] m 1 < 2.9 × 10 −7 eV. (2.13) Ifm 1 is larger than this number, there will be no era of intermediate RHN matter domination and consequently the decays of the N 1 will not produce enough entropy to lead to an appreciable dilution of the inflationary tensor mode background (see the following discussion in section 3.1). This bound implies together with (2.9) that the lightest active neutrino mass has to be smaller than 2.9 × 10 −7 eV meaning that for normal-ordering (NO) we consider the following neutrino spectrum [103] m 1 ≃ 0, m 2 ≃ ∆m 2 sol. ≃ 8.6 × 10 −3 eV, m 3 ≃ ∆m 2 sol. + ∆m 2 atm. ≃ 0.05 eV. (2.14) For the inverted ordering (IO) we would instead have a quasi-degenerate spectrum [103]   where we used that during matter domination a ∼ H −2/3 together with H(T dec. ) = Γ 1 and H(T dom. ) ∼ T 2 dom. /M Pl. at the transition from radiation to matter domination. Throughout this work we assume an initial equilibrium distribution for N 1 . For small Yukawa couplings giving rise tom 1 < 10 −3 eV [102] the interactions in (2.1) do not suffice to establish equilibrium in the radiation dominated plasma after inflationary reheating at T RH . Hence our scenario precludes thermal leptogenesis and is sensitive to the initial conditions of the radiation bath. This is why we assume the initial population of RHN is produced by additional interactions such as couplings to the inflaton φ [107] like e.g.
for a production during reheating, or new gauge bosons from e.g. GUTs [108,109] or gauged B-L [70]. Concentrating on the case of a U(1) B-L gauge boson with mass m Z ′ = g B-L v B-L > T RH as an example, the scattering rate of N 1 with the SM quarks and leptons via off-shell Z ′ would read approximately This interaction freezes-out while the N 1 are still relativistic (T FO > 10M 1 ) as long as v B-L > 7 × 10 11 GeV · M 1 7.5 × 10 8 GeV (

2.22)
The impact of the underlying U(1) B-L breaking on stochastic GWs is briefly explained in section 3.2.

Non-thermal leptogenesis
We assume the inflationary reheating dynamics satisfy M 2 , M 3 > T max > M 1 so that we can focus on the decays of the lightest RHN N 1 . In this context we defined T max > T RH as the largest temperature during the epoch of inflationary reheating [110][111][112], which ends with a radiation bath of the temperature T RH . Alternatively, if one assumes only M 2 , M 3 ≳ (3 − 10) × M 1 , the population of N 2,3 will have decayed away long before N 1 decays, as a consequence of their larger Yukawa couplings needed to explain the observed neutrino masses. Further we assume there is no primordial lepton asymmetry e.g. from the decays of N 2,3 . Since the N 1 are too weakly coupled, they would not be able to erase this preexisting asymmetry [113]. However for realistic light neutrino masses the N 2,3 will be in the strong washout regimem 2,3 > 10 −3 eV, so that inverse decays LH → N 2,3 destroy a large portion of the asymmetry produced by the decays of N 2,3 . The lepton asymmetry n B-L /s, defined in terms of the number density of leptons minus anti-leptons normalized to the entropy density s, can be converted into a baryon asymmetry via the electroweak sphaleron process. For the RHN dominated scenario one finds a baryon asymmetry of [102] The parameter ε 1 denotes the CP-violating decay parameter encoding the amount of leptonic asymmetry produced per decay of N 1 . The sphaleron redistribution coefficient is found to be c ph. = 28/79 [114] and the term ω, that will be determined later in this paragraph, parameterizes the washout of the lepton asymmetry. Our analysis is different from the more commonly studied case of non-thermal leptogenesis immediately after inflationary reheating [115,116], where T dec. /M 1 would have to be replaced with T RH /m φ with m φ being the inflaton mass, because here the RHN decay takes place much later, after it had time to dominate the energy budget of the universe. The factor of T dec. /M 1 < 2% comes from n N /s, which can be obtained from energy conservation (ρ tot. = M 1 n N before the decay) leading to and can be understood as the entropy dilution from the N 1 reheating: The dimensionless dilution factor from the entropy produced by the instantaneous 1 out-of-equilibrium decay of the dominating RHN N 1 [61,69,70] reads In this context we denote the average of g * (T ) over the decay period as ⟨g * (T )⟩ and we assume that ⟨g * (T )⟩ ≃ g * (T dec. ). To obtain the second line we assumed for the initial abundance n i N /s that N 1 decoupled from the plasma while relativistic to maximize the amount of entropy produced [70], see also (2.22). For hierarchical RHN spectrum (M 3 > M 2 > M 1 ) the decay parameter from the interference between tree-level and one-loop vertex-and self-energy-corrections is found to be [118] where the upper limit (for normal ordered neutrino masses) reads [119] ε max = 3 16π It is worth mentioning that while the small required value ofm 1 in (2.12) necessitates small values of |R 1i | 2 , this does not automatically force |ε 1 | hier. to be tiny, since this quantity depends only on a ratio of squared R-matrix elements. For completeness let us mention that for a degenerate spectrum with M 3 > M 2 ≃ M 1 the self-energy graph gets resonantly enhanced and the estimate gets modified as [118] (2.29) We estimate the baryonic asymmetry for a general value of ε 1 where we chosem 1 for matter domination according to (2.13). One can compute the observed n B /s from the baryon-to-photon-ratio in (1.1) by making use of s ≃ 7.04 n γ . The required mass M 1 for the hierarchical spectrum can be obtained from (2.28) and depends intimately on the details of the active neutrino mass spectrum. Note that unlike the usual Davidson-Ibarra bound M 1 ≳ 10 9 GeV [119] our estimate depends on the parameterm 1 due to the entropy produced in the RHN decay. It is not surprising that this bound can be slightly lower than the Davidson-Ibarra limit, as the out-of-equilibrium RHN abundance at T dec. can be larger than the typically assumed relativistic thermal yield. Fitting M 1 ,m 1 to the baryon asymmetry of the universe leads to T dec. ≳ 3.3 × 10 6 GeV [120] and the condition M 1 > T dec. is always satisfied for the range ofm 1 we consider (see the discussion above (2.13)). It is important to point out that our present treatment ignores flavour effects [121][122][123][124] such as the charged lepton Yukawa interactions being fast compared to the Hubble scale at different temperatures. These effects can change the asymmetry and consequently the Davidson-Ibarra bound by order one numbers [124] and are expected to be most relevant in the strong washout regimem 1 > 10 −3 eV [122] not applicable here. Now let us take into account the washout of the asymmetry instantaneously produced at T dec. . Because the universe transitions back to a second phase of radiation domination at T dec. , we can reuse the standard estimates for washout. Since the inverse decay requires an on-shell N 1 it gets Boltzmann-suppressed and scales as [71] Γ Consequently for T < T dec. < M 1 we can neglect the washout from inverse decays. That leaves the scattering processes Here one does not include the resonant contribution from on-shell N 1 , as they are already included in the decay term of the Boltzmann equations [59] and the masses of N 2,3 are not kinematically accessible. For T ≪ M 1 the scattering term can be expressed as [71] ∆W In the above we used equations (2.14) and (2.15) for the sum of neutrino masses m ν . This process is negligible, if the absolute value of the exponent is ≲ 0.1 [125], which corresponds to the bound compatible with the findings of [59], indicating that our parameter space (see (2.32)) will be save from any kind of washout: ω ≃ 1.

Dark Matter and Dark Radiation Co-genesis
Dark Matter could be included in Seesaw models via a lightest RHN with keV-scale masses [126,127] produced via either active-to-sterile oscillations [128,129] or gauge interactions [70]. The neutrino mass mediated by a keV-scale N 1 as DM is expected to be smaller than O 10 −5 eV [126]. Since then N 2 would have to play the role of the decaying particle for leptogenesis and we would have to require the associated effective neutrino mass to be below O 10 −7 eV for matter domination (see (2.13)), we would not be able to explain both of the observed neutrino mass splittings in (2.14) and (2.15). Consequently we consider an additional particle as the DM. The out-of-equilibrium decay of a heavy N 1 to this particle might then populate the dark matter abundance. A schematic model for this purpose consists of adding a gauge singlet Majorana fermion ψ and a real singlet scalar σ, either of which (or both) could play the role of dark matter a priori. This approach was first considered in reference [67] for the context of asymmetric dark matter and later in [68] for the case of CP-conserving decays to DM. The relevant couplings are For the sake of minimality we assumed that ψ is a Majorana fermion. It might as well be a Dirac fermion, if we were to introduce a vector-like partner for it. We assume a general renormalizable scalar potential V (H, σ) for the real scalar σ and that M 1 ≫ m ψ + m σ . Additionally all portal couplings are presumed to be small enough to prevent thermal abundances of ψ, σ in the early universe. The decay width of N 1 to ψσ reads where the factor of 1/2 compared to (2.6) arises because this decay has singlets and not doublets in the final state. We define The discussion in section 2.2 assumed that Γ 1 was the leading decay mode of N 1 determining the temperature T dec. at the end of the matter dominated phase in (2.12). Generally speaking this temperature should be calculated from Max [Γ 1 , Γ ψ ] instead. In order to use the parameter region from section 2.2 we will set BR L ≥ BR ψ . In the following we will assume that ψ is the DM, because as long as σ does not receive a vev [67] it has only a suppressed decay mode to ν L σ for m ψ > m σ via ν L − N mixing, that will be discussed in a moment. Its yield is different from the typical Freeze-in approach [64,130] since the decaying RHN is not in thermal equilibrium with the rest of the bath anymore. It also differs from the super-WIMP [63], because the RHN is relativistic at decoupling unlike the non-relativistic WIMP that decays to DM. For our case one finds [66,131] from which we deduce that One can see that the DM abundance only constraints the product m ψ BR ψ and we use it as a free parameter in the upcoming sections about gravitational waves instead of just m ψ . For small branching fractions our scenario leads to light dark matter. Fermionic DM is only gravitationally bound to the DM halo of our galaxy if m ψ ≳ O(100 eV) [132]. In order to comply with bounds from structure formation, that constrain the free-streaming scale of dark matter, we have to demand that [133] m ψ ≳ O(10 keV). (2.44) Both of these constraints illustrate why we need BR ψ ≪ BR L , which translates to y 1 ≪ λ 1i . In the regime m σ < m ψ the following decay from ν L − N 1,2,3 mixing after electroweak symmetry breaking becomes kinematically allowed [134] and we assume that m σ ≪ m ψ : (2.47) We depict the allowed parameter space in figure 1. Once can see that the parameters m 1 = 10 −7 , 10 −10 eV violate the lifetime constraint, because for each constantm 1 the purple relic abundance iso-contour line is above the red line for τ ψ = 250×10 9 yr. The only viable parameter point in this plot hasm 1 = 10 −15 eV in agreement with (2.47), because here the lifetime iso-contour is above the line for the relic density and we find dark matter close to the GeV-scale. The previous limits only apply for m ψ > m σ . In general the scalar σ couples to the SM Higgs through the following terms (2.48) Figure 1. Parameter space for the dark matter mass m ψ versus the branching ratio BR ψ of the RHN decay to dark matter. Contours with (straight, dashed, dotted) lines correspond tõ m 1 = 10 −7 , 10 −10 , 10 −15 eV . The purple contours reproduce the observed dark matter relic abundance and above the contour the abundance would be too large (for fixedm 1 ). The gray regions are excluded because of unsuccessful structure formation (Lyman-α) and dark matter not being gravitationally bound (Tremaine-Gunn). On the red contours for the DM lifetime from ψ → ν L σ is equal to the observational limit and in the colored region above (for fixedm 1 ) the lifetime would be too small. This excludes the lines withm 1 = 10 −7 , 10 −10 eV, meaning that here onlym 1 = 10 −15 eV is viable for DM. Note that lifetime bound disappears for m ψ < m σ , in which case the entire purple region is allowed. The area in light orange is excluded by our assumption BR ψ ≪ BR L ≃ 1 and the orange region would be excluded, if the real scalar also produced in the RHN decay was stable and light enough to be dark radiation (see the discussion below (2.51)).
For m σ > m h it could decay to the SM Higgs. If this is kinematically forbidden, there could be decay modes lighter SM fermions such as e.g. the electron σ → e + e − e + e − via off-shell SM Higgs bosons. In case σ has no vev these decays require the κ coupling. If σ is too light to decay to SM states or the couplings λ Hσ and κ are very small, then the relic abundance of σ survives until today. In this case and assuming that λ Hσ and κ are small enough to avoid thermalization with the SM plasma, the non-thermal σ could still exist in the form of dark radiation. Its energy density is found from n σ = n ψ = BR ψ n N to be [136] ρ dec.

(2.49)
and we compute the abundance of dark radiation, conventionally parameterized as the number of additional neutrinos as [137] assuming again that M 1 ≫ m σ : We see that σ would lead to too much dark radiation compared with the current Planck bound ∆N Planck+BAO eff.
≃ 0.28 [138] unless we make the branching ratio BR ψ , which also controls the DM production, smaller than about 20% (see figure 1). However we saw previously that BR ψ can be far below a percent for heavy enough DM, which is why we do not necessarily expect observable dark radiation. BBN sets a bound of ∆N BBN eff.
≃ 0.4 [139]. The projected sensitivities of upcoming experiments read ∆N proj.  [142,143] and NASA's PICO mission [144] or ∆N eff. ≲ 0.12 for CORE [145], the South Pole Telescope [146] as well as the Simons observatory [147]. Before closing let us emphasize again that σ only counts as dark radiation when it is very light and stable or long-lived.

Distortion of the inflationary tensor mode spectrum
We assume primordial inflation ended in an epoch of reheating, creating a Standard Model plasma of radiation with an initial temperature T RH set by the reheating dynamics. Gravitational waves produced during inflation first leave the horizon and have constant amplitudes while outside the horizon. After they re-enter the horizon the amplitude becomes damped. The power spectrum of gravitational waves (GWs) today can be written as a function of the wave-number k = 2πf with f being the frequency Ω GW (k) = 1 12 where a 0 = 1 and H 0 ≃ 2.2 × 10 −4 Mpc −1 [148] are the scale factor and expansion rate today and P T denotes the spectrum of tensor modes. It is parameterized in terms of the primordial power spectrum from inflation P prim.
as well as a transfer function T 2 T (k). This transfer function describes the propagation of GWs h ij in the Friedmann-Lemaitre-Robertson-Walker background where primes denote derivatives with respect to conformal time, after the horizon re-entry at a temperature of T in that depends on the wave-number via [33] T in = 5.8 × 10 6 GeV · 106.75 g * (T in ) (3.4) The inflationary tensor power spectrum is conventionally parameterized in terms of its amplitude A T and its spectral index n T at the pivot scale k * = 0.05 Mpc −1 [149] P prim.
This amplitude is related to the scalar power spectrum P ξ (k * ) = 2.0989 × 10 −9 [149] via the tensor-to-scalar-ratio r < 0.035 [150] A T (k * ) = r P ξ (k * ). (3.6) Observations of the cosmic microwave background only constrain the scalar spectral index to be n S = 0.9649 ± 0.0042 [149], which is why we take n T as a constant free parameter. The case of n T > 0 (< 0) is known as a blue-tilted (red-tilted) spectrum. Standard single field slow-roll inflation predicts a red-tilted spectrum, as the tensor spectral index n T satisifes the so-called consistency relation n T = −r/8 [151], however this does not rule out the possibilities of a blue-tilted spectrum, which is well motivated in various scenarios including e.g. string gas cosmology [152], super-inflation models [153], G-inflation [154], non-commutative inflation [155,156], particle production during inflation [157,158], and several others [159]. Here we will also seek to investigate such scenarios from the perspective of models of the early universe and leptogenesis. An epoch of early or intermediate matter domination would change the transfer function compared to the standard case of radiation domination, and hence the expansion of the background is imprinted in the damping of the gravitational wave amplitude. References [27,33,34,[160][161][162] computed this transfer function numerically and found a compact analytical expression with a fitting function F (k) in terms of the total matter density Ω m = 0.31, the first spherical Bessel function j 1 (z k ) and z k ≡ k τ 0 with τ 0 = 2/H 0 [148] being the conformal time today. The factors of the relativistic degrees of freedom encode the expansion of the universe and we use the fitting functions of reference [162] for g * (T in ) and g * S (T in ) with the present day values g 0 * = 3.36 and g 0 * S = 3.91, whereas the Bessel function describes the damping of the gravitational wave amplitude after horizon re-entry. In the limit z k ≫ 1, which always holds for the frequencies we are interested in, we can trade the oscillatory j 1 (z k ) for 1/( √ 2z k ). Note that in references [34,162] the correct limiting behavior was mentioned for the wrong limit z k ≪ 1 (for which one would obtain j 1 (z k ) ∼ z k instead) . We employ the most recent results of [162] for the fitting function F (k). Without intermediate matter domination it reads whereas including an epoch of RHN domination leads to Here we introduce where all quantities with a subscript (superscript) "0" are evaluated today and we set h = 0.7. The entropy dilution factor ∆ was defined in (2.25) and the fit functions read Physically T 1 describes the transition from a radiation dominated phase to a matter dominated epoch and T 2 the case of going from matter domination to radiation domination. T 3 has the same physical interpretation as T 1 but allows for a better numerical fit [162]. One deduces from the wave-number k dec. = 2πf sup. at the time of RHN decay in (3.12) that the gravitational wave spectrum gets suppressed by the entropy dilution for frequencies above f sup. ≃ 2.7 × 10 −10 Hz T dec. Here Ω IMD GW was computed from (3.10) and takes the intermediate matter domination (IMD) from the RHN into account, whereas Ω standard GW from (3.9) appears in the absence of RHN domination.

Other GW sources
So far, when it comes to gravitational waves, most studies involving the Seesaw mechanism have focused on the dynamics of e.g. the U(1) B-L breaking, which underlies the RHN Majorana masses in unified gauge theories [108,109]. The dynamics of the scalar responsible for breaking this gauge symmetry can source a separate stochastic gravitational wave background by means of a first order [93,94,96,97] or second order [35,84,95] phase transition as well as via the formation of a network of cosmic strings [84,86,92,98] via the Kibble mechanism [163]. If the phase transition or the formation of topological defects happens before inflation -and the symmetry is never (non-)thermally restored -any trace of the B-L transition will be diluted away due to the exponential expansion of space-time. The symmetry is broken throughout inflation and reheating if [164] v B-L > Max H I 2π , T max. , (3.22) where the first term is the Gibbons-Hawkings temperature [165] in terms of the Hubble rate during inflation H I and the second term the maximum temperature during reheating [110][111][112], which can be drastically larger than the temperature of the radiation bath at the end of reheating T RH . Since T max. depends on the reheating scenario, the best we can do to get an estimate on v B-L is to assume that H I /(2π) > T max. and saturate the current CMB-limit on H I ≲ 2.5 × 10 14 GeV [166] leading to v B-L ≳ 4 × 10 13 GeV. (3.23) This further motivates why we consider high scale leptogenesis. Moreover this bound is compatible with the condition (2.22) for a thermalized population of N 1 from B-L gauge scatterings. Also note that one could even consider a case, where no additional degrees of freedom except the RHN are added to the SM below the Planck scale, so that there would be no source for the stochastic GW background (in this case the initial thermal RHN abundance would have to come from inflaton decays). Consequently our high scale scenario without a stochastic GW background, being essentially independent of the dynamics of the U(1) B-L transition and the associated scalar, can be viewed as complementary to the existing analyses.

Dark radiation bounds from BBN and CMB decoupling
The energy density in gravitational waves should be smaller than the limit on dark radiation encoded in ∆N eff. from Big Bang Nucleosynthesis and CMB observations (see the discussion below (2.51) for bounds and projections on ∆N eff. ) [217] f =∞ The lower limit of the integration is f min ≃ 10 −10 Hz for BBN and f min ≃ 10 −18 Hz for the CMB. In practice, when e.g. plotting many GW spectra simultaneously, and as a first estimate we neglect the frequency dependence to constrain the energy density of the peak for a given spectrum

Impact of free-streaming particles
As shown in the seminal work [218] and expanded upon in e.g. [219][220][221][222], there is a damping effect on the GW amplitude from free-streaming particles whose mean free path is larger than the Hubble scale. Free streaming particles such as the active neutrinos, the RHN, additional sources of dark radiation or gravitational waves themselves contribute to  anisotropic stress-energy tensor and can reduce the primordial GW amplitude by up to 35.6% [218]. In this work we neglect this effect to focus on the damping from the RHN induced matter dominated epoch as a first estimate, since percent level effects will only become relevant once we have actual data.

General results
In the following we fix r = 0.035 [150] and vary the reheating temperature as well as M 1 ,m 1 together with n T ≥ 0. We depict some example spectra in figures 2 and 3, where we reproduced the figures from reference [223]. We depict the constraints from LIGO/VIRGO [167][168][169][170][171][172] and NANOGRAV [204][205][206][207][208] observations, the CMB as well as BBN as shaded regions in our plots 2-6. It is important to note that the depicted projection for the sensitivity of U-DECIGO [27,[190][191][192][193][194] is optimistic, but we do not employ the most optimistic case known as U-DECIGO-corr, which assumes that the noise of the instrument is only given by the irreducible quantum noise [191] and should therefore treated as a hypothetical best case scenario. The proposal for BBO [187][188][189] is also a bit speculative, because it is supposed to eventually succeed the currently planned LISA mission [185,186]. To remind   the reader of these potential caveats we depict the sensitivities for U-DECIGO and BBO with dashed-dotted lines in the figures 2-6. The plots in figures 4 and 5 depict the case where we fix M 1 as function ofm 1 according to (2.32) in order to reproduce the observed baryon asymmetry via leptogenesis. In the aforementioned plot we also depict which values of m ψ BR ψ would be needed according to (2.43) to fit the dark matter relic abundance for a givenm 1 . The labels "no IMD" in 2, 3 and "no intermediate matter dom." in 4-5 refer to the scenario without RHN domination computed from (3.9), where the only dilution arises from inflationary reheating. One can clearly see in 4 and 5 that the primordial tensor modes get diluted by the entropy released in the RHN decay for frequencies above f sup. ≃ 0.1 Hz, see (3.20). Furthermore one can observe in 4-6 that there is second break in the spectra at frequencies larger than f sup. ∼ T dec. . This is due to the inflationary reheating at T RH and since our scenario is defined by the regime T dec. < M 1 < T RH the second break occurs at a larger frequency. The same figures also show a small subleading suppression of frequencies larger than O(10 −9 Hz), which is due to the entropy released in the QCD phase transition [40]. Irrespective of the value of n T , one can deduce from 2-6 that LiteBIRD [211] will already probe the inflationary tensor modes in the 10 −16 − 10 −18 Hz range. For n T = 0 we find that U-DECIGO [27,[190][191][192][193][194] has the best chance to distinguish our entropy suppressed spectra from the standard case without RHN domination depicted by the dashed line in 4. In case neither BBO [187][188][189] nor U-DECIGO [27,[190][191][192][193][194] detect the tensor mode background expected from inflation, this does not have to rule out pri- mordial gravitational waves and could be a tell-tale sign of scenarios with entropy dilution, such as ours. In the next section we will analyze this in terms of the SNR. The case of n T = 0.5 without RHN domination would start to be probed by the dark radiation bounds in (3.27) from BBN [139] and Planck [138] (see the dashed line in 5) and is only borderline compatible with the existing LIGO/VIRGO [167][168][169][170][171][172] observations . An attempt to explain the recent anomaly in the 12.5-year dataset [208] of the NANOGRAV collaboration [204][205][206][207] with primordial tensor modes would require an extremely large n T ≃ 0.85. The challenge is then to have enough entropy dilution to comply with the dark radiation and LIGO/VIRGO bounds. We depict a spectrum for M 1 = 10 7 GeV,m 1 = 10 −17 eV that could be the source of the anomaly in figure 6 for the case without leptogenesis. The reason for abandoning leptogenesis is simply that with such a large n T the peak of the GW energy density at the typical frequency f sup. = 0.1 Hz (before the dilution kicks in) will already be far too large to comply with the dark radiation bounds. Therefore one needs a spectrum where the damping (which is only proportional tom 1 see (3.21)) occurs at lower decay temperatures and hence lower frequencies (set by both M 1 andm 1 see (2.12)). This is why we chose a value of M 1 = 10 7 GeV below the leptogenesis bound in (2.32). On top of that we set T RH = 5 × 10 12 GeV, so that the GWs at large frequencies beyond LIGO/VIRGO do not come into tension with the dark radiation bound due to the damping from inflationary reheating. These estimates illustrate, why we would need a rather contrived scenario and we do not pursue the aforementioned anomaly further in this work.

Signal-to-noise ratio
We use the SNR defined in (3.25) to determine the region in the M 1 versusm 1 parameter space, where a detection of primordial gravitational waves can be claimed for a SNR threshold of ten over four years of observation time. For n T = 0 we find that BBO [187][188][189], µ-ARES [196] and U-DECIGO [27,[190][191][192][193][194] are the most relevant experiments that have a chance of probing the primordial GW background, as can be deduced from figure 4. For n T > 0 there are a lot more experiments that can probe our GW spectra, which is why we focus on AEDGE [177,195], BBO [187][188][189], the Einstein Telescope (ET) [181,182] and LISA [185,186]. Of course there are also other currently developed experiments, such as the radio telescope SKA [198][199][200], that become relevant for n T > 0. The parameter space for n T = 0 was displayed in 7, whereas figure 8 showcases n T = 0.1, 0.2 and 9 the cases of n T = 0.3, 0.5. The region in (2.32) that leads to the observed baryon asymmetry via leptogenesis was shaded in gray. For n T = 0.1 one can conclude from the left plot in figure 8 that the SNR threshold for ET [181,182] will start to probe the edge of the parameter space for leptogenesis in the regimem 1 ≲ 10 −11 eV. For n T > 0 we see in 8-9 that AEDGE [177,195], BBO [187][188][189] and LISA [185,186] probe the entire parameter space for leptogenesis. We impose the following constraints in figures 7-9: Successful BBN requires that the RHN decay temperature in (2.12) is at least 10 MeV [224,225], which was depicted as a brown region. RHN with masses above 10 14 GeV could destabilize the electroweak vacuum [226,227]. We do not show the bound M 1 ≲ 10 7 GeV [228][229][230][231] from the naturalness of the Higgs mass under corrections from its couplings to the RHN, as it would basically exclude our entire parameter space in (2.32). The last bound comes from the observed neutrino masses: Due to the perturbativity of the RHN Yukawa coupling λ ij < √ 4π and the need to reproduce at least one mass eigenstate with m ν = 0.05 eV we find that M 2 ≲ 3.8 × 10 15 GeV. This together with our assumption that M 2 > 3M 1 means that we have to require at least M 1 ≲ 10 15 GeV. In all plots we fixed T RH = 10 16 GeV so that even the heaviest N 1 allowed by the previous considerations would be present in the plasma. As mentioned in the previous section we find that U-DECIGO [27,[190][191][192][193][194] is the best candidate to test our setup compared to the case with no decaying RHN for n T = 0. A future non-observation of the inflationary tensor mode spectrum could be explained by a decaying N 1 withm 1 < 10 −14 eV and a mass of M 1 ≳ 10 4 GeV (the precise number depends on the BBN bound on the RHN decay temperature of at least 10 MeV). By fixing m 1 as a function of M 1 for leptogenesis via (2.32) we plot the SNR as a function of M 1 on the right side of 7. Here the SNR for µ-ARES [196] is constant because the peak of its sensitivity is situated at a frequency below f sup. ≃ 0.1 Hz and it is therefore blind to the entropy damping. For cosmologies with n T > 0 we find that the SNR for U-DECIGO [27,[190][191][192][193][194] is always larger than 10 in the depicted parameter space, which is why we focus on different detectors. BBO [187][188][189] is a promising candidate for a detection of primordial GWs with both n T = 0 and n T > 0 (compare the plots in 7 and 8, 9). For n T ≳ 0.5 the dark radiation bound becomes important again and we show the contour ∆N proj. eff = 0.06 for CMB Stage IV [142,143] computed via (3.26) on the right side of figure  9. For completeness we display the SNR as a function of M 1 (withm 1 fixed by leptogenesis (2.32)) for n T = 0.1, 0.5 in figure 10.

Conclusions and Discussions
We focused on the minimal Seesaw model, which adds only three right handed neutrinos (RHN) to the SM, and demonstrated that an epoch of right handed neutrino domination, with a Yukawa coupling corresponding tom 1 < 2.9 × 10 −7 eV, can realize baryogenesis Figure 10. SNR as a function of M 1 , wherem 1 was fixed for leptogenesis via (2.32) with n T = 0.1 (left) and n T = 0.5 (right). In both plots we fixed T RH = 10 16 GeV via leptogenesis for a mass of M 1 ≳ 2.4 × 10 8 GeV · 2 × 10 −7 eV/m 1 (see (2.32)). Since the effective mass ism 1 is too small for a thermal RHN population, we had to assume a different production channel via either inflaton decays or B-L gauge scatterings for the initial RHN abundance. Furthermore such a smallm 1 requires that one of the SM neutrinos is approximately massless compared to the other two. The amplitude of gravitational waves that re-enter the horizon before the end of the RHN matter dominated epoch is damped by a factor proportional to the entropy released in the RHN decay. We discussed the detection possibilities of primordial GWs and computed the signal-to-noise ratio for various detectors such as AEDGE [177,195], BBO [187][188][189], DECIGO [27,[190][191][192][193][194], Einstein Telescope [181,182], LISA [185,186] or µ-ARES [196] as well as for several spectral tilts n T ≥ 0 of the tensor mode spectrum. Additionally we determined the regions in the M 1 versusm 1 parameter space in which the signal-to-noise ratio (SNR) is larger than ten over a four year observational period in the figures 7-9. Our main finding is that high scale leptogenesis can have an observable imprint on the gravitational waves from inflation. Further we discussed under which conditions our scenario leads to the dominant GW signal. Since fixing M 1 as a function ofm 1 for successful leptogenesis by saturating the maximum of the CP-violating decay parameter ε 1 for a hierarchical spectrum (see (2.28)) completely determines the RHN decay temperature to be T dec. ≃ 3.3 × 10 6 GeV, we find a constant characteristic frequency of f sup. ≃ 0.1 Hz (see (3.20) and figures 4-5), above which the suppression of the GW amplitude manifests itself. The same RHN can also have a second potentially suppressed decay mode to a stable fermion ψ, that is responsible for the dark matter abundance, if the product of the DM mass and the branching fraction of the RHN decay to DM satisfies m ψ BR ψ ≃ 85 eV · 2 × 10 −7 eV/m 1 . In order for the dark matter do be heavy enough for successful structure formation (m ψ > O(10 keV)) for fixed m 1 we typically need a small branching ratio BR ψ ≪ 1. Such a small branching fraction can also suppress the amount of BSM dark radiation ∆N eff. ≃ 0.06 · (BR ψ /4%), that could potentially be generated, if the scalar produced together with ψ is very light and survives until today. This particular scenario leads to GeV-scale DM decaying to dark radiation and SM neutrinos, which necessitatesm 1 < 9.7 × 10 −15 eV in order to have DM with a large enough lifetime on cosmological scales and the right relic abundance.