Nano-Hertz gravitational waves from collapsing domain walls associated with freeze-in dark matter in light of pulsar timing array observations

Evidence for a stochastic gravitational wave background in the nHz frequency band is recently reported by four pulsar timing array collaborations NANOGrav, EPTA, CPTA, and PPTA. It can be interpreted by gravitational waves from collapsing domain walls in the early universe. We assume such domain walls arising from the spontaneous breaking of a $Z_2$ symmetry in a scalar field theory, where a tiny $Z_2$-violating potential is required to make domain walls unstable. We propose that this $Z_2$-violating potential is radiatively induced by a feeble Yukawa coupling between the scalar field and a fermion field, which is also responsible for dark matter production via the freeze-in mechanism. Combining the pulsar timing array data and the observed dark matter relic density, we find that the model parameters can be narrowed down to small ranges.

Domain walls (DWs) are two-dimensional topological defects which could be formed when a discrete symmetry of the scalar potential is spontaneously broken in the early universe [97].They are boundaries separating spatial regions with different degenerate vacua.Stable DWs are thought to be a cosmological problem [98].As the universe expands, the DW energy density decreases slower than radiation and matter, and would soon dominate the total energy density.Moreover, large-scale density fluctuations induced by DWs could easily exceed those observed in the cosmic microwave background.
Nonetheless, it is allowed if DWs collapse at a very early epoch [99][100][101].Such unstable DWs can be realized if the discrete symmetry is explicitly broken by a small potential term that gives an energy bias among the minima of the potential.The bias induces a volume pressure force acting on the DWs that leads to their collapse.Collapsing DWs significantly produce GWs [102][103][104][105], which would form a stochastic background remaining to the present time, and it could be the one probed by recent PTA experiments.
In this work, we consider a scalar field S with a spontaneously broken Z 2 -symmetric potential to be the origin of DWs.These DWs can be described by the kink solution of the equation of motion.Since the DWs should collapse before they overclose the universe, a tiny but nonzero Z 2 -violating potential needs to be added.We will present a possibility that the Z 2 -violating potential is radiatively originated from a Yukawa interaction between S and a fermionic field χ.The dominant effect comes from one-loop tadpole diagrams for S.
Our further analysis shows that the Yukawa coupling should be feeble for reproducing the observed GW data, and this is also required by the freeze-in mechanism of dark matter (DM) production.Therefore, it is possible that the fermion χ, acting as a feebly interacting massive particle (FIMP) [106,107], constitutes the DM relic of the universe.There are a lot of recent studies on exploring such FIMPs [108][109][110][111][112][113][114][115][116][117][118][119].We will explore the interplay between the PTA observations of the SGWB and the freeze-in DM 1 .
The remainder of the paper is outlined as follows.In Sec.II, we discuss unstable DWs from the spontaneous breaking of an approximate Z 2 symmetry and the resulting GWs.In Sec.III, the freeze-in DM production and the induced Z 2 -violating potential are studied.In Sec.IV, we investigate the parameter ranges simultaneously fulfilling the recent PTA GW data and the observed DM relic density.Section V gives a summary.

II. DOMAIN WALLS AND GRAVITATIONAL WAVES
In this work, we consider the following Lagrangian for scalar fields, where H is the SM Higgs field and S is a real scalar field that is a SM gauge singlet.The zero-temperature potential V 0 (H, S) = V Z 2 + V vio consists of the Z 2 -conserving terms which respect a Z 2 symmetry S → −S, and the Z 2 -violating terms If one removes the Z 2 -violating terms, then the Lagrangian has the Z 2 symmetry, which would be spontaneously broken for a negative mass parameter −µ 2 S at low temperatures.In the phase where both the electroweak and Z 2 symmetries are broken, H and S develop nonvanishing vacuum expectation values (VEVs) ⟨H⟩ = (0, v/ √ 2) T and ⟨S⟩ = ±v s with v ≈ 246 GeV and v s > 0. We assume a hierarchy of v s ≫ v, implying that the Z 2 symmetry is spontaneously broken at a scale much higher than the electroweak scale.Furthermore, we assume the Higgs mass parameter µ 2 H > 0 and the portal coupling λ HS < 0, and the reason will be explained below.Thus, when the S field is integrated out below the scale v s , the effective mass parameter for H becomes µ 2 H + λ HS v 2 s /2 < 0, leading to the spontaneous breaking of the electroweak symmetry.That is to say, the electroweak symmetry breaking is essentially induced by the large VEV of S.
At high temperatures, the electroweak and Z 2 symmetries would be restored due to thermal corrections to the scalar potential.In the high-temperature limit, the effective potential becomes where T is the temperature, and δm 2 H (T ) and δm 2 S (T ) are thermal corrections to the masses of H and S given by where g ′ and g are the U(1) Y and SU(2) L gauge couplings, and y t is the Yukawa coupling of the top quark.
At a sufficiently high temperature, because of the positive contributions from the thermal masses, both the electroweak and Z 2 symmetries are restored.As the universe expands and cools down, these symmetries become broken at some critical temperatures.Since these phase transitions happen at an era after reheating, DWs could be produced after the spontaneous breaking of the Z 2 symmetry [97].
A DW corresponds to a kink solution of the equation of motion for the scalar field S given by [98] where the direction perpendicular to the DW is assumed to lie along the z-axis.Thus, S(z) approaches the VEVs ±v s for z → ±∞.The DW locates at z = 0 with a thickness δ ≈ ( λ S /2v s ) −1 , separating two domains with S(z) > 0 and S(z) < 0. By integrating the energy density along the z-direction, the surface energy density of the DW, i.e., its tension, is given by Inside each domain with S ∼ S(±∞) ≈ ±v s , we can parametrize H and S as where h(x) and s(x) are quantum fields describing the fluctuations above the vacuum.For v s ≫ v, and the masses squared of the scalar bosons h and s are approximately given by In order to ensure that the two vacua are local minima, the quartic couplings in the Z 2 -conserving potential should satisfy Further considering the condition to obtain a negative effective mass parameter for the H field driven by the VEV of S, the viable range of λ HS expressed by the h boson mass and the VEVs is Once DWs are created, their tension σ acts to straighten them against the friction from the interaction with the cosmic plasma.If the friction effect is important, the SGWB spectrum induced by DWs could be significantly different from the case without friction [122].In this work, however, the interactions between DWs and the particles in the thermal bath are mediated by the SM Higgs field and only the sufficiently massive SM particles are relevant.Such interactions are highly suppressed by the small mixing between s and h because of the hierarchy v s ≫ v. Furthermore, at sufficiently low temperatures2 , the friction force due to the massive SM particles, whose masses arise from the electroweak symmetry breaking, is significantly damped by the exponentially suppressed number densities of these particles [123].
In addition, the Higgs profile inside the DWs is also relevant to friction.If µ 2 H < 0 and λ HS > 0, the SM Higgs boson mass is given by m . Inside a DW, the S field value is close to zero, and forces the H field to take a large value ∼ v s .Consequently, the SM particles coupled to H become very heavy inside the DW, and their reflection probability with the DW would be highly increased, leading to significant friction [122].Nonetheless, because we have assumed µ 2 H > 0 and λ HS < 0, the H field value would vanish inside the DWs, and the friction can be safely neglected in this study.
Since the friction is negligible, DWs will quickly enter the scaling regime and their energy density evolves as where A ≈ 0.8 ± 0.1 is a numerical factor given by lattice simulation [124].ρ DW ∝ t −1 implies that the DW energy density is diluted more slowly than radiation and matter.Therefore, if DWs are stable, they would soon dominate the evolution of the universe, and it conflicts with cosmological observations.This can be evaded if an explicit Z 2 -violating potential like Eq. (3) presents.
A small Z 2 -violating potential generates a small energy bias between the two minima of the total potential.It leads to a volume pressure force acting on the DWs.Thus, the walls could collapse at a very early epoch before they overclose the universe, and would not cause a cosmological problem.With the Z 2 -violating potential (3), the minima are shifted to where we have neglected the contribution from the |H| 2 S term for v ≪ v s .We define Ŝ(x) ≡ S(x) + δ and rewrite the potential with the redefined scalar field Ŝ: where The two minima of the potential are now located at Ŝ = ±v s .Ŝ = +v s corresponds to the true vacuum, while Ŝ = −v s corresponds to the false vacuum with slightly higher energy.The energy difference between them is [123] V The volume pressure force caused by this energy bias acts on the DWs and tends to make the false vacuum domains shrink.The collapse of DWs begin when the volume pressure force becomes comparable to the tension force.As a result, the annihilation temperature of DWs can be estimated by [123,124] T ann = 34.1 MeV A −1/2 ï g * (T ann ) 10 where g * represents the effective number of relativistic degrees of freedom for the energy density of the plasma and its numerical value depending on the temperature can be found in Ref. [125].
There are two lower bounds on the energy bias between the two minima V bias [123].The first one is given by the requirement that DWs should collapse before they dominate the universe.Moreover, the energetic particles produced from DW collapse could destroy the light elements generated in the big bang nucleosynthesis (BBN).Thus, we should require that DWs annihilate before the BBN epoch.This leads to a second lower bound as The stochastic GWs from collapsing DWs can be estimated by numerical simulations [104,105,124].The frequency spectrum of the SGWB is commonly characterized by where ρ GW is the GW energy density and ρ c is the critical energy density of the universe.At high frequencies, the simulations show that the GW spectrum behaves as Ω GW ∝ f −1 .At small frequencies, the spectrum scales as Ω GW ∝ f 3 because of causality [124,126].The peak of the spectrum at the DW annihilation temperature T ann can be expressed as [124] Ω peak where εGW = 0.7 ± 0.4 is derived from numerical simulation.α * represents the ratio of the GW energy density to the radiation energy density ρ rad at T ann , i.e., Taking into account the dilution of the GW energy density due to the cosmological expansion, the peak amplitude of the SGWB spectrum at the present time can be expressed as [123,124] where g * s denotes the effective number of relativistic degrees of freedom for the entropy density.The present GW peak frequency can be estimated by the Hubble rate at T ann taking into account the redshift effect, given by Thus, the present SGWB spectrum induced by collapsing DWs can be evaluated by In Fig. 1, we show the GW spectra generated by DWs for some benchmark parameters, compared with the reconstructed posterior distributions for the NANOGrav (gray violins) [2] and EPTA (orange violins) [4] signals.We find that the spectra with σ ∼ O(10 17 ) GeV 3 and V bias ∼ O(10 −3 ) GeV 4 can explain the PTA observations.As discussed above, DWs must annihilate before they dominate the energy of the universe.Thus, the time t ann when DWs annihilate should be earlier than the time t dom when DWs would dominate.According to Eqs. ( 18), ( 19) and (24), this gives upper limits on the GW spectra from collapsing DWs.In Fig. 1, the unavailable region corresponding to t ann > t dom is shaded by the brown color.

III. FREEZE-IN DARK MATTER AND THE INDUCED Z 2 -VIOLATING POTENTIAL
So far, the Z 2 -violating potential is introduced by hand.In the following, we will consider it to be generated by loops of fermionic dark matter through a Yukawa interaction with the scalar field S. To be precise, we introduce a Dirac fermion field χ, which is a singlet under all the SM FIG. 1. GW spectra generated by collapsing DWs for four sets of parameters.The gray and orange violins represent the reconstructed posterior distributions for the NANOGrav (gray violins) [2] and EPTA (orange violins) [4] observations of nHz GWs, respectively.The brown region is excluded by the requirement that DWs should annihilate before they dominate the universe.
gauge symmetries.The Lagrangian involving χ is where y χ is the Yukawa coupling constant.When S acquires a nonzero VEV, ⟨S⟩ ≈ ±v s , the mass of χ receives a correction, m In this work, we assume that m χ ≫ y χ v s , so m (±) χ ≈ m χ holds.After reheating, s bosons are in thermal equilibrium with the SM particles due to the |H| 2 S 2 interaction, while χ fermions are assumed to be out of equilibrium with nearly vanishing number density.This requires y χ to be a feeble coupling constant.In this case, χ fermions could be produced via s decays but never reach thermal equilibrium if y χ is extremely small, say, ∼ O(10 −10 ).This is the well-known freeze-in mechanism of DM production [106], and χ acts as a DM candidate.The evolution of the DM number density n χ is determined by the Boltzmann equation [127] where x i ≡ m i /T , and Γ s→χ χ is the s → χ χ partial decay width given by for m χ ≪ m s .The function K1 (x 1 , x 2 , x 3 , η 1 , η 2 , η 3 ) is defined as By solving Eq. ( 28), the χ number density at the present time t 0 can be approximated by [127] n where g * (m s /3) ≈ 108, s 0 is the present entropy density, and M Pl ≈ 2.4 × 10 18 GeV is the reduced Planck mass.The χ relic density is then given by On the other hand, the observed value of the DM relic density is Ω DM h 2 = 0.1200 ± 0.0012 [128], which implies that the Yukawa coupling y χ should be feeble.
In the potential (15), the Z 2 -violating term is characterized by the parameter ϵ, which is related to κ 1 and κ 3 via Eq.( 16).Taking λ S ∼ 0.2, σ ∼ 10 17 GeV 3 , and V bias ∼ 10 −3 GeV 4 , which lead to a GW spectrum accounting for the recent PTA observations, we obtain ϵ ∼ O(10 −26 ) from Eqs. ( 8) and (17).Note that the S χχ Yukawa interaction explicitly breaks the Z 2 symmetry even if the tree-level Z 2 -violating potential is absent.It is natural to conjecture that such an extremely tiny ϵ is originated from the feeble Yukawa interaction through χ loops.
In Fig. 2, we show some one-loop diagrams relevant to the generation of the Z 2 -violating couplings κ 1 and κ 3 .The first two diagrams contains the one-loop tadpole diagrams of S and they give the dominant contributions to κ 1 and κ 3 .Although the third diagram also contributes to κ 3 , it is negligible compared with the second diagram due to a y 3 χ suppression.Our further calculation leads to the ϵ value at the m s scale as where ϵ = 0 at some ultraviolet (UV) scale Λ UV is assumed.Below we assume Λ UV = M Pl .Thus, λ S ∼ 0.2 and v s ∼ 10 5 GeV lead to ϵ(m s ) ∼ 0.6y χ (m χ /m s ) 3 .By taking λ S = 0.2, σ = 10 17 GeV 3 , and V bias = 3.3 × 10 −3 GeV 4 , which correspond to the GW spectrum denoted by the green line in Fig. 1, we obtain v s = 6.19 × 10 5 GeV, m s = 3.91 × 10 5 GeV, ϵ = 3.58 × 10 −26 , T ann = 163 MeV, Ω peak GW h 2 = 9.44 × 10 −8 , and f peak = 2.18 × 10 −8 Hz.Then, y χ and m χ are related by Eq. (33).For this set of parameters, the predicted DM relic density as a function of y χ is shown in Fig. 3, with the upper horizontal axis corresponding to m χ .We find that both the extremely tiny ϵ ∼ O(10 −26 ) and the observed DM relic density Ω DM h 2 = 0.12 can be naturally explained by the feeble Yukawa coupling y χ ∼ O(10 −10 ).Therefore, our theory can simultaneously explain the recent PTA observations of nHz GWs and the DM relic via freeze-in production.

IV. FAVORED PARAMETER RANGES
In this section, we investigate the parameter ranges favored by both the PTA GW observations and the observed DM relic density.Our model has four free parameters, which can be chosen to be λ S , y χ , m s , and m χ .In the following analysis, we fix the quartic coupling λ S = 0.2 to reduce one free parameter.
In Fig. 4, the GW spectra for four benchmark points with m χ = 1.1-2.5 GeV and m s = (3.5-5)× 10 5 GeV are shown.For all these benchmark points, the Yukawa coupling y χ is adjusted to give the mean value of the observed DM relic density Ω χ h 2 = 0.12.We can see that the GW spectrum is quite sensitive to m χ and m s , limiting them varying within roughly one order of magnitude.
In our model, after the annihilation of DWs, most of their energy releases to the SM thermal bath.In this case, the NANOGrav collaboration has reconstructed the posterior distributions of (T ann , α * ) accounting for the observed nHz GW signal, as shown in the left panel of Fig.   in Ref. [2].We use this result to study the favored parameter regions.In Fig. 5, the Yukawa coupling is fixed as y χ = 6.4 × 10 −10 , and the deep blue and light blue regions in the m χ -m S plane correspond to the 68% and 95% Bayesian credible regions favored by the NANOGrav data.It shows a high correlation between m χ and m S , which can be understood by the behavior of the GW peak amplitude, Ω peak GW h 2 ∝ m 10 s /m 6 χ , inferred from Eqs. ( 24) and ( 33).In Fig. 5, Ω χ h 2 = 0.12 corresponds to a red line, which intersects the regions favored by the NANOGrav GW observation.Thus, both the PTA GW signal and the observed DM relic density can be well interpreted.Moreover, the brown and gray regions are excluded by the requirements that DWs should collapse before they dominate the universe and before the BBN epoch, respectively.We find that these constraints do not exclude the NANOGrav 95% Bayesian credible region.
The intersection of the Ω χ h 2 = 0.12 line and the NANOGrav favored regions sensitively depends on the value of y χ .Our calculation shows that such a intersection can only happen for 4.6 × 10 −10 ≲ y χ ≲ 8.7 × 10 −10 .In the left and right panels of Fig. 6, we demonstrate the results for y χ = 4.6 × 10 −10 and y χ = 8.7 × 10 −10 , respectively.In both cases, the Ω χ h 2 = 0.12 line can only touch the edge of the NANOGrav 95% Bayesian credible region.To sum up, for λ S = 0.2, we find that the preferred parameter ranges where our model can simultaneously explain the NANOGrav GW signal and the DM relic density are 4.6 × 10 −10 ≲ y χ ≲ 8.7 × 10 −10 , 0.17 GeV ≲ m χ ≲ 7.5 GeV, 8.1 × 10 4 GeV ≲ m s ≲ 10 6 GeV.(34) V. SUMMARY In this work, we studied the interplay between the SGWB from collapsing DWs and FIMP DM.The newly introduced real scalar field S satisfies a Z 2 symmetry in the tree-level potential, which, however, is violated by the Yukawa coupling y χ with a fermion field χ.The linear and cubic terms of S can be induced by the Yukawa coupling at one-loop level, and they explicitly break the Z 2 symmetry of the potential, leading to an energy bias between the two minima.Thus, after the spontaneous breaking of the Z 2 symmetry, unstable DWs would be formed.We considered that the Yukawa coupling is feeble, i.e., y χ ∼ O(10 −10 ), and the χ fermions become FIMPs that are produced by the freeze-in mechanism, accounting for dark matter in the universe.
On the other hand, four PTA collaborations have recently reported evidences of a SGWB at nHz frequencies, which could be produced by the collapse of DWs in the early universe.Since the tiny Z 2 -violating potential induced by the feeble Yukawa coupling leads to unstable DWs, it is possible that our model can explain the PTA data.Comparing with the posterior distributions in the GW spectrum reconstructed by the NANOGrav and EPTA data, our analysis showed that the Z 2 -violating coefficient should be as tiny as ϵ ∼ 10 −26 , which can be naturally induced by the feeble Yukawa couplings y χ ∼ O(10 −10 ) at one-loop level.Thus, our scenario is very suitable for interpreting the PTA observations of the nHz SGWB.
Moreover, we investigated the parameter regions where both the PTA GW observations and the DM relic density can be simultaneously explained.We found that the parameters should satisfy y χ ∈ (4.6 × 10 −10 , 8.7 × 10 −10 ), m χ ∈ (0.17, 7.5) GeV, and m s ∈ (8.1 × 10 4 , 10 6 ) GeV for a fixed quartic scalar couplings λ S = 0.2.The corresponding regions also fulfill the requirements that DWs should collapse before they overclose the universe and they should not affect the BBN.

I. Introduction 2 II. Domain walls and gravitational waves 3 III 7 IV
. Freeze-in dark matter and the induced Z 2 -violating potential

FIG. 2 .
FIG. 2. Examples of one-loop diagrams for generating the Z 2 -violating potential terms.

h 2 t
FIG.4.Same as in Fig.1, but for GW spectra corresponding to different values of m χ and m s with λ S = 0.2 fixed and y χ adjusted to give Ω χ h 2 = 0.12.

= 6 . 4 × 10 − 10 FIG. 5 .
FIG.5.Parameter regions favored by the NANOGrav GW signal in the m χ -m s plane for y χ = 6.4×10 −10 .The deep blue and light blue regions corresponds to the 68% and 95% Bayesian credible regions favored by the NANOGrav data, respectively.The red line denotes the mean value of the Planck observation of the DM relic density, Ω χ h 2 = 0.12.The brown and gray regions are excluded because DWs would dominate the universe and would inject energetic particles to affect BBN, respectively.