Fingerprint of GeV scale right-handed neutrinos on inflationary gravitational waves and PTA data

We show that the seesaw mechanisms that exhibit right-handed neutrino mass-dependent non-standard post-inflationary cosmology make blue-tilted inflationary gravitational waves (GW) compatible with the recent findings of nHz stochastic GW background by the pulsar-timing arrays (PTAs) for high reheating temperatures. The right-handed neutrino (RHN) mass scale has to be $\mathcal{O}(\rm GeV)$. Remarkably, such a scenario produces a correlated signature testable by the future LIGO run. In addition to contributing to the active neutrino masses, $\mathcal{O}(\rm GeV)$ RHNs generate baryon asymmetry of the universe via low-scale-leptogenesis. They can be searched for in collider experiments. Therefore, the recent detection by PTAs is not only exciting for GWs in the nHz range; it paves the way to test and constrain well-studied mechanisms, such as seesaws, with a low-frequency and a correlated measurement of high-frequency GW spectral features, complementary to particle physics searches.

fit to the old as well as the new data [1,6,[35][36][37][38][39]; though, it should be noted that such BGWs can be produced in models which, in general, do not correspond to standard slow-roll inflation, see, e.g., [40][41][42][43][44][45][46][47].The parameter space of such a fit is, however, restrictive.This is because GWs with large blue-tilt, considering the spectrum is still a power-law at higher frequencies, saturate Big Bang nucleosynthesis (BBN) bound on the effective number of neutrino species, disfavoring any post-inflationary cosmology founded on high reheating temperature (T RH ≳ 10 GeV) after inflation [9,39].Nonetheless, if a non-standard matter epoch leads to entropy production between the reheating after inflation and the most recent radiation domination before the BBN [48], BGWs get suppressed and provide a good fit to PTA data for high reheating temperatures.Now, contrary to the standard case, such a scenario allows the overall GW spectrum to span decades of frequencies with characteristic spectral features testable, e.g., by the LIGO [49][50][51].
In this letter, we show the seesaw mechanisms with GeV scale right-handed neutrinos (RHN), which are now being extensively discussed within the context of lowmass sterile neutrino searches (see a review: [52]), naturally provide an RHN mass-dependent matter epoch to fit the PTA data with BGWs (we shall see later that a high T RH is also a requirement of the model).The scenario correlates RHN masses with the amplitude and spectral features in the BGWs verifiable at high-frequency GW detectors, providing a novel synergic search for GeV scale RH neutrinos, which, besides generating active neutrino masses via seesaw, also offer successful baryogenesis via low-scale leptogenesis [53][54][55][56].
The theoretical framework is founded primarily on this question: What is the origin of small (here GeV) RHN masses?Although the electroweak naturalness condition puts an upper bound on the RHN masses: M i < 10 7 GeV [57,58], generally, in GeV scale seesaw scenarios, the origin of such small RHN masses are not addressed; they are considered to be the bare masses in the theory.Nonetheless, the seesaw Lagrangian provides all the degrees of freedom if we suppose RHN masses originate from a phase transition driven by a scalar field [59,60].In which case, the RHN field N couples to a scalar field Φ as L ∼ f N N N Φ (omitting family indices), with f N being a Yukawa coupling.After the phase transition, Φ obtains its vacuum expectation value v Φ and generates RHN mass as M = f N v Φ .On the other hand, if kinematically allowed, Φ can decay to a pair of RHNs (Φ → N N ), with the decay rate Γ ∝ f 2 N .RHN masses and the decay width (or lifetime) of Φ are now connected via f N .For a fixed v Φ (which may be large), one obtains a required small value of M for a small f N , making Φ long-lived.We shall see that the long-lived Φ dominates the universe's energy budget as a matter component before decaying.The smaller the M , the longer the lifetime of Φ.This results in a longer duration of matter domination-hence larger entropy production and more suppressed BGWs.The time and amount of entropy production correlate the amplitude and spectral features of BGWs with RHN mass scale M .For more technical details, please see Ref. [60] where the idea explained above was introduced.

II. RHN MASS-DEPENDENT MATTER DOMINATION
Obtaining RHN masses via a phase transition is not ad-hoc; the coupling f N N N Φ can appear as U (1) B−L symmetric coupling following the breaking of a grand unification group [61][62][63][64].In that case, N and Φ have B − L charge 1 and −2, respectively.Therefore, for concreteness, we shall consider a B − L phase transition, i.e., as the temperature drops, the scalar field rolls from Φ = 0 towards Φ = v Φ , breaking the B − L symmetry.The finite temperature potential that restores the symmetry at higher temperatures is given by [65,66] where D, E and T 0 are functions of gauge coupling g ′ , the self-interaction coupling λ, and . The last term in Eq.( 1) generates a potential barrier causing a secondary minimum at Φ ̸ = 0, which at T = T c degenerates with the Φ = 0 one.The potential barrier vanishes at T 0 (≲ T c ), making the minimum at Φ = 0 a maximum.The critical temperature T c and the field value Φ c ≡ Φ (T c ) are given by [60] T The transition dynamics can be described as a rolling of Φ if it smoothly transits to Φ = v Φ , which can be quantified roughly with the order parameter Φ c /T c ≪ 1.In this case, the field rolls because the potential barrier disappears very quickly [59,60].The condition Φ c /T c ≪ 1 can be fulfilled for λ ≃ g ′ 3 and g ′ ≲ 10 −2 , which correspond to Φ c /T c ≲ 0.08 [60].Once it rolls down, the field oscillates around v Φ .For V (Φ) = αΦ β , the equation of state of such an oscillation phase is computed as ω = (β − 2)(β + 2) −1 [60].Assuming the oscillation of the scalar field is driven by the dominant quadratic term in the potential and expanding the zero temperature potential around v Φ , we obtain α = λv 2 Φ and β = 2. Therefore, the scalar field behaves like matter (w = 0).One can also compute the angular frequency of oscillation as m Φ = √ 2λv Φ .The decay channels determine the lifetime of Φ.For λ ≃ g ′ 3 and g ′ ≪ 1, Φ → z ′ z ′ is not allowed from kinematic consideration.Here z ′ is the U (1) B−L gauge boson.Another decay mode Φ → hh, where h is the SM model Higgs, does not contribute if Φ is sequestered from SM Higgs at the tree level (even at the radiative level, it is not efficient due to the discussed small values of f N ).Then, the competing decay channels are Φ → N N and Φ → f f V ( The corresponding two-body decay; Φ → f f is suppressed due to chirality flip, e.g., [67]), where f and V are SM fermions and vector bosons.The former corresponds to a tree-level process, whereas the latter is a one-loop (z ′ z ′ f ) triangle process.The strengths of these two processes are determined by the couplings f N and g ′ .To control the duration of matter domination with small RHN mass via f N , the process Φ In which case, the entropy produced (κ) by the decay of Φ is given by (which amounts typically equivalent to the ratio of the temperatures corresponding to the time when Φ dominates and when it decays, see, e.g., the supplementary material) [60] where MP l = 2.4 × 10 18 GeV is the reduced Planck constant, The analytical formula for κ is very precise, and we shall use it to compute the GW spectrum.Let's mention the following conditions that the model complies with.I) As mentioned, Γ Φ N ≳ Γ f f V , II) Φ decays before BBN (T ∼ 5 MeV), III) the transition happens following reheating after inflation, i.e., T c ≲ T RH , IV) the vacuum energy does not dominate the radiation at T c (ρ Φ (T c ) < ρ R (T c )); violation of which leads to a second period of inflation.We shall see shortly that there are three additional constraints, excluding the PTA data, and the model has six free parameters.Therefore, the model withstands with 6:8 capacity against the constraints (including PTA data), leading to extremely robust predictions.Specifically, we shall see that the recent PTA data can be explained only for a constrained range of RHN masses, which can be probed by the fu-ture LIGO run and is also interesting for future collider FCC-ee.

III. FIT TO THE NANOGRAV DATA AND PREDICTIONS
GWs are described with the perturbed FLRW line element: where τ is the conformal time, a(τ ) is the scale factor.The transverse and traceless part of the 3 × 3 symmetric matrix h ij ; ∂ i h ij = 0 and δ ij h ij = 0, characterizes the GWs.Following the linearized evolution equation considering subdominant contribution from anisotropy stress tensor π ij [68], and solving the Fourier space propagation equation for h ij , one obtains the gravitational wave energy density as [69] Ω where H 0 ≃ 2.2 × 10 −4 Mpc −1 and τ 0 = 1.4 × 10 4 Mpc.The quantity P T (k) represents the primordial power spectrum connecting to the inflation models: where r ≲ 0.06 [70] (constraint V) is the tensor-to-scalarratio, k = | ⃗ k| = 2πf with f being the frequency of the GWs at the present time at a 0 = 1, A s ≃ 2 × 10 −9 is the scalar perturbation amplitude at the pivot scale k * = 0.01Mpc −1 and n T is the tensor spectral index.The simplest single-field slow-roll inflation models satisfy a consistency relation: n T = −r/8 [71].We shall treat n T > 0 as constant, ignoring scale dependence due to higher-order corrections [72].The most important quantity in the discussion is the transfer function T T (τ 0 , k) given by [73-78] where F (k) reads with j 1 (kτ 0 ) being the spherical Bessel function, Ω m = 0.31, g * 0 = 3.36, g * 0s = 3.91.We use the scale-dependent g * 0(s) (T k,in ) in Eq.( 9) from [78][79][80], where T k,in is the temperature corresponding to the horizon entry of kth mode.The T i (ζ)s are given by where ζ i ≡ k/k i , with k i s being the modes: k Φ = 1.7 × 10 14 g * s (T Φ ) 106.75 k ΦR = 1.7 × 10 14 κ 2/3 g * s (T Φ ) 106.75 k R = 1.7 × 10 14 κ −1/3 g * s (T RH ) 106.75 crossing the horizon at standard matter-radiation equality temperature T eq , at T Φ ≃ 90 Γ Φ N MP l when Φ decays, at T ΦR when Φ dominates the energy density, and at T RH , respectively.Given the above set of equations and using κ from Eq.( 3), we evaluate Eq.( 6) to obtain the GW spectrum and to fit the NANOGrav data, considering two more constraints.VI) LIGO bound on SGWB, which roughly reads Ω GW (35 Hz)h 2 ≤ 6.8×10 −9 [49], and VII) BBN bound on the effective number of neutrino species: where ∆N eff ≲ 0.2 [81].f low corresponds to the frequency entering the horizon during BBN epoch.On the other hand, the Hubble rate at the end of inflation determines the upper limit: f high = a end H end /2π.For numerical computations, f high ≃ 10 5 Hz would suffice because the spectrum falls and the integration saturates.We follow the NANOGrav parametrization for the GW energy density to perform a power-law fit to the new data within the frequency range f ∈ 10 −9 Hz, f yr .The parametrization reads with Ω yr = 2π 2 3H 2 0 A 2 f 2 yr and f yr = 1yr −1 ≃ 32 nHz.Fitting the data requires comparing Eq.( 6) and Eq.( 17), then extracting the values of the amplitude A and the spectral index γ that lie within the 1, 2, 3σ contours (cf.middle panel of Fig. 1) reported by the NANOGrav [1].In the left panel of Fig. 1, we show a GW spectrum consistent with all the constraints.To produce the figure, the following benchmark values for the model parameters have been used: T RH = 10 13 GeV, n T = 0.9, r = 0.06, v Φ = 10 14 GeV, g ′ = 10 −2.9 , and M = 16 GeV.The Left: The spectrum has been generated for M = 16 GeV.For other benchmarks, please see the text.Besides NANOGrav and LIGO, sensitivities of SKA [82], LISA [83], DECIGO [84], and ET [85] are shown.The vertical dashed blue line represents the frequency corresponding to Eq.( 18).Middle: The red ⋆ represents the fit-point corresponding to the GW spectrum on the left.Right: An allowed parameter space on the g ′ −M plane (white).All the colored regions are excluded.The NANOGrav signal cannot be reproduced in the brown region at 2σ. BBN constraint on ∆N eff disfavors the grey region.The LIGO-O3 bound on SGWB excludes the purple region.The region right to the sky-blue line (indicated by arrows) represents the future ALIGO sensitivity.In the pink region (top left corner), the three-body decay of Φ is more dominant than the RH neutrino pair production.In the red region on the top, one has Tc > TRH.
corresponding values of A and γ are shown in the middle panel with the red '⋆'.The fit is not very different from the standard BGW-fit without intermediate matter domination [39].This is because, within the NANOGrav frequency range, the transfer function T 1 (ζ) determines the spectral shape, which results in γ ∼ 5 − n T ≃ 4. Note also that the first peak of the spectrum occurs at a frequency f Φ > f yr so that the NANOGrav band can be fitted with a pure power-law Ω GW (f ≲ f yr ) ∼ f n T .An analytical expression for f Φ can be derived from Eq.( 14), which is given by One can do a more exhaustive fit by varying the parameters.For instance, in the right panel, we varied g ′ and M , keeping the rest fixed to their benchmark values.Nonetheless, there is no qualitative difference − allowed values of the other parameters lie near the benchmarks.This does not change significantly even though one varies all the parameters.As stated earlier, this is because of a bunch of constraints that the model complies with.Among all the constraints, the first, i.e., Γ Φ N ≳ Γ f f V , the third, i.e., T RH ≳ T c , and the sixth, i.e., LIGO bound on SGWB, are much stronger, making the allowed parameter space consistent with the NANOGrav data, stringent.This model fits the NANOGrav data at 2σ with M ∈ [1,47] GeV for random values of other parameters around the benchmark and without violating any constraints.Remarkably, consistency with the NANOGrav data makes the model extremely predictive, not only in terms of RHN masses; the infrared tail of the second peak in the GW spectrum passes through the sensitivity region of advanced LIGO (ALIGO), for a significant portion of  13 TeV run [52,86,87].Future sensitivities of SHiP [88], MATHUSLA [89], NA62 [90], and FCC-ee [91]  allowed parameter space (see the right panel in Fig. 1).Therefore, any non-observation of SGWB by ALIGO would potentially rule out a large parameter space of the model, provided that the model fits NANOGrav data.Not only that, the latter also motivates us to combine, for the first time, the particle physics sensitivity curves for GeV scale RHN searches with the ALIGO projection, shown in Fig. 2 with the vertical sky-blue band representing the RHN mass range M ∈ [2.5,47] GeV.Particle physics experiments are sensitive to RHN masses and their mixing (|U α | 2 ∼ m ν /M ) to active neutrinos, where m ν is the active neutrino mass scale (≃ 0.05 eV).On the contrary, predictions of this model depend on the former, allowing us to identify the region independent on |U µ | in Fig. 2. We conclude with the following remarks covering some additional aspects of the work.1) The framework can be extended to other variants, such as inverse-seesaw, which offers additional phenomenology.2) One may straightforwardly obtain analytical expressions of the peak and dip frequencies in terms of RHN masses using Eq.( 14)-( 16) (via T Φ ). 3) We do not present explicit computation of leptogenesis.It would be interesting to reproduce calculations, e.g., of [56] (including entropy production), leading to a parameter space sensitive to |U α | and M .4) This one is perhaps the most interesting: the possibility of obtaining GWs from cosmic strings.The occurrence of B − L phase transition would naturally produce cosmic gauge strings.We, however, work with unconventional values of λ and g ′ , which in general are taken O(1) in the Nambu-Goto simulations [92][93][94].Additionally, for the parameter space consistent with the NANOGrav data, the cosmic string width δ w ∼ 1/ √ λv Φ constitutes a considerable fraction of the horizon H(T c ) −1 (relatively thick strings).Claiming GWs from cosmic strings in this model thus requires a straightforward assumption (which we are less confident about): results of the numerical simulations also hold for our preferred parameter range.An existence of GWs from cosmic strings, nonetheless, would produce further spectral distortion to the BGWs shown in Fig. 1, making a combined peak-plateau-peak spectrum instead of a peak-dip-peak one (supplementary material can be seen).This distinguishes the model from any other matter domination+BGW scenario, even at the level of the GW spectrum.

IV. SUMMARY
We discuss a novel framework to probe seesaw models with GeV scale right-handed neutrinos (RHN) with the recent pulsar timing (PTA) data interpreted as stochastic gravitational waves background (SGWB) from inflation.A fit to the PTA data with inflationary GWs predicts the mass scale of RHN to be O(GeV) and a PTA-LIGO correlation on SGWB.While any non-observation of SGWB by advanced LIGO (ALIGO) would rule out a large parameter space of the model, the recent PTA data motivates us to combine the particle physics sensitivity curves for low mass RHN searches with the future LIGO projection for the mass range M ∼ [2.5,47] GeV.We performed the fit with the NANOGrav 15 yrs data.We do not expect a significant qualitative change in our results if the fit is performed by combining the data from all the PTAs, because the A − γ global contours reported by the IPTA collaboration are similar [95] to the NANOGrav ones.where z = T c /T , and from the third of Eq.( 19), the time-temperature relation has been derived as The amount of entropy production is computed numerically by solving and then computing the ratio of S ∼ a 3 /z 3 after and before the scalar field decays.The evolution of the relevant quantities with z = T c /T is presented in Fig. 4.

Combined gravitational waves spectrum including cosmic strings
In this scenario, cosmic strings appear after the spontaneous breaking of U (1) B−L symmetry.Following the formation, the strings get randomly distributed in space, forming close loops plus a network of horizon-size long strings.When two segments of long strings cross, they inter-commute and form loops.Long strings are characterised by a correlation length L = µ/ρ ∞ , with ρ ∞ being the long string energy density and µ is the string tension defined as µ = πv 2 Φ h λ 2g ′2 .The quantity h is a slowly varying function with . Therefore, in our analysis, the string tension is given by µ Generally, a string network oscillates to enter a scaling evolution phase, characterized by stretching the correlation length due to the Hubble expansion and the fragmentation of the long strings into loops.Numerical simulations also support this phase.These loops oscillate independently and produce GWs.An attractor solution exists between these two competing dynamics: the scaling regime.In this regime, L scales as cosmic time t, which corresponds to ρ ∞ ∝ t −2 .Therefore, the network tracks cosmological background energy density ρ bg ∝ t −2 with a small constant proportional to Gµ, and does not dominate the energy density of the Universe.The evolution of a radiating loop of initial size l i = αt i is described as l(t) = αt i − ΓGµ(t − t i ), where Γ ≃ 50 and α ≃ 0.1.The total energy loss from a loop can be decomposed into a set of normal-mode oscillations with frequencies f j = 2j/l j = a(t 0 )/a(t)f , where j = 1, 2, 3...j max (j max → ∞).The jth mode GW density parameter is given by Ω where n (t, l j ) is the scaling loop number density, which can be computed analytically as well as with numerical simulations.Without any intermediate matter-dominated epoch, the GWs arising from Eq.( 23) can be described with a peak at a low frequency owing to the GW radiation from the loops formed in the radiation era and decay in the standard matter era, plus a plateau at high frequency: that arises due to the loop formation and decay only in the radiation era.In Eq.( 24), ϵ r = α/ΓGµ ≫ 1, A r ≃ 5.4 and Ω r ∼ 9 × 10 −5 .
In the presence of a new matter epoch before the most recent radiation era, the above description remains the same; barring, the spectrum becomes red at a high frequency f CS Φ : Ω 1 GW (f > f CS Φ ) ∝ f −1 .The frequency f CS Φ can be calculated analytically and is given by where z eq ≃ 3387 is the red-shift at the standard matter-radiation equality taking place at time t eq , and t Φ < t BBN is the time when the matter domination ends.In Fig. 5, we show the resulting GW spectrum produced by cosmic strings ( for j = 1 ) for the same benchmark used to produce the BGW spectrum.Notice that when combined with the BGWs, it creates a plateau in the LISA sensitivity region.This signature is exclusive to the model under consideration.The falling GW spectrum from cosmic strings can not affect the overall spectrum because, in that frequency range, the BGWs dominate.
FIG. 1.Left: The spectrum has been generated for M = 16 GeV.For other benchmarks, please see the text.Besides NANOGrav and LIGO, sensitivities of SKA[82], LISA[83], DECIGO[84], and ET[85] are shown.The vertical dashed blue line represents the frequency corresponding to Eq.(18).Middle: The red ⋆ represents the fit-point corresponding to the GW spectrum on the left.Right: An allowed parameter space on the g ′ −M plane (white).All the colored regions are excluded.The NANOGrav signal cannot be reproduced in the brown region at 2σ. BBN constraint on ∆N eff disfavors the grey region.The LIGO-O3 bound on SGWB excludes the purple region.The region right to the sky-blue line (indicated by arrows) represents the future ALIGO sensitivity.In the pink region (top left corner), the three-body decay of Φ is more dominant than the RH neutrino pair production.In the red region on the top, one has Tc > TRH.

FIG. 2 .
FIG. 2. Particle physics exclusion and projections for RHNmixing with muon flavor.The grey region and region above the brown curve are excluded from the previous experiments and CMS13 TeV run[52,86,87].Future sensitivities of SHiP[88], MATHUSLA[89], NA62[90], and FCC-ee[91] are shown with green, orange, blue, and pink curves.The region named BBN is excluded; otherwise, the decay product of RHNs would contradict BBN prediction.The dashed line represents the mixing in the canonical seesaw |Uµ| 2 ∼ mν /M .The vertical red band represents the RHN mass range M ∼[1,47] GeV consistent with NANOGrav 2σ data.The vertical skyblue band represents the RHN mass range M ∼ [2.5, 47] GeV consistent with NANOGrav 2σ data and testable with SGWB searches by advanced LIGO.
FIG. 2. Particle physics exclusion and projections for RHNmixing with muon flavor.The grey region and region above the brown curve are excluded from the previous experiments and CMS13 TeV run[52,86,87].Future sensitivities of SHiP[88], MATHUSLA[89], NA62[90], and FCC-ee[91] are shown with green, orange, blue, and pink curves.The region named BBN is excluded; otherwise, the decay product of RHNs would contradict BBN prediction.The dashed line represents the mixing in the canonical seesaw |Uµ| 2 ∼ mν /M .The vertical red band represents the RHN mass range M ∼[1,47] GeV consistent with NANOGrav 2σ data.The vertical skyblue band represents the RHN mass range M ∼ [2.5, 47] GeV consistent with NANOGrav 2σ data and testable with SGWB searches by advanced LIGO.

ACKNOWLEDGMENTR.
Samanta is supported by the project International Mobility MSCA-IF IV FZU -CZ.02.2.69/0.0/0.0/20079/0017754 and acknowledges European Structural and Investment Fund and the Czech Ministry of Education, Youth and Sports.

FIG. 4 .
FIG.4.Evolution of the normalized energy densities of the radiation, the scalar field, and total entropy (normalized to its initial value).The horizontal dashed line represents an excellent match of the analytical approximation for the entropy production given in Eq.(3) (in the letter) with the numerical one.
FIG. 5. Blue-tilted GW spectrum (blue) combined with GWs from gauge cosmic strings (red) originating due to the (B − L) symmetry breaking in this model. t