NANOGrav spectral index $\gamma=3$ from melting domain walls

We discuss cosmic domain walls described by a tension red-shifting with the expansion of the Universe. These melting domain walls emit gravitational waves with the low-frequency spectral shape $\Omega_{gw}\propto f^{2}$ corresponding to the spectral index $\gamma=3$ favoured by the recent NANOGrav 15 yrs data. We discuss a concrete high-energy physics scenario leading to such a melting domain wall network in the early Universe. This scenario involves a feebly coupled scalar field, which can serve as a promising dark matter candidate. We identify parameters of the model matching the gravitational wave characteristics observed in the NANOGrav data. The dark matter mass is pushed to the ultra-light range below $10^{-11}-10^{-12}\,\text{eV}$ which is accessible through planned observations thanks to the effects of superradiance of rotating black holes.

Introduction and Summary.Recently several pulsar timing arrays (PTAs) such as NANOGrav [1,2], EPTA (including InPTA) [3][4][5], PPTA [6,7], and CPTA [8] reported evidence for a common-spectrum signal in each dataset, with inter-pulsar angular correlations described by the Hellings-Downs (HD) curve [9], pointing to a breakthrough discovery of nHz stochastic gravitational wave (GW) background.Signals from all the PTAs are in a good agreement, and in this article we shall focus on the NANOGrav 15yrs data with the largest statistical significance.The source of nHz GWs remains unknown, but the best fit power law GW spectrum Ω gw ∝ f1.2−2.4 at 68% CL [10] reported by NANOGrav 1 disfavors simple GW-driven models of supermassive black hole binaries (SMBHBs) predicting Ω gw ∝ f 2/3 at 95% CL [10,15,16].While statistical and environmental effects may alleviate the tension [10,15,17,18], the latter motivates to investigate a possible cosmological origin of the observed GW signal.We explore this possibility in the present Letter and focus on GWs from cosmic domain walls (DWs) [19].
These melting DWs serve as a source of GWs with the spectrum distinguishable from the one provided by constant tension DWs.Most notably, the low-frequency GW spectrum from melting DWs behaves as 2 Ω gw ∝ f 2 [34], which should be compared with Ω gw ∝ f 3 [37] in the case of constant tension walls.The larger signal at small frequencies stems from the fact that the network of melting DWs efficiently emits GWs over an extended period of time: while the most energetic GWs are produced at the network formation, later emission from somewhat melted DWs feeds into the low energy tail of the spectrum.This contrasts sharply with the constant tension case, where GWs are mainly emitted at the end of wall evolution right before dissolving, e.g., due to slight breaking of Z 2 -symmetry.Note that there is no contradiction with causality considerations [38][39][40][41][42], which typically lead to Ω gw ∝ f 3 .Indeed, the standard steeper shape assumes a finite operation of the GW source, shorter than the Hubble time.In contrast, in the scenario [34,35] we discuss here, GWs are efficiently produced by the time-extended source over many Hubble time intervals.
Remarkably, the behaviour Ω gw ∝ f 2 better fits NANOGrav 15 yrs data compared to 3 Ω gw ∝ f 3 .It is conventional to parameterise the PTA GW signal as with γ being the spectral index and f yr = 1 yr −1 ≃ 32 nHz and Ω yr = 5.8 × 10 −8 .The NANOGrav best-fit value of the spectral index reads γ = 3.2 ± 0.6 [1], which is automatically recovered in the scenario with melting DWs.We demonstrate that matching to other characteristics of the NANOGrav 15 yr signal, i.e., maximal frequency f yr and energy density Ω yr , allows us to unambiguously define coupling constants of the scalar field constituting melting DWs.In particular, this scalar field should be extremely weakly coupled making it a suitable dark matter candidate, provided that its mass is confined to the ultra-light range.For such low masses, superradiance [44][45][46] plays an important role by triggering instability of rotating black holes with astrophysical masses [47].This leads to potentially observable spindown of rotating black holes and to stochastic GW background due to gravitational radiation of the bosonic condensate forming around black holes, see, e.g., Ref. [48].
Overview of melting domain walls scenario.The scenario giving rise to melting DWs involves the Z 2symmetric model of a real scalar field χ, which interacts through the portal coupling with a scalar multiplet ϕ from the primordial thermal bath: where M χ , λ χ , and g 2 are the bare mass, quartic selfinteraction constant of the field χ, and portal coupling constant, respectively [34,35].We assume that particles ϕ are relativistic at the times of interest, which fixes the variance of the field ϕ to be ⟨|ϕ| 2 ⟩ = N T 2 12 , where N counts the number of degrees of freedom associated with ϕ.The sign of the portal coupling in the scenario is fixed as g 2 > 0. This induces tachyon instability in the twofield system, which is tamed, provided that the following condition is obeyed: where λ ϕ is the quartic self-interaction constant of the multiplet ϕ (The ratio of the coupling constants β we introduced above will appear in the relations below).As a result, the effective potential characterizing the field χ exhibits spontaneous symmetry breaking leading to the non-zero temperature-dependent expectation value: In the expanding Universe, this temperature-dependence induces time-dependence, which is crucial for our further discussions.At some (lower) temperature the bare mass term becomes relevant, and the symmetry restores with ⟨χ⟩ = 0, i.e., the inverse phase transition happens.However, at the times of the cosmological evolution we are interested, the bare mass M χ is negligible; it will be included only when considering dark matter implications of the model.
Spontaneous breaking of Z 2 -symmetry leads to the formation of DWs, provided that the background field χ is set to zero, i.e., χ = 0, prior to falling into the minima of symmetry breaking potential; see Ref. [34] for details.DWs are often unwelcome in cosmology because they quickly begin to dominate the evolution of the Universe, in contradiction with observational data.This problem is absent in our case, exactly due to the time dependence of the expectation value ⟨χ⟩, as it will become clear shortly.
The Universe temperature at the time of DW formation is defined by the balance of the Hubble friction and the tachyonic thermal mass; it is estimated as where g * (T ) counts the number of relativistic degrees of freedom at the temperature T , and M P l ≈ 2.44•10 18 GeV is the reduced Planck mass.The constant B here takes into account the finite duration of the field χ roll to the minimum of its potential; B ≃ 1 for the infinitely fast roll, but generically it takes values in the range 1 ≲ B ≲ 10 3 , see Ref. [34].
The DW tension (mass per unit area) is given by σ = 2 2λ χ • ⟨χ⟩ 3 /3.The energy density of DWs in the scaling regime with one (a few) DWs per "horizon" volume H −3 , where H is the Hubble rate, is estimated as ρ wall ∼ σH.Using Eq. ( 4), where we neglect the bare mass, one finds that the energy density of melting DWs redshifts as ρ wall ∝ 1/a 5 at the radiation-dominated stage, which is in contrast to the scenario with constant tension DWs yielding ρ wall ∝ 1/a 2 .Hence, the energy density of melting DWs drops faster than the energy density of radiation, and there is no DW problem in the Universe.
GWs from melting domain walls.In both scenarios with constant tension and melting DWs, most energetic GWs are emitted within a short time interval (of the order of the Hubble time): close to the moment of the DW formation t i [34] in the case of melting DWs and near the time of the network dissolution for constant tension DWs [37].In this time interval, which is crucial for defining peak properties of GWs, one can approximate the tension of melting DWs as constant, and use some of results of analysis of Ref. [37] supported by numerical simulations performed assuming constant tension DWs.In particular, the peak frequency is estimated by the Hubble rate at the time t i , i.e., H i .Taking into account the redshift, the present-day peak frequency reads [34,35] which gives upon substituting of Eq. ( 5): Similarly, a simple estimate based on dimensional analysis gives a rather accurate prediction for the fractional spectral energy density of GWs at peak [37].Including correction factors taken from simulations of Ref. [37], one can write where the coefficients ϵ gw and A account for the efficiency of GW emission and scaling, correspondingly; one has ϵ gw A 2 ≈ 0.5.In essence, this value of ϵ gw A 2 is the only input from lattice simulations with constant tension walls, while the rest is fixed on dimensional grounds.Unless simulations with melting DWs show a considerably lower value of ϵ gw A 2 , our further discussion is unaffected.Combining Eqs. ( 4), ( 5), ( 8), using definition (3), and taking into account the redshift of GW fractional energy density during matter-dominated stage, we obtain [34,35] where h 0 = 0.67 is the reduced Hubble constant [49].Note that the peak frequency and energy density of GWs are largely determined by the coupling constants of the field χ; this property will be exploited later to recover these constants using PTA observations.To discriminate between GWs emitted by melting domain walls and other sources, one should obtain the GW spectrum.For this purpose we observe that GW emission coming from the later times t > t i occurs at the characteristic frequency f ≃ H(t)a(t)/a 0 ∝ T (t), -in full analogy with Eq. ( 6); it is obvious that f < f p .Furthermore, using the same considerations, which led to Eq. ( 8), one obtains that the fractional energy density of GWs emitted at the times t behaves as Ω gw (t) ∝ T 2 (t).Assuming that these GWs sourced at the times t give the main contribution to the spectral energy density at frequency f , we obtain the low frequency part of the spectrum [34]: For more details see Appendix, where we prove that GW emission coming from the times earlier and later than ∼ t leaves the behaviour in Eq. (10) intact 4 .Note that Eq. ( 10) is in contrast to the result obtained in the case of constant tension DWs and many other sources, i.e., first-order phase transitions and (stable) cosmic strings, giving Ω gw h 2 0 ∝ f 3 .This does not imply violation of causality: indeed, according to the discussion above, lowfrequency GW emission still fulfills f ≃ H(t)a(t)/a 0 and hence follows from on-horizon dynamics of melting DWs.Concerning the high-frequency part of the spectrum, we do not study it in detail, since in our scenario it is outside of the domain probed by NANOGrav (see below).We may safely assume that it is not different from the case of constant tension walls, i.e., there is a power law decrease Ω gw ∝ 1/f at f > f p followed by the exponential suppression at frequencies corresponding to the inverse width of DWs [37].
Figure 1 demonstrates that the predicted GW signal is compatible with the NANOGrav signal for the set of consistent values of model parameters.Below, we explain the notations used and the assumed choice of model constants.GW spectral energy density associated with the NANOGrav, or more generally, with the PTA signal is conventionally parameterised by the amplitude A related to Ω yr in Eq. ( 1) by where H 0 is the Hubble constant; recall that f yr = 1 yr −1 ≃ 32 nHz.The best fit to the NANOGrav signal is provided by the values A ≃ 6.4 +4.2 −2.7 × 10 −15 and γ = 3.2 ± 0.6 [1].The latter automatically agrees well with the model prediction (10).
To achieve the best fit value of A, which corresponds to rather large GW energy density Ω yr ≃ 5.8×10 −8 , we first set f p ≃ f yr , so that Ω gw,p ≃ Ω yr .The reason for this choice will become clear a posteriori.Combining Eqs. ( 5) and (7), we obtain T i ≃ 1.3 GeV and g * (T i ) ≃ 75.Now we fix the constants entering the GW energy density (9): β = 1 , B = 1 , N = 24.Finally, using this and Eq. ( 7), we can fix the constant g, i.e., g = 10 −18 .This implies a tiny portal coupling g 2 = 10 −36 , while β = 1 translates into the self-interaction constant λ χ = 10 −72 .This demonstrates the point made in Ref. [35] that GWs from melting DWs can be used to identify otherwise inaccessible extremely weak interactions.Such tiny coupling constants are not unfamiliar in physics; they are characteristic for axion-like particles [51].However, the two setups have different symmetry properties, and thus should not be confused.Note also that we chose the constants β and B to be close to the lower bounds of allowed values, see Eq. (3), in order to achieve the observed value Ω yr .This also explains the choice f p ≃ f yr , because for f p ≫ f yr , one would need to assume too large Ω gw,p ≫ Ω yr .It is important to stress that one can accommodate larger values of parameters β and B by a moderate increase of the number of degrees of freedom N .Indeed, increasing β by factor ξ requires only an increase of N by smaller factor ξ 1/2 .On the other hand, a change of parameter B by factor ζ requires a corresponding increase of N by a FIG. 1. Top-Left: Spectral shape of GWs emitted by melting domain walls is shown with the solid red line versus sensitivity curves of various current and planned PTAs and GW interferometers.The plot has been produced for g = 10 −18 , β = 1, B = 1, N = 24 and g * (Ti) = 75.The spectrum is cut off at fp ≃ fyr (vertical red dashed line).The shaded regions show the sensitivity of SKA [52], GAIA and THEIA [53], µARES [54], LISA [55], DECIGO [56], and BBO [57].Top-Right: 68% and 95% CL regions for the amplitude A and spectral index γ of a power-law fit to the observed GW signal (blue and green contours, respectively).At the reference frequency f = fyr, the NANOGrav best-fit values are A ≃ 6.4 +4.2 −2.7 × 10 −15 and γ = 3.2 ± 0.6.The amplitude and spectral index, A ≃ 6.3 × 10 −15 and γ = 3, predicted for the same set of model parameters as in the left plot, are marked with the red 'star.'For comparison, we show 68% and 95% CL regions for A and γ predicted for GW-driven supermassive black hole binary populations with circular orbits (purple and orange contours, respectively), and the best-fit value γ = 13/3 with the red dashed line [16].Bottom: Zoomed-in plot of the GW spectrum against the NANOGrav 15yrs data (dashed violins) [1] within the frequency range f ∈ 10 −9 , fyr Hz. much smaller factor ζ 1/4 .Let us comment on the properties of the field ϕ.Our interpretation of PTA signal in terms of GWs from melting DWs bounds the mass m ϕ of the field ϕ as m ϕ ≪ 1 GeV, which is enforced by the requirement that ϕ is relativistic at the times of DWs formation.As a result, one runs the risk of spoiling a well established picture of BBN.There are two ways of avoiding this.One is to assume that the particles ϕ decoupled from primor-dial plasma at very early times, and thus contribute insignificantly to the effective number of neutrino species N ef f .In that case, however, the effective temperature T ϕ describing the system of particles ϕ is lower than the Universe temperature.This tends to decrease GW energy density according to Eq. ( 8), but the decrease can be (partially) compensated by the sharp change of degrees of freedom number g * (T ) around QCD phase transition.Another way to handle the problem is to assume that the particles ϕ have mass m ϕ in the MeV range, i.e., 1 MeV ≪ m ϕ ≪ 1 GeV.That is, the particles ϕ become non-relativistic sometime before BBN and then decay into SM species in one or another way.In that case, one can also consider the scenario with the effective temperature T ϕ higher than the Universe temperature T .
Implications for dark matter.The field χ being very feebly coupled to the primordial thermal bath is a proper dark matter candidate.It should be stressed that for such a tiny portal coupling, g 2 ≃ 10 −36 , neither freeze-out nor freeze-in production mechanisms are efficient.Yet it is possible to generate the right dark matter abundance even with this tiny coupling constant.We briefly comment on two production mechanisms below and identify the mass M χ as a function of GW parameters assuming that the field χ constitutes all dark matter.i) Dark matter production via the direct phase transition [34].Oscillations of the field χ around the minima of its potential naturally feed into dark matter.These oscillations start at the times t ≃ t i , when the DW network is created, and continue till present unless the particles χ are unstable.In that case, the observed dark matter abundance is achieved for extremely small M χ : ) Dark matter via inverse phase transition [34], cf.Refs.[58,59].Dark matter production occurs also in the case, when there is an efficient decay channel for the aforementioned oscillations, and the field χ settles to the minimum of its potential.Yet coherent oscillations are produced at the inverse phase transition because symmetry restoration is a non-adiabatic process.In that case, one has where T sym ≃ m ϕ is the Universe temperature at the inverse phase transition.
We observe that in both cases GW parameters favoured by NANOGrav data imply ultra-light dark matter masses M χ .Notably, with these values of M χ , our scenario predicts superradiance instability of rotating black holes with astrophysical masses [47,48].This suggests a complementary way of testing the model, in particular, the future LISA observations will probe the masses of dark matter particles corresponding to the direct phase transition, while the LIGO data may be used to test the masses involved in the inverse phase transition [60,61].
Discussions.We have shown that the properties of GWs emitted by the network of melting DWs are consistent with the signal detected at PTAs.Keeping in mind that melting DWs do not overclose the Universe and the constituent field χ serves as a suitable dark matter candidate, this makes them interesting objects deserving further investigation.Perhaps the most important prospect for future studies is the numerical study of melting DWs evolution and eventually more precise determination of GW parameters, i.e., peak frequency, energy density, and the spectral shape including the high-frequency range.In particular, the details of the formation of melting walls and settling them into the scaling regime are yet to be better understood.The precise characteristics of those are crucial given that the most energetic GWs signals are coming from the earliest stages of the wall network evolution.With the current estimates of GW parameters, the NANOGrav signal is fitted in a very narrow range of model constants.Therefore, with more detailed information on the signal/improved predictions of GW properties, one will have a chance to rule out the proposed interpretation of the GW signal or establish it on firmer grounds.On a more theoretical side, it is interesting to embed the field ϕ into a realistic particle physics scenario.While in the present work we assumed that ϕ is in equilibrium with primordial plasma, it is worth investigating situations, where ϕ decouples from plasma prior to DWs formation or has never reached thermal equilibrium.
Acknowledgments.EB acknowledges support of ANR grant StronG (ANR-22-CE31-0015-01).DG acknowledges support of the scientific program of the National Center for Physics and Mathematics, section 5 "Particle Physics and Cosmology", stage 2023-2025.SR acknowledges the European Structural and Investment Funds and the Czech Ministry of Education, Youth and Sports (Project CoGraDS -CZ.02.1.01/0.0/0.0/15003/0000437).RS acknowledges the project MSCA-IF IV FZU -CZ.02.2.69/0.0/0.0/20079/0017754, European Structural and Investment Fund, and the Czech Ministry of Education, Youth and Sports.AV was supported by the Czech Science Foundation (GA ČR), project 20-28525S and is thankful to Enrico Barausse for discussions.

APPENDIX
In the main text, to derive the low frequency behaviour of GW spectrum, i.e., Ω gw (f ) ∝ f 2 , we accounted only for emission at the times t defined from f ∼ H(t)a(t)/a 0 .Let us ensure that the contribution due to GWs emitted at earlier and later times than t, does not affect this behaviour.For this purpose, we split the overall time range of GW emission in small intervals ∆t k .Then, the contributions of GWs at frequency f emitted in these time intervals sum up to where Ωgw,p (∆t k ) ≡ Ωgw,p (t k + ∆t k ) − Ωgw,p (t k ) is the fractional energy density of GWs emitted in the interval ∆t k , and fp,k ∼ H(t k )a(t k )/a 0 is the peak frequency of these GWs; the function S encoding spectral shape of GWs emitted in the interval ∆t k will be concretized shortly.The 'tilde' notation is introduced to avoid confusion with the analogous quantities characterizing emission at the absolute peak, i.e., f p and Ω gw,p .For infinitely narrow intervals ∆t k → dt k , one replaces Ωgw,p (∆t k ) → d Ωgw,p .
Consequently, one casts Eq. ( 14) in the integral form: or equivalently where we used that Ωgw,p = Ω gw,p • fp Note that the lower integration limit f min corresponds to the moment of time, when the network of melting DWs disappears at the inverse phase transition.This will not play a profound role in what follows.
For concreteness, we use the following spectral shape where we keep the numbers p > 0 and q > 0 generic for the time being.Substituting this into Eq.( 17) and switching to the integration constant x = f / fp , we obtain For p > 2, the integral here is convergent and it is saturated at x ∼ 1, which means that the peaks with fp ∼ f give the main contribution.Being interested in the regime f min ≪ f ≪ f p , we replace in which case one can take the integral explicitly: Hence, for p > 2 and q > 0 we recover our Ω gw ∝ f 2 behaviour.
To complete the proof, recall that the exponent p describes the low-frequency tail of GWs emitted during the short time interval dt k at radiation-domination.Under these conditions, causality arguments fix p = 3 [35,36,37], which is sufficient for convergence of the integral in Eq. (20).The value of q is fixed by simulations with constant tension domain walls to be q = 1; however, we do not need to assume this, as convergence of the integral ( 22) is warranted for arbitrary q > 0. Having this said, we confirm the low frequency behaviour Ω gw ∝ f 2 .