Gravitational waves from phase transitions and cosmic strings in neutrino mass models with multiple majorons

We explore the origin of Majorana masses within the majoron model and how this can lead to the generation of a distinguishable primordial stochastic background of gravitational waves. We first show how in the simplest majoron model only a contribution from cosmic string can be within the reach of planned experiments. We then consider extensions containing multiple complex scalars, demonstrating how in this case a spectrum comprising contributions from both a strong first order phase transition and cosmic strings can naturally emerge. We show that the interplay between multiple scalar fields can amplify the phase transition signal, potentially leading to double peaks over the wideband sloped spectrum from cosmic strings. We also underscore the possibility of observing such a gravitational wave background to provide insights into the reheating temperature of the universe. We conclude highlighting how the model can be naturally combined with scenarios addressing the origin of matter of the universe, where baryogenesis occurs via leptogenesis and a right-handed neutrino plays the role of dark matter.


Introduction
The discovery of gravitational waves (GWs) [1] opens new opportunities to test physics beyond the standard model (BSM).This is particularly interesting for those models currently evading constraints from colliders and, more generally, from laboratory experiments.Even though GWs have so far been detected only from astrophysical sources, there are many different processes in the early universe that could lead to the production of detectable primordial stochastic GW backgrounds.In particular, a production from the vibration of cosmic strings [2] and from strong first order phase transitions [3][4][5] provide quite realistic and testable mechanisms within various extensions of the standard model (SM) [6][7][8].
These two GW production mechanisms are usually studied separately.In this paper we show how within the majoron model [9], an extension of the SM explaining neutrino masses and mixing, a GW spectrum is produced where both sources can give a non-negligible contribution and fall within the sensitivity of planned experiments.In the majoron model a type-I seesaw [10][11][12][13][14] Lagrangian results as the outcome of a global U (1) L spontaneous symmetry breaking and Majorana masses are generated by the vacuum expectation value (VEV) of a single complex scalar field.The massless Goldstone boson, identified as the imaginary part of the complex scalar, is dubbed as majoron.The model can also nicely embed leptogenesis for the explanation of the matter-antimatter asymmetry of the universe [15].
The idea that a strong first order electroweak phase transition associated to the lepton number symmetry breaking can generate a stochastic GW background has been explored in [16].In this case a coupling of the complex scalar to the SM Higgs field was considered.A phase transition within the dark sector of the majoron model, disconnected from the electroweak phase transition, was considered in [17], where non-renormalisable operators and explicit symmetry breaking terms have been included in order to enhance the signal.Moreover, a low-scale phase transition, in the keV-MeV range was also considered in order to reproduce the NANOGrav putative signal at very low frequencies (∼ 10 −9 Hz) [18].
A first order phase transition from U (1) L -symmetry breaking in the dark sector, with no coupling of the complex scalar field to the SM Higgs field, was also considered in [19] without resorting either to explicit symmetry breaking terms or to non-renormalizable operators.Both the case of low and high scale phase transition were explored.It was found that at low scales the NANOGrav result cannot be explained, unless one invokes some enhancement from some unaccounted new effect.On the other hand, it was found that at high energy scales the signal can be sufficiently large to fall within the sensitivity of future experiments such as µAres [20], DECIGO [21], AEDGE [22], AION [23], LISA [24], Einstein Telescope (ET) [25], BBO [26] and CE [27].However, this result relied on the introduction of an external auxiliary real scalar field undergoing its own phase transition occurring prior to the complex scalar field phase transition.Once the auxiliary scalar gets a VEV, its mixing with the complex scalar field generates a zero-temperature barrier described by a cubic term in the effective potential of the latter, leading to a strong first order phase transition and detectable GW spectrum.
In this paper, we show how the role of the auxiliary field can be nicely played by a second complex scalar in a multiple majoron model.We discuss neutrino mass models with spontaneous breaking of multiple global lepton number symmetries, typically with hierarchical scales.The three right-handed (RH) neutrino masses are then generated by different complex scalars each undergoing its own independent phase transition occurring, in general, at different energy scales and breaking lepton number along a specific direction in flavour space.We have, then, what could be referred to as a (RH neutrino) flavoured majoron model.Importantly, we show that a contribution from the vibration of cosmic strings generated from the spontaneous breaking of the global lepton number symmetry has also to be taken into account to derive the GW spectrum of these models.The overall spectrum then is the sum of contributions from both production mechanisms: a contribution from strong first order phase transitions and a contribution from the vibration of cosmic strings.For sufficiently strong phase transitions, the resultant signal looks like one or more peaks (from phase transition) over a slanted plateau (from cosmic string).
The paper is organised as follows.In Section 2 we review the traditional single majoron model where the right-right Majorana mass term, with three RH neutrinos, is generated by a single complex scalar field breaking total lepton number symmetry.The differences in the Majorana masses are then to be ascribed to different couplings.Even in this traditional setup we point out that a GW production from the vibration of cosmic strings, not accounted for in previous works, should be considered and can give a detectable signal.In Section 3 we extend the model with an additional complex scalar with its respective global lepton number symmetry, whose spontaneous breaking gives mass to the two lighter RH neutrinos.In this way only two distinguished phase transitions occur with hierarchical energy scales.We show that the resulting GW spectrum is, in general, the sum of two contributions, one from the lower scale phase transition and one from the vibration of cosmic strings created at the highest scale symmetry breaking.The corresponding phase transition does not produce a sizeable contribution to the GW spectrum, but it results into a VEV of the complex scalar field that generates a term entering the effective potential describing the second phase transition at a lower scale.This term strongly enhances the production of GWs during the second phase transition.In this way the high scale complex scalar associated with the majoron field provides the external auxiliary scalar that had to be assumed in [19], so that the model is self-contained and does not rely on external assumptions.Finally, in Section 4 we consider the case when all three RH neutrino masses are associated with different complex scalars, each charged under a different global lepton number symmetry.At high temperatures one has the restoration of a U (1) While the temperature decreases, a sequential breaking of each U (1) L I symmetry occurs at a different scale accompanied by a different phase transition.In this case we show that the GW spectrum now can receive a contribution from both the two lower scale phase transitions and still from the vibration of cosmic strings at the highest scale symmetry breaking.We show that such a spectrum may have twin peaks from phase transition signals over a slightly sloped plateau of the cosmic string signal.We draw conclusions in Section 5 and point out that the GW spectrum of the model can provide us important information about the reheating temperature of the universe, and that the model fits naturally within a unified framework of solving the puzzles of baryon asymmetry and dark matter.

Primordial GW stochastic background in the single majoron model
In this section, we first review the main features of the single majoron model and then discuss the generation of a stochastic background of primordial GWs.

The single majoron model
The traditional single majoron model is a simple extension of the SM [9], where the spontaneous breaking of a global U L (1) symmetry generates a Majorana mass term for the RH neutrinos.The SM field content is then augmented with N RH neutrino fields N I (I = 1, 2, . . ., N ) and a complex scalar singlet, where the real component is CP -even and the imaginary component is CP -odd.The new scalar ϕ has a tree level potential V 0 (ϕ).For definiteness, we consider the well motivated case N = 3.The tree-level extension of the SM Lagrangian is then given by where Φ is the dual Higgs doublet.In the early universe, above a critical temperature T c , one has ⟨ϕ⟩ = 0 so that the RH neutrinos are massless.Moreover, since the lepton doublets L α and the RH neutrinos N I have L = 1, and ϕ has L = −2, lepton number is conserved.Below T c , the U L (1) symmetry is broken and the scalar ϕ acquires a vacuum expectation value ⟨ϕ⟩ = v 0 e iθ 0 / √ 2. In this way the RH neutrinos become massive with Majorana masses M I = v 0 λ I / √ 2. This leads to lepton number violation and small Majorana masses for the SM neutrinos via type-I seesaw mechanism.We assume that T c ≫ T ew ∼ 100 GeV, so that the majoron phase transition occurs prior to the electroweak phase transition and, therefore, the Majorana mass term is generated before the Dirac mass term.
Let us consider the simple tree level potential where λ is real and positive, in a way that the potential is bounded from below, and µ 2 is real and positive to ensure the existence of degenerate nontrivial stable minima with ⟨ϕ⟩ = v 0 e iθ / √ 2 with 0 ≤ θ < 2π and where v 0 ≡ µ 2 /λ.After spontaneous symmetry breaking, we can rewrite ϕ as where S is a massive field with m 2 S = 2λv 2 0 and J is the majoron, a massless Goldstone field.Moreover, RH neutrino masses M I = λ I v 0 / √ 2 are generated by the VEV of ϕ and these lead to a light neutrino mass matrix given by the (type-I) seesaw formula where v ew = 246 GeV is the standard Higgs VEV.Notice that the potential in Eq. (2.3) corresponds to a minimal choice where we are neglecting possible mixing terms between the new complex scalar field ϕ and the standard Higgs boson.In this way, the phase transition involves only the dark sector, consisting only of ϕ and the three RH neutrinos.Moreover, we are not considering non-renormalisable terms, so that the model is UV-complete.Since all minima are equivalent, one can always redefine θ in a way that the symmetry is broken along the direction θ 0 = 0, without loss of generality.The minimum of the potential lies along the real axis and, for all purposes, one can consider the potential as a function of φ, so that one has: Let us now discuss the generation of a primordial stochastic background of GWs.There are two possible sources in the majoron model.The first is an associated strong first order phase transition [19] that we discuss in the subsection 2.2.The second is the network of cosmic strings generated by the breaking of the global U (1) L symmetry that we discuss in the subsection 2.3.The latter has not been discussed before within a majoron model, though it is analogous to the U (1) B−L spontaneous symmetry breaking discussed, for example, in [28][29][30][31][32].
1 Notice that J is not the imaginary part of ϕ but rather related to the phase fluctuation around the vev ⟨ϕ⟩.Analogously, S is not the real part of ϕ but rather related to the radial fluctuation.This is easy to see if one rewrites the fluctuation δϕ about the vev ⟨ϕ⟩ starting from ϕ = [(⟨φ⟩ + δφ)/ √ 2] e i (θ 0 +δθ) .One can then easily identify v0 = ⟨φ⟩, S = δφ and J = v0 δθ.This shows that Eq. (2.4) is equivalent to

Stochastic GW background from first order phase transition
The scalar field and the three RH neutrinos form what we refer to as the dark sector.
The dark sector interacts with the SM sector only via the Yukawa interactions.In the early universe finite temperature effects need to be taken into account.They will drive a phase transition, occurring in the dark sector, from the metastable vacuum at ϕ = 0, where lepton number is conserved and RH neutrinos are massless, to the true stable vacuum at ϕ = v 0 / √ 2, where lepton number is non-conserved and RH neutrino are massive (for a recent review on early universe phase transitions and GWs see [33]).They are described in terms of a finite-temperature effective potential V T eff (ϕ).At temperatures above a critical temperature T c , finite temperature effects will induce symmetry restoration [34].When temperature drops down the critical temperature, the phase transition occurs and, in the zero temperature limit, the tree-level potential V 0 (ϕ) is recovered, in the broken symmetry phase. 2  The finite-temperature effective potential can be calculated perturbatively at one-loop [35] and is given by the sum of three terms, where the zero-temperature one-loop contribution V 0 1 (ϕ) is given by the Coleman-Weinberg potential.This can be written, using cut-off regularization, as [35][36][37][38] The pre-factor of two in the second line accounts for two degrees of freedom for each RH neutrino species.The one-loop thermal potential is given by [36-38] where the thermal functions are The functions m 2 ϕ (ϕ) and M 2 I (ϕ) are the shifted masses given by 2 Notice that the reheating temperature of the universe TRH needs to be higher than Tc for both symmetry restoration and symmetry breaking to occur.If it is lower, the universe history starts directly in the broken phase and there is no phase transition.For this reason, finding evidence for a phase transition and establishing the value of Tc would straightforwardly place a lower bound on TRH. and where we specialized their dependence as a function of φ since, even when thermal effects are included, all the study of the dynamics can be done along the real axis of ϕ without loss of generality.This time the sum over the RH neutrino species, the only fermions coupling to ϕ, should only include those that are fully thermalised prior to the phase transition, while we can neglect the contribution from those that are not.RH neutrinos thermalise at a temperature [39,40] T eq I ≃ 0.2 where m eq ≡ [16π 5/2 g ⋆ ρ /(3 is the usual effective equilibrium neutrino mass and we simply defined vew = v ew / √ 2 ≃ 174 GeV.Note that the quantity M = v2 ew /m eq = (3 GeV g SM ρ /g ⋆ ρ is independent of the electroweak scale v ew .The condition for the thermalisation of the RH neutrino species N I prior to the phase transition can then be written as The equilibration temperature T eq and the condition Eq. (2.14) can also be conveniently expressed in terms of the dimensionless RH neutrino decay parameters obtaining, respectively,

.16)
Taking into account the measured values of the solar and atmospheric neutrino mass scales, from the seesaw formula it can be shown that all three RH neutrino species can satisfy the condition of thermalisation and this is what we assume for simplicity following [19]. 3Of course we also assume T RH ≫ T c for the phase transition to occur (as noticed in footnote 1).Another important thermal effect to be taken into account is that the tree-level shifted mass have to be replaced by resummed thermal masses [41] where the Debye mass Π ϕ is given by In this expression one has d scalar = 2 for the case of a complex scalar we are considering.The quantity M denotes either the mass of the heaviest RH neutrino in the case of hierarchical RH neutrino mass spectrum (in which case N = 1), or a common mass in the case of quasidegenerate RH neutrinos (in which case N is the number of RH neutrinos).This allows us to reduce the number of parameters while spanning the space between N = 1 (hierarchical RH neutrinos) and N = 3 (quasi-degenerate RH neutrinos).
With this replacement and neglecting O((M I /T ) 6 ) terms in the high temperature expansion of the thermal functions, one obtains the dressed effective potential [19,42,43] In this expression we introduced where T 0 is the destabilisation temperature defined by The dimensionless constant coefficients D and A are given by Finally, the dimensionless temperature dependent coefficient λ T is given by Notice that one has to impose M 1 < m S in order for the massive scalar S to decay into RH neutrinos in a way that its thermal abundance does not overclose the universe.However, this condition is easily satisfied, since the scalar and RH neutrino masses are roughly of the same order-of-magnitude as v 0 .At very high temperatures the cubic term in the effective potential (2.19) is negligible and one has symmetry restoration.However, while temperature drops down, there is a particular time when a second minimum at a nonzero value of φ forms.When temperature further decreases, a barrier separates the two coexisting minima.The critical temperature T c is defined as that special temperature when the two minima become degenerate.Until this time, the probability that a bubble of the false vacuum nucleates vanishes but below the critical temperature it is nonzero.The nucleation probability per unit time and per unit volume can be expressed in terms of the Euclidean action S E as [44]: (2.24) At finite temperatures one has where the quantity S 3 is the spatial Euclidean action given by The physical solution for φ minimizing S 3 (φ, T ) can be found solving the EoM with boundary conditions (dφ/dr) r=0 = 0 and φ(r → ∞) = 0. Since for T ≥ T c the nucleation probability vanishes, one has lim T →T − c S E → ∞, while on the other hand lim T →T 0 S E → 0, so that at T 0 all space will be in the true vacuum and the phase transition comes to its end. 4If the phase transition is quick enough, then one can describe the phase transition as occurring within a narrow interval of temperatures about a particular value T ⋆ such that T c > T ⋆ > T 0 .The temperature T ⋆ is referred to as the phase transition temperature and it is usually identified with the percolation temperature, defined as the temperature at which the fraction of space still in the false vacuum is 1/e.The fraction of space filled by the false vacuum at time t is given by [46,47] where a(t) is the scale factor and v w is the bubble wall velocity.Therefore, P (t ⋆ ) = 1/e corresponds to I(t ⋆ ) = 1, where t ⋆ ≡ t(T ⋆ ).It can be shown [48] that at T ⋆ the Euclidean action has to satisfy where H ⋆ = H(t ⋆ ).This equation allows to calculate T ⋆ and S E (T ⋆ ) having derived S E (T ) from the solution of the EoM.The calculation of the GW spectrum produced during the phase transition is characterised by two quantities.The first is β ≡ Γ/Γ, the rate of variation of the nucleation rate.Its inverse, β −1 , gives the time scale of the phase transition.In our case, we are interested in the scenario of fast phase transition, for β −1 ≪ H −1 , so that, with a first order expansion of the Euclidean action about t This provides a sufficiently good approximation for β/H ⋆ ≳ 100 [48].The second quantity characterising the phase transition is the strength of the phase transition α defined as where ε(T ⋆ ) is the latent heat released during the phase transition and ρ(T ⋆ ) is the total energy density of the plasma, including both SM and dark sector degrees of freedom.The latent heat can be calculated using where ), and in the first relation, from thermodynamics, ∆s is the entropy density variation and the free energy of the system has been identified with the effective potential.Notice that in our case, V T⋆ eff (ϕ false 1 ) = 0. Also notice that the constraint β/H ⋆ ≫ 1 for the validity of Eq. (2.30) implies a constraint α ≪ 1, since the two quantities are not completely independent of each other with β/H ⋆ ∝ α −2 [49].For definiteness, we will then impose α ≤ 0.3, corresponding typically to β/H ⋆ ≳ 100.The total energy density of the plasma can be expressed, as usual, as The number of the total ultrarelativistic degrees of freedom g ρ (T ) is in this case given by the sum of two contributions, one from the SM and one from the dark sector, explicitly, one has g ρ (T ) = g SM ρ (T ) + g dark ρ (T ), where g SM ρ (T ⋆ ) = 106.75 and g dark ρ Let us now calculate the GW spectrum defined as where ρ c0 is the critical energy density and ρ GW0 is the energy density of GW, produced during the phase transition, both calculated at the present time.We assume that the phase transition occurs in the detonation regime, i.e., with supersonic bubble wall velocities, v w ≥ c s = 1/ √ 3, that is typically verified in the regime α ≤ 0.3 we are considering.Moreover, the dominant contribution to the GW spectrum typically comes from sound waves in the plasma, with a sub-dominant contribution from magnetohydrodynamic (MHD) turbulence, so that [24].
A numerical fit to the the GW spectrum that is the result of semi-analytical methods and at the same time takes into account the results of numerical simulations, quite reliable in the regime α ≤ 0.3 we are considering, yields [24,50,51] where the redshift factor r gw (t ⋆ , t 0 ), evolving Ω gw⋆ ≡ ρ gw⋆ /ρ c⋆ into Ω gw0 ≡ ρ gw0 /ρ c0 , is given by [52] r gw (t and in the numerical expression we used: g γ = 2, g S⋆ = g ρ⋆ , g S0 = 43/11 ≃ 3.91 , Ω γ0 = 0.537 × 10 −4 (0.6875/h) 2 .Replacing the expression for the mean bubble separation R ⋆ = (8π) 1/3 v w /β, valid in the detonation regime we are assuming, we obtain the numerical expression (2.37) The normalised spectral shape function is given by S sw (f ) ≃ 0.687 S sw (f ) with where f sw is the peak frequency given by (2.39) Notice that we have normalized the number of degrees of freedom to the SM value since we are discussing phase transitions at or above the electroweak scale.The efficiency factor κ(α) measures how much of the vacuum energy is converted to bulk kinetic energy.We adopt Jouguet detonation solutions since we assume that the plasma velocity behind the bubble wall is equal to the speed of sound.Then, the efficiency factor is [53,54] κ(α) ≃ α 0.73 + 0.083 and the bubble wall velocity is v w (α) = v J (α), where Jouguet solutions provide a simple prescription but a rigorous description would require numerical solutions of the Boltzmann equations [54].The prefactor Ω gw in Eq. (2.35) is calculated from numerical simulations and a recent analysis shows that in the regime we are considering, for α ≤ 0.3 and v w = v J ≳ c s , it takes values approximately in the range Ω gw = 10 −3 -10 −2 [51], with the exact value depending on additional parameters necessary to simulate the GW production from sound waves, such as friction, that we do not describe in our analysis.For this reason we show in all results bands of GW spectra corresponding to this range of values for Ω gw rather than a single curve.This should also account for the use of simple Jouguet solutions for v w rather than solutions of Boltzmann equations, also depending on friction as additional parameter.Finally, notice that in Eq. (2.35) there is also a suppression factor Υ(α, β/H ⋆ ) < 1 which decreases with the strength of the phase transition and is given by [55,56]: where the product of the lifetime of the sound waves τ sw with the Hubble expansion parameter at the time of the phase transition can, in turn, be expressed in terms of α and β/H ⋆ as For the MHD turbulence contribution, the GW spectrum is given by [24] h 2 Ω tb0 (f ) = 3.28 × 10 −4 106.75 where κ tb (α) = ϵκ(α), with ϵ ≈ 0.05 representing the fraction of bulk motion of the plasma which is turbulent.The other quantities are given by (2.47) Let us now calculate the GW spectrum within the majoron model.If we consider the minimal tree level potential in Eq. (2.3), there is a simple solution of the EoM for the Euclidean action given by [19,37] where we defined the dimensionless parameter In Fig. 1 we show, with blue bands, the GW spectra corresponding to the three benchmark choices for the values of v 0 , λ and M in Table 1.The thin dashed and dot-dashed lines correspond to sound wave and MHD turbulence contributions to the GW spectra, whereas the thick lines represent the combined spectra.The contribution from sound wave is dominant at the peak amplitude, while the MHD turbulence modifies the high-frequency tail.We also show the sensitivity regions of LIGO [57,58] and some planned/proposed  experiments, µAres [20], LISA [24], BBO [26], DECIGO [21], AEDGE [22], AION [23], ET [25] and CE [27].
Considering that these three choices are those found in a scan that maximise the signal in respective peak frequencies, it should be clear that the contribution from phase transitions in the case of the minimal model is far below the experimental sensitivity.In this way we confirm the conclusions found in [19].In the next subsection we point out, however, that at least for large values of v 0 , the contribution from cosmic strings could be detectable in future experiments even in this minimal model.
Before concluding we should also mention that one could think to add an explicit symmetry breaking cubic term in the tree level potential.However, as noticed in [19], its coefficient is upper bounded by the observation that it unavoidably generates also a linear term in the effective potential.This tends to remove the barrier between the two vacua so that, if the coefficient is too large, there is no first order phase transition and, therefore, no GW production.For this reason, we do not pursue this scenario.

GW from global cosmic strings
Spontaneously breaking the U (1) L symmetry at high energies by the complex scalar ϕ generates a global cosmic string network, which dominantly radiates Goldstone bosons, and sub-dominantly emits gravitational waves [59].Compared to the Nambu-Goto stringinduced almost flat gravitational wave spectrum associated with a gauged symmetry breaking, the global cosmic string-induced gravitational waves are typically suppressed, and their amplitude mildly falls off with frequency for most of the spectrum of interest.This makes their detection at interferometers more challenging, unless the symmetry breaking scale v 0 is above 10 14 GeV.In this section we briefly review the dynamics of global cosmic strings using the velocity-dependent one scale (VOS) model [60][61][62][63][64][65] and the associated gravitational wave spectrum following [66]. 5 The global cosmic string network consists of horizon sized long strings that randomly intersect and form sub-horizon sized loops at the intersections.The network shrinks and loses energy with time, but eventually enters a scaling regime where the average inter-string separation scale L, and the ratio of the energy density of the network to the total background energy density remain constant.For Nambu-Goto strings, this energy is radiation from the string loops predominantly in the form of gravitational waves.However, for global strings the leading mode of energy radiation is from the emission of Goldstone particles, and only a fraction of the energy is radiated as gravitational waves.
The energy density of the global string network can be expressed as where µ(t) is the energy per unit length of the long strings, and ξ(t) is a dimensionless parameter which represents the number of long strings per horizon volume.While for Nambu-Goto strings, µ is a constant, for global strings it has a logarithmic dependence on the ratio of two scales, a macroscopic scale L(t) close to the Hubble scale, and a microscopic scale δ(t) ∼ 1/(λv 0 ) representing the width of the string core, (2.52) 5 We assume that symmetry breaking and the consequent formation of the cosmic string network occur after inflation.In the case that the symmetry breaking happens during inflation, the GW spectrum gets modified by inflation, see for example [67] and references therein.We also assume that the cosmic strings are stable.This is always the case if a matter parity remains unbroken after U (1)L symmetry breaking (for example, see discussion in [29]).
The evolution of the inter-string separation scale L(t) and the average long string velocity v are given by a system of coupled differential equations, ) The first term on the RHS of Eq. (2.54) represents dilution effect from Hubble expansion.
The second term gives a negligible thermal friction effect with a characteristic scale ℓ f ∝ µ/T 3 .The third term stands for the loop chopping effect, where c is the rate of loop chopping.The fourth term represents the backreaction due to the Goldstone emission.The quantity q is a momentum parameter.In analogy with Nambu-Goto strings, the solution of Eqs.(2.54) and (2.55) can be expressed as where ∆ ≡ κ/(N (q + c)), κ ≡ sv 5 /(1 − ∆) 5/2 , and n = 3, 4 corresponds to matter and radiation domination, respectively.Fitting data extracted from the simulation results in Refs.[68][69][70], the VOS model parameters can be approximated as [66] {c, q, κ} ≃ {0.497, 0.284, 5.827}. (2.58) Since ξ(t) = t 2 /L 2 (t) from Eq. (2.51), Eq. (2.56) can be used to express ξ as a function of N (t).Eq. (2.53) then expresses N (t) as a function of t.Similarly, v can be expressed as function of t from Eq. (2.57).
Assuming that the loop size during formation of the string network is given by ℓ i ∼ αt i , where α is a dimensionless O(1) parameter not to be confused with the strength of the phase transition, and the fraction of energy density of the strings contributing to gravitational wave F α ∼ 0.1 (for α ∼ 0.1) [71,72], the formation rate of string loops is given by where F α ∼ 1 is the loop size distribution function, and E loop ≡ c v ξ 3/2 is the loop emission parameter.
After formation, the string loop rapidly oscillates and radiates energy in the form of Goldstone particles and gravitational waves until disappearing completely with a rate [59] where we assume the benchmark values Γ ∼ 50 [71,[73][74][75] and Γ a ∼ 65 [76,77].The size of a loop initial length ℓ i = α t i at a later time can be expressed as where the second and third terms represent the decrease in loop size for gravitational wave emission and Goldstone emission, respectively.It is useful to decompose the radiation into a set of normal modes fk = 2k/ l, where k = 1, 2, 3, . .., and l ≡ ℓ( t) is the instantaneous size of a loop when it radiates at t. Accordingly, the radiation parameters can be decomposed as a , where and the normalization factor is approximately ∞ j=1 j −4/3 ≃ 3.60.Taking redshift into account, the observed frequency at today's interferometers is where t 0 is present time and the scale factor today is a(t 0 ) ≡ 1.The relic gravitational wave amplitude is summed over all normal modes From Eqs. (2.59) and (2.61), the contribution from an individual k mode can be expressed as where t f is the formation time of the string network.Heaviside theta functions ensure causality and energy conservation.t (k) i represents the time when a loop is formed, which emits gravitational wave at the time t, and is given by where the loop size can be written as l = 2 k a( t)/f .The frequency spectrum of the gravitational wave amplitude is calculated by numerically evaluating Eq. (2.65), and summing from k = 1 to a large value to ensure convergence.Evidently, the spectrum can be divided into three regions.The first region corresponds to very high frequencies starting from a cutoff value f v 0 , where the signal falls off, and the exact shape depends on the initial conditions and very early stages of the string network evolution not fully captured by the VOS model.The cutoff value f v 0 is related to the time when the Goldstone radiation becomes significant.In the intermediate radiation dominated region f eq < f < f v 0 , the spectrum gradually declines as log 3 ( l−1 /f ).In the matter dominated region f 0 < f < f eq , the spectrum behaves as f −1/3 .The frequency f eq is related to the time of matter-radiation equality, while f 0 to the emission at the present time.These characteristic frequencies are given by ∼ 10 10 Hz, (2.67) f eq ∼ 1.8 × 10 −7 Hz. (2.69) Although the gravitational wave spectrum from global cosmic strings span over a very wide frequency range, for our purposes we will be concerned mostly in the µ-Hz to kilo-Hz range, where some of the planned interferometers are sensitive.This range falls under f eq < f < f v 0 .The gravitational wave spectrum can be approximately expressed in this regime by [66] where z eq ≃ 8000 [78], and ∆ R (f ) represents the effect of varying number of relativistic degrees of freedom over time: . (2.71) We note that we have focused on the GW spectrum from decaying cosmic string loops.Another approach, the Abelian-Higgs model, that takes into account the contribution from long string network is expected to be subdominant [28].
For our numerical calculations we have set α = 0.1 as the peak value of the loop sizes at the time of their formation inspired by results from Nambu-Goto string simulations [71,72].The resulting GW spectrum is modified by up to an order of magnitude if we deviate from this choice [66].In the standard radiation-dominated cosmology, α ≲ 0.1 leads to smaller lifetime of the loops, higher string tension, and larger available string energy density to produce GWs.Furthermore, smaller α implies that loops emit GW at higher frequency, hence the amplitude for a given frequency observed today is higher in the frequency range we are interested in.For α ≫ 0.1, the loops are long-lived and resemble the scenario of Nambu-Goto strings, where the GW amplitude rises with α as Ω GW ∝ α 1/2 .On the other hand, a recent simulation of global cosmic string [68] suggests a log-normal distribution of α, in which case the GW amplitude is also enhanced by a factor of few [66].Our choice of α ∼ 0.1 can therefore be treated as a conservative choice as far as the GW amplitude is concerned.
There are several constraints on the global cosmic string formation scale ∼ v 0 .The dominant radiation mode from global strings is emission of Goldstone bosons.Assuming they remain massless, the upper limit on the total relic radiation energy density from CMB ∆N eff ≲ 0.2 [78] implies v 0 ≲ 3.5 × 10 15 GeV [66].If we assume standard cosmology, non-observation of gravitational waves at Parkes Pulsar Timing Array (PPTA) [72,79,80] gives an upper bound v 0 < 2 × 10 15 GeV.Other constraints from inflation scale and CMB anisotropy bound require v 0 ≲ O(10 15 ) GeV [81,82].Hence we consider the global lepton number symmetry violation at scales ≲ 10 15 GeV.Furthermore, we require T RH ≳ v 0 to ensure that the lepton number symmetry is restored in the early universe and symmetry breaking can take place at the scale ∼ v 0 .
We show the global cosmic string induced GW signals for v 0 = 10 14 and 10 15 GeV in Fig. 1 with red curves.The former is within the sensitivity of upcoming interferometers µAres, DECIGO and BBO, whereas the latter might be probed at LISA, AEDGE and Einstein Telescope as well.The phase transition signals for v 0 = 10 14 and 10 15 GeV remain buried under their respective cosmic string signals. 6We therefore conclude that the single majoron model can still be probed in GW interferometers through its cosmic string signal as long as the global lepton number symmetry is spontaneously broken in between 10 14 and 10 15 GeV.

GW from Majorana mass genesis in a two-majoron model
As we discussed, the GW contribution to the stochastic background from a phase transition in the single majoron model is by far below the sensitivity of planned experiments.The reason for the suppressed signal amplitude can be traced back to the fact that for a single scalar, the cubic term is strictly temperature dependent and vanishes at zero temperature.It was noticed in [19] that the signal can be strongly enhanced if an auxiliary scalar field is introduced.This would undergo its own phase transition getting its final VEV prior to the phase transition of the original scalar.In this way a bi-quadratic mixing term could be added to the tree level potential.This term generates a zero temperature barrier in the thermal effective potential able to enhance the strength of the phase transition α and, consequently, the GW spectrum. 7The nature of the auxiliary scalar field was not specified in [19].Here we propose a model with two majorons where the auxiliary scalar field is identified as a complex scalar field charged under a new global lepton number symmetry.
For definiteness, we call the complex scalars ϕ 1 , ϕ 3 , and their respective global lepton number symmetries U (1) L 1 , U (1) L 3 .The Lagrangian can be written as (I = 1, 2, 3) As before, we ignore any mixing between the SM Higgs doublet Φ with the complex scalars ϕ 1 , ϕ 3 .Here ϕ 3 couples only to the RH neutrino N 3 , whereas ϕ 1 couples to both N 1 and N 2 . 8This can be ensured by giving nonzero U (1) L 1 charges to N 1 and N 2 and half of their complementary charge to ϕ 1 , whereas N 3 and ϕ 3 have similar complementary charges under U (1) L 3 only.Furthermore, we have chosen a basis where ϕ 1 and ϕ 3 only couple to the diagonal elements of the RH neutrino mass matrix.
Analogously to the single majoron model, we write the complex fields as ϕ 1 = φ 1 e iθ 1 / √ 2 and ϕ 3 = φ 3 e iθ 3 / √ 2 and assume that the vacuum expectation values are along the real axis, After spontaneous breaking of both U (1) symmetries, J 1 = v 1 δθ 1 and J 3 = v 3 δθ 3 are identified as two majorons.We further assume the hierarchy v 3 ≫ v 1 , so that the RH neutrino mass spectrum is hierarchical The U (1) L 1 × U (1) L 3 symmetry allows the usual quadratic and quartic terms for both ϕ 1 and ϕ 2 .It also allows a quartic mixing between the two scalars, so that the tree level potential can now be written as At sufficiently high temperatures, both symmetries are restored.At temperatures T ∼ v 3 , spontaneous breaking of U (1) L 3 generates the massless majoron field J 3 .From this first phase transition we can expect a negligible contribution to the GW spectrum at observable frequencies, as we have seen in the previous section.After the ϕ 3 phase transition has completed, and ϕ 3 has settled down to its VEV v 3 , the ϕ 1 phase transition starts.Interestingly, as we are going to show, the nonzero mixing of ϕ 1 with ϕ 3 implies that the tree-level zerotemperature effective potential of ϕ 1 now gains a cubic term, which can make the phase transition of ϕ 1 strong enough to produce an observable GW spectrum.
Let us then now focus on the phase transition of ϕ 1 at a lower scale.Writing the potential Eq. (3.2) in terms of the real fields φ 1 , φ 3 , the minimization conditions yield Because of the mixing term in the potential, the scalar mass matrix has non-vanishing off-diagonal terms.Since φ 3 has already completed the phase transition, we can write φ 3 = v 3 + δφ 3 .On the other hand, we have to use the unshifted field φ 1 , since we want to describe its phase transition.Following [83], the mass matrix can be diagonalized by rotating the basis vectors, so that φ 1 and φ 3 can be expressed in terms of the new mass eigenstates φ1 and δ φ3 where the rotation angle can be determined, assuming In order to see the impact of the mixing term on the phase transition of φ 1 , we expand φ 3 as given by Eq. (3.6) in the potential in Eq. (3.2).The reason for this is that, since the phase transition of φ 3 has completed by the time φ 1 undergoes a phase transition, we can expand φ 3 around its VEV v 3 .Notice that from Eq. (3.6), one can see that the phase transition of φ 1 will induce a small shift of the φ 3 VEV given by φ1 sin θ.However, since θ ∝ v 1 /v 3 is tiny, this has no effect on the φ 1 phase transition.Although the mass eigenstates above are for the tree-level, zero-temperature, shifted fields, we will treat the aforementioned expansions simply as change of basis.In this basis, expanding the quartic mixing term in Eq. (3.2) in terms of the mass eigenstates yields a cubic term for φ1 , The dots denote the presence of additional terms ∝ (v 1 /v 3 ) 2 that can be neglected.Furthermore, to leading order, since the mixing angle θ is very small, the mass eigenstate φ1 almost coincides with φ 1 . 9As anticipated, the net effect is that a non-vanishing zero temperature cubic term appears in the thermal effective potential of φ 1 that can now be written as where C = ζ 2 v 1 /(4λ 3 ).Comparing Eq. (3.10) to Eq. (2.19), the expressions for M T , A and λ T are obtained from Eqs. (2.20)-(2.23)with the replacement λ → λ 1 , v 0 → v 1 .Since the heaviest RH neutrino mass M 3 ∼ v 3 ≫ v 1 , we can assume that the N 3 's have fully decayed at the onset of the ϕ 1 phase transition.On the other hand, we can assume that both lighter RH neutrinos are fully thermalised and, therefore, take N = 2. 10  The cubic term at zero temperature helps to strengthen the phase transition of ϕ 1 .To illustrate this, we perform a random scan over the model parameters (λ 1 , v 1 , C) 11 in the 9 There are other subleading terms that are suppressed by the small ratio v1/v3 and can be neglected.
For a full derivation of all terms one can start from Eq. (3.2) and rewrite it in terms of φ1 and φ3 as One can then rewrite φ1 and φ3 in terms of φ1 and δ φ3 using Eqs.(3.5) and (3.6).It is easy to see that, neglecting O(θ) and O(δ φ3) subleading terms, one obtains the thermal effective potential Eq. (3.10) describing the dynamics of φ1. 10 On the other hand, notice that in the Coleman-Weinberg potential in Eq. (2.8) one still has three RH neutrinos.This mismatch would produce a logarithmic term in the effective thermal potential, as pointed out in [19].However, as it has been shown there and we verified, neglecting this term is a very good approximation in the calculation of the GW spectrum. 11Strictly speaking, the model parameters are the coefficients that appear in the potential Eq. (3.2), namely, µ1, λ1, µ2, λ2 and ζ. v1, v3 and C can be expressed in terms of these parameters.We choose λ1, v1, λ3, v3 and C as free parameters.To get a better understanding of how T ⋆ , α and β/H ⋆ depend on the model parameters, we look at a two-dimensional slice of the parameter space in terms of {λ 1 , v 1 } setting M = 0.15v 1 and C = 0.002v 1 .The results are shown in Fig. 3.We find that T ⋆ ∼ O(v 1 ) and is nearly independent of λ 1 .On the other hand α is essentially determined by λ 1 , and peaks near λ 1 ∼ O(10 −4 ).Finally, β/H ⋆ depends on both λ 1 and v 1 .
We now look at the gravitational wave spectrum for the three benchmark points listed in Table 2. 12 The benchmark points have been chosen to maximize the GW amplitude from first order phase transition in their respective peak frequencies.The resulting signals are shown in Fig. 4, along with the GW spectrum from the global cosmic strings for v 3 = 10 15 , 5 × 10 14 , 2 × 10 14 and 10 14 GeV.As expected, the dominant contribution is from sound waves (thin dashed lines), whereas the MHD turbulence (thin dot-dashed lines) contributes to the high-frequency tail of the spectra.The peak amplitude of these signals are consistent with what one would expect from the range of α and β/H ⋆ where our calculation of GW spectrum is valid, as discussed in Appendix A.
The peak amplitude of the benchmark points A and B are sensitive to DECIGO, BBO, B.P. AEDGE, and ET, CE, respectively, while point C peaks at a higher frequency.In all cases, the peak amplitude is larger than the global cosmic string induced spectrum for v 3 ≲ 10 15 GeV.For any benchmark point and a given v 3 , the combined gravitational wave spectrum would look like a peak towering above the slightly tilted plateau. 13While the wideband nature of the global cosmic string induced signal offers detection possibility at multiple interferometers, the larger peak from first order phase transition provides better visibility.
Combining the two features, a unique gravitational wave signal emerges for models with two scalars, one breaking a global U (1) symmetry at ultraviolet scales and the other undergoing a strong first order phase transition at lower scales.14

GW from Majorana mass genesis in a three-majoron
A straightforward generalization of the model is to include three complex scalars with hierarchical VEVs, so that each scalar gives mass to one of the RH neutrinos, The Lagrangian has a U (1) L 1 × U (1) L 2 × U (1) L 3 symmetry, with each U (1) corresponding to each scalar. 15We denote the VEVs as ⟨ϕ I ⟩ ≡ v I and without loss of generality assume The tree-level scalar potential is given by Table 3. Benchmark point for gravitational wave signal from order phase transition of φ 1 , denoted by D (I = 1) and φ 2 , denoted by E (I = 2).Taken together, the parameters given in this table constitute one benchmark point for the three-majoron model.
After spontaneous breaking of the global U (1) symmetries, the three RH neutrinos get nonzero Majorana mass from the VEV of ϕ 1 , ϕ 2 , ϕ 3 .Assuming these VEVs are along the radial axis, one can identify the phase part of the complex scalars as massless majorons.
The mixing terms ζ IJ in Eq. (4.2) introduce a zero temperature cubic term to the effective potential of a scalar with smaller VEV.As before, the phase transition of ϕ 3 occurring at around the scale v 3 is not expected to generate any strong gravitational wave signal, since there is no zero temperature cubic term in its effective potential.However, the spontaneous breaking of U (1) L 3 at this scale would generate global cosmic string induced gravitational waves, which can be probed if v 3 ≳ 10 14 GeV.Suppose the phase transition of ϕ 3 is completed before the universe cools down to the scale v 2 , when ϕ 2 undergoes a phase transition.The quartic mixing of ϕ 2 with ϕ 3 now introduces a zero temperature cubic term to the thermal effective potential of ϕ 2 , resulting in a strong first order phase transition and associated gravitational wave from the sound waves.At this stage ϕ 1 does not play any role in the phase transition of ϕ 2 .Then, during the phase transition of ϕ 1 at around the scale v 1 , the other two scalars have already completed their phase transition and together they would introduce an effective zero temperature cubic term from their mixing with ϕ 1 , resulting in a strong phase transition and subsequent gravitational wave signal.
Because of the assumed hierarchy v 3 ≫ v 2 ≫ v 1 , the cubic terms for the phase transition of ϕ 2 and ϕ 1 depend predominantly on the mixing parameters ζ 23 and ζ 12 , respectively, C 2 ≃ v 2 ζ 23 /(4λ 3 ) and C 1 ≃ v 1 ζ 12 /(4λ 2 ).The mass parameters µ I can be expressed from minimization of the zero-temperature tree-level potential as As before, we will treat the symmetry breaking scales v I and the cubic terms C I as model parameters instead of µ I and ζ IJ , since the latter can be determined from the former in conjugation with the rest of the parameters λ I .Typically the percolation temperature T ⋆ is proportional to the VEV of the corresponding scalar undergoing the phase transition.From Eq. (2.39), this implies that the combined effect of the phase transition of the three scalars may yield a double peaked gravitational wave spectrum, with one peak at a lower frequency due to the phase transition of ϕ 1 , and another peak at a higher frequency due to the phase transition of ϕ 2 .Together with a global cosmic string induced gravitational wave spectrum from U (1) L 3 breaking, the combined amplitude of the gravitational wave signal may resemble twin peaks over a slightly slanted plateau, if the phase transition signals are sufficiently strong. 16In Table 3, we show a benchmark point consisting of the phase transition of ϕ 1 , denoted by D and the phase transition of ϕ 2 , denoted by E , that together with v 3 = 2×10 14 GeV generate the combined gravitational wave signal shown in Fig. 4 17 .Notice that in this case we have assumed that at each phase transition N = 1, corresponding to a situation where only the RH neutrino species N I , coupling to its associated scalar field ϕ I undergoing the phase transition, is fully thermalised, while the other two either have fully decayed or have not yet thermalised.This assumption is quite natural because of the strong hierarchy we are assuming for the v I 's, implying that a strong hierarchy of the RH neutrino mass spectrum and in turn of the equilibration temperatures (see Eq. (2.13).

Conclusion
We have investigated the gravitational wave signatures of the majoron model of neutrino mass generation and have identified two sources of gravitational waves.In the simplest 16 This is, of course, assuming TRH > v3, otherwise the signals that can be generated only above a given TRH would not appear. 17The corresponding model parameters are ζ12 = 0.0032, ζ23 = 0.0051, µ1 = 4.47 × 10 12 GeV, µ2 = 1.01 × 10 13 GeV, µ3 = 6.32 × 10 12 GeV, taking λ3 = 0.001 and ζ13 = 0.001.
single majoron model, a complex scalar couples to the RH neutrinos and generates their Majorana mass after spontaneously breaking the global lepton number symmetry.The breaking of a global symmetry creates global cosmic strings which can produce gravitational waves, with a different spectrum as compared to that from local Nambu-Goto strings.In the observable frequency window, the amplitude of this signal mildly declines as log 3 ], and still remains sensitive to upcoming GW interferometers if the symmetry is broken at a scale in between 10 14 − 10 15 GeV.However, there is a possible additional source of GWs in this model, since the complex scalar gets a nonzero vacuum expectation value and in the process might undergo a first order phase transition.If such a phase transition is sufficiently strong, it could generate a peaked GW signal which may tower over the global cosmic string signal.
For the simplest model with just one complex scalar coupling to all three RH neutrinos, we confirm the result of Ref. [19] that the phase transition signal is too feeble to be detected.However, we point out that the global cosmic string signal even in this model can be detected if the lepton number symmetry is broken at around 10 14 − 10 15 GeV.
We then considered an extended majoron model, introducing two complex scalars with hierarchical vacuum expectation values, one giving mass to the heaviest RH neutrino and the other to the remaining two lighter ones.Assuming the scalars are charged under separate lepton number symmetries and have a quartic mixing between them, we explored the global cosmic string induced GW spectrum which is generated when the heaviest RH neutrino gets a mass.We showed that, while the phase transition of the associated scalar remains weak, its mixing with the other scalar introduces a zero-temperature cubic term to the potential of the latter, and greatly enhances the GW signal from its phase transition.We have discussed examples where the combined GW spectrum of the model may have an observable bump or peak due to the phase transition signal, visible in the slanted plateau region from the cosmic string signal, where such a bump may appear anywhere over the whole range of observable frequencies.
Finally, we have discussed an interesting possibility of a double peaked spectrum which may occur over the global cosmic string plateau region, where such a spectrum may arise from an extension of the majoron model to include three complex scalars.This rather plausible model is easily implemented for a hierarchical RH neutrino mass spectrum, where each RH Neutrino gets its mass from the spontaneous breaking of its respective lepton number symmetry.Such a double peaked spectrum provides a characteristic signature of the three majoron model of neutrino mass generation.
We have also noticed how the observation of such a GW spectrum would give us a precious information on the cosmological history and in particular on the reheating temperature of the universe.We have implicitly assumed that this was higher than all vacuum expectation values and critical temperatures so that the GW spectra are produced through the entire range of corresponding frequencies.However, if the reheating temperature is below the vacuum expectation value of one of the complex scalar fields, then the phase transition would not take place and the signal would be absent.At the same time it should be mentioned that the model we have presented can be clearly combined with (minimal) leptogenesis [15] since the decays of the RH neutrinos would produce a B − L asymmetry that can then be partly converted into a baryon asymmetry.Therefore, the observation of the GW spectra in this model would also provide a strong test of leptogenesis.Moreover, as proposed in [84], a phase transition of the complex scalar field can be also associated to the production of a dark RH neutrino playing the role of dark matter [85].Future GW experiments have then the potential to shed light on neutrino mass genesis, cosmological history and origin of matter of the universe.

Figure 1 .
Figure1.The blue bands denote the contribution to the predicted GW spectrum from first order phase transition for v 0 = 10 15 GeV (solid line), 10 14 GeV (dashed line) and 1 TeV (dotted line).The red lines denote the contribution from cosmic strings for v 0 = 10 15 GeV (solid line) and 10 14 GeV (dashed line), the signal for v 0 = 1 TeV is too suppressed to show here.The (green) shadowed regions show the sensitivity curves of the indicated experiments.

Figure 2 .
Figure 2. GW parameters α and β/H ⋆ from scan over the model parameters λ 1 , v 1 , M and C. The color bar represents log 10 T ⋆ /GeV for each point.

Figure 3 .
Figure 3. Dependence of the GW parameters T ⋆ , α and β/H ⋆ on the model parameters λ 1 and v 1 , setting M = 0.15 v 1 and C = 0.002 v 1 .

Figure 4 .
Figure 4. Gravitational wave spectrum from first order phase transition of v 1 for three benchmark points shown in Table 2 and from global cosmic string formed by spontaneous breaking of the global U (1) L3 symmetry by ϕ 3 for four representative cases.The shaded region in the phase transition signals represent uncertainties in the calculation of the gravitational wave amplitude.Sensitivities and upper bounds from various upcoming and present interferometers are also shown.See main text for more details.

Figure 5 .
Figure 5. Gravitational wave spectrum from first order phase transition of ϕ 1 and ϕ 2 , with corresponding parameters shown in rows D and E , respectively, in Table 3, and from cosmic string formed by spontaneous breaking of the global U (1) L3 symmetry by ϕ 3 at v 3 = 2 × 10 14 GeV.Individual GW contributions are shown with dotted, dashed and dotdashed while the combined spectrum is shown with a solid curve.Band in GW spectrum from phase transition represent the possible O(0.1) suppression in the parameter Ω gw .

Figure 6 .
Figure 6.Contours of the peak GW amplitude from sound waves as a function of α and β/H ⋆ .
Using this expression for the Euclidean action, for a given choice of the model parameters v 0 , λ and M , one can calculate the critical temperature using Eq.(2.29).From this one can calculate the parameters α and β/H ⋆ and then finally derive the GW spectrum from Eq. (2.35).