Towards distinguishing Dirac from Majorana neutrino mass with gravitational waves

We propose a new method towards distinguishing the Dirac versus Majorana nature of neutrino masses from the spectrum of gravitational waves (GWs) associated with neutrino mass genesis. Motivated by the principle of generating small neutrino masses without tiny Yukawa couplings, we assume generic seesaw mechanisms for both Majorana and Dirac neutrino masses. For Majorana neutrinos, we further assume a spontaneously broken gauged $U(1)_{B-L}$ symmetry, independently of the type of Majorana seesaw mechanism, which gives a cosmic string induced GW signal flat over a wide range of frequencies. For Dirac neutrinos, we assume the spontaneous breaking of a $Z_2$ symmetry, the minimal symmetry choice associated with all Dirac seesaw mechanisms, which is softly broken, generating a peaked GW spectrum from the annihilation of the resulting domain walls. In fact, the GW spectra for all types of Dirac seesaws with such a broken $Z_2$ symmetry are identical, subject to a mild caveat. As an illustrative example, we study the simplest respective type-I seesaw mechanisms, and show that the striking difference in the shapes of the GW spectra can help differentiate between these Dirac and Majorana seesaws, complementing results of neutrinoless double beta decay experiments. We also discuss detailed implications of the recent NANOGrav data for Majorana and Dirac seesaw models.


Introduction
Charged fermion masses in the standard model (SM) are necessarily of the Dirac type because of electric charge conservation.Neutrino mass, on the other hand, may be of two types: Dirac or Majorana, where the latter possibility arises due to the fact that neutrinos are electrically neutral.If neutrinos were massless particles, as originally envisioned in the SM, their nature, i.e.Dirac vs Majorana, would not be distinguishable in weak interactions.However, as oscillation experiments have confirmed, they possess nonzero, albeit tiny, mass [1][2][3].
Dirac and Majorana neutrino mass are traditionally distinguished experimentally by neutrinoless double beta decay [4,5].This process is allowed only in the former case, where the first entry of the Majorana neutrino mass matrix m ββ is model dependent.Extensive experimental efforts are currently underway for detecting neutrinoless double beta decay, achieving upper bounds |m ββ | ≲ O(10−100) meV [6], and are expected to gain another order of sensitivity in the next decade [7].However, the non-observation of neutrinoless double beta decay will not be decisive about the Majorana or Dirac nature of neutrinos.At this point, it remains interesting to seek astrophysical or cosmological probes of distinguishing the nature of the neutrino mass (see [8][9][10][11] for recent studies), and this work is motivated by such considerations.
We begin by recalling that neutrino mass is associated with the breaking of separate lepton numbers L e , L µ , L τ .Dirac neutrino mass preserves total lepton number L = L e + L µ + L τ , while Majorana mass breaks it.In the latter case, the small mass of the neutrinos may originate from the dimension-five Weinberg operator lℓHH, where ℓ represents the lepton doublets and H is a Higgs doublet, breaking the lepton number by two units.This could be associated with the spontaneous breaking of an Abelian U (1) L symmetry which may be global, or, when combined with baryon number, gauged U (1) B−L .The occurrence of cosmic strings in the early universe is a consequence in both scenarios, and their subsequent decay can produce detectable gravitational wave (GW) signatures [12][13][14].This offers a well known possible observable indication of Majorana neutrino mass.
A convincing ultraviolet completion of the Weinberg operator is achieved by introducing right-handed neutrinos that get large Majorana masses after spontaneous breaking of a U (1) L or U (1) B−L symmetry, the latter case opening up the possibility of a gauged Abelian symmetry which may be anomaly free if there are precisely three right-handed neutrinos.The type-I seesaw mechanism [15][16][17][18][19] then provides an elegant explanation for the generation of light neutrino masses, avoiding the need for extremely small Yukawa couplings.Other realizations of the Weinberg operator, such as the type-II [20][21][22][23][24] and type-III [25] seesaw mechanisms that involve different intermediate particles, can also generate small neutrino masses without tiny Yukawa couplings.All these mechanisms may yield cosmic string induced GW signals if the U (1) L or U (1) B−L symmetry is spontaneously broken.This provides a generic signature for all seesaw mechanisms responsible for Majorana masses.
If neutrinos are Dirac particles, the most minimal extension of the SM would be to add two or three right-handed SM singlet neutrinos ν R with tiny tree-level Yukawa couplings y D , defined by y D lν R H, together with a conserved U (1) L or U (1) B−L symmetry.However such an approach involves tiny Yukawa couplings y D a million times smaller that that of the electron.There have been many attempts which yield Dirac neutrinos without relying on such tiny Yukawa couplings .Each of these mechanisms has its own experimental implications, but most studies have not considered GW signatures.One particularly attractive idea is to invoke a Dirac seesaw mechanism as an ultraviolet completion of the dimension-five operator 1  Λ lν R Hσ, with both tree-level and one-loop realizations, where σ is a scalar and Λ denotes the scale of heavy intermediate particles.1 Although there are three types of Dirac seesaw mechanisms, they all have one thing in common in their minimal formulations, namely they rely on a spontaneously broken Z 2 symmetry under which σ and ν R are odd [48].This leads to domain wall formation, and, including a soft Z 2 breaking term, domain wall annihilation and GW production, providing a distinctive signature generic to all Dirac seesaw mechanisms.Intriguingly, the GW spectrum is determined only by the 1 The Dirac seesaw mechanism, whose minimal version requires a Z2 symmetry, has the phenomenological advantage over just the inclusion of Dirac mass terms in that it can facilitate leptogenesis in its minimal realization [47].Furthermore, the particles mediating the one-loop realizations of the Dirac seesaw operator can be dark matter candidates [48], analogous to the scotogenic model [49] in the case of a Majorana seesaw.These are in addition to the theoretical motivation of explaining the smallness of the neutrino Dirac Yukawa couplings, which is not accounted for in the usual Dirac mass model (without a Z2 symmetry).potential of σ, thus yielding identical signals for all Dirac seesaw mechanisms provided σ does not couple to a scalar mediator.Motivated by the above considerations, and generic assumptions, it seems possible that the GW spectrum can distinguish the Dirac from Majorana seesaw mechanisms, the former yielding a peaked spectrum (from the Z 2 domain walls) and the latter a flat spectrum (from the U (1) L or U (1) B−L cosmic strings).For definiteness, we focus on the type-I seesaw mechanism for the Majorana mass generation, whereas for Dirac mass generation, without loss of generality we focus on the type-I Dirac seesaw mechanism, since the GW signatures are identical for all Dirac seesaw types.2An illustrative example of this idea is to analyze the spectral shape of the GW signal at nHz frequencies in the two types of neutrino mass models, where several pulsar timing arrays (PTAs) including NANOGrav [50][51][52][53][54][55][56][57], EPTA [58][59][60][61][62][63], PPTA [64][65][66], and CPTA [67] have reported convincing evidence of a stochastic gravitational wave background which cannot arise from a population of inspiraling supermassive black hole binaries.We show that the Majorana mass model which yields a cosmic string induced GW signal is unlikely to produce such a signal, whereas the domain wall induced GW signal from the Dirac mass model remarkably matches the PTA result.
The organization of the paper is as follows.In section 2 we discuss models of neutrino mass generation in the Majorana and Dirac cases.Sections 3 and 4 focus on the production of gravitational waves, specifically from cosmic strings in the context of Majorana mass generation, and from domain walls in relation to Dirac mass generation.The resulting signals are examined in section 5, followed by concluding remarks in section 6.

Neutrino masses via the seesaw mechanism
In this section we explore the generation of Majorana and Dirac neutrino masses via dimension-five operators that realize the respective seesaw mechanism.Both operators produce small neutrino masses without assuming tiny Yukawa couplings.In both cases a family of tree-level models emerges which share a common theme.For the Majorana seeesaw, the B − L symmetry must be broken, which we assume to be spontaneously broken from an ultraviolet U (1) B−L symmetry.For the Dirac seesaw, U (1) B−L must remain unbroken, but naturally preventing tree-level Dirac mass necessitates a Z 2 symmetry, which must be spontaneously broken in the process of generating the Dirac mass of the SM neutrinos.A key observation is that the respective breaking of symmetries in the Majorana and Dirac seesaws, regardless of the specific tree-level model, yields strikingly different GW signatures.

Majorana seesaw
As an example of a seesaw model with a spontaneously broken U (1) B−L symmetry, we consider a type-I seesaw in which the SM is extended with three right-handed neutrinos Ni and a scalar ϕ, both singlet under the SM gauge groups.However, the model has an anomaly-free gauged B − L symmetry, under which ϕ has two units of charge and Ni have a single unit of charge.The Lagrangian of the model is given by which yields the diagram, The right-handed neutrinos acquire heavy Majorana masses after the B − L symmetry is spontaneously broken when ϕ gets a nonzero vacuum expectation value (VEV).Light neutrino masses are generated by integrating out the heavy right-handed neutrinos, and their mass matrix is given by where M N is the mass matrix of the right-handed neutrinos.
The most elegant type-I seesaw models invoke a gauged U (1) B−L instead of a bare mass term for the right-handed neutrinos explicitly breaking the B − L symmetry.This is motivated by anomaly cancellation, separating the scale of the right-handed neutrino masses from the grand unification scale, where the breaking chains may contain an intermediate U (1) B−L in grand unified theories such as SO (10), Pati-Salam and left-right symmetric models.
The main characteristic of this scenario is the breaking of the U (1) B−L symmetry, which creates a cosmic string network that eventually decays and produces a stochastic gravitational wave background.As we will discuss in section 3, such GW signals are nearly flat for a vast range of observable frequencies in gravitational wave interferometers.We note that there could be a secondary contribution to the GW spectrum if the scalar ϕ undergoes a first order phase transition (FOPT) when it spontaneously breaks the U (1) B−L symmetry.However, it is well known that the signal from FOPT of a single scalar is typically suppressed compared to the cosmic string signal, particularly when the U (1) symmetry is broken at a sufficiently high scale.Hence, we will not consider the FOPT signal for this model.
. Three tree-level realizations of the dimension-five Dirac seesaw operator lν R Hσ.A Z 2 symmetry is required in all three cases to forbid a tree-level Dirac mass term lν R H.

Dirac seesaw
Analogous to the Majorana seesaw, there are three tree-level realizations of the Dirac seesaw operator lHν R σ depending on the mediator field.In Fig. 1 we show the Feynman diagrams corresponding to these cases.
The intermediate Dirac fermion ∆ is a singlet under SM SU (2) L , and the intermediate scalar ρ and fermion Σ transform as doublets.In all three cases the scalar σ and the right-handed neutrino ν R are gauge singlets.
The Dirac-ness of the neutrino mass can be ensured by imposing a U (1) L or U (1) B−L symmetry which remains unbroken.It was noted in Ref. [48] that the Dirac mass generated via any of the three diagrams is guaranteed to be the leading contribution by imposing a Z 2 symmetry under which ν R and σ are odd.This statement holds even for one-loop realizations of the Dirac seesaw operator [48].The essential point is that all minimal ultraviolet completions of the operator lν R Hσ at the tree and one-loop level requires a Z 2 symmetry, which is spontaneously broken when σ gets a nonzero VEV.
As a concrete example, we will discuss the type-I scenario for the remainder of the paper, but we emphasize that our main argument, i.e. spontaneous breaking of the Z 2 symmetry, would be a common feature of all relevant scenarios.The Lagrangian of the type-I model is given by We assume that the Dirac fermion ∆ is heavy.After integrating it out, when the SM Higgs and the new scalar get VEVs v and u, respectively, an effective Dirac mass term M D Lν R for the light neutrinos is generated, where is the Dirac mass matrix suppressed by the large eigenvalues of the mass matrix M ∆ of the heavy fermions ∆.
The scalar σ spontaneously breaks the Z 2 symmetry when it acquires a nonzero VEV, necessary for Dirac mass generation.This leads to the creation of domain walls.Long-lived domain walls are dangerous for cosmology if they dominate the energy density of the Universe.However, they can be made to annihilate into gravitational waves by softly breaking the Z 2 symmetry, which lifts the degeneracy between the two Z 2 vacua.This leads to characteristic GW signals peaked at a single frequency.Since the global lepton number symmetry remains unbroken due to the Dirac nature of the neutrinos, this setup does not lead to the generation of cosmic strings and an associated flat GW spectrum as in the Majorana case.

GWs from cosmic strings in Majorana seesaw model
To make the discussion self-contained, we briefly recount how GWs are produced from the spontaneous breaking of a gauged U (1) symmetry.The breaking of the U (1) B−L symmetry leads to the creation of a horizon-length string network [68].We specifically focus on Nambu-Goto cosmic strings that lose energy primarily through loop formation and emission of gravitational radiation.The energy density in the string network is diluted by producing closed string loops [69], about 10% of which are large loops and the remaining are highly boosted smaller loops [70][71][72][73].The formation of the loops from long string networks can be described using the velocity-dependent one-scale model [74,75].The loop formation rate is assumed to be equal to the rate of energy loss of the evolving long string network in a cosmological background, and is given by with the parameter values α ≃ F α ≃ 0.1, C eff ≃ 0.5 and 5.7 during matter and radiation domination, respectively, are found from lattice simulations [76].
While the kinetic energy of the smaller loops are diluted by simple redshifting, the larger loops oscillate and emit energy in the form of gravitational waves at a constant rate, where Γ ≃ 50 is a dimensionless constant [77], G is the Newton's constant and µ is the tension in the strings.Typically µ ∼ O(Λ), where Λ is the scale of the U (1) symmetry breaking.As a consequence of emitting gravitational radiation, the initial length of a large loop created by the network at time t i , given by l i = αt i , decreases as until the loop completely disappears.The total energy loss from a loop can be decomposed into normal modes with frequency fk = 2k/ℓ at a time t, where k = 1, 2, 3, . . . is the mode number.Accounting for redshift evolution, the frequency today becomes f = [a( t)/a(t 0 )] fk , where t 0 is the current time.The relative emission rate per mode is given by Combining Eqs.(3.1), (3.2) and (3.3), and integrating over the emission time, the gravitational wave amplitude of the k-th mode is given by where ρ c = 3H 2 0 /(8πG) is the critical energy density, t is the formation time of loops contributing to the k-th mode and is given by Summing over all modes, we get the total amplitude of the gravitational waves where the sum can be easily evaluated using 4 GWs from domain walls in Dirac seesaw model We now focus on the GW spectrum generated in the type-I realization of the Dirac seesaw operator.We will argue that the GW spectrum is identical for any tree-level or one-loop realization of the Dirac seesaw operator in which the scalar σ does not mix with the SM Higgs or other scalars.We assume a simple potential for σ: This potential has two degenerate minima at σ = ±u and is symmetric under a Z 2 transformation σ → −σ.Domain walls are formed around the boundaries of these two minima.The symmetry is spontaneously broken when the scalar chooses one of the two vacua.This choice depends on random fluctuations of the field and is made independently at spatially distant regions in space, creating the so-called 'domains'.Domain walls are formed around the boundaries of these domains.We assume that the domain walls have a static planar configuration perpendicular to the z direction.Introducing a kinetic term 1 2 (∂ µ σ) 2 , the field equation for σ(z) is given by which yields the solution, for the boundary condition σ(z → ±∞) → ±u.The surface energy density (also called tension) of the wall can be derived from integrating the 00 component of the stress-energy tensor T µν = (dσ/dz) 2 diag(+1, −1, −1, 0), and is given by Domain walls can be very long-lived and may dominate the energy density of the Universe, alter its equation of state and lead to rapid expansion inconsistent with standard cosmology.Even if their energy density is subdominant today, domain walls may produce excessive density perturbations observable in the CMB epoch if their surface energy density is above O(MeV 3 ) [78].
An interesting solution to the domain wall problem is to softly break the discrete symmetry that lifts the degeneracy between the vacua.We introduce an explicit breaking term in the potential, where ϵ is a dimensionless constant.The overall potential V (σ) + ∆V (σ) still has two minima at σ = ±u, but with a difference in the potential at these points: The probability p − of a domain ending up in the −u vacuum ('false' vacuum) is smaller compared to p + of it being in the +u vacuum ('true' vacuum), their ratio being related to the potential difference where is the potential difference between the maximum and the +u minimum.Treating the system as a three-dimensional lattice, percolation theory predicts that an infinite cluster of the false vacuum appears in space if the corresponding probability is above the threshold p c = 0.311 [79].This yields an upper bound on the bias potential for the generation of domain walls, V bias < V 0 log 1−pc pc = 0.795V 0 , which can be written as As long as the bias potential is below this limit, domain walls are created and their dynamics is mostly controlled by the surface energy density.The energy density of the wall in this regime is given by a scaling solution [80], where A ≃ 0.8 ± 0.1 is the so-called area parameter [81].A volume pressure p v ∼ V bias tends to shrink the false vacuum region.Domain walls collapse when the volume pressure overcomes the pressure from surface energy density, which happens at where C ann = 5 for Z 2 breaking [82].
Assuming that the domain walls annihilate during the radiation dominated era and annihilation happens instantaneously at t = t ann , the peak amplitude of the generated gravitational waves at present time t 0 can be expressed as [83] and the peak frequency is given by where the parameter εGW is estimated to be εGW ≃ 0.7±0.4[81].Here, g * s is the relativistic degrees of freedom for the entropy density at temperature T .Numerical simulations suggest that the gravitational wave amplitude rises as Ω GW ∝ f 3 for f < f peak and falls off as Ω GW ∝ f −1 for f > f peak .Joined by a smoothing function [84], the gravitational wave spectrum can be written as Here the low-frequency slope of the signal is required by causality to be a = 3 [84], and numerical simulations suggest that the high-frequency slope b and the width near the spectral peak c are close to unity [81,85].Uncertainties in the numerical simulations allow 0.5 ≤ b ≤ 1 and 0.3 ≤ c ≤ 3 [86].
It is evident from Eqs. (4.12)-(4.15) that the GW spectrum depends on E and V bias , both of which are determined by the Z 2 -symmetric and explicit-breaking part of the potential of the scalar σ.The specific realization of the Dirac seesaw operator does not impact the GW signal as long as σ does not have a mixing term with another scalar in its potential, although the Dirac masses of the neutrinos do depend on the underlying model [48].
It should be noted that another possibility for obtaining GWs in this model is through a first order phase transition induced by the scalar σ.The bias potential in Eq. (4.5) contains a cubic term which creates a barrier between the true and false vacua at zero temperature.However, we have explicitly checked for 10 −6 < λ < 1, 10 3 < u < 10 15 GeV and 10 −9 < ϵ < 1 that a FOPT either does not occur or is extremely weak.This is because the linear term in the potential, also controlled by the parameter ϵ, tends to remove the barrier between the two vacua.It is conceivable that in a variation of the model, a FOPT would occur, which would produce a peaked GW spectrum, more akin to that from domain walls than cosmic strings.

Results
In this section we discuss the gravitational wave signatures from cosmic strings and domain walls in the context of Majorana and Dirac seesaw models, respectively.Existing and planned interferometers probe frequencies from 10 −9 to 10 4 Hz range.In the nanoHz range (10 −9 − 10 −7 Hz), currently operating pulsar timing arrays (PTA) EPTA [87] and NANOGrav [88] have set upper bounds on the stochastic GW background, and the upcoming SKA [89] and IPTA [90] interferometers will have much greater sensitivity.µAres [91] will be sensitive to the µ-Hz to Hz band.The mHz to Hz band will be further probed by future laser interferometers LISA [92], BBO [93] and DECIGO [94,95], as well as by atomic interferometers AION [96] and AEDGE [97].Around the 100 Hz, Advanced LIGO+VIRGO [98] have set an upper limit [99,100] and their future upgrades will improve on the sensitivity by at least an order [99].Einstein Telescope (ET) [101] and Cosmic Explorer (CE) [102] are planned to operate in the same band with three orders of magnitude greater sensitivity.
For the Nambu-Goto cosmic string network, the only free parameter in Eq. (3.5) is µ of O(Λ 2 ), where Λ is the scale of the U (1) B−L symmetry breaking that generates the Majorana masses of the right-handed neutrinos.In Fig. 2 we show the GW spectrum for Λ = 10 14 , 10 13 , 10 12 , 10 11 and 10 10 GeV, corresponding to a high scale of the righthanded neutrino masses.For comparison, we also show the sensitivity and upper bounds of various interferometers spanning a large range of frequencies from nano-Hz to kilo-Hz.The characteristic shape of the cosmic string induced GW signal is a rising spectrum at low frequencies which plateaus at higher frequencies.The height of this plateau is proportional to the symmetry breaking scale.The signals for Λ ≳ 10 14 GeV are ruled out by EPTA, whereas signals for smaller scales are within the sensitivity of several interferometers.
For the Dirac seesaw, the parameter space is subject to various physical constraints that impact the formation and stability of the domain walls [83,103].If the bias potential is sufficiently small, domain walls collapse too late and may dominate the energy density of the Universe.The time at which domain walls become dominant is given by (5.1) Requiring t ann < t dom yields a lower bound on the bias potential, V bias > 4C ann A 2 E 2 /(3M 2 Pl ), which can be written as (5.2) Even if the domain walls decay before they overclose the Universe, their decay products may destroy the light element abundances created by Big Bang nucleosynthesis (BBN).
Assuming that a significant fraction of the energy density of the domain walls is converted into energetic particles, constraints on energy injection at the epoch of BBN require t ann ≲ t BBN ≃ 0.01 sec [104,105], which can be written as (5.3) Equations (5.2) and (5.3), together with Eq. (4.9), constrain the parameter space for annihilation of domain walls and subsequent gravitational wave production.In terms of the scalar VEV u and bias potential V bias , and choosing A = 0.8, C ann = 5, these constraints can be expressed as t ann < t dom : V bias GeV 4 > 6.42 × 10 −37 λ u GeV V bias < 0.795V 0 : (5.6) Finally, from Eq. (2.4), if we assume that the mediator fermion mass is below the Planck scale O(10 19 ) GeV, and the heaviest light neutrino mass is around O(0.1) eV, we find Here we have assumed a single mediator responsible for the O(0.1) eV neutrino mass, and a single coupling y = y L = y R associated with it.For Yukawa couplings y ≳ O(10 −2 ), this implies an upper bound on the scale of Z 2 symmetry breaking, u ≲ O(10 11 ) GeV.The constraints of Eqs.(5.4)-(5.6)are depicted by the gray shaded regions in Fig. 3 for (a) λ = 1 and (b) λ = 10 −3 .For smaller λ, the upper left region expands while the other two regions shrink, as expected from Eqs. (5.4)- (5.6).The peak frequencies of the gravitational waves f peak = 10 −8 , 10 −6 , . . . 10 2 Hz, are marked by dots.Their colors represent the amplitude of the GW signal at the corresponding peak frequency.We find that amplitudes above Ω peak GW h 2 ∼ 10 −6 are ruled out by Eq. (5.4), while peak frequencies below f peak ∼ 10 −8 Hz are ruled out by Eq. (5.5).Interestingly, the allowed region can still generate GW signals within the sensitivity of various interferometers.
The four benchmark points listed in Table 1 are chosen from the allowed parameter space.The last column of the table gives the upper bound on the Yukawa coupling, assum-  1. ing that the mediator mass lies below the Planck scale.We note that the values of such Yukawa couplings cover the range of the third family charged fermion Yukawa couplings in the SM, and only exceed this range by less than an order of magnitude.The gravitational wave spectra for these benchmark points are shown in Fig. 4. Benchmark point 1 can be probed by SKA, while 2 and 3 can be probed by µAres, LISA, AEDGE, DECIGO, BBO, AION, and 4 by Advanced LIGO+VIRGO, ET and CE, among others.We set b = c = 1 for the spectral shape, but slightly different values still yield a peaked spectrum for a = 3.In Fig. 5, we show the signal-to-noise ratio [106]  of the generated GW spectrum in the parameter space of the model for various interferometers operating from nHz to kHz.Here Ω exp (f )h 2 is the effective strain noise power spectral density [107], τ = 4 years is the typical observation time, and f min and f max define the range of frequencies in which an interferometer is sensitive.We only show the parameter space for SNR ≥ 10, which is the threshold for detection in individual interferometers.The results show that detectable signals arise if the domain walls annihilate when they are close to dominating the energy density of the Universe.
The main difference between the signals for the Majorana seesaw model and the Dirac seesaw model is their spectral shape.While cosmic string signals for the former are mostly flat for observable frequencies, domain wall signals for the latter are peaked.We expect that cosmic string signals should be detected at multiple interferometers at different frequency bands, whereas domain wall signals are likely to be detected in only a narrow frequency range.Such a detection will provide valuable information about the nature of neutrino mass genesis and will complement results from neutrinoless double beta experiments.The blue curves correspond to the GW signal from the seesaw mechanism alone, and the red curves correspond to the combined signal from the seesaw mechanism and supermassive black hole binaries (SMBHBs) with A BHB = 10 −15.6 and γ BHB = 4.7.

Implications of NANOGrav data
Finally, we comment on the implications of the recent PTA signals on our Majorana and Dirac seesaw models.For concreteness, we compare the predictions of the models to the NANOGrav 15-year dataset [50][51][52][53][54][55][56][57] using the PTArcade software [108,109].In Fig. 6 we show the posterior distributions and allowed regions of the model parameters.
For the Majorana seesaw model, we show the 68% (darker blue) and 95% (lighter blue) Bayesian credible regions for the single parameter Λ.For the Dirac seesaw model, we show both the reconstructed 1D marginalized distributions (diagonal plots) and the 2D distribution (off-diagonal), along with 68% (darker blue) and 95% (lighter blue) Bayesian credible regions.In the same plots, we also show the fit to the combined contribution (red) from the respective models and from inspiraling supermassive black hole binaries (SMBHBs) [51] where H 0 = h×100 kms −1 Mpc −1 is the Hubble constant, using the best fit parameters from numerical simulations A BHB ≃ 10 −15.6 and γ BHB ≃ 4.7 [53].We notice that the inclusion of the SMBHB contribution does not noticeably enlarge the parameter space of either model.
For the Majorana seesaw model, the U (1) B−L breaking scale Λ is in the range 10 14 GeV-10 14. 16 GeV at the 68% credible level.For the Dirac seesaw model, the Z 2 breaking scale u and bias potential V bias are between 10 5.5 GeV-10 −5.76 GeV and 10 −2.8 GeV 4 -10 −1.8 GeV 4 , respectively, at the 68% credible level.
We choose the best fit points of the respective models and plot the GW spectrum in Fig. 7.The Majorana seesaw model with B − L symmetry broken at Λ = 1.2 × 10 14 GeV generates the maroon curve, and the Dirac seesaw model with benchmark point 1 yields the dark blue curve at PTA frequencies.The blue violins show the NANOGrav 15yr data and the orange violins show the second data release from EPTA [58][59][60][61][62][63].For comparison, we also show the signal (black dashed line) from SMBHBs with the best fit values of A BHB and γ BHB .The spectral slope of the cosmic string signal from the Majorana mass model provides an even poorer fit to the PTA signal than the SMBHB signal.On the other hand, the domain wall signal from the Dirac mass model fits the PTA result remarkably well.This suggests that Majorana mass generation from the spontaneous breaking of a gauged B − L symmetry at around 10 14 GeV is not favored by PTA data, while Dirac mass generation from the spontaneous breaking of the Z 2 symmetry in Eq. ( 2.3) at around 10 5−6 GeV is favored.While this statement relies on several assumptions about the astrophysical background in PTA data, and is not necessarily an inevitable outcome of Dirac or Majorana neutrino mass generation, it provides an example of the complementary information one may expect from GW signatures of the respective seesaw models.In the absence of any direct evidence of the nature of neutrino mass, such indirect probes are important and encourage further exploration in this direction.

Discussion and outlook
We have considered a novel possibility to distinguish between Dirac and Majorana neutrino mass through the difference in the gravitational wave spectra, assuming a seesaw mechanism in each case.For definiteness, we considered two simple classes of the respective type-I seesaw models, both inspired by the assumption that small neutrino masses be generated without tiny Yukawa couplings.The Majorana seesaw is assumed to involve spontaneous breaking of lepton number, generating cosmic strings, while the Dirac seesaw is assumed to involve spontaneous breaking of a discrete Z 2 symmetry (the minimal choice) leading to domain walls.The resulting very different shapes of the GW spectra, after the cosmological relics decay, flat in the Majorana seesaw case, and peaked in the Dirac seesaw case may help distinguish between the two seesaw mechanisms.While the method by itself is not conclusive, it may provide valuable additional information about the mass generation mechanism of neutrinos, which is combinable with other indirect evidence to determine the nature of neutrino mass.
We emphasize that the primary difference between Majorana and Dirac neutrino masses is whether or not lepton number is broken.In the case of Majorana masses, it is broken by two units.In a well motivated scenario, the lepton number symmetry, or equivalently, the B − L symmetry, is exact at ultraviolet scales, but is spontaneously broken when a scalar charged under the gauged U (1) B−L symmetry gets a nonzero VEV.This gives Majorana mass to the right-handed neutrinos, which further generates small masses for the SM neutrinos via the type-I seesaw mechanism.The breaking of the U (1) B−L symmetry triggers the creation of cosmic strings in the early Universe.The string network loses energy via the production of string loops, some of which emit gravitational waves.The GWs have a flat spectrum over a wide range of frequencies, and their amplitude is related to the scale of symmetry breaking.Hence, the detection of a flat spectrum of stochastic gravitational wave background may imply a Majorana nature of the neutrinos and shed light on the scale at which such masses are generated.
On the other hand, the Dirac seesaw mechanisms to generate a small Dirac mass for the SM neutrinos utilize the effective operator Lν R Hσ.To keep lepton number symmetry unbroken, and to prohibit a tree-level Dirac mass term for the SM neutrinos requires ν R and σ to be non-trivially charged under a Z 2 symmetry, which is spontaneously broken by the VEV of σ.The breaking of a discrete symmetry creates a domain wall network, which poses a threat to the standard cosmology if it is long-lived and/or dominates the energy density of Universe.Domain walls can be made to annihilate by softly breaking the discrete symmetry, thereby lifting the degeneracy between the Z 2 symmetric vacua and creating a bias potential that tends to collapse the walls.This leads to a characteristic GW signal peaked at a frequency determined by the scale of spontaneous and soft symmetry breaking.
Interestingly, there are several ultraviolet realizations of the Dirac seesaw operator at the tree-level and at one-loop level.The characteristics of the GW signals produced in the realizations of the Dirac seesaw operator are identical, since these are solely determined by the potential of the scalar σ, as long as σ does not mix with other scalars.We have shown that depending on the parameter space, such signals may be probed by various terrestrial and satellite-based interferometers.We compared the predictions of the two classes of mass generation models to the NANOGrav dataset and found that Dirac seesaw models are favored.
The scope of our proposal is limited in various respects.While our setup is wellmotivated, it is possible that the B−L symmetry is explicitly broken in the case of Majorana neutrinos, or that the scale of breaking is small enough that the signal is undetectable.Since our discussion is focused on a subclass of Dirac and Majorana mass generation models, our results are representative of the type of models, and not of the specific nature of the neutrinos (which can be probed by neutrinoless double beta decay experiments).It is also arguable that other models, unrelated to neutrino mass generation, can produce similar signals from cosmic strings or domain walls.Furthermore, we have assumed a typical radiation-dominated Universe after inflation, and modifications to the standard cosmic evolution may affect the GW spectra.Nevertheless, the model classes we have investigated are simple, minimalistic, and highly predictive under the set of assumptions made.Hence, our work should be viewed as an attempt to extract whatever meaningful information one can from GW signals in the absence of a positive result from neutrinoless double beta decay experiments.In particular, given that a nonobservation of neutrinoless double beta decay does not necessarily imply that neutrinos are Dirac particles, the peaked domain wall signal may provide indirect evidence of the Dirac nature of neutrinos when combined with data from other terrestrial experiments.

Figure 2 .
Figure 2. Gravitational wave spectrum induced by cosmic strings generated via the spontaneous breaking of the gauged U (1) B−L symmetry responsible for Majorana mass of the neutrinos.Λ denotes the scale of symmetry breaking.

6 ,Figure 3 .
Figure 3. Parameter space of the Dirac seesaw model for creation and subsequent annihilation of the domain wall network.Gray shaded regions show the parameter space ruled out by physical constraints.Dotted lines in the allowed parameter space represent contours of peak frequency, and the colors of the dots indicate the amplitude of the gravitational wave signal at that frequency.

Figure 4 .
Figure 4. Gravitational wave spectrum from annihilation of domain walls created by soft-breaking of the Z 2 symmetry in a Dirac neutrino mass model.Benchmark points 1 -4 are listed in Table1.

Figure 5 .
Figure 5. Signal-to-noise ratio of the gravitational wave signal in the Dirac mass model (setting λ = 1) at some interferometers.Gray regions are in conflict with the constraints shown Fig. 3a.Colored regions represent SNR ≥ 10.

4 ]Figure 6 .
Figure 6.Posterior distributions for the Majorana seesaw (MS) model parameter Λ (left) and the 1σ and 2σ credible regions for the Dirac seesaw (DS) model parameters u and V bias (right).The blue curves correspond to the GW signal from the seesaw mechanism alone, and the red curves correspond to the combined signal from the seesaw mechanism and supermassive black hole binaries (SMBHBs) with A BHB = 10 −15.6 and γ BHB = 4.7.

Figure 7 .
Figure 7. Gravitational wave spectrum at PTA frequencies from the best fit parameters of the Dirac mass model (benchmark point 1 in Table 1) labeled DW, and Majorana mass model (B − L breaking scale Λ = 1.2 × 10 14 GeV) labeled CS.Recent results from NANOGrav 15yr data (blue violins) and EPTA data (orange violins) are also shown.The dashed line is the signal from SMBHBs with A BHB = 10 −15.6 and γ BHB = 4.7 [51].

Table 1 .
Benchmark points for gravitational wave signals from domain walls with λ = 1.