Probing Leptogenesis and Pre-BBN Universe with Gravitational Waves Spectral Shapes

On the frequency-amplitude plane, Gravitational Waves (GWs) from cosmic strings show a flat plateau at higher frequencies due to the string loop dynamics in standard radiation dominated post-inflationary epoch. The spectrum may show an abrupt upward or a downward trend beyond a turning point frequency $f_*$, if the primordial dark age prior to the Big Bang Nucleosynthesis (BBN), exhibits non-standard cosmic histories. We argue that such a spectral break followed by a rising GW amplitude which is a consequence of a post-inflationary equation of state ($\omega>1/3$) stiffer than the radiation ($\omega=1/3$), could also be a strong hint of a leptogenesis in the seesaw model of neutrino masses. Dynamical generation of the right handed (RH) neutrino masses by a gauged $U(1)$ symmetry breaking leads to the formation of a network of cosmic strings which emits stochastic GWs. A gravitational interaction of the lepton current by an operator of the form $\partial_\mu R j^\mu$--which can be generated in the seesaw model at the two-loop level through RH neutrino mediation, naturally seeks a stiffer equation of state to efficiently produce baryon asymmetry proportional to $1-3\omega$. We discuss how GWs with reasonably strong amplitudes complemented by a neutrino-less double beta decay signal could probe the onset of the most recent radiation domination and lightest RH neutrino mass at the intermediate scales.

On the frequency-amplitude plane, Gravitational Waves (GWs) from cosmic strings show a flat plateau at higher frequencies due to the string loop dynamics in standard radiation dominated post-inflationary epoch. The spectrum may show an abrupt upward or a downward trend beyond a turning point frequency f * , if the primordial dark age prior to the Big Bang Nucleosynthesis (BBN), exhibits non-standard cosmic histories. We argue that such a spectral break followed by a rising GW amplitude which is a consequence of a postinflationary equation of state (ω > 1/3) stiffer than the radiation (ω = 1/3), could also be a strong hint of a leptogenesis in the seesaw model of neutrino masses. Dynamical generation of the right handed (RH) neutrino masses by a gauged U (1) symmetry breaking leads to the formation of a network of cosmic strings which emits stochastic GWs. A gravitational interaction of the lepton current by an operator of the form ∂ µ Rj µ -which can be generated in the seesaw model at the two-loop level through RH neutrino mediation, naturally seeks a stiffer equation of state to efficiently produce baryon asymmetry proportional to 1 − 3ω. We discuss how GWs with reasonably strong amplitudes complemented by a neutrino-less double beta decay signal could probe the onset of the most recent radiation domination and lightest RH neutrino mass at the intermediate scales. Leptogenesis [1][2][3][4][5][6][7] is a simple mechanism to explain the observed baryon asymmetry of the universe [8]. The right handed (RH) heavy neutrinos which are introduced in the Standard Model (SM) to generate light neutrino masses (Type-I seesaw), decay CP asymmetrically to create lepton asymmetry which is then converted to baryon asymmetry via Sphaleron transition [9]. When it comes to the testability of leptogenesis, there are subtleties. If the heavy neutrino masses are not protected by any symmetry [10], it is quite natural to assume that they are hierarchical in nature like any other family of SM fermions. In that case, the lightest RH mass scale is bounded from below M ≳ 10 9 GeV [11]which is beyond the reach of the present collider experiments. Nonetheless, still the colliders and other low energy neutrino experiments can probe leptogenesis mechanisms that do not constitute hierarchical RH neutrinos-starting from O(TeV) to O(MeV) scale heavy neutrinos [12][13][14][15]. A shift of attention from the collider experiments to the Gravitational Waves (GWs) physics is not less interesting in terms of testing leptogenesis. Particularly, this new cosmic frontier, in which after the discovery of GWs from black hole mergers by LIGO and Virgo collaboration [16,17], plenty of efforts are being made to detect primordial GWs from the Early Universe (EU) within a wide range of frequencies-starting from the Pulsar Timing Arrays (PTAs, ∼ nHz) to the LIGO(∼ 25Hz). A network of cosmic strings [18][19][20] which is a generic consequence of breaking symmetries such as U (1), is one of the prominent sources of strong stochastic primordial gravitational waves which can be tested in a complementary way in most of the planned GW detectors. Numerical simulations based on the Nambu-Goto action [21,22] indicate that cosmic string loops loose energy dominantly via GW radiation, if the underlying broken symmetry corresponds to a local gauge symmetry. In the context of seesaw, this sounds music to the ears, since such a gauge symmetry is U (1) B−L [23][24][25], breaking of which could be responsible for the dynamical generation of the heavy RH masses and hence the lepton number violation as well as creation of a network of cosmic strings. Having this set-up, there could be two categories to look for the GWs as a probe of leptogenesis. Category A: A scale separation between the RH masses and the typical Grand Unified Theory (GUT) scale (∼ 10 16 GeV), imposed by seesaw perturbativity condition and the neutrino oscillation data [26] implies that residual symmetries like U (1) B−L protects the RH neutrinos to get mass at the GUT scale. Therefore, breaking of that symmetry at a later stage and consequent emission of GWs from cosmic strings are natural probes of the scale of leptogenesis. In this case, it is the amplitude (GW energy density normalised by the critical energy density) of the GWs that matters as a probe and this approach has been taken in Refs. [27,28]. Category B: To make the testability more robust, along with the amplitudes, one can associate leptogenesis also to the spectral shapes of the GWs [29,30]. Cosmic string loops that originate and decay in the radiation domination, exhibit a flat plateau on the amplitude-frequency plane at the higher frequencies. This spectral shape may show an upward or a downward trend if something other than radiation dominates the energy density of the EU before the onset (T * ) of most recent radiation domination prior to the BBN (T ∼ 5 MeV) [31][32][33]. Such a non-standard cosmic history that is responsible for this spectral break which along with the GW amplitude, one aims to claim also as a probe, should therefore be a natural/well-motivated call from the perspective of leptogenesis. Two well-known scenarios in this context can be opted for. Category B1: A matter domination (ω = 0 < 1/3) [34,35]. Category B2: Scenarios such as kination (ω = 1 > 1/3) [36,37]. For the former (latter), one finds a spectral break followed by a downward (upward) going GW amplitude [38][39][40][41]. Two leptogenesis mechanisms in the Category B1-a low-scale leptogenesis and a leptogenesis from ultralight primordial black holes (M P BH ≲ 13g) have been studied in Ref. [29] and Ref. [30] respectively. In this article, we discuss a scenario that falls in the Category B2, i.e., interpreting a flat then a spectral break followed by a rising GW amplitude as a signature of leptogenesis.
Note that, two crucial ingredients for this typical signal are of course cosmic string network itself and then a non-standard equation of state (ω = 1 in our discussion). In the context of leptogenesis from decays [1], though the former is a natural consequence in the sense of Category A [27], a stiffer equation of state is not an indispensable criterion. However, in seesaw models, even when the Lagrangian is minimally coupled to gravity, through massive RH neutrino mediation one can generate an operator of the form ∂ µ Rj µ /M 2 at two-loop level [42][43][44] (see also Ref. [28,45] for a flavour generalisation and Ref. [46] for a recent review), where R is the Ricci scalar and j µ is the lepton current. This operator is a well-studied operator [47][48][49][50] with the corresponding mechanism dubbed as "gravitational lepto/baryogenesis" and produces final baryon asymmetry proportional toṘ ∝ (1 − 3ω). Interestingly, note now that two primary ingredients of the GW signal are also natural requirements to obtain non-zero lepton asymmetry, i.e., the symmetry breaking which gives rise to massive RH neutrinos (mediate in the loops [44]) as well as cosmic strings and then an equation of state ω ≠ 1/3 [51]. We shall discuss later on, that indeed a stiffer equation of state is needed to efficiently produce lepton asymmetry. Plateau amplitudes corresponding to Gµ ≲ 10 −12 with G being the Newton constant and µ being the string tension, with a post LISA spectral break supplemented by a potential test in neutrino-less double beta decay experiments, make the scenario generally robust. The above introduction summarises the basic idea and the main results of this paper. The next sections are dedicated to a more detailed description and technicalities.

II. GRAVITATIONAL WAVES FROM COSMIC STRINGS
Cosmic strings may originate as the fundamental or composite objects in string theory [52,53] as well as topological defects from spontaneous symmetry breaking (SSB) when the vacuum manifold M has a non-trivial first homotopy group π 1 (M). A theory with spontaneous breaking of a U (1) symmetry exhibits string solution [19,20], since π 1 (M) = Z. An example of a field theory containing string like solution is a theory of U (1)-charged complex scalar field φ that in the context of seesaw could be a SM scalar singlet φ B−L which is responsible for the dynamical generation of RH neutrino masses. After the formation, strings get randomly distributed in space and form a network of horizon-size long strings [54,55] characterised by a correlation length L = µ/ρ ∞ , where µ-the string tension or energy per unit length is in general constant (however, e.g., in case of global strings [56] and recently introduced melting strings[57] µ ∼ f (T )) and typically taken to be the square of the symmetry breaking scale Λ CS and ρ ∞ is the long string energy density. When two segments of long strings cross each other they inter-commute and form loops with a probability P = 1 [58] (exceptions [59]). A string network may interact strongly with thermal plasma and thereby its motion gets damped [60]. After the damping stops, the strings oscillate and enter a phase of scaling evolution that constitute two competing dynamics namely the stretching of the correlation length due to the cosmic expansion and fragmentation of the long strings into loops which oscillate independently and produce particle radiation or gravitational waves [61][62][63]. Out of these two competing dynamics, there is an attractor solution called the scaling regime [64][65][66] in which the characteristic length scales as L ∼ t. This implies, for constant string tension, ρ ∞ ∝ t −2 .
Therefore, the network tracks any cosmological background energy density with the same equation of state and hence cosmic strings do not dominate the energy density of the universe like any other defects. The loops radiate GWs at a constant rate which sets up the time evolution of a loop of initial size [61,63] and the initial loops size parameter α ≃ 0.1-a value preferred by numerical simulations [67,68]. The total energy loss from a loop is decomposed into a set of normal-mode oscillations with frequencies f k = 2k/l = a(t 0 )/a(t)f , where k = 1, 2, 3...k max (k max is for numerical purpose, otherwise ∞) and f is the frequency observed today. Given the loop number density n t , l k , the present time gravitational wave density parameter is given by The quantity Γ k depends on the small scale structures of the loop and is given by Γ (k) = Γk −δ ζ(δ) , e.g., δ = 4/3 and 5/3 for cusps and kinks [69]. The integration in Eq.II.1 is subjected to a Heaviside func- , end of damping(t fric )] and l cric is the critical length below which massive particle radiation dominates over GWs [70,71]. Both these Θ functions set a high-frequency cut-off in the spectrum (a systematic analysis can be found in Ref. [40]).
The most important aspect to obtain the GW spectrum is the computation of the loop number density n t , l k which we calculate from the Velocity-dependent-One-Scale (VOS) model [72][73][74] which assumes the loop production function to be a delta function, i.e. all the loops are created with the same fraction of the horizon size with a fixed value of α. Given a general equation of state parameter ω, the number density n ω t , l k is computed as where β = 2/3(1 + ω) and we assume A β = 29.6 (ω = 1), 5.4 (w = 1/3) and 0.39 (ω = 0) [74] is a step-function while changing the cosmological epochs. The most interesting feature of GWs from cosmic string is that the amplitude increases with the symmetry breaking scale Λ CS . This can be seen by computing the Ω GW , considering loop production as well as decay in the radiation domination which gives an expression for a flat plateau at higher frequencies (see AUX A for an exact formula) where r = α/ΓGµ and Ω r ≃ 9 × 10 −5 . Such strong GWs as a consequence of a very high scale symmetry breaking thus serves as an outstanding probe of particle physics models [75][76][77][78][79][80][81]. Possibly the most important recent development is the finding of a stochastic common spectrum process across 45 pulsars by NANOGrav PTA [82], which if interpreted as GWs, corresponds to a strong amplitude and is better fitted with cosmic strings [28,83,84] than the single value power spectral density as predicted by supermassive black hole models. Let's also mention that a very recent analysis by PPTA [85] is in agreement with the NANOGrav result. In presence of an additional early kination era, the entire GW spectrum is determined by four dynamics. I) A peak at a lower frequency-caused by the loops which are produced in the radiation era and decay in the standard matter era. II) The flat plateau, Ω plt GW , as mention while describing Eq.II.3. III) A spectral break at f * = eq -so called the turning point frequency [30,39,40], followed by a rising GW amplitude Ω , caused by modified redshifting of the GWs during kination era V) a second turning point frequency f ∆ after which the GWs amplitude falls, e.g, due to particle productions below l < l cric = β m µ −1/2 (ΓGµ) m , with β m ∼ O(1) and m = 1, 2 for loops with kinks or cusps [70,71]. If the falling is caused due to thermal friction, then one needs to consider the damping of the smaller loops along with the long-string network for t < t f ric , discarding any GWs production by the smaller loops, i.e, the entire dynamics is completely frozen until t f ric [60]. In fact, in our computation we do not take into account any GWs produced from smaller loops prior to t f ric and consider that the falling is due to particle production which sets the high-frequency cut-off that is much more stronger (appears at lower frequencies) than the friction cut-off [40]. Note also that if the two turning-point frequencies are close to each other, potentially the GW detectors could see a small bump after the flat plateau with a peak amplitude ≃ Ω plt GW (f ∆ /f * ). Nevertheless, as we show in the next section that given a successful leptogenesis, the second turning point frequency as well as small bumps are most likely to be outside the frequency range of the GW detectors.
Before concluding the section, we note two important points. Firstly, the VOS model overestimates the number density of the loops by an order of magnitude compared to the numerical simulations [67]. This is due to the fact that VOS model considers all the loops are of same size at the production. However, there could be a distribution of α. Numerical simulation finds that only 10% of the energy of the long-string network goes to the large loops (α ≃ 0.1) while the rest 90% goes to the highly boosted smaller loops that do not contribute to the GWs. This fact is taken into account by including a normalisation factor F α ∼ 0.1 in Eq.II.2 [74]. Secondly, the amplitude beyond f * goes as f 1 even after taking into account high-k modes (see AUX A) unlike the case of an early matter domination where the same changes from f −1 → f −1/3 for cusps like structures [29,30,40].

III. GRAVITATIONAL LEPTOGENESIS, RESULTS AND DISCUSSION
The idea behind gravitational leptogenesis [48] is, a C and CP-violating operator L CP V ∼ b∂ µ Rj µ ∼ b∂ µ R¯ γ µ with b as a real effective coupling, corresponds to a chemical potential µ = bṘ for the lepton number in the theory. Therefore, the normalised (by photon density n γ ∼ T 3 ) equilibrium lepton asymmetry (using standard Fermi-Dirac statistics with energies E ± = E ±µ) is given by N eq B−L ∼ bṘ T . Interestingly, L CP V can be generated in a UV framework using the seesaw Lagrangian even when it is minimally coupled to gravity (see e.g., Ref. [43] for an in-depth discussion, sec.II of Ref. [28] for a brief summary). As a computational insight, one calculates an effective h vertex corresponding to the operator L CP V using a conformally flat metric g µν = (1 + h)η µν with R = −3∂ 2 h, capitalising the fact that the coupling 'b' is independent of the choice of background.
In seesaw model, a similar h vertex that manifests the L CP V operator, can be constructed at two-loop level, where the Higgs and the RH masses mediate the loops. Then simply comparing the coefficients of both the vertices up to linear order in h, the coupling b can be calculated in terms of the Yukawa coupling f (where, f αi¯ LαH N Ri is the Yukawa interaction in seesaw, with Lα , H and N R being the lepton doublet, Higgs and RH fields respectively) and RH neutrino masses M i . The expression for the equilibrium asymmetry then reads The above expression could be modulated by a factor (M 2 j /M 2 i ) γ , where γ = 0, 1. However, γ = 0 appears to be the most natural solution which can be calculated exactly [43,44]. In any case, even if one considers γ = 1 or the 'hierarchical enhancement', tuning the complex part in k 2 ij , correct baryon asymmetry can always be reproduced. The most important part is, N B−L ∝Ṙ ∝ 1 − 3ω which is still vanishing in radiation domination at high temperatures with SM-QCD thermodynamic potential [51]. Therefore, a general cosmological background other than radiation that is quite a natural call now, always stems a non-vanishing equilibrium asymmetry unless the Yukawa couplings are real or purely imaginary. In the EU, any dynamically produced lepton asymmetry tracks the N eq B−L if the interaction that causes the asymmetry production is strong enough. When the interaction rate becomes weaker (compared to the Hubble expansion), the asymmetry freezes out with the potential to reproduce correct baryon asymmetry N B−L ∼ 6 × 10 −8 [8]. In seesaw model, ∆L = 2 interactions [86] play this role. The general evolution equation that governs the entire dynamics is given by , where z * = M 1 /T * and m i is the i-th light neutrino mass with i = 1, 2, 3. All the exact expressions can be found in AUX B. Before proceeding further, let us mention that we do not include the charged lepton flavour effects in this analysis for simplicity. Nonetheless, a systematic description with flavour issues can be found in Ref. [45] along with a more finer description in Ref. [28]. To proceed further, the process consists of two distinct temperature regimes. At a higher temperature T in ∼ Λ CS , as soon as the symmetry breaks, the RH neutrinos become massive and Eq.III.2 starts acting without W ID which is negligible at this regime. In this gravitational leptogenesis scenario, typically, z in (= M 1 /T in ) can be constrained with so called weak field condition as z in ≥ M 1 /M P l , whereM P l is the reduced Planck constant. Once the asymmetry freezes out, at the lower temperatures, it faces a washout by the inverse decays which are strongly active at T ∼ M 1 . The final asymmetry is therefore of the form N where the dimensionless washout/decay parameter K 1 is a function of Yukawa couplings. Eq.III.3 that matches with the numerical solutions of the Eq.III.2 with quite a high accuracy, is the master equation which we use to present all the results. Prior to the explanation of Fig.1, let's introduce a parametrisation of the Yukawa matrix as m D = U mΩ M , where m D = f v with v = 174 GeV, U is the leptonic mixing matrix and Ω is a 3 × 3 complex orthogonal matrix with a standard parametrisation in terms of three complex rotation matrices [28] with complex angles θ ij = x ij +iy ij . In general, Ω is a completely 'free' matrix unless one invokes additional symmetries to fix the flavour structure of the theory. A plethora of works is dedicated in this direction [10]. With this orthogonal parametrisation it is easy to show that the equilibrium asymmetry is independent of U . Therefore, as far as the seesaw parameters are concerned, the light, heavy neutrino masses and the orthogonal matrix take part in the process. The decay parameter can also be expressed in terms of these parameters as with m * ≃ 10 −3 being the equilibrium neutrino mass [4]. In Fig.1, we show the variation of the produced asymmetry with the lightest neutrino mass for three benchmark values; M 1 = 10 6,7,8 GeV with a fixed orthogonal matrix and different values of z * . The basic nature of the curves is quite interesting. Let's focus on the z * = 10 3 curve (yellow) for M 1 = 10 8 GeV. It shows a plateau until m 1 ≃ 10 −2 eV, then an increase followed by a downfall at large m 1 values. First of all, for w = 1, the parameter κ ∼ z −1 * and therefore for large values of z * , the strength of the ∆L = 2 process becomes so weak that the asymmetry instantly freezes out without tracking the equilibrium number density. The coefficient f κ does not change much until m 1 ∼ 10 −2 eV and then increases for m 1 ≳ 10 −2 eV [28]. This increase in f κ pushes the asymmetry more towards the equilibrium and hence the overall magnitude of N B−L increases for m 1 ≳ 10 −2 eV. A downfall at large m 1 is caused by the exponential term in Eq.III.3. The washout is in fact modulated by two parameters, K 1 and z * . However, for large values of m 1 , the parameter K 1 becomes huge and therefore, even if one has a large z * , the frozen out asymmetry is completely washed out. On the other hand, when z * is small, e.g., z * = 10 2 , the overall magnitude of N B−L decreases since β ∼ z 3 * . In this case however, z −1 * suppression in κ is not that significant compared to the previous one. Until m 1 ∼ 10 −2 eV, it shows the constant behaviour due to the mentioned nature of f κ , however, at large m 1 values, it becomes strong enough to maintain the asymmetry in equilibrium for a period of time. The downfall is mostly dominated due to this equilibrium asymmetry tracking and not due to the late time washout. Note that for ω < 1/3, for a fixed value of z * , κ increases (causes delayed freeze out and hence dilution of the asymmetry N G0 B−L ) and β decreases (causes a decrease in N eq B−L ). A concrete example is a matter domination, i.e., ω = 0, where κ ∼ z * and β ∼ z −3/2 * . Moreover, these kind of scenarios are inclusive of a late time entropy production which dilutes the produced asymmetry significantly [30,35]. Therefore, ω < 1/3 scenarios are utterly inefficient. This possibly strengthens the claim that in the future, should the GW detectors find a flat and then a rising signal, RH neutrino induced gravitational leptogenesis with a stiffer equation of state is a natural mechanism to associate with, since both of them, successful leptogenesis and the GW signal, are triggered by common theoretical ingredients.
In Fig.2, we show the future sensitivities of the GW detectors such as LISA [87], BBO [88], CE [89], ET [90] on the Gµ − T * plane. In the case of strong GW amplitudes, the most stringent constraint comes from the effective number of neutrino species which reads ∫ df f −1 Ω GW (f )h 2 < 5.6 × 10 −6 ∆N ef f . Considering ∆N ef f ≤ 0.2, the peak of the spectrum at f ∆ , and taking into account contributions from the infinite number of modes that give a factor of ζ(7/3) amplification , where we consider particle production from cusps [71]. The corresponding region has been shaded in red. We have ignored the variation of the effective relativistic degrees of freedom even when T * is below the QCD phase transition. Proper temperature dependence of the same, would include a factor of 1.5-3 modification. Since we are entirely onto the gravitational leptogenesis (to motivate ω ≠ 1/3), we take M max 1 ∼ 10 8 GeV so that the contribution from the decays are negligible. This gives an upper bound on the T in (Λ CS ) that corresponds to Gµ ≲ 10 −12 . Therefore, the mechanism can be tested with reasonably strong GW amplitudes even for the flat part (Eq.II.3). For strong amplitudes, the spectral breaks are likely to happen at high-frequency GW detectors like CE and ET plus the bump like signals (f * and f ∆ are close to each other) in general lie outside those detectors. In Fig.2, the black point represented by ♠ (♣), should (not) be a signal (see a supplementary Fig.3).
We shall end the discussion with a 'Neutrino-Gravitational Waves Complementarity (NGWC)' or more generally, how this type of GW signal could be supplemented by low energy neutrino experiments. NGWC depends on the z * and flavour structure of the theory or more precisely, on the orthogonal matrix. From Fig.1, it can be seen that, depending on the RH neutrino mass (hence Gµ), various z * values are sensitive to the neutrinoless double beta decay experiments (the N B−L curves intersect with the N Obs B−L at the same time falls within the vertical green region). For the parameter set in Fig.1, the NGWC points fall unfortunately in the red region as well as they are well outside the GW detectors (showed by the ♡ and ♢ points in Fig.2). However, if one decreases the y ij , to produce correct N B−L , for a fixed value of Gµ one needs larger values of z * , meaning the NGWC points would move towards the left side, i.e., towards the smaller values of T * . The entire picture can be encapsulated within the triangle drawn on the test parameter space shaded in green in Fig.2. The red horizontal arm represents the constant Gµ line along which the entries of Ω decrease as one goes from larger to smaller T * . The yellow arm represents the constant z * line, as one goes along the line towards smaller Gµ values, entries of Ω increase and the . A fall at a high frequency is due to the particle production from cusps for l < l cric = µ −1/2 (ΓGµ) 2 [40,71]. We have shown the spectrum only for the fundamental mode. blue arm represents the constant (already predicted) orthogonal matrix and as one goes towards the higher values of Gµ, z * decreases or in other words, T * increases. The blue arm is of great interest. If one has a completely determined orthogonal matrix, from Fig.1 the NGWC points can be determined with the sets of M 1 and T * . This means the blue arm is a line of predictions from the GW experiments, i.e., we can predict at which amplitude and at which frequency the spectral break would occur. The triangle as a whole can be pushed towards the larger T * values increasing y ij . This implies, seesaw models which exhibit an orthogonal matrix with large imaginary part entries, would likely to show the spectral break at higher frequencies and therefore may not be tested with the planned detectors. These models are dubbed as 'boosted' seesaw models where the light neutrino basis vectors and heavy neutrino basis vectors are strongly misaligned [91]. On the other hand, models with flavour structures close to 'form dominance' [92] that typically predicts a real orthogonal matrix (Ω = P , where P is a permutation matrix), would show a spectral break within the frequency range of the current or planned GW detectors.
The frequency derivative of ρ GW is given by the quantity P GW (t, f k ) represents the power emitted by the loops and is given by (see e.g., Ref. [68]) Integrating Eq.III.6 over the loop lengths gives From Eq.III.7 and Eq.III.5 one gets and therefore the energy density corresponding to the mode 'k' is given by Using the VOS equations and considering the loop production function as a delta function (see, e.g., Ref. [74]), it is easy to obtain the most general formula for the number density in an expanding background that scales as a ∼ t β . The expression is given by . (III.10) The Eq.III.9 can be expressed in the conventional form that are used in many papers (e.g., Ref. [39,40]) using the time dependence of the loop length which gives initial time t Now using Eq.III.11, the number density in Eq.III.10 can be re-expressed as (III. 12) time t * at which the most recent radiation domination begins, an approximate frequency up to which the spectrum shows a flat plateau is given by Similarly, using the critical length l cric = µ −1/2 (ΓGµ) 2 for cusp like structures, the second turning point frequency can be computed as where we have assumed that f ∆ /f * ≃ t * /t ∆ (cf. Eq.III.21), and for simplicity we consider g * (T ) = g * ≃ 106 throughout. Therefore, for post QCD phase transition T ≲ 200 MeV, the formula is bit errorful.
To observe both the frequencies distinctively, one should have f ∆ > f * . This gives the following restriction on the parameter space which is shown by the red region in Fig.2.
A4. The BBN limit: To be consistent with the the number of effective neutrino species, the GW energy density has to comply with with ∆N ef f < 0.2. Considering the dominant contribution from the non-flat part after the first turning point frequency f * , the following constraint on the parameter space can be obtained On the other hand considering the analytic approach, i.e., using Velocity dependent One Scale (VOS) the same is obtained as where A r = 5.4. As mentioned before, the VOS model assumes all the loops are of same length at creation. However, at the moment of creation, the loops may follow a distribution depending on α. If so, the above formula should be modified as n(t, l k (t)) = A r ∫ w(α)N α dα l k (t) + ΓGµt 5/2t 3/2 .

(III.27)
Therefore, to the make VOS formula in Eq.III.26 consistent with the numerical result, one has to normalise Eq.III.26, i.e., n(t, l k (t)) = F α A r N α l k (t) + ΓGµt . (III.38) The most general Boltzmann equations (BEs) for leptogenesis with seesaw Lagrangian minimally coupled to gravity are where the first equation governs the production of RH neutrinos and the first term in the second equation represents the contribution to the lepton asymmetry from RH neutrino decays. Since we are neglecting the contribution from decays, only the second equation with the terms in 'bold' is relevant. Note that recently in Ref. [44], another curvature-induced evolution term that modulates of the asymmetry production dynamics at ultra-high temperatures has been introduced. We neglect that term in our computation. However, that will not change the qualitative features of our final results. To obtain a simpler form of the Boltzmann equation it is convenient to simplify the expression of the equilibrium asymmetry and the lepton number violating processes. Using the orthogonal parametrisation of the Dirac neutrino mass matrix m D = U mΩ M , the equilibrium asymmetry can be expressed as a power law in z as