The split majoron model confronts the NANOGrav signal and cosmological tensions

Abstract

In the light of the evidence of a gravitational wave background from the NANOGrav 15yr data set, we reconsider the split majoron model as a new physics extension of the standard model able to generate a needed contribution to solve the current tension between the data and the standard interpretation in terms of inspiraling supermassive black hole massive binaries. In the split majoron model the seesaw right-handed neutrinos acquire Majorana masses from spontaneous symmetry breaking of global $U(1)_{B-L}$ in a strong first order phase transition of a complex scalar field occurring above the electroweak scale. The final vacuum expectation value couples to a second complex scalar field undergoing a low scale phase transition occurring after neutrino decoupling. Such a coupling enhances the strength of this second low scale first order phase transition and can generate a sizeable primordial gravitational wave background contributing to the NANOGrav 15yr signal. Some amount of extra-radiation is generated after neutron-to-proton ration freeze-out but prior to nucleosynthesis. This can be either made compatible with current upper bound from primordial deuterium measurements or even be used to solve a potential deuterium problem. Moreover, the free streaming length of light neutrinos can be suppressed by their interactions with the resulting majoron background and this mildly ameliorates existing cosmological tensions. Thus cosmological observations nicely provide independent motivations for the model.

1 Introduction

The NANOGrav collaboration has found evidence for a gravitational wave (GW) background at $\sim$ nHZ frequencies in the 15-year data set [1, 2, 3, 4]. This strongly relies on the observed correlations among 67 pulsars following an expected Hellings-Downs pattern for a stochastic GW background [5]. A simple baseline model is provided by a standard interpretation in terms of inspiraling supermassive black hole binaries (SMBHBs) with a fiducial $f^{-2/3}$ characteristic strain spectrum. Such a baseline model provides a poor fit to the data and some deviation is currently favoured. In particular, the collaboration finds that models where in addition to SMBHBs one also has a contribution from new physics, provide a better fit to the NANOGrav data than the baseline model, resulting in Bayes factors between 10 and 100 [4].

First order phase transitions at low scales could potentially provide such an additional contribution. For temperatures of the phase transition in the range 1 MeV – 1 GeV, the resulting GW background may explain the entire NANOGrav signal [1, 4]. However, when a realistic model is considered, one needs also to take into account the cosmological constraints on the amount of extra radiation from big bang nucleosynthesis (BBN) and CMB anisotropies. A phase transition associated to the spontaneous breaking of a $U(1)_{L^{\prime}}$ symmetry, where a Majorana mass term is generated, has been previously discussed [64] as a potential origin for the NANOGrav signal from 12.5-year data set [2, 3]. In this case a complex scalar field gets a non-vanishing vacuum expectation value at the end of the phase transition and a right-handed neutrino, typically the lightest, coupling to it acquires a Majorana mass. The phase transition involves only a few additional degrees of freedom forming a dark sector, and some of them can decay into ordinary neutrinos potentially producing extra radiation so that cosmological constraints need to be considered. It has been shown that these can be respected if the phase transition occurs after neutrino decoupling and if the dark sector (re-)thermalises only with decoupled ordinary neutrinos. In this case the amount of extra radiation does not exceed upper bounds from big bang nucleosynthesis and CMB temperature anisotropies. However, in [64] it was concluded that the amplitude of the NANOGrav signal was too high to be explained by such a phase transition since the peak of the predicted spectrum was two orders of magnitude below the signal. This conclusion was based on the 12.5-year data, and on a way to calculate the sound wave contribution to the GW spectrum valid for values of the strength of the phase transition $\alpha\lesssim 0.1$ that is now outdated [65]. In this paper we reexamine this conclusion in the light of the 15 year data set and adopting an improved description of the sound wave contribution, applicable for larger values $\alpha\leq 0.6$ [66]. We introduce different improvements in the description of the phase transition in the dark sector coupled with neutrinos, distinguishing the different temperature and strength parameter compared to the visible sector, calculating the ultra-relativistic degrees of freedom at the phase transition occurring during electron-positron annihilations, including a suppression factor taking into account that the duration of gravitational wave production is in general shorter than the duration of the phase transition. We confirm that such a phase transition can hardly reproduce the whole signal but can be combined with the contribution from the SMBHB baseline model to improve the fit of the signal. On the other hand, we notice that the split majoron model receives independent motivations, since it can address different cosmological tensions. Not only it can ameliorate the well known Hubble tension, and more generally it improves the fit of cosmological observations compared to the $\Lambda$CDM model, but we notice that it also provides a solution to a potential deuterium problem that is suggested by latest measurement and reanalysis of relevant nuclear rates (D(d,n)³He and D(d,p)³H).

The paper is structured as follows. In Section 2 we discuss the split majoron model. In Section 3 we discuss the cosmological constraints deriving by the presence of extra-radiation in the model. We also discuss how the model can address a potential deuterium problem. In Section 4 we review the calculation of the GW spectrum and show the results we obtain confronting the NANOGrav 15 year-data set signal. In Section 5 we draw our conclusions and discuss future developments .

2 The split majoron model

The split majoron model was sketched in [64]. It can be regarded as an extension at low energies of the multiple majoron model proposed in [67], albeit with important distinctions and phenomenological implications. Compared to the traditional majoron model [68], we have two complex scalar fields each undergoing its own first order phase transition, one at high scale, above the electroweak scale, and one at much lower scale, dictated by the possibility to address the NANOGrav signal. If we denote by $\phi$ and $\phi^{\prime}$ the two complex scalar fields, respectively, we can write the Lagrangian as ($I=1,\dots,N$ and $I^{\prime}=N+1,\dots,N+N^{\prime}$):

$$-{\cal L}_{N_{I}+N_{I^{\prime}}+\phi+\phi^{\prime}} = \overline{L_{\alpha}}\,h_{\alpha I}\,N_{I}\,\widetilde{\Phi}+{\lambda_{I}\over 2}\,\phi\,\overline{N_{I}^{c}}\,N_{I}$$

$$+\overline{L_{\alpha}}\,h_{\alpha I^{\prime}}\,N_{I^{\prime}}\,\widetilde{\Phi}+{\lambda_{I^{\prime}}\over 2}\,\phi^{\prime}\,\overline{N_{I^{\prime}}^{c}}\,N_{I^{\prime}}+V_{0}(\phi,\phi^{\prime})+{\rm h.c.}\,,$$

where $\Phi$ is the SM Higgs doublet, $\widetilde{\Phi}$ its dual and the $N_{I},N_{I^{\prime}}$ are the RH neutrinos coupling, respectively, to $\phi$ and $\phi^{\prime}$. Imposing that the Lagrangian (2) obeys a $U(1)_{\sum_{I}L_{I}}\times U(1)_{\sum_{I^{\prime}}L_{I^{\prime}}}$ symmetry, we can take as (renormalisable) tree level potential (with no $\phi-\Phi$ and $\phi^{\prime}-\Phi$ couplings)

$$V_{0}(\phi,\phi^{\prime})=-{\mu^{2}}\,|\phi|^{2}+\lambda\,|\phi|^{4}-{\mu^{\prime}}^{2}|\phi^{\prime}|^{2}+\lambda^{\prime}|\phi^{\prime}|^{4}+\zeta|\phi|^{2}\,|\phi^{\prime}|^{2}\,.$$ (2)

We will assume that $\phi$ undergoes a phase transition, breaking a $U(1)_{\sum_{I=1}^{N}L_{I}}$ global symmetry, at some scale above the electroweak scale. In the broken phase we can rewrite $\phi$ as

$$\phi={e^{i\theta}\over\sqrt{2}}\,\left(v_{0}+S+i\,J\right)\,,$$ (3)

where $v_{0}$ is the $\phi$ vacuum expectation value, $S$ is a massive boson field with mass $m_{S}=\sqrt{2\lambda}\,v_{0}$ and $J$ is a majoron, a massless Goldstone field. The vacuum expectation value of $\phi$ generates RH neutrino masses $M_{I}=\lambda_{I}\,v_{0}/\sqrt{2}$. After electroweak symmetry breaking, the vacuum expectation value of the Higgs generates Dirac neutrino masses $m_{D\alpha I}=v_{\rm ew}\,h_{\alpha I}/\sqrt{2}$ and $m_{D\alpha I^{\prime}}=v_{\rm ew}\,h_{\alpha I^{\prime}}/\sqrt{2}$, where $v_{\rm ew}=246\,{\rm GeV}$ is the standard Higgs vacuum expectation value. In the case of the RH neutrinos $N_{I}$, their Majorana masses lead, via type-I seesaw mechanism, to a light neutrino mass matrix given by the seesaw formula

$$(m_{\nu})_{\alpha\beta}=-{v_{\rm ew}^{2}\over 2}\,{h_{\alpha I}h_{\beta I}\over M_{I}}\,.$$ (4)

Notice that we are writing the neutrino Yukawa matrices in the flavour basis where both charged leptons and Majorana mass matrices are diagonal. The Yukawa couplings $h_{\alpha I^{\prime}}$ have to be taken much smaller than usual massive fermions Yukawa couplings or even vanishing, as we will point out.

Eventually, at a scale much below the electroweak scale, $\phi^{\prime}$ also undergoes a first order phase transition breaking the $U(1)_{\sum_{I^{\prime}=1}^{N^{\prime}}L_{I^{\prime}}}$ symmetry. In the broken phase we can rewrite $\phi^{\prime}$ as

$$\phi^{\prime}={e^{i\theta}\over\sqrt{2}}\,\left(v^{\prime}_{0}+S^{\prime}+i\,J^{\prime}\right)\,,$$ (5)

where $v^{\prime}_{0}$ is the $\phi^{\prime}$ vacuum expectation value, $S^{\prime}$ is a massive boson field with mass $m_{S^{\prime}}=\sqrt{2\lambda^{\prime}}\,v^{\prime}_{0}$ and $J^{\prime}$ is a (second) massless majoron. The vacuum expectation value of $\phi^{\prime}$ generates RH neutrino masses $M_{I^{\prime}}=\lambda_{I^{\prime}}\,v^{\prime}_{0}/\sqrt{2}$. In the following, for the description of the phase transition, it will also prove convenient to introduce the real scalar field $\varphi^{\prime}$, such that $\phi^{\prime}=(\varphi^{\prime}/\sqrt{2})\,e^{i\,\theta^{\prime}}$.

Let us now discuss two different cases we will consider. First, we can have a minimal case with $N=2$ and $N^{\prime}=1$. The seesaw formula generates the atmospheric and solar neutrino mass scales while the lightest neutrino would be massless. However, after the electroweak symmetry breaking and before the $\phi^{\prime}$ phase transition, the small Yukawa couplings $h_{\alpha 3}$ generate a small Dirac neutrino mass for the lightest neutrino in a way to have a hybrid case where two neutrino mass eigenstates are Majorana neutrinos and the lightest is a Dirac neutrino. Finally, at the $\phi^{\prime}$ phase transition a Majorana mass $M_{3}$ is generated and one has a second low scale seesaw mechanism (‘mini-seesaw’) giving rise to a lightest neutrino mass $m_{1}=\sum_{\alpha}|m_{D\alpha 3}|^{2}/M_{3}$.

In a second case one has $N=3$ and a generic $N^{\prime}$. In this case the Yukawa couplings $h_{\alpha I^{\prime}}$ can even vanish. The RH neutrinos $N_{I^{\prime}}$ acquire a Majorana mass at the $\phi^{\prime}$ phase transition but they do not contribute to the ordinary neutrino masses. They can be regarded as massive neutral leptons in the dark sector, with no interactions with the visible sector (including the seesaw neutrinos).

As we will better explain in Section 4, the mixing between the two complex scalar fields $\phi$ and $\phi^{\prime}$ significantly increases the strength of the $\phi^{\prime}$ phase transition. This will be crucial in enhancing the amplitude of the generated GW spectrum observable in the NANOGrav frequencies. Before delving into the details of the GW production, we discuss some cosmological constraints on our setup and its role in potentially alleviating cosmological tensions.

3 Cosmological constraints and connection to cosmological tensions

Let us now consider the impact of cosmological constraints coming from big bang nucleosynthesis and CMB anisotropies on the model from the amount of extra-radiation (also sometimes referred to as dark radiation). To this end, we first carefully calculate the evolution of the number of degrees of freedom in the model.

3.1 Evolution of the ultra-relativistic degrees of freedom in the SM and dark sector

The number of energy density ultra-relativistic degrees of freedom $g_{\rho}(T)$ is defined as usual by $\rho_{R}(T)\equiv g_{\rho}(T)(\pi^{2}/30)\,T^{4}$, where $\rho_{R}(T)$ is the energy density in radiation. In our case it receives contributions from the SM sector and from the dark sector, so that we can write $g_{\rho}(T)=g_{\rho}^{\rm SM}(T)+g_{\rho}^{\rm D}(T)$. At the $\phi$ phase transition, occurring at a phase transition temperature $T_{\star}$ above the electroweak scale, one has for the SM contribution $g_{\rho}^{\rm SM}(T_{\star})=106.75$ and for the dark sector contribution

$$g_{\rho}^{\rm D}(T_{\star})=g_{J+S}+{7\over 4}\,N\,,$$ (6)

where $g_{J+S}=2$. Notice here we are assuming that the $N$ seesaw neutrinos all thermalise at the $\phi$ phase transition. This is something that can always be realised since all their decay parameters, defined as $K_{I}\equiv(h^{\dagger}\,h)_{II}\,\bar{v}_{\rm ew}^{2}/(M_{I}\,m_{\rm eq})$ with the effective equilibrium neutrino mass $m_{\rm eq}=[16\pi^{5/2}\sqrt{g^{\star}_{\rho}}/(3\sqrt{5})]\,(v_{\rm eq}/M_{\rm P})$ and $\bar{v}_{\rm ew}=v_{\rm ew}/\sqrt{2}=174\,{\rm GeV}$, can be larger than unity in agreement with neutrino oscillation experiments. Therefore, at the high scale phase transition the dark sector is in thermal equilibrium with the SM sector thanks to the seesaw neutrino Yukawa couplings.

After the $\phi$ phase transition, all massive particles in the dark sector, $S$ plus the $N$ seesaw neutrinos will quickly decay, while the massless majoron $J$ will contribute to dark radiation. We can then track the evolution of $g_{\rho}(T)$ at temperatures below $T_{\star}$ and prior to the low scale phase transition occurring at a temperature $T_{\star}^{\prime}$ and also prior to any potential process of rethermalisation of the dark sector that we will discuss later.

In particular, we can focus on temperatures $T\ll m_{\mu}\sim 100\,{\rm MeV}$. In this case the SM contribution can be written as [69]

$$g^{\rm SM}_{\rho}(T\ll 100\,{\rm MeV})=g^{\gamma+e^{\pm}+3\nu}_{\rho}(T)=2+{7\over 8}\,\left[4\,g_{\rho}^{e}(T)+6\,r_{\nu}^{4}(T)\right]\,,$$ (7)

where the number of energy density ultra-relativistic degrees of freedom of electrons per single spin degree is given by

$$g_{\rho}^{e}(T)={120\over 7\,\pi^{4}}\,\int_{0}^{\infty}\,dx\,{x^{2}\,\sqrt{x^{2}+z^{2}}\over e^{\sqrt{x^{2}+z^{2}}}+1}\,,$$ (8)

with $z\equiv m_{e}/T$. Above the electron mass one has $g_{\rho}^{e}(T\gg m_{e})=1$, while of course $g_{\rho}^{e}(T)\rightarrow 0$ for $T/m_{e}\rightarrow 0$. The neutrino-to-photon temperature ratio $r_{\nu}(T)\equiv T_{\nu}(T)/T$ can, as usual, be calculated using entropy conservation,

$$r_{\nu}(T)=\left({2\over 11}\right)^{1\over 3}\,\left[g_{s}^{\gamma+e^{\pm}}(T)\right]^{1\over 3}\,,$$ (9)

where

$$g_{s}^{\gamma+e^{\pm}}(T)=2+{7\over 2}\,g^{e}_{s}(T)\,,$$ (10)

having defined the contribution to the number of entropy density ultra-relativistic degrees of freedom of electrons (per single spin degree of freedom) as

$$g_{s}^{e}(T)={8\over 7}{45\over 4\pi^{4}}\,\int_{0}^{\infty}\,dx\,{x^{2}\,\sqrt{x^{2}+z^{2}}+{1\over 3}\,{x^{4}\over\sqrt{x^{2}+z^{2}}}\over e^{\sqrt{x^{2}+z^{2}}}+1}\,.$$ (11)

One can again verify that $g_{s}^{e}(T\gg m_{e})=1$ and $g_{s}^{e}(T)\rightarrow 0$ for $T/m_{e}\rightarrow 0$, so that one recovers the well known result $r_{\nu}(T\ll m_{e})=(4/11)^{1/3}$. With this function one can write the SM number of entropy density ultra-relativistic degrees of freedom as

$$g^{\rm SM}_{s}(T\ll m_{\mu})=g^{\gamma+e^{\pm}+3\nu}_{s}(T)=2+{7\over 8}\,\left[4\,g_{s}^{e}(T)+6\,r^{3}_{\nu}(T)\right]\,.$$ (12)

For $T\ll m_{e}$ one recovers the well known results $g^{\rm SM}_{s}(T\ll m_{e})=43/11\simeq 3.91$ and $g^{\rm SM}_{\rho}(T\ll m_{e})\simeq 3.36$.

Let us now focus on the dark sector contribution. This is very easy to calculate since one has simply $g_{\rho}^{\rm D}(T)=g_{J}\,[r_{\rm D}(T)]^{4}$ and $g_{s}^{\rm D}(T)=g_{J}\,[r_{\rm D}(T)]^{3}$, where $g_{J}=1$ and where the dark sector-to-photon temperature ratio $r_{\rm D}(T)$ can again be calculated from entropy conservation as

$$r_{\rm D}(T)=\left[{g^{\rm SM}_{s}(T)\over g^{\rm SM}_{s}(T_{\star})}\right]^{1/3}\,.$$ (13)

For example, for $m_{\mu}\gg T\gg m_{e}$, one finds $r_{\rm D}(T)=(43/427)^{1/3}\simeq 0.465$. We can also rewrite $g_{\rho}^{\rm D}(T)$ in terms of the extra-effective number of neutrino species $\Delta N_{\nu}(T)$ defined by

$$g_{\rho}^{\rm D}(T)\equiv{7\over 4}\,\Delta N_{\nu}(T)\,[r_{\nu}(T)]^{4}\,.$$ (14)

Again, in the particular case $m_{\mu}\gg T\gg m_{e}$, one finds

$$\Delta N_{\nu}(T)={4\over 7}\,g_{J}\,[r_{\rm D}(T)]^{4}={4\over 7}\,\left({43\over 427}\right)^{4\over 3}\simeq 0.027\,.$$ (15)

Such a small amount of extra radiation is in agreement with all cosmological constraints that we summarise here:

• Primordial helium-4 abundance measurements combined with the baryon abundance extracted from cosmic microwave background (CMB) anisotropies place a constraint on $\Delta N_{\nu}^{\rm eff}(t)$ at $t=t_{\rm f}\sim 1\,{\rm s}$, the time of freeze-out of the neutron-to-proton ratio [70]:

  $$\Delta N_{\nu}(t_{\rm f})\simeq-0.1\pm 0.3\ \Rightarrow\ \Delta N_{\nu}(t_{\rm f})\lesssim 0.5\,\ \ (95\%\,{\rm C.L.})\,.$$ (16)

• From measurements of the primordial deuterium abundance at the time of nucleosynthesis, $t_{\rm nuc}\simeq 310\,{\rm s}$, corresponding to $T_{\rm nuc}\simeq 65\,{\rm keV}$ [71]:

  $$\Delta N_{\nu}(t_{\rm nuc})=-0.05\pm 0.22\ \Rightarrow\ \Delta N_{\nu}(t_{\rm nuc})\lesssim 0.4\,\ \ (95\%\,{\rm C.L.})\,.$$ (17)

• CMB temperature and polarization anisotropies constrain $\Delta N_{\nu}(t)$ at recombination, when $T\simeq T_{\rm rec}\simeq 0.26\,{\rm eV}$. The Planck collaboration finds [72]

  $$\Delta N_{\nu}(t_{\rm rec})=-0.05\pm 0.17\ \Rightarrow\ \Delta N_{\nu}(t_{\rm rec})\lesssim 0.3\,\ \ (95\%\,{\rm C.L.})\,.$$ (18)

Let us now consider the low scale phase transition, assuming first that this occurs at a temperature $T^{\prime}_{\star}$ above neutrino decoupling temperature $T_{\nu}^{\rm dec}\sim 1\,{\rm MeV}$, so that $r_{\nu}(T^{\prime}_{\star})=1$, but below 1 GeV. At such low temperatures, Yukawa couplings are ineffective to rethermalise the dark sector [64]. On the other hand, the coupling term $\zeta\,J^{2}\,|\phi^{\prime}|^{2}$ can thermalise $\phi^{\prime}$ and the $N^{\prime}$ RH neutrinos with $J$ at a common temperature $T_{\rm D}$. Therefore, at the $\phi^{\prime}$ phase transition, the dark sector will have a temperature $T^{\prime}_{{\rm D}\star}\simeq 0.465\,T^{\prime}_{\star}$ and with such a small temperature one would obtain a GW production much below the NANOGrav signal. Notice that after the phase transition the second majoron $J^{\prime}$ would give a contribution to $\Delta N_{\nu}(T)$ equal to that one from $J$, in a way that one would obtain $\Delta N_{\nu}(T)\simeq 0.05$.

One could envisage some interaction able to rethermalise the dark sector so that $r_{\rm D}(T^{\prime}_{\star})=1$. However, in this case the thermalised $J^{\prime}$ abundance would correspond to have $\Delta N_{\nu}(T)\simeq 8/14\simeq 0.6$ throughout BBN and recombination, in disagreement with the cosmological constraints we have just reviewed. For this reason we now consider, as in [64], the case when the low scale phase transition occurs well below neutrino decoupling (i.e., $T^{\prime}_{\star}\lesssim 1\,{\rm MeV}$).

In this case a rethermalisation between the dark sector and just the decoupled ordinary neutrino background can occur without violating the cosmological constraints. Prior to the $\phi^{\prime}$ phase transition one has the interactions

$$-{\cal L}_{\nu{\rm D}}={i\over 2}\,\sum_{i=2,3}\widetilde{\lambda}_{i}\,\overline{\nu_{i}}\,\gamma^{5}\,\nu_{i}\,J$$ (19)

that can thermalise the majoron $J$ with the ordinary neutrinos, and also the complex scalar field $\phi^{\prime}$ via the coupling $\zeta\,J^{2}|\phi^{\prime}|^{2}$. This interacts with the $N_{I^{\prime}}$’s that also thermalise prior to the phase transition. The lightest neutrino $\nu_{1}$ thermalises after the $\phi^{\prime}$ phase transition interacting with $J^{\prime}$ via an interaction term analogous to the one in Eq. (19). In this way ordinary neutrinos would lose part of their energy that is transferred to the dark sector, so that they reach a common temperature $T_{\nu{\rm D}}$ given by [75, 76, 77]

$$r_{\nu{\rm D}}\equiv{T_{\nu{\rm D}}\over T}=r_{\nu}(T)\,\left({N_{\nu}^{\rm SM}(T)\over N_{\nu}^{\rm SM}(T)+N^{\prime}+12/7+4\,\Delta g/7}\right)^{1\over 4}\,,$$ (21)

where $T_{\nu}(T)$, given by Eq. (9), is the standard neutrino temperature (i.e., in the absence of the dark sector) so that one simply has $r_{\nu{\rm D}}(T)=r_{\nu}(T)\,[T_{\nu{\rm D}}(T)/T_{\nu}(T)]$. Notice that $N_{\nu}^{\rm SM}(T\gg m_{e}/2)=3$ and $N_{\nu}^{\rm SM}(T\ll m_{e}/2)=3.043$ for the predicted SM value of the effective neutrino species [78].

3.2 Hubble tension

The minimal content of the dark sector is given by $J,\phi^{\prime}$ and $N^{\prime}$ RH neutrinos. However, we can also account for the possibility of the existence of $\Delta g$ extra massless degrees of freedom. For example, for $N^{\prime}=1$ and $\Delta g=0,1,2,3$, one finds respectively $T_{\nu{\rm D}}/T_{\nu}=0.815,0.784,0.76$. One can also calculate the amount of extra-radiation at temperatures much below the electron mass, obtaining

$$\Delta N_{\nu}\simeq 3.043\left[\,\left({3.043+N^{\prime}+12/7+4\Delta g/7\over 3.043+N^{\prime}+12/7+4\Delta g/7-N_{\rm h}}\right)^{1\over 3}-1\right]\,,$$ (22)

where $N_{\rm h}$ is the number of massive states that decay after the phase transition and produce the excess radiation. In our case these states are given by $S^{\prime}$ and the $N^{\prime}$ RH neutrinos so that $N_{h}=N^{\prime}+1$. For $N^{\prime}=1$ and $\Delta g=0$ one obtains $\Delta N_{\nu}\simeq 0.465$. In this case, such an amount of extra-radiation can actually even be beneficial in order to ameliorate the Hubble tension [76, 79] compared to the $\Lambda$CDM model since one has a simultaneous injection of extra radiation together with a reduction of the neutrino free streaming length due to the interactions between the ordinary neutrinos and the majorons. For this reason a low energy scale majoron model of this kind is a leading candidate to resolve the cosmological tensions within the $\Lambda$CDM model [80]. Recently a new analysis of this model has been presented in [81] where the authors find an improvement at the level of $1\sigma$ compared to the $\Lambda$CDM model. It is then interesting that this kind of model can link the NANOGrav signal, that we are going to discuss in the next section, to the cosmological tensions.

3.3 Deuterium problem

If the rethermalisation occurs at a temperature above 65 keV, one should worry about the constraint Eq. (17) from deuterium. In this case, one can reduce the amount of extra-radiation increasing the number of massless degrees of freedom in the dark sector considering $\Delta g\neq 0$. For example, for $\Delta g=1,2,3$ one obtains, respectively, $\Delta N_{\nu}=0.41,0.37,0.33$. Therefore, an increase of the degrees of freedom in the dark sector actually produces a reduction of the amount of extra-radiation making it compatible also with deuterium constraints. However, notice that it is actually interesting that the model predicts some increase of the deuterium abundance compared to standard big bang nucleosynthesis (SBBN). There is indeed a potential tension with the current measurement of primordial deuterium abundance within SBBN. The experimental value is found to be [82] ${\rm D/H}=(2.527\pm 0.030)\times 10^{-5}$. Using a calculation of D(d,n)³He and D(d,p)³H nuclear rates based on theoretical ab-initio energy dependencies the authors of [83] find, as SBBN prediction, ${\rm D/H}=(2.439\pm 0.037)\times 10^{-5}$, showing a $\sim 2\sigma$ tension with the experimental value. Since the primordial deuterium abundance scales with $\Delta N_{\nu}$ approximately as [84] $({\rm D/H})(\Delta N_{\nu})=({\rm D/H})^{\rm SBBN}(1+0.135\,\Delta N_{\nu})^{0.8}$, one finds that $\Delta N_{\nu}(t_{\rm nuc})\simeq 0.3$ would solve the tension. However, using a polynomial expansion of the S-factors of the above-mentioned nuclear rates the authors of [71] find ${\rm D/H}=(2.54\pm 0.07)\times 10^{-5}$, a predicted value that would be essentially in agreement with the experimental value and that places the upper bound on $\Delta N_{\nu}(t_{\rm nuc})$ given in Eq. (17). New and more accurate data on the nuclear rates should be able to establish which one of the two descriptions is more reliable, thus confirming or ruling out the tension [85]. In case it will be confirmed, the split majoron model would be not only a natural candidate to explain the tension but, very importantly, it would also offer a simultaneous solution to the other cosmological tensions and, as we are now going to discuss, realise an intriguing connection with the NANOGrav signal.

4 GW spectrum predictions confronting the NANOGrav signal

We first briefly review how the first order phase transition parameters relevant for the production of GW spectrum in the split majoron model are calculated and refer the interested reader to Ref. [67] for a broader discussion. The finite-temperature effective potential for the (real) scalar $\varphi^{\prime}$ can be calculated perturbatively at one-loop level and is the summation of zero temperature tree-level, one-loop Coleman-Weinberg potential and one-loop thermal potential. Using thermal expansion of the one-loop thermal potential, this can be converted in a dressed effective potential given by

$$V^{T_{\nu{\rm D}}}_{\rm eff}(\varphi^{\prime})\simeq{1\over 2}\,\widetilde{M}_{T_{\nu{\rm D}}}^{2}\,\varphi^{{}^{\prime}2}-(A\,T_{\nu{\rm D}}+C)\,\varphi^{{}^{\prime}3}+\frac{1}{4}\lambda_{T_{\nu{\rm D}}}\,\varphi^{{}^{\prime}4}\,,$$ (24)

where notice that the common neutrino-dark sector temperature $T_{\nu{\rm D}}$ replaces the photon temperature $T$. However, once the calculations are done in terms of $T_{\nu{\rm D}}$, everything can then be more conveniently expressed in terms of the standard $T$ simply using $T_{\nu{\rm D}}(T)=r_{\nu{\rm D}}(T)\,T$. Here, a zero-temperature cubic term $C=\zeta^{2}v_{0}^{\prime}/(2\lambda)$ is introduced due to the presence of the scalar $\phi$ with a high scale vacuum expectation value during the phase transition of $\phi^{\prime}$ at a lower scale. This term significantly enhances the strength of the phase transition. The other parameters in Eq. (24) are given by

$$\widetilde{M}_{T_{\nu{\rm D}}}^{2}\equiv 2\,D\,(T_{\nu{\rm D}}^{2}-\overline{T}_{\nu{\rm D}}^{2})\,,$$ (25)

where the destabilisation temperature $\overline{T}_{\nu{\rm D}}$ is given by

$$2\,D\,\overline{T}_{\nu{\rm D}}^{2}=\lambda^{\prime}\,v_{0}^{{}^{\prime}2}+{N^{\prime}\over 8\,\pi^{2}}\,{M^{{}^{\prime}4}\over v_{0}^{{}^{\prime}2}}-{3\over 8\,\pi^{2}}\lambda^{{}^{\prime}2}\,v_{0}^{{}^{\prime}2}\,,$$ (26)

and the dimensionless constant coefficients $D$ and $A$ are expressed as

$$D={{\lambda^{\prime}\over 8}+{N^{\prime}\over 24}\,{M^{{}^{\prime}2}\over v_{0}^{{}^{\prime}2}}}\,\;\;\;\;\mbox{\rm and}\;\;\;\;A={(3\,\lambda^{\prime})^{3/2}\over 12\pi}\,\,.$$ (27)

The dimensionless temperature dependent quartic coefficient $\lambda_{T_{\nu{\rm D}}}$ is given by

$$\lambda_{T_{\nu{\rm D}}}=\lambda^{\prime}-\frac{N^{\prime}\,M^{{}^{\prime}4}}{8\,\pi^{2}\,v_{0}^{{}^{\prime}4}}\,\log{a_{F}\,T_{\nu{\rm D}}^{2}\over e^{3/2}\,M^{{}^{\prime}2}}+{9\lambda^{2}\over 16\pi^{2}}\,\log{a_{B}\,T_{\nu{\rm D}}^{2}\over e^{3/2}\,m_{S}^{{}^{\prime}2}}\,.$$ (28)

The cubic term is negligible at very high temperatures and the potential is symmetric with respect to $\phi^{\prime}$. However, at lower temperatures it becomes important and a stable second minimum forms at a nonzero $\phi^{\prime}$. At the critical temperature the two minima are degenerate and below the critical temperature bubbles can nucleate from the false vacuum to the true vacuum with nonzero probability. We refer to $T^{\prime}_{{\nu{\rm D}}\star}$ as the characteristic phase transition temperature and identify it with the percolation temperature, when $1/e$ fraction of space is still in the false vacuum. It is related to the corresponding temperature of the SM sector simply by $T^{\prime}_{{\nu{\rm D}}\star}=T^{\prime}_{\star}\,r_{\nu{\rm D}}(T^{\prime}_{\star})$ (we always use $T$ as independent variable and from this we calculate $T_{\nu{\rm D}}$). Two other parameters relevant for the calculation of the GW spectrum from first order phase transitions are $\alpha$ and $\beta/H_{\star}$, where the first denotes the strength of the phase transition and the latter describes the inverse of the duration of the phase transition. These parameters are defined as

$${\beta\over H_{\star}}\simeq T^{\prime}_{\star}\left.{d(S_{3}/T_{\nu{\rm D}})\over dT}\right|_{T^{\prime}_{\star}}\,,\qquad\text{and}\qquad\alpha\equiv{\varepsilon(T^{\prime}_{{\nu{\rm D}}{\star}})\over\rho(T^{\prime}_{\star})}\,,$$ (29)

where $S_{3}$ is the spatial Euclidean action, $\varepsilon(T^{\prime}_{{\nu{\rm D}}\star})$ is the latent heat released during the phase transition and $\rho(T^{\prime}_{\star})$ is the total energy density of the plasma, including both SM and dark sector degrees of freedom. An approximate analytical estimate for calculating $S_{3}/T_{\nu{\rm D}}$, and from this $T^{\prime}_{{\nu{\rm D}}\star}$, in terms of the model parameters can be found in Ref. [67]. In calculating $\alpha$ for phase transition at low temperatures, one must be careful about various cosmological constraints, as outlined in section 3.

We now proceed to calculate the GW spectrum of the model relevant for nanoHZ frequencies. Assuming that first order phase transition occurs in the detonation regime where bubble wall velocity $v_{w}>c_{s}=1/\sqrt{3}$, the dominant contribution to the GW spectrum mainly comes from sound waves in the plasma. Numerical simulations confirm for $\alpha\lesssim 0.1$ [65, 88, 89] the analytical result from the sound shell model [90]:

$$h^{2}\Omega_{\rm sw0}(f)=3\,h^{2}\,r_{\rm gw}(t_{\star},t_{0})\,\widetilde{\Omega}_{\rm gw}\,H_{\star}\,R_{\star}\,\left[\frac{\kappa(\alpha_{\nu{\rm D}})\,\alpha}{1+\alpha}\right]^{2}\,\widetilde{S}_{\rm sw}(f)\,\Upsilon(\alpha,\alpha_{\nu{\rm D}},\beta/H_{\star})\,,$$ (30)

where $R_{\star}$ is the mean bubble separation and a standard relation is $R_{\star}=(8\pi)^{1/3}\,v_{\rm w}/\beta$. Notice that the parameter $\alpha_{\nu{\rm D}}\equiv\varepsilon(T^{\prime}_{\star})/\rho_{\nu{\rm D}}(T^{\prime}_{\star})$ replaces $\alpha$ inside $\kappa$ and $\Upsilon$ [73, 74], where we simple defined $\rho_{\nu{\rm D}}\equiv\rho_{\nu}+\rho_{\rm D}$. We adopt Jouguet detonation solution for which the efficiency factor is given by  [91, 92]

$$\kappa(\alpha_{\nu{\rm D}})\simeq{\alpha_{\nu{\rm D}}\over 0.73+0.083\sqrt{\alpha_{\nu{\rm D}}}+\alpha_{\nu{\rm D}}}\,,$$ (31)

and the bubble wall velocity $v_{\rm w}(\alpha_{\rm D})=v_{\rm J}(\alpha_{\rm D})$, where

$$v_{\rm J}(\alpha_{\nu{\rm D}})\equiv\frac{\sqrt{1/3}+\sqrt{\alpha_{\nu{\rm D}}^{2}+2\alpha_{\nu{\rm D}}/3}}{1+\alpha_{\nu{\rm D}}}\,.$$ (32)

The suppression factor $\Upsilon(\alpha,\alpha_{\nu{\rm D}},\beta/H_{\star})\leq 1$ takes into account the finite lifetime of the soundwaves and is given by [93, 94]:

$$\Upsilon(\alpha,\alpha_{\nu{\rm D}},\beta/H_{\star})=1-{1\over\sqrt{1+2\,H_{\star}\tau_{\rm sw}}}\,,$$ (33)

where we can write

$$H_{\star}\tau_{\rm sw}=(8\,\pi)^{1\over 3}{v_{\rm w}\over\beta/H_{\star}}\left[{1+\alpha\over\kappa(\alpha_{\nu{\rm D}})\,\alpha}\right]^{1/2}\,.$$ (34)

Only in the ideal asymptotic limit $H_{\star}\tau_{\rm sw}\rightarrow\infty$ one has $\Upsilon=1$. The prefactor $\widetilde{\Omega}_{\rm gw}$ is a dimensionless number given by an integral over all wave numbers $k$ [95]

$$\widetilde{\Omega}_{\rm gw}=\int{dk\over k}\,{(k\,L_{\rm f})^{3}\over 2\pi^{2}}\,\widetilde{P}_{\rm GW}(k\,L_{\rm f})\ ,$$ (35)

where $L_{\rm f}$ is a characteristic length scale in the velocity field, $\widetilde{P}_{\rm GW}$ is the GW power spectrum and it is found [95, 88]

$$\widetilde{\Omega}_{\rm gw}={(0.8\pm 0.1)\over 2\pi^{3}}\sim 10^{-2}\,.$$ (36)

The redshift factor $r_{\rm gw}(t_{\star},t_{0})$, evolving $\Omega_{\rm gw\star}\equiv\rho_{\rm gw\star}/\rho_{{\rm c}\star}$ to $\Omega_{\rm gw0}\equiv\rho_{\rm gw0}/\rho_{\rm c0}$, is given by [96]

$$r_{\rm gw}(t_{\star},t_{0})=\left({a_{\star}\over a_{0}}\right)^{4}\,\left({H_{\star}\over H_{0}}\right)^{2}=\left({g_{s0}\over g_{s\star}}\right)^{4\over 3}\,{g_{\rho\star}\over g_{\gamma}}\,\Omega_{\gamma 0}\,,$$ (37)

where $\Omega_{\gamma 0}\equiv\rho_{\gamma 0}/\rho_{{\rm c}0}$. The normalised spectral shape function $\widetilde{S}_{\rm sw}(f)$ is given by $\widetilde{S}_{\rm sw}(f)\simeq 0.687\,S_{\rm sw}(f)$ with

$$S_{\rm sw}(f)=\left(\frac{f}{f_{\rm sw}}\right)^{3}\left[\frac{7}{4+3({f/f_{\rm sw}})^{2}}\right]^{7/2}\,,$$ (38)

where $f_{{\rm sw}}$ is the peak frequency at the present time. This is simply obtained redshifting the peak frequency at the time of the phase transition: $f_{{\rm sw}}=f_{{\rm sw}\star}a_{\star}$. The peak frequency at the phase transition is given, in terms of $v_{\rm w}$ and $\beta/H_{\star}$, by [88]

$$f_{{\rm sw}\star}=\kappa\,{\beta/H_{\star}\over v_{\rm w}}\,H_{\star}\,,$$ (39)

with $\kappa\simeq 0.54$. From entropy conservation one can write $a_{\star}=T_{0}\,g_{\rm s}^{1/3}(T\ll m_{e})/(T_{\star}\,g^{1/3}_{\rm s\star})$ and from the Friedmann equation $H_{\star}\simeq 1.66\,T_{\star}^{2}\,g_{\rho\star}^{1/2}/M_{\rm Pl}$, where $T_{0}\simeq 2.35\times 10^{-4}\,{\rm eV}$ is the CMB temperature and $g_{s}(T\ll m_{e})=g_{s}^{\rm SM}(T\ll m_{e})\simeq 3.91$. In this way one obtains for the peak frequency at the present time

$$f_{{\rm sw}}\simeq 1.66\,\kappa\,T_{0}\,g_{{\rm S}0}^{1/3}\,{\beta/H_{\star}\over v_{\rm w}}\,{T_{\star}\,g_{\rho\star}^{1/2}\over M_{\rm Pl}\,g^{1/3}_{\rm S\star}}\simeq 4.1\,\mu{\rm Hz}\,\frac{1}{v_{\rm w}}\frac{\beta}{H_{\star}}\frac{T_{\star}}{\rm 100\,GeV}\frac{g_{\rho\star}^{1/2}}{g_{S\star}^{1/3}}$$ (40)

For phase transitions above the electroweak scale the peak frequency one has $g_{\rho\star}=g_{s\star}=g_{\star}$. In this way Eq. (40) specialises into

$$f_{\rm sw}=8.9\,\mu{\rm Hz}\,\frac{1}{v_{\rm w}}\frac{\beta}{H_{\star}}\frac{T_{\star}}{\rm 100\,GeV}\,\left(\frac{g_{\rho}^{\star}}{106.75}\right)^{1/6}\,.$$ (41)

We can also write a numerical expression for the redshift factor

$$r_{\rm gw}(t_{\star},t_{0})\simeq 3.5\times 10^{-5}\,\left({106.75\over g_{\star}}\right)^{1\over 3}\,\left({0.68\over h}\right)^{2}\,$$ (42)

and, finally, for the GW spectrum from sound waves

$$h^{2}\Omega_{\rm sw0}(f)=0.97\times 10^{-6}\,{\widetilde{\Omega}_{\rm gw}\over 10^{-2}}\,\frac{v_{\rm w}(\alpha)}{\beta/H_{\star}}\left[\frac{\kappa(\alpha_{\nu{\rm D}})\,\alpha}{1+\alpha}\right]^{2}\left(\frac{106.75}{g_{\star}}\right)^{1/3}\,S_{\rm sw}(f)\,\Upsilon(\alpha,\alpha_{\nu{\rm D}},\beta/H_{\star})\,.$$ (43)

Let us now turn to the case of our interest, a low scale phase transition for $T^{\prime}_{\star}\lesssim 1\,{\rm MeV}$. In this case one has $g_{\rho}(T^{\prime}_{\star})\neq g_{s}(T^{\prime}_{\star})$, specifically:

$$g_{\rho\star}^{\prime}\equiv g_{\rho}(T^{\prime}_{\star}) = g^{\gamma+e^{\pm}}_{\rho}(T^{\prime}_{\star})+g^{3\nu}_{\rho}(T^{\prime}_{\star})+g_{\rho}^{\rm D}(T^{\prime}_{\star})$$

$$= 2+{7\over 2}\,g_{\rho}^{e}(T^{\prime}_{\star})+\left({21\over 4}\,+g_{J}+g_{S^{\prime}}+g_{J^{\prime}}+{7\over 4}\,N^{\prime}+\Delta g\right)\,r_{\nu{\rm D}}^{4}(T^{\prime}_{\star})\,,$$

and

$$g_{s\star}^{\prime}\equiv g_{s}(T^{\prime}_{\star}) = g^{\gamma+e^{\pm}}_{s}(T^{\prime}_{\star})+g^{3\nu}_{s}(T^{\prime}_{\star})+g_{\rho}^{\rm D}(T^{\prime}_{\star})$$

$$= 2+{7\over 2}\,g_{s}^{e}(T^{\prime}_{\star})+\left({21\over 4}+g_{J}+g_{S^{\prime}}+g_{J^{\prime}}+{7\over 4}\,N^{\prime}+\Delta g\right)\,r_{\nu{\rm D}}^{3}(T^{\prime}_{\star})\,.$$

In the limit $T^{\prime}_{\star}\gg m_{e}/2$ one has $r_{\nu{\rm D}}\simeq g_{\rho}^{e}\simeq 1$. If, for definiteness, we consider the minimal case with $\Delta g=0$ and $N^{\prime}=1$, one has then, for $T^{\prime}_{\star}\gg m_{e}/2$, $g_{\rho\star}^{\prime}\simeq g_{s\star}^{\prime}\simeq 62/4$. We can then conveniently rewrite numerically:

$$r_{\rm gw}(t^{\prime}_{\star},t_{0})\simeq 6.6\times 10^{-5}\,\left({0.68\over h}\right)^{2}\,\left(\frac{15.5}{g^{\prime}_{s\star}}\right)^{4\over 3}\,\frac{g^{\prime}_{\rho\star}}{15.5}\,,$$ (46)

$$f_{{\rm sw}}\simeq 6.47\,{\rm nHz}\,\frac{1}{v_{\rm w}}\frac{\beta/H_{\star}}{100}\frac{T_{\star}}{\rm 1\,MeV}\frac{(g_{\rho\star}/15.5)^{1/2}}{(g_{s\star}/15.5)^{1/3}}$$ (47)

and

$$h^{2}\Omega_{\rm sw0}(f)=1.845\times 10^{-6}\,{\widetilde{\Omega}_{\rm gw}\over 10^{-2}}\,\frac{v_{\rm w}(\alpha)}{\beta/H_{\star}}\left[\frac{\kappa(\alpha_{\nu{\rm D}})\,\alpha}{1+\alpha}\right]^{2}\left(\frac{15.5}{g^{\prime}_{s\star}}\right)^{4/3}\left(\frac{g^{\prime}_{\rho\star}}{15.5}\right)\,S_{\rm sw}(f)\,\Upsilon(\alpha,\alpha_{\nu{\rm D}},\beta/H_{\star})\,.$$ (48)

For $\alpha\gtrsim 0.1$ one expects strong deviation from (30) that can be expressed in terms of a function $\xi(f;\alpha,v_{\rm w},\beta/H_{\star},\dots)$. This function is currently undetermined. Here we mention a few effects that have been studied and that contribute to $\xi(f;\alpha,v_{\rm w},\beta/H_{\star},\dots)$.

• The expression (30) neglects a contribution from turbulent motion of the dark sector plasma after the phase transition. This contribution is certainly subdominant for $\alpha\lesssim 0.1$ but it might become sizeable for $\alpha\gtrsim 0.1$, though its determination requires a better theoretical understanding [97].

• In [66] it was found in numerical simulations that for the integral on the whole range of frequencies, i.e., for the sound wave contribution to the GW energy density parameter, one has a suppression by a factor $0.1$–$1$ for values $\alpha\lesssim 0.6$ and $v_{\rm w}\simeq c_{\rm s}$ compared to the expected result one obtains integrating Eq. (30). There are currently no well established results on the spectral shape function deviation for $\alpha>0.1$ compared to the broken power law Eq. (38). In order to account for such an indetermination, we show the GW spectrum for bands corresponding to $\xi=0.1$–$1$ in our plots in Figs. 2 and 3.

  We refer the interested reader to Ref. [67] for more details about the known issues and caveats in using the above expressions for calculating the GW spectrum. The GW spectrum plots are obtained for a set of benchmark points given in Table 1 and 2 in Figs. 2 and 3, respectively.

• Recently, it has been found in [99] that, when the bubbles expand as deflagrations, the heating of the fluid in front of the phase boundary suppresses the nucleation rate increasing the mean bubble separation and enhancing the gravitational wave signal by a factor of up to order ten. This enhancement increases for increasing values of $\alpha$ and low values of $v_{\rm w}$, so that it is sizeable only in the case of deflagrations ($v_{\rm w}<c_{\rm s}$), while it is negligible in the case of detonations ($v_{\rm w}>c_{\rm s}$), the case we have considered. In any case this effect can only partially compensate the suppression effect mentioned in the previous point.

• Another possible effect leading to a strong enhancement can come from density fluctuations if $\delta T/\bar{T}\gtrsim 1/(\beta/H_{\star})$ [100]. From the reported results, the enhancement might be up to an order of magnitude. However, there are no specific calculations and at the moment such an effect should be regarded as potential.

We can conclude that, within current knowledge, Eq. (30) should be regarded as an upper bound of the GW spectrum from first order phase transitions in the dark sector, likely for values $\alpha<0.6$ from existing numerical simulations. Even for higher values of $\alpha$ there is currently no real reason to think there can be a strong enhancement, rather a suppression, except for the hope that turbulence might become dominant and produce $\xi\gg 1$. For this reason in the following we will show results using Eq. (30) as a plausible upper bound. We will comment again in the final section about the possibility to evade such an upper bound.

4.1 Results

First of all we have produced scatter plots in the plane $\beta/H_{\star}$ versus $\alpha$ over the four parameters $v^{\prime}_{0}$, $M^{\prime}$, $\lambda^{\prime}$, $C$ and for the three values $N^{\prime}=1,3,5$. The results are shown in the three panels in Fig. 1. The shadowed regions for $\alpha>0.6$ indicate that in this regime there are no firm results from numerical simulations and for this reason we do not show benchmark points for such high values of $\alpha$. We also highlight benchmark points for which we show the GW spectrum in Fig. 2 and Fig. 3 for values of the parameters showed in the two tables.

Figure 1: Scatter plot in the plane $\beta/H_{\star}$ versus $\alpha$ over the four parameters $v^{\prime}_{0}$, $M^{\prime}$, $\lambda^{\prime}$, $C$ and for the three values $N^{\prime}=1,3,5$ corresponding to the three panels. The shadowed region indicates that for $\alpha\gtrsim 0.6$ we do not have a reliable expression for the GW spectrum. In the first panel, for $N^{\prime}=1$, the diamond, lower and upper triangles indicate the three benchmark points in Table 1. The diamond in the first panel, the square in the second panel and the circle in the third panel indicate the three benchmark points in Table 2.

B.P.	$N^{\prime}$	$\lambda^{\prime}$	$v^{\prime}_{0}$/keV	$M^{\prime}$/keV	$C$/keV	$\alpha$	$\alpha_{\nu{\rm D}}$	$\kappa_{\nu{\rm D}}$	$\beta/H_{\star}$	$T_{\star}$/keV	$v_{\rm w}$	$\Upsilon$
A	1	$0.0013$	$54.85$	$16.08$	$0.96$	$0.45$	$2.06$	$0.74$	$423.93$	$276.70$	$0.96$	$0.014$
B	1	$0.001$	$71.0$	$20.0$	$0.75$	$0.52$	$2.40$	$0.74$	$424.0$	$240.58$	$0.97$	$0.013$
C	1	$0.001$	$83.0$	$23.0$	$1.70$	$0.60$	$2.62$	$0.75$	$399.73$	$515.11$	$0.97$	$0.013$
D	1	$0.001$	$144.0$	$40.0$	$3.0$	$0.59$	$2.56$	$0.75$	$393.63$	$888.35$	$0.97$	$0.013$

Table 1: Benchmark points for GW signals from first order phase transition of $\phi^{\prime}$ for $N^{\prime}=1$.

Figure 2: Left: GW spectrum at NANOGrav for $N^{\prime}=1$ and different $\alpha$. Right: Strain spectrum compared to best fit from NANOGrav 15-yr data. Benchmark points are given in Table 1.

B.P.	$N^{\prime}$	$\lambda^{\prime}$	$v^{\prime}_{0}$/keV	$M^{\prime}$/keV	$C$/keV	$\alpha$	$\alpha_{\nu{\rm D}}$	$\kappa_{\nu{\rm D}}$	$\beta/H_{\star}$	$T_{\star}$/keV	$v_{\rm w}$	$\Upsilon$
B	1	$0.001$	$71.0$	$20.0$	$0.75$	$0.52$	$2.40$	$0.74$	$424.0$	$240.58$	$0.97$	$0.013$
E	3	$0.001$	$129.0$	$29.0$	$1.20$	$0.59$	$2.09$	$0.71$	$420.93$	$262.61$	$0.96$	$0.013$
F	5	$0.0013$	$86.7$	$14.18$	$0.87$	$0.59$	$1.87$	$0.69$	$463.13$	$218.65$	$0.96$	$0.012$

Table 2: Benchmark points for phase transition of $\phi^{\prime}$ with $N^{\prime}=1,3,5$, respectively.

Figure 3: Left: GW spectrum at NANOGrav with different $N^{\prime}$. Right: Strain spectrum compared to best fit from NANOGrav 15-yr data. Benchmark points are given in Table 2.

The left panel of Fig. 2 shows four GW spectra, corresponding to the four benchmark points in Table 1, peaking in the frequency range probed by NANOGrav for $N^{\prime}=1$. The peak amplitude of the signals are comparable, while the peak frequency shifts. In Fig. 3, we show two more benchmark points E and F, for $N^{\prime}=3$ and $N^{\prime}=5$ respectively. The resulting spectra is very similar to the benchmark point B, showing that the maximum GW signal that can be achieved in this model in the NANOGrav frequencies is somewhat independent of $N^{\prime}$. In these plots we have shown the sensitivity of some interferometers (SKA [101], THEIA [102]) in the relevant frequency range with green shaded regions, and the recent NANOGrav [1, 4] and EPTA [103, 104] results with blue and orange violins. These represent the symmetrical representations of the 1D marginalized posterior probability density distributions of the GW energy density at each sampling frequency of the NANOGrav 15-yr and EPTA [104] data, respectively. We have also shown the baseline signal expected from supermassive black hole binaries (SMBHB), modeling their GW spectrum as a power-law fit following Ref. [4], where the dashed line shows the best fit and the bands correspond to $2\sigma$ deviations. Our model predicts larger amplitude than the worst case scenario of the baseline SMBHB model.

In the right panel of Figs. 2 and 3 we show the dimensionless strain $h_{c}(f)$ of the GW signals, given by

$$h_{c}(f)=\sqrt{\frac{3\mathcal{H}_{0}^{2}\ \Omega_{\rm GW}(f)}{2\pi^{2}f^{2}}},$$ (49)

where $\mathcal{H}_{0}\approx 68$ km/s/Mpc is today’s Hubble rate. We compare the results with the spectral slope $\beta={d\log h_{c}(f)}/{d\log f}$ modeling the NANOGrav strain spectrum with a simple power law of the form $h_{c}(f)=A_{\rm GW}(f/f_{\rm PTA})^{\beta}$. Expressing $\beta$ in terms of another parameter $\gamma_{\rm GW}=3-2\beta$, the $1\sigma$ fit to NANOGrav 15-yr data gives $\gamma_{\rm GW}\simeq 3.2\pm 0.6$ around $f\sim 1/(10\rm yr)$ [1]. This favorable range is shown with bands superimposed on our strain plot. We find that the spectral tilt of the phase transition signal is in tension in some range of the frequency band probed by NANOGrav.

5 Final remarks

Let us draw some final remarks on the results we obtained and how these can be further extended and improved.

• Our results are compatible with those presented in [64]. The differences can be mainly understood in terms of the different expression we are using to describe the GW spectrum from sound waves, the Eq. (30). This supersedes the expression used in [64] based on [65]. The suppression factor taking into account the shorter duration of the stage of GW production compared to the duration of the phase transition is somehow compensated by the fact that the new expression we are using is extended to higher values of $\alpha$. However, our description of bubble velocity in terms of Jouguet solutions should be clearly replaced by a more advanced one taking into account friction though we expect slight changes since the GW spectrum scales just linearly with $v_{\rm w}$. Another important difference is that compared to [65] the peak frequency is more than halved for the same values of all relevant parameters such as $T_{\star}$. This explains why we obtain higher values of $T^{\prime}_{\star}\sim 100\,{\rm keV}$ for the peak value to be in the nHz range spanned by the NANOGrav signal. Also, notice that we have improved the calculation of the GW spectrum taking properly into account the different temperature of the dark sector and calculating the efficiency in terms of $\alpha_{\nu{\rm D}}$ rather than $\alpha$.

• The peak amplitude we find is at most $h^{2}\Omega_{\rm gw}(f)\sim 10^{-11}$ at the NANOGrav frequencies and cannot reproduce the whole NANOGrav signal. However, it can help the contribution from SMBHBs to improve the fit, one of the two options for the presence of new physics suggested by the NANOGrav collaboration analysis. This is certainly sufficient to make our results interesting, also considering that the model we studied is independently motivated by the cosmological tensions. Clearly, it would be interesting to perform a statistical analysis to find the best fit parameters in our model and to quantify the statistical significance.

• The possibility to have a higher peak amplitude, corresponding to $\xi\gg 1$, cannot be excluded but from current results from numerical simulations it seems unlikely. We just notice that increasing the value of $N^{\prime}$, values of $\alpha$ higher than 0.6 are possible and since firm predictions are missing for such strong first order transitions, one cannot exclude large enhancement coming, for example, from not yet understood contribution from turbulence. A specific account of the effect of small fluctuations within our model might also offer potentially a way to obtain $\xi\gg 1$.

• The values $T^{\prime}_{\star}\sim 100\,{\rm keV}$ that we found in our solutions that enter the NANOGrav frequency range, imply $\Delta N_{\nu}\simeq 0.4$ at the time of nucleosynthesis, for $N^{\prime}=1$. This is in marginal agreement with the constraint Eq. (17) from primordial deuterium measurements but it can be fully reconciled just simply assuming extra degrees of freedom in the dark sector ($\Delta g\neq 0$). On the other hand, such an amount of dark radiation at the time of nucleosynthesis might be even beneficial to solve a potential deuterium problem. Actually if such a deviation from SBBN should be confirmed, this would provide quite a strong support to the model. Moreover, as discussed, the same amount of dark radiation at recombination can ameliorate the cosmological tensions. In this respect a dedicated analysis within our model would be certainly desirable.

• One could think to explore a scenario with $T^{\prime}_{\star}\gg 1\,{\rm MeV}$, with a massive majoron $J^{\prime}$ quickly decaying before big bang nucleosynthesis thus avoiding all cosmological constraints [73, 74]. This would require an extension of the model introducing explicit symmetry breaking terms giving mass to $J^{\prime}$. However, it has been noticed that the introduction of these terms leading to a majoron mass larger than about 1 eV would actually jeopardise the occurrence of a first order phase transition [64]. For this reason this possibility does not seem viable, since the decay rate of extremely ultra-relativistic particles is strongly suppressed.

• Finally, let us comment on the possibility to add to the tree level potential in Eq. (2) a usual Higgs portal interaction of the form $\lambda_{\Phi\phi^{\prime}}|\Phi|^{2}\,|\phi^{\prime}|^{2}$. This is not forbidden by the $U(1)_{L}$ symmetry and it is potentially interesting since one could directly consider the Higgs as the auxiliary scalar field needed to enhance the strength of the $\phi^{\prime}$ phase transition instead of $\phi$, making the model more minimal. However, it is easy to see that such a possibility is excluded by the constraints on the Higgs invisible decay width [105, 106] that place an upper bound $\lambda_{\Phi\phi^{\prime}}\lesssim 0.03$ [107, 108]. If we write the Higgs potential as $V_{0}^{\rm SM}(\Phi)=-\mu^{2}_{\rm ew}\,|\Phi|^{2}+\lambda_{\rm ew}|\Phi|^{4}$, then $v_{\rm ew}=-\mu^{2}_{\rm ew}/(2\lambda_{\rm ew})\simeq 174\,{\rm GeV}$. Expanding $\Phi$ about the electroweak VEV, one has $|\Phi|=(0,v_{\rm ew}+h/\sqrt{2})^{T}$, where $h$ is the Higgs boson field with mass $m_{h}=2\sqrt{\lambda_{\rm ew}}\,v_{\rm eq}\simeq 125\,{\rm GeV}$, so that one has $\lambda_{\rm ew}=m^{2}_{h}/(4\,v^{2}_{\rm ew})\simeq 0.13$. After electroweak symmetry breaking, the Higgs portal term would give an additional contribution to the zero temperature cubic term $C$ given by

  $$C_{\Phi}={\lambda_{\Phi\phi^{\prime}}^{2}\,v_{0}^{\prime}\over 2\lambda_{\rm ew}}\lesssim 0.35\,{\rm keV}\,{v_{0}^{\prime}\over 100\,{\rm keV}}\,.$$ (50)

  If this is compared to the values of $C$ obtained for the benchmark points, it seems that this contribution is sub-dominant. However, since it is only marginally sub-dominant, it would be certainly interesting to explore in more detail the very attractive possibility that the NANOGrav signal might have some connection with a potential contribution of majorons to the Higgs invisible decay width that should be discovered at colliders.

In conclusion the split majoron model is an appealing possibility to address part of the NANOGrav signal and the cosmological tensions, including, potentially, the deuterium problem. Moreover, it can have connections with other different phenomenologies. In any case it is certainly a clear example of how, with the evidence of a GW cosmological background from NANAOGrav, GWs have opened a new era in our quest of new physics. This should certainly alleviate the regret for the non-evidence of new physics at the LHC (at least so far).

Acknowledgments

We acknowledge financial support from the STFC Consolidated Grant ST/T000775/1. We also would like to thank David Weir for many useful discussions.