Imprints of an early matter-dominated era arising from dark matter dilution mechanism on cosmic string dynamics and gravitational wave signatures

Modern cosmology originates from the study of the universe expansion and element formation [1]. Further observations give rise to the standard $ \Lambda\mathrm{CDM} $ cosmological model, which suggests that after inflation and reheating, the universe successively experiences a radiation-dominated (RD) era, a matter-dominated era (MD), and a dark-energy-dominated era. Main evidences involve the observations of the cosmic expansion [2], the big bang nucleosynthesis (BBN) [3, 4], the cosmic microwave background [5−8], and the accelerating expansion at late times [9, 10].

Nevertheless, these observational results can hardly date back to pre-BBN eras. Therefore, it is essential to maintain an open perspective regarding the early history of the universe before BBN [11]. Various hypotheses beyond the standard cosmic history have been put forward, such as an early matter-dominated (EMD) era [12−19], a kination-dominated era [20−25], and an intermediate inflationary era [26−29]. Verification of such novel hypotheses about the cosmic history necessitates observational probes capable of accessing epochs prior to BBN, when the universe was opaque to photons and traditional electromagnetic detection methods are ineffective. Nonetheless, this challenge can be addressed through observations of gravitational waves (GWs), a new messenger first detected in 2015 [30]. Unlike electromagnetic radiation, GWs can propagate freely through space, preserving information from the early universe and reaching us in the present day.

To explore the pre-BBN history, the stochastic gravitational wave background (SGWB) originating from cosmic strings (CSs) provides a compelling avenue of investigation [31−42]. CSs are one-dimensional topological defects predicted by a variety of new physics theories beyond the standard model (SM), particularly those involving a spontaneously broken $ \mathrm{U}(1) $ symmetry [43, 44]. CS loops are expected to persist as long-lasting sources, emitting GWs from their formation epoch to the present day. Once formed, the CS network rapidly evolve into a scaling regime [45−48], where their correlation length scales proportionally with the Hubble radius. As a result, the resulting SGWB spectrum encodes a wealth of information about the cosmic history. Furthermore, the SGWB generated by CSs spans an exceptionally broad frequency range, making it detectable across a diverse array of GW experiments. These include pulsar timing arrays (PTAs) operating in $ 10^{-9}\text{–}10^{-7}\; \text{Hz} $ [49−51], space-borne interferometers sensitive to $ 10^{-4}\text{–}10^{-1}\; \text{Hz} $ [52−54], and ground-based interferometers covering $ 10\text{–}10^3\; \text{Hz} $ [55−58]. Thus, future GW experiments hold immense potential to unravel the mysteries of the early universe.

In this work, we explore how the SGWB spectrum originated from a preexisting CS network is modified by an EMD era that arises within the framework of the dark matter (DM) dilution mechanism [59−68]. This mechanism addresses the DM overproduction problem by relying on significant entropy production from the decays of a long-lived particle, referred to as the "dilutor". For the dilution mechanism to be effective, the dilutor must dominate the energy density of the universe for a finite period, thereby introducing an EMD era that interrupts the conventional RD era. To quantify the impact of this EMD era, we will trace the evolution of the total energy density of the universe by solving the relevant Boltzmann equations. This allows us to derive the resulting SGWB spectrum, which reflects the influence of the EMD era on the cosmic string dynamics and the GW emissions.

This paper is outlined as follows. In Section II, we review the formation of CS loops and analyze their number densities during RD and MD eras, which are essential for assessing the SGWB spectrum induced by CS loops. In Section III, we investigate an EMD era emerging from the DM dilution mechanism and examine its effects on the evolution of the CS network, the loop number densities, and the resulting SGWB spectrum. Finally, in Section Ⅳ, we summarize our findings.

In a scalar field theory with a global or gauge $ \mathrm{U}(1) $ symmetry, CSs could be formed after the spontaneous breaking of the $ \mathrm{U}(1) $ symmetry in the early universe. They are one-dimensional topological defects concentrating energy of the scalar field (and the gauge field for the gauged case) [69]. In the Nambu-Goto approximation, CSs are described as infinitely thin objects with tension μ, which is the energy per unit length. The dimensionless quantity $ G\mu $ is commonly used to represent the CS tension. Considering the dynamics of CSs, these combined factors establish CSs as promising sources of GWs [35, 70].

Cosmic strings are generated randomly in the early universe, leading to the formation of a CS network. Long strings with super-horizon lengths could intersect with each other to form CS loops, whose relativistic oscillations can effectively emit GWs. Although small structures in long strings can also produce GWs, their contributions are generally negligible, compared to those generated by loops [37, 71, 72]. In practice, it is sufficient to consider only the stochastic GWs originating from CS loops.

The GW emission power of CS loops is given by [72, 73]

where the coefficient Γ is estimated to be about 50 [74], and G represents the Newtonian gravitational constant. The frequencies of GWs emitted by a CS loop of length l are [70]

where n denotes the harmonic modes of the loop oscillation. Thus, the expression for the power can be rewritten as

where $ P_n $ is the dimensionless emission power in units of $ G\mu^2 $ for a single mode n, which can be estimated by numerical simulations for GW emissions in RD and MD eras [74].

Introducing $ n_\mathrm{CS}(l,t) \,\mathrm{d}l $ as the number density of CS loops of length l at cosmic time t, the energy density of GWs emitted from CS loops per unit time at the emission time $ t_\mathrm{e} $ can be expressed as

Using Eq. (2), we derive

Note that this result corresponds to the GW emission time $ t_\mathrm{e} $, but GW experiments only receive signals at the present time $ t_0 $. Hence, the effect of the cosmological redshift must be accounted for, and the present GW frequency is given by $ f = a(t_\mathrm{e})f_\mathrm{e} $, while the present GW energy density is $ \rho_\mathrm{GW} = \rho_\mathrm{GW}(t_\mathrm{e}) a^4(t_\mathrm{e}) $, where $ a(t) $ is the scale factor normalized such that $ a(t_0)=1 $. Integrating Eq. (5) over time, we arrive at

where

with $ t_\ast $ denoting the cosmic time when CS loops start to radiate GWs.

The frequency spectrum of the SGWB induced by CS loops is commonly characterized by a dimensionless quantity

where $ \rho_\mathrm{c} = 3H_0^2 /(8\pi G) $ is the present critical density and $ H_0 = 100h\; \mathrm{km}\; \mathrm{s}^{-1}\; \mathrm{Mpc}^{-1} $ is the Hubble constant with $ h = 0.674\pm 0.005 $ [8]. In order to calculate Eq. (8), we need to know the evolution of the scalar factor $ a(t) $, which can be computed from

where $ H(t) $ is the Hubble expansion rate.

In the standard $ \Lambda\mathrm{CDM} $ model, the Hubble rate can be expressed as [75]

where $ z = a^{-1} -1 $ is the cosmological redshift. $ \Omega_\mathrm{r} = 1.68\times (5.38\pm 0.15)\times 10^{-5} $, $ \Omega_\mathrm{m} = 0.315\pm 0.007 $, and $ \Omega_\Lambda = 0.685\pm 0.007 $ [76] are cosmological constants representing the energy fraction of radiation, matter, and dark energy, respectively.

is a function of the redshift z introduced to account for the changes in relativistic degrees of freedom, where $ g_\star(z) $ and $ g_{\star s}(z) $ are the effective numbers of relativistic degrees of freedom for the energy and entropy densities, respectively. Considering the evolution of $ g_\star(z) $ and $ g_{\star s}(z) $, we can approximate $ \mathcal{G}(z) $ as a piecewise function which changes its value at the epochs of electron-positron annihilation $ (z\simeq 10^9) $ and QCD phase transition $ (z\simeq2\times10^{12}) $: [75]

However, if an EMD era exists, the evolution of the Hubble rate $ H(t) $ would be altered, necessitating modifications to the expressions above. This situation will be discussed in detail in Section III.

For calculating the SGWB spectrum induced by CS loops, it is crucial to estimate the loop number density distribution $ n_\mathrm{CS}(l,t) $. There are two models for evaluating $ n_\mathrm{CS}(l,t) $ based on numerical simulations and the scaling nature of the CS network. The first is the BOS model, introduced by Blanco-Pillado, Olum, and Shlaer [77]. This model uses the horizon distance as the sole kinematic scale to describe the scaling behavior and extrapolates the loop production functions from numerical simulations to derive $ n_\mathrm{CS}(l,t) $ in both RD and MD eras. The second is the LRS model, proposed by Lorenz, Ringeval, and Sakellariadou [78]. In this model, the number density distribution of produced loops per unit time is assumed to follow a power law in the scaling regime, and the gravitational backreaction effect is included to reduce loop production below a certain scale. Compared to the BOS model, the SGWB spectrum predicted by the LRS model has a significantly higher amplitude. Consequently, the constraint on the CS tension derived from the LIGO-Virgo O3 dataset is much stringent for the LRS model than for the BOS model [79]. Additionally, future space-borne interferometers will more easily detect the SGWB predicted by the LRS model [80].

Nevertheless, these two models rely on the scaling behavior of the CS network, which may be violated by the inclusion of an EMD era. To account for this nonscaling effect, we instead adopt the velocity-dependent one-scale (VOS) model [72, 81, 82], which describes the dynamical evolution of the CS network in a generic manner. This model has also been used to interpret the positive evidence of the nHz SGWB reported by PTA experiments [83−85]. We will demonstrate below that the predictions of the VOS model align with those of the BOS model when the scaling regime is achieved and appropriate loop production functions are assumed. Furthermore, the VOS model can be directly applied to the cosmic history with an EMD era, as will be discussed in Subsection ⅢC.

In the VOS model, the dynamics of the CS network is characterized by two fundamental parameters, the correlation length L and the root-mean-square (RMS) velocity v of string segments. Thus, the energy density of long strings can be expressed as $ \rho = \mu / L^2 $, and its evolution in the universe is governed by [81]

where the friction effect is neglected. The term $ -2H \rho $ accounts for the dilution and stretching of the strings due to the cosmic expansion, while the term $ -2 H v^2 \rho $ captures the energy loss caused by the redshifting of velocities. The term $ -\tilde{c}v\rho/L $ represents the energy removed through loop formation, where $ \tilde{c}\simeq 0.23 $ [86] is the loop chopping efficiency. Additionally, the evolution of the RMS velocity v is determined by [81]

with [86]

Introducing a dimensionless quantity $ \xi \equiv L/t $, which represents the correlation length L normalized by the cosmic time t, the evolution equations above can be rewritten as

By numerically solving these equations, one finds that the solutions rapidly converge to constant values of ξ and v, independent of the initial conditions. This indicates that the CS network quickly evolves into a linear scaling regime characterized by $ L \propto t $. By setting $ \dot{\xi} = \dot{v} = 0 $, the scaling solutions for RD and MD eras can be derived as [87]

The increase in the energy density of CS loops $ \rho_\circ $ is driven by the energy transfer from long strings, expressed as [37, 87, 88]

where the Lorentz factor $ \gamma_v \equiv (1-v^2)^{-1/2} $ incorporates the energy loss caused by the redshifting of the loop velocity, and $ \mathcal{F} $ is a coefficient that characterizes the fraction of slow, large loops responsible for the dominant contribution to GW emissions. The values of $ \mathcal{F} $ differ between RD and MD eras, denoted as $ \mathcal{F}_\mathrm{r} $ and $ \mathcal{F}_\mathrm{m} $, respectively. These values can be determined by comparison with numerical simulations in the scaling regime [77, 89].

The energy in a CS loop of length l is expressed as $ \mu l $. Thus, the loop production function, defined by [37, 87]

quantifies the rate of increase in the number density of CS loops per unit length per unit time. Because of GW emission, a CS loop of length $ l' $ at time $ t' $ loses energy $ \Gamma G\mu^2 (t-t') $ by time t and consequently shrinks to a length $ l = l' - \Gamma G\mu (t-t') $. Therefore, the loop number density per unit loop length at time t is given by

with $ l'(t) = l+\Gamma G\mu (t-t') $. $ t_\mathrm{ini} $ denotes the initial time of loop production, and the scale factors account for the dilution because of cosmic expansion. To determine $ n_\mathrm{CS}(l,t) $ for calculating the energy density of emitted GWs according to Eq. (4), knowledge of the loop production function $ \mathcal{P}(l,t) $ is required.

During an RD era, numerical simulations [77] reveal that most large loops are generated with a length l characterized by a fraction $ \alpha_\mathrm{r} = l/L $ of the correlation length L. $ \alpha_\mathrm{r} $ satisfies $ \alpha_\mathrm{r} \xi_\star \simeq 0.1 $, where a quantity with a subscript $ \star $ is evaluated at the loop production time $ t_\star $. Consequently, the loop production function can be accurately approximated by

where $ \delta(x) $ denotes the Dirac δ function, which fixes the loop production time $ t_\star $ through the relation $ \alpha_\mathrm{r}\xi_\star = {l_\star}/{t_\star} $. Substituting the expression for $ \mathcal{P}_\mathrm{r}(l,t) $ into Eq. (22) with $ a \propto t^{1/2} $ in the RD era, we derive an analytical expression for the loop number density [87, 88, 90],

where $ \Theta(x) $ represents the Heaviside step function. By accounting for the decrease in loop length due to GW emission, we have $ \alpha_\mathrm{r}\xi_\star {t_\star} = {l_\star} = l + \Gamma G \mu(t - t_\star) $, leading to

In the scaling regime of the RD era, we have $ \xi_\star = \xi_\mathrm{r} $, $ v_\star = v_\mathrm{r} $, and $ \dot{\xi}_\star = 0 $. Adopting

the loop number density formed in an RD era for $ \alpha_\mathrm{r} \xi_\star = 0.1 \gg \Gamma G \mu $ reduces to

which agrees with the result from numerical simulations presented in Ref. [77].

On the other hand, numerical simulations for an MD era suggest that loop production occurs on a wider range of scales compared to that in an RD era and can be effectively described by a power-law distribution of $ l^{-1.69} $ with a cutoff [77]. This motivate us to assume the loop production function as

where the theta function gives a cutoff at $ l = \alpha_\mathrm{m} L $ with $ \alpha_\mathrm{m} $ satisfying $ \alpha_\mathrm{m} \xi_\star \simeq 0.18 $ [77]. Using $ a \propto t^{2/3} $ in the MD era, we derive the loop number density

where $ t_\mathrm{m} $ is the initial time of the MD era. Under the assumption that the $ t' $-dependence of $ \xi' $ and $ v' $ is negligible and that $ l' = l + \Gamma G \mu(t - t') \simeq l $ for $ \Gamma G \mu \ll l/(t - t') $, we perform integration by parts to obtain

where $ t_\star $ is determined by $ \alpha_\mathrm{m}\xi_\star {t_\star} = l + \Gamma G \mu(t - t_\star) $. For $ t \gg t_\mathrm{m} $, the last term in the above expression can be omitted, leading to

Moreover, for small $ \Gamma G\mu $, we can use the approximations $ t_\star \simeq (\alpha_\mathrm{m}\xi_\star)^{-1}(l + \Gamma G\mu) $ and $ t' \simeq (l/t')^{-1}(l + \Gamma G\mu) $ to express the loop number density in the MD era as

In the scaling regime, where $ \xi = \xi_\star = \xi_\mathrm{m} $ and $ v = v_\star = v_\mathrm{m} $, taking $ \alpha_\mathrm{m} \xi_\star = 0.18 $ and

the loop number density formed in the MD era becomes

which is also consistent with the simulation result in Ref. [77].

We now analyze the SGWB spectrum generated by CS loops within the standard $ \Lambda\mathrm{CDM} $ cosmological model, utilizing the expressions (27) and (34) for the loop number density in the scaling regime. In the $ \Lambda\mathrm{CDM} $ model, the MD era succeeds the RD era at the time of matter-radiation equality, $ t_\mathrm{eq} = 51.1\pm 0.8 \; \mathrm{kyr} $ [76]. Notably, CS loops formed during the RD era and surviving into the MD era also contribute to the loop number density in the MD era, as given by [77]

By inserting the scale factor $ a(t) $ and the loop number density $ n_\mathrm{CS}(l,t) $ into Eqs. (7) and (8), the SGWB spectrum can be obtained.

In Fig. 1, we demonstrate the SGWB spectra in the $ \Lambda\mathrm{CDM} $ model for the CS tension parameter $ G\mu = 10^{-9}, 10^{-11}, 5\times 10^{-13}, 10^{-15} $, and $ 5\times 10^{-18} $. At sufficiently high frequencies, the spectra tend to be flat, primarily due to GW emissions during the RD era [77]. For smaller values of $ G\mu $, the total GW emission power is reduced, allowing CS loops to survive longer. At a specific cosmic time, the average length of loops would be smaller, leading to GWs emitted at higher frequencies. As a result, the SGWB spectrum shifts downward in amplitude and rightward in frequency as $ G\mu $ decreases [74].

Figure 1.  (color online) SGWB spectra for various values of $ G\mu $ in the standard $ \Lambda\mathrm{CDM} $ cosmological model.

In the following section, we will examine the effects of an additional EMD era inserted into the RD era, focusing on its impact on the scale factor, the loop number density, and the SGWB spectrum. This modification may violate the scaling behavior, necessitating the use of the more general expressions for the loop number density as in Eqs. (24) and (32).

In this section, we discuss an EMD era motivated by the DM dilution mechanism and study its impact on the SGWB spectrum generated by cosmic strings.

An EMD era is an MD period spanning from cosmic time $ t_1 $ to $ t_2 $, embedded within the conventional RD era. Thus, there exists an initial RD era before $ t_1 $, and a subsequent RD era after $ t_2 $. The second RD era ends at the time of matter-radiation equality $ t_\mathrm{eq} $, after which the universe transitions into the final MD era, corresponding to the traditional cosmological evolution. This timeline is depicted in Fig. 2.

Figure 2.  (color online) Sketch of the cosmic timeline featuring an EMD era. In chronological order, the universe evolves through the original RD era, the EMD era, the second RD era, and the final MD era.

We propose the following assumption for the origin of such an EMD era. A massive, long-lived particle Y decouples from thermal equilibrium and transitions into nonrelativistic matter, dominating the universe as the temperature decreases and initiating an EMD era at time $ t_1 $. The slow decays of Y into SM particles release entropy into the plasma. In the context of the DM dilution mechanism, this entropy injection dilutes the DM abundance to the observed value. Once the majority of Y particles have decayed, the EMD era concludes at time $ t_2 $.

Now we discuss the evolution of the related energy densities. The radiation and matter energy densities, $ \rho_\mathrm{r} $ and $ \rho_\mathrm{m} $, scale with the scale factor a as [91]

The scale factor evolves as $ a \propto t^{1/2} $ in an RD era and $ a \propto t^{2/3} $ in an MD era. Starting from an initial temperature $ T_\mathrm{ini} $ in the original RD era, the Y particles constitute a matter component with energy density $ \rho_\mathrm{Y} $, while the plasma of SM particles has energy density $ \rho_\mathrm{SM} $. Consequently, the energy densities evolve as

during the original RD era, and

during the EMD era. Notably, the energy density of Y decreases more slowly than that of SM particles during the original RD era, leading to the onset of the EMD era.

In this subsection, we discuss the DM dilution mechanism as the origin of the EMD era. As severe constraints from DM detection experiments have reduced interest in weakly interacting massive particles with masses around the electroweak scale as DM candidates, a lighter DM candidate X has gained prominence due to its ability to evade direct detection bounds. Nevertheless, thermal production of X particles with low annihilation cross sections typically results in an overproduction problem [92, 93]. Within the DM dilution mechanism [60−64, 66, 67], the overproduction of DM particles is mitigated by entropy injection from the decays of a dilutor particle Y, which dominates the universe for a period, thereby inducing an EMD era.

Following Ref. [68], we discuss the DM dilution mechanism within the minimal left-right symmetric model (LRSM) as a concrete and illustrative example. In this model, the lightest right-handed neutrino $ N_1 $ serves as the light DM candidate X, and the next-to-lightest right-handed neutrino $ N_2 $ or a neutral Higgs boson Δ acts as the dilutor Y. As a Majorana fermion, the DM candidate X has two degrees of freedom and its number density in thermal equilibrium and in the relativistic limit is given by [91]

where T denotes the temperature and $ \zeta(x) $ represents the Riemann zeta function. For sufficiently weak interactions, X particles freezes out from the plasma at a temperature $ T_\mathrm{FO}\gg m_X $, leaving a number density in the comoving volume of

where $ s = 2\pi^2 g_{\star s}(T) T^3/45 $ is the entropy density.

If X particles evolve without interference, $ Y_X $ remains conserved, leading to a present-day relic energy density fraction of [68]

where $ s_0 = 2891.2\; \mathrm{cm}^{-3} $ and $ \rho_\mathrm{c} = 1.05\times 10^{-5} h^2\; \mathrm{GeV}/\mathrm{cm}^3 $ [76] are the entropy density and the critical energy density at present, respectively. Thus, for $ m_X \lesssim \text{keV} $, $ \Omega^0_X $ exceeds the observed cold DM relic density $ \Omega_\mathrm{CDM} = 0.265 $ [76] by at least one order of magnitude, causing an overproduction problem.

Introducing a long-lived dilutor Y with a mass $ m_Y $ much larger than $ m_X $ can effectively address the overproduction problem of X particles. First, during the RD era, both Y and X particles decouple relativistically at a similar temperature, resulting in comparable yields, $ Y_Y \simeq Y_X $. Second, because of $ m_Y\gg m_X $, Y particles become nonrelativistic at a relatively high temperature, while X particles remain relativistic. Consequently, Y particles quickly dominate the energy density of the universe, initiating an EMD era. Finally, when the lifetime of Y particles comes to an end, they decay into SM particles and X particles, injecting entropy and consequently diluting the energy density $ \rho_X $ of X particles.

This situation can be analyzed using the Boltzmann equations [67]

where $ \Gamma_Y $ is the decay width of the dilutor Y, y represents the energy fraction carried away by X particles among all the decay products of Y, and $ B_X $ is the branching ratio of the decay channel $ Y \rightarrow X $. Assuming that X particles remain ultra-relativistic throughout the aforementioned process, and using

we can solve the Boltzmann equations.

In the minimal LRSM with the $ \mathrm{SU}(3)_\mathrm{C} \times \mathrm{SU}(2)_\mathrm{L} \times \mathrm{SU}(2)_\mathrm{R} \times \mathrm{U}(1)_{B-L} $ gauge symmetry [94, 95] considered in Ref. [68], there are two scenarios for the DM dilution mechanism, where the DM candidate X is the lightest right-handed neutrino $ N_1 $.

● Scenario 1: the dilutor Y is the next-to-lightest right-handed neutrino $ N_2 $, which undergoes a three-body decay mediated by a right-handed gauge boson $ W^\pm_\mathrm{R} $ into two charged leptons $ \ell \ell' $ and one $ N_1 $. The related right-handed charged current interactions are described by the Lagrangian

where g is the unified $ \mathrm{SU}(2)_\mathrm{L}\times \mathrm{SU}(2)_\mathrm{R} $ gauge coupling, and $ V_\mathrm{PMNS}^\mathrm{R} $ and $ V_\mathrm{CKM}^\mathrm{R} $ are the right-handed Pontecorvo-Maki-Nakagawa-Sakata (PMNS) and Cabibbo-Kobayashi-Maskawa (CKM) matrices that connect the flavor and mass bases. The decay channels of Y include $ N_2 \rightarrow N_1\ell\ell' $, $ N_2 \rightarrow \ell q\bar{q}' $, and $ N_2 \rightarrow \ell W $.

● Scenario 2: the dilutor Y is the neutral Higgs boson Δ, which is associated with the origin of the Majorana neutrino masses and can decay into two $ N_1 $ particles. The relevant Yukawa couplings are given by the Lagrangian

where the family indices are suppressed. Φ is an $ \mathrm{SU}(2)_\mathrm{L}\times \mathrm{SU}(2)_\mathrm{R} $ bidoublet scalar, with $ \tilde{\Phi} = \mathrm{i}\sigma_2 \Phi^\ast \mathrm{i}\sigma_2 $. $ \Delta_\mathrm{L} $ and $ \Delta_\mathrm{R} $ are $ \mathrm{SU}(2)_\mathrm{L} $ and $ \mathrm{SU}(2)_\mathrm{R} $ triplet scalars, respectively. $ Q_{\mathrm{L},\mathrm{R}} $ and $ L_{\mathrm{L},\mathrm{R}} $ denote left-handed and right-handed quarks and leptons. $ Y_q $, $ \tilde{Y}_q $, $ Y_l $, $ \tilde{Y}_l $, $ Y_{\Delta_\mathrm{L}} $, and $ Y_{\Delta_\mathrm{R}} $ represent the Yukawa coupling matrices. The dilutor Δ is the Higgs boson arising from the electrically neutral component of $ \Delta_\mathrm{R} $. Its decay channels involve $ \Delta \rightarrow N_1 N_1 $, $ \Delta \rightarrow W_\mathrm{R}^\ast W_\mathrm{R}^\ast \to \text{4 fermions} $, $ \Delta \rightarrow \gamma \gamma $, and so on.

We will consider both scenarios in the following analysis.

We solve the Boltzmann equations for the benchmark points (BPs) specified in Tables 1 and 2 for Scenarios 1 and 2, respectively. The values of y, $ m_{N_1} $, $ m_{N_2} $, $ m_\Delta $, $ m_{W_\mathrm{R}} $, the ratio $ \tan\beta $ of two vacuum expectation values from the bidoublet, and the mixing angle $ \theta_{\Delta h} $ between Δ and the $ 125\; \text{GeV} $ Higgs boson h are adopted based on the results presented in Ref. [68]. The dilutor decay width $ \Gamma_Y $ and the branching ratio $ B_X $ are calculated using the formulas provided therein. All BPs are selected to ensure the correct DM relic abundance through the DM dilution mechanism. BP1a and BP2a will be the primary focus of the following analysis, while the remaining BPs will provide additional results for comparison.

Scenario 1	BP1a	BP1b	BP1c
Common	$ m_{N_2}= 200 \mathrm{GeV} $, $ y = 0.35 $, $ \tan\beta = 0.5 $
$ m_{N_1} $	$ 6.5 \mathrm{keV} $	$ 10 \mathrm{keV} $	$ 30 \mathrm{keV} $
$ m_{W_\mathrm{R}} $	$ 5\times 10^7 \mathrm{GeV} $	$ 6\times 10^7 \mathrm{GeV} $	$ 7\times 10^7 \mathrm{GeV} $
$ \Gamma_Y $	$ 2.22\times 10^{-23} \mathrm{GeV} $	$ 1.07\times 10^{-23} \mathrm{GeV} $	$ 5.77\times 10^{-24} \mathrm{GeV} $
$ B_X $	$ 4.41\times 10^{-3} $	$ 4.41\times 10^{-3} $	$ 4.41\times 10^{-3} $

Table 1.  Benchmark points in Scenario 1 where $ N_2 $ serves as the dilutor Y.

Scenario 2	BP2a	BP2b	BP2c
Common	$ m_{N_1}= 6.5 \mathrm{keV} $, $ y = 1 $, $ \theta_{\Delta h} = 0 $
$ m_\Delta $	$ 1 \mathrm{TeV} $	$ 10^3 \mathrm{TeV} $	$ 10^6 \mathrm{TeV} $
$ m_{W_\mathrm{R}} $	$ 10^{11} \mathrm{GeV} $	$ 3.16\times 10^{12} \mathrm{GeV} $	$ 10^{14} \mathrm{GeV} $
$ \Gamma_Y $	$ 3.51\times 10^{-21} \mathrm{GeV} $	$ 3.51\times 10^{-15} \mathrm{GeV} $	$ 3.51\times 10^{-9} \mathrm{GeV} $
$ B_X $	$ 5.11\times 10^{-12} $	$ 5.11\times 10^{-18} $	$ 2.42\times 10^{-20} $

Table 2.  Benchmark points in Scenario 2 where Δ serves as the dilutor Y.

We adopt the initial temperature as $ T_\mathrm{ini} = m_Y/10 $, at which Y particles have decoupled with $ Y_Y \simeq Y_X $ and become nonrelativistic, while X particles are relativistic. The initial conditions for the energy densities of Y, X, and SM particles are

The obtained energy densities $ \rho_Y $, $ \rho_X $, and $ \rho_\mathrm{SM} $ as functions of cosmic time t are illustrated in Fig. 3 for BP1a in Scenario 1 and BP2a in Scenario 2.

Figure 3.  (color online) Evolution of the energy densities $ \rho_Y $, $ \rho_X $, and $ \rho_\mathrm{SM} $ for BP1a in Scenario 1 with the dilutor $ N_2 $ (a) and BP2a in Scenario 2 with the dilutor Δ (b).

As we can see from Figs. 3(a) and 3(b), the initial energy density of the dilutor Y in both BPs is lower than $ \rho_\mathrm{SM} $ by one order of magnitude. Nonetheless, the decrease in $ \rho_Y $ is slower than that of $ \rho_\mathrm{SM} $ as the universe expands, leading to $ \rho_Y = \rho_\mathrm{SM} $ at $ t_1 \simeq 3\times {{10}^{17}}\; \text{GeV}^{-1} $ and $ t_1 \simeq 2\times {{10}^{16}}\; \text{GeV}^{-1} $ in BP1a and BP2a, respectively. This indicates the beginning of the EMD era. Subsequently, Y particles dominate the universe until their decays become effective, leading to the end of the EMD era at $ t_2 \simeq 5\times {{10}^{22}}\; \text{GeV}^{-1} $ for BP1a and $ t_2 \simeq 3\times {{10}^{20}}\; \text{GeV}^{-1} $ for BP2a. Since $ \Gamma_Y $ in BP1a is lower than that in BP2a, as shown in Tables 1 and 2, the decays of Y particles occur at later times, leading to a longer duration of the EMD era. the Y decays inject entropy into the plasma, increasing $ \rho_\mathrm{SM} $, the plasma temperature, and the entropy density s. Consequently, the X comoving number density $ Y_X = n_X/s $ is significantly diluted, ensuring that the relic abundance of X particles being consistent with the observation in both BPs.

By substituting the obtained $ \rho_Y $, $ \rho_X $, and $ \rho_\mathrm{SM} $ into Eq. (45) and solving Eq. (9) with $ a(t_0) = 1 $, we determine the evolution of the scale factor $ a(t) $ in both BPs, as illustrated in Fig. 4. Compared with the standard $ \Lambda\mathrm{CDM} $ cosmological model, the presence of the EMD era reduces the scale factor prior to $ t_2 $. This occurs because the scale factor evolves as $ a \propto t^{2/3} $ during an MD era, increasing more rapidly with time than $ a \propto t^{1/2} $ during an RD era. Therefore, and a smaller a is required at the onset of the EMD era to ensure $ a(t_0) = 1 $ at the present time $ t_0 $.

Figure 4.  (color online) Evolution of the scale factor $ a(t) $ (blue solid lines) for BP1a in Scenario 1 with the dilutor $ N_2 $ (a) and BP2a in Scenario 2 with the dilutor Δ (b). For comparison, the result from the standard $ \Lambda\mathrm{CDM} $ model, which does not include the EMD era, is represented by red dotted lines.

Since the insertion of the EMD era affects the evolution of the CS network in the VOS model, it is necessary to examine how the normalized correlation length ξ and the RMS velocity v change with time. We should consider these variations when calculating the number density of CS loops during the EMD era. To achieve this, we solve Eqs. (16) and (17) with the Hubble rate $ H(t) $ modified by the presence of the EMD era, and obtain ξ and v being functions of cosmic time t for BP1a and BP2a, as demonstrated by blue solid lines in Fig. 5. For comparison, the results of the standard $ \Lambda\mathrm{CDM} $ model without the EMD era are also presented as red dotted lines.

Figure 5.  (color online) Evolution of the normalized correlation length ξ (left panels) and the RMS velocity v (right panels) of the CS network in the presence of the EMD era (blue solid lines) for BP1a (upper panels) and BP2a (lower panels). The wheat-colored bands denote the duration of the EMD era, while the vertical dashed lines indicate $ t_\mathrm{eq} $. For comparison, the results from the standard $ \Lambda\mathrm{CDM} $ model are also plotted as red dotted lines.

During the RD era in the $ \Lambda\mathrm{CDM} $ model, ξ and v basically maintain the scaling values $ \xi_\mathrm{r} = 0.271 $ and $ v_\mathrm{r} = 0.662 $, as stated in Eq. (18). Some deviations around $ t \sim 10^{20}\; \text{GeV}^{-1} $ and $ t \sim 10^{26}\; \text{GeV}^{-1} $ are attributed to changes in relativistic degrees of freedom shown in Eq. (12). For $ t \gtrsim t_\mathrm{eq} $, the universe transits to the MD era, resulting in an increase in ξ and a decrease in v. It is important to note that the time axis in Fig. 5 is on a logarithmic scale, indicating that the changes in ξ and v occur relatively slowly. Nonetheless, this represents a mild nonscaling effect.

When the EMD era occurs, ξ increases, while v decreases, and both quantities gradually return to their scaling values after the EMD era ends. Therefore, the EMD era also introduces a nonscaling effect, which will be incorporated into the calculation of the CS loop number density. Furthermore, we estimate the magnitudes of the time derivatives $ \dot{\xi} $ and $ \dot{v} $ by calculating the evolution of their ratios to $ \xi/t $ and $ v/t $, respectively, as shown in Fig. 6. We find that, during the EMD era for both BPs, $ |t\dot{\xi}/\xi| $ and $ |t\dot{v}/v| $ are smaller than $ 0.15 $ and $ 0.03 $, respectively. Thus, $ \dot{\xi} $ and $ \dot{v} $ are negligible compared to $ \xi/t $ and $ v/t $, and this will simplify the calculations below.

Figure 6.  (color online) Evolution of $ t\dot{\xi}/\xi $ (left panels) and $ t\dot{v}/v $ (right panels) in the presence of the EMD era for BP1a (upper panels) and BP2a (lower panels).

Following the derivation from Eq. (29) to Eq. (32), we evaluate the CS loop number density formed during the EMD era as

where $ \Gamma G \mu \ll \alpha_\mathrm{m}\xi_\star=0.18 $ has been used in the third step.

In addition, similar to Eq. (35), we must account for the number densities of CS loops produced in a preceding era and surviving into subsequent eras, with respect to the timeline in Fig. 2. For clarity, we define the notations for these loop number densities as follows.

● $ n_\mathrm{CS}^\mathrm{rm\ast} $: CS loops formed in the original RD era and surviving into the EMD era.

● $ n_\mathrm{CS}^\mathrm{mr} $: CS loops formed in the EMD era and surviving into the second RD era.

● $ n_\mathrm{CS}^\mathrm{rm} $: CS loops formed in the second RD era and surviving into the final MD era.

● $ n_\mathrm{CS}^\mathrm{rmr} $: CS loops formed in the original RD era and surviving into the second RD era.

● $ n_\mathrm{CS}^\mathrm{mrm} $: CS loops formed in the EMD era and surviving into the final MD era.

● $ n_\mathrm{CS}^\mathrm{rmrm} $: CS loops formed in the original RD era and surviving into the final MD era.

Following the approach in Ref. [77] for accounting for the cosmological expansion, $ n_\mathrm{CS}^\mathrm{rm} $ can be determined by

where $ l_\mathrm{eq} = l + \Gamma G\mu (t-t_\mathrm{eq}) $. This leads to

which is just Eq. (35). Similarly, by introducing $ l_{1,2} = l + \Gamma G\mu (t-t_{1,2}) $, we obtain

Moreover, the number densities of CS loops produced during the EMD era and surviving into the second RD and final MD eras are given by

We now analyze the SGWB spectrum generated by cosmic strings, incorporating modifications to the scale factor $ a(t) $ and the CS loop number density $ n_\mathrm{CS}(l,t) $ induced by the EMD era. Based on Eqs. (7) and (8), the GW spectra for $ G\mu = 10^{-11} $ are presented in Figs. 7(a) and 7(b) for BP1a and BP2a, respectively. The contributions from CS loops formed during different eras are demonstrated separately. Note that, in earlier eras, the CS loops are smaller in length and generate GWs at higher frequencies.

Figure 7.  (color online) GW spectra from cosmic strings with $ G\mu = 10^{-11} $ modified by the EMD era for BP1a in Scenario 1 (a) and BP2a in Scenario 2 (b). The red solid lines represent the total GW spectra, which include contributions from CS loops across all eras, while the blue, green, olive, and cyan dotted lines correspond to the contributions from CS loops formed during the original RD, EMD, second RD, and final MD eras, respectively. For comparison, various constraints and sensitivity curves from several GW detection experiments are also shown.

In order to assess the sensitivity of GW detection experiments, we also show the constraints from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) [96], the European Pulsar Timing Array (EPTA) [97], and the Parkes Pulsar Timing Array (PPTA) [98], as well as the sensitivity curves of the International Pulsar Timing Array (IPTA) [99], the Square Kilometer Array (SKA) [100], the Laser Interferometer Space Antenna (LISA) [52], TianQin [101], Taiji [102], LIGO [103], the Cosmic Explorer (CE) [104], and the Einstein Telescope (ET) [105].

Compared to the GW spectrum for $ G\mu = 10^{-11} $ in the $ \Lambda\mathrm{CDM} $ model illustrated in Fig. 1, the spectra in both BPs with the EMD era display a suppression at high frequencies $ f \gtrsim 10\; \text{Hz} $, which corresponds to the contributions from CS loops formed during the original RD and EMD eras. The primary reason for this suppression is as follows. Since the scaling behavior of the CS network is only slightly violated, the correlation length L remains approximately proportional to the cosmic time t, suggesting that the lengths of the generated loops are positively correlated with the scale factor $ a(t) $. As shown in Fig. 4, the EMD era reduces the scale factor before $ t_2 $. This means that CS loops with a given initial length l, which is inversely proportional to the frequencies of the emitted GWs through $ f_\mathrm{e} = 2n/l $, are formed at a later time, when the energy densities of both CS loops and emitted GWs are reduced. As a result, the GW spectrum at sufficiently high frequencies is suppressed, and the end time of the EMD era, $ t_2 $, is related to the GW frequency at which this suppression effect becomes significant. For the parameters listed in Table 1, the dilutor $ N_2 $ in BP1a has a longer lifetime than the dilutor Δ in BP2a, resulting in a longer EMD era and, consequently, a stronger suppression.

We proceed to investigate how the GW spectrum changes with different parameters. In Fig. 8, the GW spectra for BP1a and BP2a are displayed with the CS tension parameter varying as $ G\mu = 10^{-9} $, $ 10^{-11} $, $ 3\times 10^{-13} $, $ 10^{-15} $, and $ 5\times 10^{-18} $. For a smaller CS tension, the GW emission power is reduced, and the lifetimes of CS loops are extended. This suggests that the loops existing at $ t_2 $ originated from earlier times with smaller initial lengths, reducing the average loop length at $ t_2 $ and increasing the frequencies of the emitted GWs. Consequently, the suppression effect caused by the EMD era begins at a higher GW frequency, as shown in Fig. 8.

Figure 8.  (color online) GW spectra modified by the EMD era for BP1a in Scenario 1 (a) and BP2a in Scenario 2 (b) assuming different values of the CS tension parameter $ G\mu $.

In addition, we consider the results for BP1b and BP1c in Scenario 1 and for BP2b and BP2c in Scenario 2, whose parameters as tabulated in Tables 1 and 2, respectively. The obtained GW spectra are illustrated in Fig. 9. For BP1a, BP1b, and BP1c in Scenario 1, the dilutor $ N_2 $ has the same mass, but the masses of the DM candidate $ N_1 $ differ, causing variations in the decay width of $ N_2 $. A smaller decay width $ \Gamma_{N_2} $ corresponds to a longer duration of the EMD era, leading to stronger suppression effects at high frequencies, as illustrated in Fig. 9(a). On the other hand, the masses of the dilutor Δ are different for BP2a, BP2b, and BP2c in Scenario 2. A heavier Δ implies that the EMD era occurs earlier, and hence a higher frequency at which the suppression of the GW spectrum commences. This is clearly demonstrated in Fig. 9(b).

Figure 9.  (color online) GW spectra with $ G\mu=1\times 10^{-11} $ modified by the EMD era for various BPs in Scenarios 1 (a) and 2 (b). The blue solid lines represents the corresponding GW spectrum in the standard $ \Lambda\mathrm{CDM} $ model without the EMD era.

Furthermore, we assess the influence of the suppression effect caused by the EMD era on the experimental sensitivity of ground-based interferometers. Given a GW spectrum $ \Omega_\mathrm{noise}(f) $ converted from the strain noise of a interferometer system, along with a signal spectrum $ \Omega_\mathrm{signal}(f) $, the signal-to-noise ratio (SNR) for a practical observational time $ t_\mathrm{obs} $ can be evaluated as [106, 107]

where $ [f_\mathrm{min}, f_\mathrm{max}] $ is the accessible frequency band. If the SNR reaches a threshold $ \varrho_\mathrm{thr} = 10 $, it is probable that the GW signal will be detected.

Using the sensitivity curves of LIGO, CE, and ET, we compute the SNRs for $ t_\mathrm{obs}=1\; \text{yr} $ as functions of the CS tension parameter $ G\mu $. The results with the EMD era for BP1a in Scenario 1 and BP2a in Scenario 2 are presented by solid lines in Figs. 10(a) and 10(b), respectively, while those for the $ \Lambda\mathrm{CDM} $ model are indicated by dotted lines. In addition, the upper limits $ G\mu < 5.3\times 10^{-11} $ and $ G\mu < 2.88\times 10^{-11} $ at 95% confidence level given by NANOGrav [96] and PPTA [108] searches for a SGWB generated by cosmic strings are indicated by the regions shaded in wheat and magenta colors. For the BPs we choose, the EMD era does not influence the GW spectrum in the $ \sim \text{nHz} $ frequency band, where the PTA experiments are sensitive, and thus these upper limits remain valid.

Figure 10.  (color online) Estimated SNRs of the LIGO (blue solid lines), CE (gold solid lines), and ET (cyan solid lines) experiments affected by the EMD era for BP1a in Scenario 1 (a) and BP2a in Scenario 2 (b). For comparison, the corresponding SNRs for the $ \Lambda\mathrm{CDM} $ models are represented by dotted lines in similar colors.

For $ G\mu \lesssim 10^{-12} $, we find that the inclusion of the EMD era has insignificant impact on the SNRs, because the suppression effect occurs at frequencies outside the sensitive bands of LIGO, CE, and ET (cf. Fig. 8). In contrast, for $ G\mu \gtrsim 10^{-12} $, the SNRs are obviously reduced. Compared to BP2a, the dilutor in BP1a has a longer lifetime, leading to a stronger suppression effect and, consequently, smaller SNRs.

In this work, we explore the influence of an EMD era on the dynamics of a preexisting CS network and the resulting SGWB. We first review the GWs generated by CS loops, highlighting the explicit dependence of the SGWB spectrum on the CS loop number density $ n_\mathrm{CS} $ and the scale factor $ a(t) $. We then analyze the evolution of the CS network based on the VOS model and discuss the scaling behavior. By assuming appropriate forms for the loop production functions, we derive expressions for the CS loop number densities in RD and MD eras, which are consistent with numerical simulations in the scaling regime. Notably, these expressions are expected to remain valid even in cases where the scaling behavior is violated.

Next, we consider a cosmic timeline that incorporates an EMD era within the conventional RD era. As an illustrative example, this EMD era is attributed to a massive, long-lived dilutor in the DM dilution mechanism, particularly within the framework of the minimal LRSM, where the dilutor is either $ N_2 $ in Scenario 1 or Δ in Scenario 2. In this context, the overproduction of the DM particle $ N_1 $ is resolved by entropy injection from the dilutor decays. By solving the Boltzmann equations, we obtain the evolution of the number densities for the dilutor, DM, and SM particles, as well as the time dependence of the scale factor. These solutions allow us to identify the onset and conclusion of the EMD era. Compared to the standard $ \Lambda\mathrm{CDM} $ cosmological model, the scale factor before the end time $ t_2 $ of the EMD era is smaller.

Furthermore, we show that the evolution of the normalized correlation length ξ and the RMS velocity v of the CS network is modified by the presence of the EMD era, exhibiting a nonscaling behavior. Subsequently, we calculate the number density of CS loops formed during the EMD era, as well as the number densities of loops generated in earlier eras and surviving into later eras. Using these results, the influence of the EMD era on the SGWB spectrum arising from CS loops is demonstrated.

We find that the EMD era induces a suppression in the SGWB spectrum at sufficiently high frequencies, which is linked to the reduction of the scale factor prior to $ t_2 $. Additionally, for a smaller CS tension, the average length of CS loops at $ t_2 $ is smaller, causing the suppression effect to begin at a higher frequency. Moreover, a smaller decay width of the dilutor implies a longer duration of the EMD era, leading to a stronger suppression effect. Furthermore, a heavier dilutor causes an EMD era to occur earlier, resulting in suppression at higher frequencies.

Finally, we estimate the SNRs of the ground-based interferometers LIGO, CE, and ET as functions of $ G\mu $. The results show that for BP1a and BP1b with $ G\mu \gtrsim 10^{-12} $, the SNRs are significantly reduced, as the suppression effect falls within the sensitive frequency bands of these experiments. This study highlights that changes in the cosmic history, such as the presence of an EMD era, can affect the dynamics of the CS network and produce observable signatures in the resulting SGWB spectrum.

The authors acknowledge Ye-Ling Zhou for helpful discussions.