Gravitational waves from a dark $U(1)_D$ phase transition in the light of NANOGrav 12.5 yr data

We study a possibility of a strong first-order phase transition (FOPT) taking place below the electroweak scale in the context of $U(1)_D$ gauge extension of the standard model. As pointed out recently by the NANOGrav collaboration, gravitational waves from such a phase transition with appropriate strength and nucleation temperature can explain their 12.5 yr data. We first find the parameter space of this minimal model consistent with NANOGrav findings considering only a complex singlet scalar and $U(1)_D$ vector boson. Existence of a singlet fermion charged under $U(1)_D$ can give rise to dark matter in this model, preferably of non-thermal type, while incorporating additional fields can also generate light neutrino masses through typical low scale seesaw mechanisms like radiative or inverse seesaw.


Introduction:
The NANOGrav collaboration has recently released their results for gravitational wave (GW) background produced from a first order phase transition (FOPT) in 45 pulsars from their 12.5 year data [1]. According to their analysis, the 12.5 yr data can be explained in terms of a FOPT occurring at a temperature below the electroweak (EW) scale although there exists degeneracy with similar signals generated by supermassive black hole binary (SMBHB) mergers. Last year, the same collaboration also reported evidence for a stochastic GW background with a power law spectrum having frequencies around the nano-Hz regime [2] which led to several interesting new physics explanations; For example, [3][4][5] studied cosmic string origins and [5][6][7][8][9] studied FOPT origin. The pulsar timing arrays (PTAs) like NANOGrav sensitive to GW of extremely low frequencies offer a complementary probe of GW background to future space-based interferometers like eLISA [10,11].
Inspired by the results from NANOGrav explained in terms of a FOPT characterized by the preferred ranges for strength (α * ) as well as temperature (T * ) of the phase transition as shown in [1], we propose a simple model to achieve such a low scale strong FOPT. For our purpose, we introduce a dark U (1) D gauge symmetry under which only a complex singlet scalar Φ and a vector like singlet fermion Ψ are charged while all the standard model (SM) particles are neutral. Since the SM particles are neutral under this gauge symmetry, one can evade strong bounds from experiments on the corresponding gauge coupling g D and gauge boson mass m Z D . We further impose a classical conformal invariance so that U (1) D symmetry breaking occurs only via radiative effects on scalar potential, naturally leading to a vacuum below EW scale. Then, a strong FOPT can take place at bubble nucleation temperature much below electroweak scale. For * Electronic address: dborah@iitg.ac.in † Electronic address: arnabdasgupta@protonmail.ch ‡ Electronic address: skkang@seoultech.ac.kr earlier works on FOPT within such Abelian gauge extended scenarios, please refer to [12][13][14][15][16][17][18] and references therein. While such dark phase transition of strongly first order and resulting gravitational waves have been discussed earlier as well, we study this possibility for the first time after NANOGrav collaboration analysed their 12.5 year data in the context of gravitational waves from the FOPT at a low temperature below EW scale [1]. In addition, we note that the dark U (1) D symmetry can also be motivated from tiny neutrino mass and dark matter (DM) which the SM fails to address. In this work, we examine how tiny neutrino masses can be generated through low scale seesaw mechanism like radiative or inverse seesaw, and a singlet fermion charged under U (1) D can be a good dark matter candidate while keeping the model parameters consistent with the results from NANOGrav.
The Model: As mentioned above, we consider a U (1) D extension of the SM. The newly introduced fields are a complex scalar Φ and a vector like fermion Ψ with U (1) D charges 2n 1 and n 1 , respectively. All the SM fields are neutral under this new gauge symmetry. The zerotemperature scalar potential at tree level is given by where H is the SM Higgs doublet. Note the absence of bare mass squared terms due to the conformal invariance imposed. The vacuum expectation value (VEV) of the singlet scalar, Φ = M/ √ 2, acquired via the running of the quartic coupling λ breaks the gauge symmetry leading to a massive gauge boson m Z D = 2g D M . In order to realize the electroweak vacuum, the coupling λ needs to be suppressed. So in our analysis we neglect the coupling λ . We also consider the Yukawa coupling (y) of the scalar singlet with fermion Ψ to be negligible compared to gauge coupling, g D y for simplicity. This assumption is for simplicity and also to make sure that the SM Higgs VEV does not affect the light singlet scalar mass. Furthermore, the Yukawa coupling (y) is taken to be negligible as to suppress its role in the renormalisation group evolution (RGE) of the singlet scalar quartic coupling.
The total effective potential can be schematically divided into following form: where V tree , V CW and V th denote the tree level scalar potential, the one-loop Coleman-Weinberg potential, and the thermal effective potential, respectively. In finitetemperature field theory, the effective potential, V CW and V thermal , are calculated by using the standard background field method [19,20]. We consider Landau gauge throughout this work. Issues related to gauge dependence in such conformal models have been discussed recently by the authors of [18]. Denoting the singlet scalar as Φ = (φ + M + iA)/ √ 2, the zero temperature effective potential up to one-loop can be written as [12] where t = log(φ/µ) with µ being the scale of renormalisation. G(t) is given by where we have ignored φ couplings with Ψ as well as the SM Higgs H for simplicity. In the above equation a 2 = 24. The gauge coupling g D (t) and quartic coupling λ(t) at renormalisation scale can be calculated by solving the corresponding RGE equations. In terms of α D = g 2 D /4π and α λ = λ/4π, the RGEs are where b = 8/3, a 1 = 10, and a 3 = 48. Taking the renormalisation scale µ to be M , the condition dV dφ | φ=M = 0 leads us to the relation, which makes α λ (0) determined by α D (0). Since running of the coupling can be solved analytically, the scalar potential can be given by [12] V 0 (φ, t) = πα λ (t) where and the coefficient C is determined by Eq. (7). Thermal contributions to the effective potential are given by where n Bi and n Fi denote the degrees of freedom (dof) of the bosonic and fermionic particles, respectively. In this expressions, J B and J F functions are defined as follows: On calculating V th , we include a contribution from daisy diagram to improve the perturbative expansion during the phase transition [21][22][23]. The daisy improved effective potential can be calculated by inserting thermal masses into the zero-temperature field dependent masses. The author of Ref. [22] proposed the thermal resummation prescription in which the thermal corrected field dependent masses are used for the calculation in V CW and V th (Parwani method). In comparison to this prescription, authors of Ref. [23] proposed alternative prescription for the thermal resummation (Arnold-Espinosa method). They include the effect of daisy diagram only for Matsubara zero-modes inside J B function defined above. In our work, we use the Arnold-Espinosa method.
As mentioned before, we ignore singlet scalar coupling to fermion and the SM Higgs and hence calculate the field dependent and thermal masses as well as the daisy diagram contribution for vector boson only.
As the evolution has two scales, φ and T , where T is the temperature of the universe, we consider the renormalisation scale parameter u instead of t as Note that Λ represents the typical scale of the theory. Now, the one-loop level effective potential is given as: where wherein, V B T is the thermal correction and V daisy is the daisy subtraction [21][22][23].
First order phase transition: The first order phase transitions proceed via tunnelling, and the corresponding spherical symmetric field configurations called bubbles are nucleated followed by expansion and coalescence 1 . The tunnelling rate per unit time per unit volume is given as where A(T ) ∼ T 4 and S 3 (T ) are determined by the dimensional analysis and given by the classical configurations, called bounce, respectively. At finite temperature, the O(3) symmetric bounce solution [26] is obtained by solving the following equation 1 For recent reviews of cosmological phase transitions, refer to [24,25].
The boundary conditions for the above differential equation are where φ false denotes the position of the false vacuum. Using φ governed by the above equation and boundary conditions, the bounce action can be written as The temperature at which the bubbles are nucleated is called the nucleation temperature T * . This can be calculated by comparing the tunnelling rate to the Hubble expansion rate as Here, assuming the usual radiation dominated universe, the Hubble parameter is given by H(T ) 1.66 √ g * T 2 /M Pl with g * being the dof of the radiation component. Thus, the rate comparison equation above leads to for g * ∼ 100 and T * ∼ 100 GeV while for lower temperature near MeV, it comes down to g * ∼ 10. If φ(T * )/T * > 1 is satisfied, where φ(T * ) is the singlet scalar VEV at the nucleation temperature, T = T * , the corresponding phase transition is conventionally called strong first order. By choosing benchmark values as α * = 0.68, T * = 2.25 MeV, g D = 0.32, m Z D = 12.6 MeV, we can portray the curves of the potential in terms of φ/M at critical and nucleation temperatures as shown in Fig. 1. Clearly, we see that φ = 0 becomes a false vacuum below the critical temperature T c and the existence of the barrier at T c indicates a strong first order phase transition driven by the singlet scalar sector, which triggers bubble production and subsequent production of GW.
The phase transition completes via the percolation of the growing bubbles. To see when the phase transition finishes, we need to estimate the percolation temperature T p at which significant volume of the Universe has been converted from the symmetric to the broken phase. Following [27,28], T p is obtained from the probability of finding a point still in the false vacuum given by The percolation temperature is then calculated by using I(T p ) = 0.34 [27] (implying that at least 34% of the comoving volume is occupied by the true vacuum). It is also necessarily required that the physical volume of the false vacuum should be decreased around percolation for successful completion of the phase transition. This requirement reads Confirming that this condition is satisfied at the percolation temperature T p , one can ensure that the phase transition successfully completes. For the same benchmark values as taken in Fig. 1, we have calculated the percolation temperature T p and checked that the condition eq.(24) is satisfied. The results and values of some parameters are presented in table I.
The amplitudes of GW depend upon the ratio of the amount of vacuum energy released by the phase transition to the radiation energy density of the universe, ρ rad = g * π 2 T 4 /30, given by where ∆V tot ≡ V tot (φ false , T ) − V tot (φ true , T ) is the free energy difference between the false and true vacuum. * is related to the change in the trace of the energymomentum tensor across the bubble wall [11,44]. The amplitude of GW is also dictated by the duration of the FOPT, denoted by the parameter β, defined as [10] β Here, α * and β/H(T ) are evaluated at T = T * . While S 3 can be evaluated using eq. (20), the effective potential at sufficiently low temperatures i.e T M can be safely approximated as In such a scenario the action can be approximated to be In our estimation for the gravitational wave amplitude we have used the above expressions eq.(28) and eq. (29) in calculating α, β and the percolation temperature T p . We note that NANOGrav collaboration has estimated the required FOPT parameters using thin shell approximation for bubble walls (envelope approximation) [46], semi-analytic approximation [47] as well as full lattice results. Here, we present the predictions of our model against the backdrop of their estimates in Fig. 2.
During a FOPT, there are three sources producing GWs: bubble collisions, sound wave of the plasma, and turbulence of the plasma [10,36,37,42,46,48]. These three contributions together give the resultant gravitational wave power spectrum given as [1]: In general, each contribution has its own peak frequency and each GW spectrum can be parametrised in the following way [1] h where the pre-factor R 7.69 × 10 −5 g −1/3 * takes in account the red-shift of the GW energy density, S(f /f 0 * ) parametrises the shape of the spectrum and ∆(v w ) is the normalization factor which depends on the bubble wall velocity v w . The Hubble parameter at T = T * is denoted by H * . Finally the peak frequency today, f 0 * , is related to the value of the peak frequency at the time of emission, f * , as follows: where g * denotes the number of relativistic degrees of freedom at the time of the phase transition. The values of the peak frequency at the time of emission, the normalisation factor, the spectral shape, and the exponents p and q are given in Table I of [1]. The efficiency factors namely, κ φ is discussed in [12,28] and κ sw is taken from [44,49]. On the other hand, the remaining efficiency factor κ turb is taken to be approximately 0.1 × κ sw [1]. The bubble wall velocities are given in [50][51][52][53][54].  Based on the formulae presented above and by choosing a benchmark choice of model as well as FOPT parameters shown in table I consistent with NANOGrav data at 95% CL, we calculate the individual contributions to GW energy density spectrum Ωh 2 (f ) from bubble collisions, sound wave of the plasma, and turbulence of the plasma as well as the total contribution to Ωh 2 (f ). In Fig. 3, the red, orange, cyan and black curves correspond to the individual contribution from turbulence of the plasma, sound wave of the plasma, bubble collisions, and the total contribution to Ωh 2 (f ), respectively. Due to the small value of FOPT strength parameter α * , as anticipated from earlier studies [55,56], the contribution from bubble collision is suppressed as can be seen in Fig. 3.
Neutrino mass: Dark Abelian gauge extension of SM can also be related to the origin of neutrino mass. Neutrino oscillation data suggest tiny but non-vanishing light neutrino masses with two large mixing [57]. Since nonzero neutrino mass and mixing can not be explained in SM, there have been several beyond standard model (BSM) proposals. It turns out that the simplest U (1) D extension like the one discussed above augmented with additional discrete symmetries of fields can explain the origin of light neutrino mass. Here we briefly mention two such possibilities for neutrino mass origin.
First we discuss a radiative origin of light neutrino masses, a natural origin of low scale seesaw. In addition to the singlet scalar Φ and the dark fermion Ψ in the minimal model discussed above, we need an additional scalar doublet χ and a scalar singlet η to realise a radiative seesaw. The required field content and their charges under U (1) D are shown in table II. The relevant terms of the leptonic Lagrangian are given by The relevant part of the scalar potential is The singlet scalar η, neutral under U (1) D is introduced in order to avoid terms in the Lagrangian breaking conformal invariance [58]. The U (1) D symmetry is broken by a nonzero VEV of Φ to a remnant Z 2 symmetry under which Ψ L,R , χ 1 , χ 2 are odd while all other fields are even. While light neutrino mass can be realised at one loop level with these Z 2 odd particles going inside the loop, the lightest Z 2 odd particle can be a stable DM  candidate. A possible one loop diagram for light neutrino mass is shown in Fig. 4. Since Z 2 odd particles take part in the loop, the origin of light neutrino masses is similar to the scotogenic mechanism [59]. The contribution from the diagram shown in Fig. 4 can be estimated as where (M Ψ ) k is the mass of pseudo-Dirac fermion states going inside the loop and I ν is the corresponding loop function written in terms of r χ1 = M 2 χ1 /M 2 χ2 , and r k = M 2 Ψ k /M 2 χ2 with M 2 χ1 = (m 2 χr1 + m 2 χi1 )/2 and M 2 χ2 = (m 2 χr2 + m 2 χi2 )/2. Use of r, i in subscripts denotes real and imaginary neutral parts of the corresponding complex scalar fields.
We now consider the realisation of another low scale seesaw, namely inverse seesaw with U (1) D symmetry. It turns out that a minimal U (1) D gauge symmetry is not sufficient to ensure the required structure of inverse seesaw mass matrix. To have a minimal possibility we consider an additional Z 4 discrete symmetry. The new fields and their transformations under the imposed symmetries are shown in table III. The relevant part of the Yukawa Lagrangian is Clearly, the lepton number violating term involves Φ which also breaks the U (1) D symmetry. Therefore, a low scale U (1) D naturally leads to a tiny lepton number violating term in the inverse seesaw mass matrix. After symmetry breaking, the light neutrino mass is given by Thus, in both the examples discussed here, the low scale U (1) D symmetry can play non-trivial role in light neutrino mass generation even though all the SM fields are neutral under this symmetry.
Dark matter and cosmological constraints: Evidences from astrophysics and cosmology suggest the presence of a non-baryonic form of matter giving rise to approximately 26% of the present universe's energy density [57]. The simplest possibility is to consider a vector like fermion Ψ having charge n ψ under U (1) D . Depending on the strength of gauge interactions, the relic abundance of DM can be realised either via thermal or non-thermal mechanisms. While the U (1) D gauge coupling was kept large in the analysis for FOPT and GW above, DM interactions with the SM can still be suppressed due to small kinetic mixing between U (1) D and U (1) Y . However, in the discussion on neutrino mass, we have introduced additional fields charged under both SM and U (1) D gauge symmetries. This will keep the one loop kinetic mixing between U (1) D and U (1) Y suppressed but still large enough to produce Z D in equilibrium. Thus, a light gauge boson with not too small kinetic mixing with U (1) Y can decay into SM leptons at late epochs (compared to neutrino decoupling temperature T ν dec ∼ O(MeV) increasing the effective relativistic degrees of freedom which is tightly constrained by Planck 2018 data as N eff = 2.99 +0. 34 −0.33 [60]. Such constraints can be satisfied if m Z D O(10 MeV) [61,62] which agrees with the benchmark value chosen in our FOPT and GW analysis. On the other hand, taking m Z D to much higher regime will not explain the NANOGrav data. Therefore, we keep its benchmark at minimum allowed value. Similar bound also exists for thermal DM masses in this regime which can annihilate into leptons. As shown by the authors of [63], such constraints from the big bang nucleosynthesis (BBN) as well as the cosmic microwave background (CMB) measurements can be satisfied if M DM O(1 MeV). On the other hand, constraints from CMB measurements disfavour such light sub-GeV thermal DM production in the early universe through s-channel annihilations into SM fermions [60]. Since fermion singlet DM in our model primarily annihilates via s-channel annihilations mediated by Z D only, cosmological constraints are severe for thermal DM mass around or below 10 MeV.
Due to the tight cosmological constraints on thermal DM with mass below 10 MeV as discussed above, we consider a non-thermal DM scenario, also known as the feebly interacting massive particle (FIMP) paradigm [64]. While we can not make g D very small, in order to satisfy the FOPT and GW criteria, we choose U (1) D charge of DM n ψ to be very small 2 . For DM mass above Z D , it can be produced in the early universe via annihilation of SM bath particles into DM, mediated by Z D . On the left panel of Fig. 5, we show the evolution of comoving DM number density Y for DM mass M DM = 1 GeV and its U (1) D charge n ψ = 3.9 × 10 −10 . The kinetic mixing of Z D with U (1) Y of the SM is taken to be approximately ∼ g D g /(16π 2 ), similar to one-loop mixing. Clearly, DM with negligible initial abundance freezes in and gets saturated at lower temperatures, giving rise to the required relic density. On the right panel of Fig. 5, we show the parameter space in terms of M DM − n ψ giving rise to correct DM relic while keeping U (1) D sector parameters fixed at g D = 0.32 and m Z D = 12.6 MeV. Since DM mass is varied all the way upto 1 MeV for the right panel plot of Fig.  5, which is below the Z D mass threshold, we consider both annihilation and decay contributions to DM relic. Clearly, smaller values of n ψ requires larger DM mass to satisfy the relic criteria. This is because, smaller DM coupling leads to smaller non-thermal abundance and hence larger mass is required to generate the observed relic abundance. While we skip other phenomenological signatures of such DM, such sub-GeV DM can have very interesting phenomenology in the context of latest experiments like XENON1T [66] as has been discussed by [67][68][69][70] among others. Such Dirac fermion DM, upon receiving a tiny Majorana mass contribution from singlet scalar, as discussed in the context of radiative neutrino mass above can give rise to inelastic DM [71,72] with interesting DM phenomenology [73].
Conclusion: Motivated by the recent NANOGrav collaboration's analysis of their 12.5 yr data implying a possible origin of stochastic GW spectrum from a first order phase transition below EW scale, we revisit the simplest possibility of a dark Abelian gauge extension of the SM. While the SM fields are neutral under this gauge symmetry, a complex scalar singlet with non-vanishing gauge charge can lead to the necessary symmetry breaking. We further consider a classical conformal invariance such that the symmetry breaking occurs through radiative corrections to the scalar potential, keeping the model minimal. While additional dark fermions can be introduced in order to explain the origin of dark matter, for the phase transition details we confine ourselves to only the singlet scalar -vector boson interactions, ignoring other scalar portal or Yukawa interactions for simplicity. We perform a numerical scan to show how a light gauge boson Z D in sub-GeV scale can explain the FOPT parameters given by [1] in order explain their data. We have also commented on the possibility of connecting such U (1) D models to light neutrino mass and dark matter in a common setup. Due to tight cosmological constraints on such light vector boson Z D as well as DM whose interactions with the SM sector are mediated by Z D via kinetic mixing, we consider a non-thermal DM scenario. By choosing U (1) D sector parameters in a way which satisfies NANOGrav data, we perform a numerical scan over DM mass in sub-GeV range and its U (1) D charge which can give rise to the correct non-thermal DM relic. Due to complementary nature of observable signatures of such minimal models, specially in the context of GW from FOPT as well as typical sub-GeV dark matter signatures, near future experiments should be able to do more scrutiny of such predictive scenarios. Additionally, while PTAs like NANOGrav offer a complementary GW window to proposed space-based interferometers, more data are need to confirm whether this is a clear detection of GW and whether it is due to FOPT or astrophysical sources like SMBHB mergers (Possible ways of distinguishing cosmological backgrounds from astrophysical foregrounds have been discussed recently in [74]).