Gravitational Wave Signals From Early Matter Domination: Interpolating Between Fast and Slow Transitions

The recent detection of gravitational wave (GW) signals has opened a new epoch in observational astronomy. While current observations have mostly expanded our knowledge of black hole and neutron star populations [1, 2, 3, 4], future experiments are expected to also probe early universe cosmology [5, 6, 7, 8], including early matter-dominated (MD) epochs [9, 10, 11, 12, 13].

Primordial curvature perturbations, such as those produced by inflation, produce GWs at second order in perturbation theory [11]. The signal strength of these perturbations depends on the cosmological history of the universe, particularly through the gravitational potential. During an MD epoch, the gravitational potential does not decay, even on subhorizon scales, in contrast to a radiation-dominated (RD) epoch. Consequently, the metric perturbations sourced at second order by the primordial curvature perturbations are enhanced as compared to the RD scenario.

Such an early MD epoch can occur if the energy density of the universe is dominated by either heavy quasi-stable objects or a rapidly oscillating, but decaying, scalar field. Heavy objects can be motivated in extensions of the Standard Model (for example, moduli fields [14, 15, 16, 17]), as part of a dark sector (for example, [18, 19, 20]), or as extended objects, including primordial black holes [21, 22] and non-topological solitons [23, 24]. Currently, there are no constraints on the existence or duration of such an epoch, provided that it ended before nucleosynthesis [25].

Since an early MD epoch can arise in such a varied array of models [26, 27, 28, 29, 30, 24, 31, 32, 33, 34, 22, 35, 23, 36, 37, 38, 39, 40, 41, 21] (for a review, see ref. [42]), detailed studies of the associated GW signal are necessary. In particular, scenarios may differ in how rapidly the early MD epoch ends. Two limits have been well-studied in the literature, that of a slow transition (as compared to the Hubble time) [13] and that of an instantaneous transition [12]. The resulting GW signals are quite different: in the fast transition limit, there is a resonant-like enhancement (known at the poltergeist mechanism [21]); not only is this absent in the slow transition, but it is in fact replaced with a suppression. This naturally leads to the question of how to interpolate between these two regimes. This has been explored in the literature by smoothing the resonant enhancement with a Gaussian distribution of primordial black holes or Q-balls [21, 24], but a preferable approach would be to fully account for the finite duration of the reheating transition. We present such an analysis in this work.

Our paper is arranged as follows. In section 2 we outline our approach, beginning with a brief review of the formalism for calculating the induced GWs. We continue by outlining how we achieve an end to the early MD era that interpolates between the gradual and sudden regimes using a time dependent decay rate. We conclude the section by explaining our treatment of the scale factor, Green’s function and perturbations, highlighting the need for numerical calculations. In section 3 we present our calculation of the GW spectrum for varying transition duration, highlighting the disappearance of the resonant-like peak as the transition becomes more gradual as well as demonstrating the unphysical nature of oscillations in the spectrum seen in previous analysis of the gradual case [13]. In section 4 we present our conclusions. We also give some technical details of the numerical methods we used to achieve accurate results in appendix B.

In this section we outline the procedure for computing the power spectrum of GWs induced from curvature perturbations. We then explain how we model a period of early matter domination that ends in a variable time frame via a parameterised, time dependent decay rate.

We work in the conformal Newtonian gauge with the perturbed metric

$$ds^{2}=a^{2}\left[-(1+2\Phi)d\eta^{2}+\left((1-2\Psi)\delta_{ij}+\frac{h_{ij}}{2}\right)dx^{i}dx^{j}\right]$$

(2.1)

where $a$ is the scale factor of the background FLRW metric, the conformal time $\eta=\int dt/a(t)$, $\Phi$ and $\Psi$ are the first order scalar perturbations and $h_{ij}$ are the tensor perturbations induced at second order. We neglect anisotropic stress, which implies $\Psi=\Phi$, which we refer to as the gravitational potential. Throughout the paper we will use primes to denote differentiation with respect to conformal time $\eta$.

Assuming a Gaussian distribution for the curvature perturbations, the power spectrum of GWs arising from scalar perturbations is given by [11, 43, 44, 45]

$$\overline{\mathcal{P}_{h}(\eta,k)}=4\int_{0}^{\infty}dv\int_{|1-v|}^{1+v}du\left(\frac{4v^{2}-(1+v^{2}-u^{2})^{2}}{4vu}\right)^{2}\overline{I^{2}(u,v,k,\eta)}\mathcal{P}_{\zeta}(uk)\mathcal{P}_{\zeta}(vk)$$

(2.2)

where the overline denotes an oscillation average over time. The power spectrum of curvature perturbations is given by

$${\mathcal{P}}_{\zeta}(k)=\Theta(k_{\rm max}-k)A_{s}\left(\frac{k}{k_{\ast}}\right)^{n_{s}-1}$$

(2.3)

with $A_{s}=2.1\times 10^{-9}$ being the amplitude at the pivot scale, $n_{s}=0.97$ the spectral tilt, and $k_{\ast}=0.05{\rm Mpc}^{-1}$ the pivot scale, where we take all of these values from ref. [46]. We will also at times consider the scale invariant case with $n_{s}=1$. The cutoff $k_{\rm max}$ should be taken to be the wavenumber of the mode that re-enters at the start of the MD era or the non-linear scale, whichever is smaller. The application of linear perturbation theory is valid for $k_{\rm max}\leq 470/\eta_{R}$ [9, 13] where $\eta_{R}$ is the time at which the universe transitions from the MD to RD era. (The specific definition of $\eta_{R}$ in the case of a gradual transition will be clarified below; we note that the time when the energy densities are equal may generally be different.) To remain consistent with previous literature [13, 12] we conservatively take $k_{\rm max}=450/\eta_{R}$. Further, the time dependence of the GWs is encapsulated in

$$I(u,v,k,\eta)=k\int_{0}^{\eta}d\bar{\eta}\;\frac{a(\bar{\eta})}{a(\eta)}kG_{k}(\eta,\bar{\eta})f(u,v,k,\bar{\eta})$$

(2.4)

where the Green’s function $G_{k}$ is the solution to the GW equation of motion

$$G^{\prime\prime}_{k}(\eta,\bar{\eta})+\left(k^{2}-\frac{a^{\prime\prime}(\eta)}{a(\eta)}\right)G_{k}(\eta,\bar{\eta})=\delta(\eta-\bar{\eta}).$$

(2.5)

The source function is given by the expression

	$\displaystyle f(u,v,k,\bar{\eta})=\frac{3}{25(1+w)}$	$\displaystyle\Bigl{(}2(5+3w)\Phi(uk,\bar{\eta})\Phi(vk,\bar{\eta})$
		$\displaystyle+4\mathcal{H}^{-1}(\Phi^{\prime}(uk,\bar{\eta})\Phi(vk,\bar{\eta})+\Phi(uk,\bar{\eta})\Phi^{\prime}(vk,\bar{\eta}))$
		$\displaystyle+4\mathcal{H}^{-2}\Phi^{\prime}(uk,\bar{\eta})\Phi^{\prime}(vk,\bar{\eta})\Bigr{)}$		(2.6)

where the equation of state is denoted by $w$, $\Phi$ is the transfer function of the gravitational potential, normalised so that $\Phi(x\rightarrow 0)=1$, and ${\cal H}=a^{\prime}/a$ is the conformal Hubble parameter.¹¹1We note that this convention is consistent with eq. (2.6). Other references, such as ref. [11], normalise $\Phi$ differently but also correspondingly alter the $w$ dependence in the source (the only place where $\Phi$ enters the analysis). Finally, we note that if the initial condition is set in the RD epoch that precedes the early MD epoch, then $\Phi_{0}=10/9$, so that $\Phi=1$ at the start of the MD era; this is done in e.g., ref. [21]. Finally, from the power spectrum, we can derive the GW abundance

	$\displaystyle\Omega_{\rm GW}(\eta,k)$	$\displaystyle=\frac{\rho_{\rm GW}(\eta,k)}{\rho_{\rm tot}(\eta)}$
		$\displaystyle=\frac{1}{24}\left(\frac{k}{\mathcal{H}(\eta)}\right)^{2}\overline{\mathcal{P}_{h}(\eta,k)}.$		(2.7)

Because previous analyses, whether of a rapid or slow transition, used piecewise approximations about $\eta_{R}$ for the scale factor, Green’s function and gravitational potential [13, 12, 21], they decomposed the GW signal into three parts: one sourced from the RD era, one sourced from the MD era, and a cross term. In the rapid scenario, the spectra is dominated by contributions sourced in the RD era, as the gravitational potential is matched at the reheating time with no decay. In the case of a gradual transition, the contributions from the RD and MD eras are very similar, so the cancellation due to the cross term is severe. We are interested in intermediate cases between the two extremes. Although in our analysis we will avoid any piecewise approximations by numerically evolving the scale factor, Green’s function and gravitational potential, we will still decompose our results into these components in Section 3 when we compare our results to previous work.

To model a reheating transition in the intermediate regime between gradual and rapid, we propose a time dependent decay rate with a parameter $\beta$ that controls the speed at which the matter evaporates. Specifically we choose a decay rate of the form,

$$\Gamma(\eta)=\Gamma_{max}(\tanh(\beta(\eta-\eta_{\Gamma}))+1)/2.$$

(2.8)

While this decay rate is centred at the time $\eta_{\Gamma}$, the reheating transition will commence at time $\eta_{R}$ when the decay becomes efficient compared to the Hubble rate $\Gamma\sim\mathcal{H}$. Consequently, to have a rapid decay with $\Gamma\gg\mathcal{H}$, it is necessary to have a time-evolving decay rate which rapidly increases.

In more detail, we note that $\eta_{R}$ depends not only on $\eta_{\Gamma}$ but also $\beta$ and $\Gamma_{\rm max}$. For small enough $\beta$, the transition will occur during the rise of $\tanh$, in a region where the decay rate is not changing rapidly, resulting in a gradual transition. Therefore, for fixed $\eta_{\Gamma}$, decreasing $\beta$ shifts the transition to earlier times. The extreme limit of this is $\beta=0/\eta_{R}$ for which the decay rate becomes constant. In the other limit, as $\beta\rightarrow\infty$, the decay rate approaches a step function which, given a large enough $\Gamma_{\rm max}$ such that $\Gamma(\eta_{\Gamma})\gg\mathcal{H}$, describes a rapid transition where all of the matter evaporates approximately instantaneously at $\eta_{\Gamma}$.

The matter energy density $\rho_{m}$, the radiation energy density $\rho_{r}$ and the scale factor $a$ will evolve according to the coupled Friedmann and continuity equations with our variable decay rate inserted as follows [47]

	$\displaystyle\rho_{m}^{\prime}$	$\displaystyle=-(3\mathcal{H}+a\Gamma(\eta))\rho_{m}$		(2.9)
	$\displaystyle\rho_{r}^{\prime}$	$\displaystyle=-4\mathcal{H}\rho_{r}+a\Gamma(\eta)\rho_{m}$		(2.10)
	$\displaystyle a^{\prime}$	$\displaystyle=\sqrt{\frac{8\pi}{3M_{\rm Pl}^{2}}(\rho_{m}+\rho_{r})}a^{2}$		(2.11)

where $M_{\rm Pl}=1/\sqrt{G}$ is the Planck mass.

In the solid black curve of figure 1 we show the scale factor obtained from the Friedmann equations for a gradual transition with constant decay rate. We see that it continuously evolves between the MD and RD regimes. In contrast, we introduce an analytic fit $a_{f}$, valid for an exactly instantaneous transition

$$\frac{a_{f}(\eta)}{a_{f}(\eta_{R})}=\begin{cases}\left(\tfrac{\eta}{\eta_{R}}\right)^{2}&\eta<\eta_{R}\\ 2\tfrac{\eta}{\eta_{R}}-1&\eta\geq\eta_{R}\end{cases}$$

(2.12)

which is constructed so that $a_{f}$ and its first derivative are continuous at $\eta_{R}$. Fits like this one have often been used in studies of GWs at second order, even in the case of gradual transitions (e.g. in ref. [13, 12, 21]). The fit has two free parameters $a_{f}(\eta_{R})$ and $\eta_{R}$ which are to be determined by fitting onto two points of our numerical solution. Note that $a_{f}(\eta_{R})$ need not be the same as $a(\eta_{R})$. In fact enforcing this would make one of the points we fit to be in the middle of the transition which we find results in a poorer overall fit.

Since a definite transition time $\eta_{R}$ is a property of an instantaneous transition, when fitting to our smooth numerical calculation, different prescriptions for choosing $\eta_{R}$ can impact the error in the fit. The two different methods we use are as follows.

In the red dashed curve of figure 1 we compute $\eta_{R}$ and $a_{f}(\eta_{R})$ by fitting to one point early in the MD era and another late in the RD era. This captures the asymptotic behaviour in both eras, and as a result the calculation of $\eta_{R}$ is a good estimate for the time at which the universe transitions from being MD to RD. Henceforth, when we refer to $\eta_{R}$ we use this method to compute it. However, due to the gradual nature of the transition we see that this approach has as much as an eight percent error during the transition and the fit is poorer during the RD era.

Alternatively we perform the fit using two points in the RD era, where the fit is purely linear. Specifically the first point is taken when $\rho_{r}=10^{4}\rho_{m}$ and the second much later, at the point where we stop integrating. This is shown in the dotted blue curve of figure 1. We see that the fit is very good during the RD era, with the actual size of the error being of order $10^{-7}$, while being more than 10% during the MD era. This fit will be useful for calculating an analytic expression for the perturbation $\Phi$ during the RD era. From here on we use $\eta_{\rm RD}$ to indicate the value of $\eta_{R}$ in eq. (2.12) resulting from this approach.

We highlight that while we make some use of the fits, we do so in a way such that the errors introduced by their use are negligible. The fit in the red dashed curved of figure 1 is only used to estimate the time of the transition $\eta_{R}$. This only enters the calculation when defining the non-linear cut-off scale $k_{\rm max}$ and when breaking the result down into contributions sourced from the MD era and those from the RD era. The other fit is only used for the analytic expression for the gravitational potential late in the RD era after its relative error has dropped below $10^{-7}$. During the MD era and the transition itself, when the error in the fits is largest, our calculation is fully numeric. We note that this discussion has all been in the context of the gradual transition limit for which these issues are the most pronounced. As the transition becomes more rapid, eq. (2.12) agrees better with the numerical solution.

Calculating the signal also requires the Green’s function. Equation (2.5) can be solved analytically in the MD and RD regimes, and it is common in the literature to define a piecewise Green’s function, changing between the two expressions at $\eta_{R}$. This is done even in the analysis of gradual transitions, e.g. ref. [13]. We improve on this by instead using a numerical solution for the Green’s function which accurately captures the behaviour during the transition. As we discuss below, we found that this eliminated unphysical oscillations (in frequency) in the resulting signal, which are present in the results of ref. [13]. In our approach, we numerically evolve eq. (2.5) in $\eta$, from $\bar{\eta}$, with initial conditions

$$G_{k}(\eta\leq\bar{\eta})=0,\;\;G_{k}^{\prime}(\bar{\eta},\bar{\eta})=1,$$

(2.13)

since we are looking for the retarded Green’s function. At sufficiently late times, we match onto the solution

$$G_{k}(\eta_{,}\bar{\eta})=E(\bar{\eta})\cos{(k\eta)}+F(\bar{\eta})\sin{(k\eta)}$$

(2.14)

since in the RD limit $a^{\prime\prime}$ vanishes. (Note that we always take $\eta\gg\eta_{R}$; that is, we assume the GWs are observed significantly after the MD epoch has ended.) This decomposition is primarily useful for analytically performing the oscillation average in eq. (2.2). We perform the matching at the time $\eta_{G}$ defined as

$$k^{2}=10^{-12}\frac{a^{\prime\prime}(\eta_{G})}{a(\eta_{G})},$$

(2.15)

which is motivated by the observation from eq. (2.5) that the RD limit is valid when $k^{2}\gg a^{\prime\prime}(\eta)/a(\eta)$. Thus, when $\bar{\eta}>\eta_{G}$ we simply use the RD Green’s function

$$G_{k}^{\rm RD}(\eta,\bar{\eta})=\frac{1}{k}\sin{(k\eta-k\bar{\eta})}.$$

(2.16)

Our results below indicate that the abrupt change in the piecewise definitions of the Green’s function produces unphysical oscillations in the signal, and we note that this matching necessarily introduces similar abrupt changes. We have defined $\eta_{G}$ in eq. (2.15) so that the resulting oscillations induced in the spectra are negligible.

For each value of $k$ we sample and interpolate $G$ over $\bar{\eta}$ densely enough so that the integral in eq. (2.4) can be performed accurately. Since the Green’s function depends only on $k$, and not $u$ or $v$, we can reuse the calculation for each point sampled in the integrals of eq. (2.2).

The predominant computational complexity in calculating the induced GWs comes from evolving the gravitational potential. This involves solving a coupled differential equation for the gravitational potential, energy density perturbations $\delta_{m/r}$, and velocity divergences $\theta_{m/r}$ of matter and radiation respectively. To allow for as direct a comparison as possible with references [13, 12] we choose to evolve the perturbations in the Newtonian gauge. This causes complications for a time dependent decay rate, as the presence of $\Phi$ in (2.1) causes time to not be synchronised with the matter. Reference [21] deals with this by solving for the perturbations in the synchronous gauge, with coordinates comoving with the matter, before transforming to the Newtonian gauge to calculate the induced gravitational wave spectrum. We instead transform the synchronous gauge equations in [47] to the Newtonian gauge while including the time dependence of the decay rate. Neglecting anisotropic stress of radiation, the perturbations for a given wavenumber $k$ evolve as

	$\displaystyle\delta_{m}^{\prime}$	$\displaystyle=-\theta_{m}+3\Phi^{\prime}-a\Gamma(\eta)\Phi-a\Gamma^{\prime}(\eta)\frac{\theta_{m}}{k^{2}}$		(2.17)
	$\displaystyle\theta_{m}^{\prime}$	$\displaystyle=-\mathcal{H}\theta_{m}+k^{2}\Phi$		(2.18)
	$\displaystyle\delta_{r}^{\prime}$	$\displaystyle=-\frac{4}{3}(\theta_{r}-3\Phi^{\prime})+a\Gamma(\eta)\frac{\rho_{m}}{\rho_{r}}(\delta_{m}-\delta_{r}+\Phi)+a\Gamma^{\prime}(\eta)\frac{\rho_{m}}{\rho_{r}}\frac{\theta_{m}}{k^{2}}$		(2.19)
	$\displaystyle\theta_{r}^{\prime}$	$\displaystyle=\frac{k^{2}}{4}\delta_{r}+k^{2}\Phi-a\Gamma(\eta)\frac{3\rho_{m}}{4\rho_{r}}\left(\frac{4}{3}\theta_{r}-\theta_{m}\right)$		(2.20)
	$\displaystyle\Phi^{\prime}$	$\displaystyle=-\frac{k^{2}\Phi+3\mathcal{H}^{2}\Phi+\frac{3}{2}\mathcal{H}^{2}\left(\displaystyle{\frac{\rho_{m}\delta_{m}+\rho_{r}\delta_{r}}{\rho_{\rm tot}}}\right)}{3\mathcal{H}}$		(2.21)

where the additional $\Gamma^{\prime}$ terms arise in the transition to Newtonian gauge, due to the non-trivial time-dependence of $\Gamma$. These terms were not included in [12], although as noted above a time-dependent decay rate is necessary to have a transition with $\Gamma\gg\mathcal{H}$. The derivation of these equations may be found in appendix A. We evolve these equations forward in time from the initial conditions

$$\delta_{m,0}=-2\Phi_{0},\;\;\delta_{r,0}=\frac{4}{3}\delta_{m,0},\;\;\theta_{m,0}=\theta_{r,0}=\frac{k^{2}\eta}{3}\Phi_{0}.$$

(2.22)

By taking $\Phi_{0}=1$ we effectively solve for the transfer function.

We compare solving this system of equations numerically to the approximations used in the literature. First, in the gradual case [13], the oscillatory period was neglected as the gravitational potential decayed significantly. In place of a numerical solution, the following fit, which is $k$ independent, was used

$$\Phi_{\rm fit}=\begin{cases}\exp\left(-\tfrac{2}{3}\left(\tfrac{\eta}{\eta_{R}}\right)^{3}\right)&\eta<\eta_{R}\\ \exp\left(-2\left(\left(\tfrac{\eta}{\eta_{R}}\right)^{2}-\tfrac{\eta}{\eta_{R}}+\tfrac{1}{3}\right)\right)&\eta\geq\eta_{R}.\end{cases}$$

(2.23)

Because we include the oscillating tail, our analysis is not $k$ independent and must be solved for each wave number, which increases the computational time. This is especially pronounced because if the $k$ dependence is neglected, then $f$ and $I$ become independent of $u$ and $v$. Consequently, $I$ factors out of the integral in eq. (2.2).

For the sudden case [12], conversely, the gravitational potential experiences no decay, instantaneously switching from a constant solution during the MD era to an oscillating solution with large amplitude at $\eta_{R}$. In this case the spectrum is almost entirely sourced by the RD component and we will see below that an analytic solution exists both for $\Phi$ and $I$ during the RD era. This means that their analysis completely avoids the expensive step of evolving the perturbations and the remaining integral for the power spectrum eq. (2.2) depends only on analytic functions of $u$ and $v$.

Our approach, as we describe below, is to numerically solve for the gravitational potential $\Phi$ before matching onto a late time analytic solution. In figure 2 we show the time dependence of the gravitational potential, calculated numerically, for wavenumber $k=100/\eta_{R}$ in the case of a rapid transition in blue and a gradual transition in dashed orange. As expected, during the MD era $\Phi$ remains constant before experiencing some period of decay and then beginning to oscillate. The amount of decay depends on how rapid the transition is, as we see that in the rapid limit the potential experiences no decay, whereas in the gradual limit it decays by roughly ten orders of magnitude. In both cases, the potential continues to oscillate and decays at a slower rate during the RD era.

The oscillations during the RD era complicate the numerical integration over $\bar{\eta}$ in eq. (2.4). These oscillations are crucial in the rapid scenario, as they give rise to the resonant poltergeist peak [12]. To capture this behaviour we numerically evolve the perturbations during the transition, but once $\delta_{r}$ comes to dominate the last term in eq. (2.21) we match onto an analytic solution. We define the time $\eta_{\rm osc}$ at which we match onto the analytic solution to be when

$$\frac{\rho_{m}\delta_{m}}{\rho_{r}\delta_{r}}=10^{-4}$$

(2.24)

and as such it will depend both on $k$ and the suddenness of the transition $\beta$.

We now describe the analytic solution that we match onto. It makes use of the fact that, away from the transition, the evolution of $\Phi$ in a single component fluid is described by [48]

$$\Phi^{\prime\prime}+3(1+w)\mathcal{H}\Phi^{\prime}+wk^{2}\Phi=0,$$

(2.25)

where, in this case, $\mathcal{H}$ is evaluated from the analytic approximation in eq. (2.12), as we are away from the transition. During the RD era this has the solution

$$\Phi_{A}(k,\eta)=A(k,\eta_{\rm osc})\mathcal{J}(k\eta)+B(k,\eta_{\rm osc})\mathcal{Y}(k\eta)$$

(2.26)

where $\mathcal{J}$ and $\mathcal{Y}$ are defined in terms of the spherical Bessel functions, $j_{1}(x)$ and $y_{1}(x)$,

	$\displaystyle\mathcal{J}(k\eta)=\frac{3\sqrt{3}\;j_{1}\left(\frac{k\eta-k\eta_{\rm RD}/2}{\sqrt{3}}\right)}{k\eta-k\eta_{\rm RD}/2}$		(2.27)
	$\displaystyle\mathcal{Y}(k\eta)=\frac{3\sqrt{3}\;y_{1}\left(\frac{k\eta-k\eta_{\rm RD}/2}{\sqrt{3}}\right)}{k\eta-k\eta_{\rm RD}/2},$		(2.28)

where we have defined $\eta_{RD}$ above. The coefficients $A(k,\eta_{\rm osc})$ and $B(k,\eta_{\rm osc})$ are given by imposing continuity of the solution and its derivative at $\eta_{\rm osc}$.

We see that during the RD era the amplitude of the potential decays as

$$\Phi_{\rm amp}=\frac{9\sqrt{A(k,\eta_{\rm osc})^{2}+B(k,\eta_{\rm osc})^{2}}}{(k(\eta-\eta_{\rm RD}/2))^{2}}$$

(2.29)

which we use to normalise the error in the analytic solution as $|\Phi-\Phi_{A}|/\Phi_{\rm amp}$, which is shown in the lower panel of figure 2. We only show the error for $\eta>\eta_{\rm osc}$ when we actually use the analytic solution, where $\eta_{\rm osc}\approx 4.4\eta_{R}$ in the gradual case and $\eta_{\rm osc}\approx\eta_{R}$ in the rapid. While we enforce a relative error of $10^{-11}$ in the ODE solver, we find that the error in the approximation is significantly larger than this, for the gradual case being of order $10^{-4}$. This is because we neglect the term $\rho_{m}\delta_{m}$ in equation (2.21) which is of order $10^{-4}$ times the dominant term $\rho_{r}\delta_{r}$ at $\eta_{\rm osc}$ by definition. The error is much better in the rapid case as the $\rho_{m}\delta_{m}$ term simply decays away faster due to the sharp turn on of the decay rate.

There is a balance between setting $\eta_{\rm osc}$ late enough that the $\rho_{m}\delta_{m}$ term is truly negligible but also early enough to avoid spending too much time accurately integrating the ODE for an oscillating solution. Since the main purpose for transitioning to an analytic expression is evaluating the integral for $I$ in eq. (2.4), the error in the amplitude translates directly to the error in $I$. The reduced accuracy in the gradual case is offset by the more than ten orders of magnitude suppression by the time oscillations begin, making their contribution to the integral negligible.

The reasoning for using the analytic expression at late times is two-fold. Firstly, as discussed above, finding an accurate numerical solution for $\Phi$ during the oscillations is computationally expensive. Utilising the analytic expression takes roughly a quarter of the time it would take to numerically integrate the ODE over the time ranges shown in figure 2. Secondly, due to the oscillations, it would be difficult to perform the integral for $I$, in eq. (2.4), numerically. The analytic expression $\Phi_{A}$ admits an analytic integral for $I$. We evaluated this integral (using Mathematica) and, after transforming the integration variable to $\bar{x}=k(\bar{\eta}-\eta_{\rm RD}/2)$ we found our result is consistent with that of appendix A of ref. [11] making the substitutions

	$\displaystyle A(u)$	$\displaystyle=\frac{9}{10}A(uk,\eta_{\rm osc})$		(2.30)
	$\displaystyle B(u)$	$\displaystyle=\frac{9}{10}B(uk,\eta_{\rm osc})$		(2.31)
	$\displaystyle C$	$\displaystyle=-\cos\left(k(\bar{\eta}-\eta_{\rm RD}/2)\right)$		(2.32)
	$\displaystyle D$	$\displaystyle=\sin\left(k(\bar{\eta}-\eta_{\rm RD}/2)\right)$		(2.33)
	$\displaystyle x_{1}$	$\displaystyle=k(\eta_{\rm osc}-\eta_{\rm RD}/2)$		(2.34)

and taking the limit $x_{2}\rightarrow\infty$. The factors of $9/10$ in $A$ and $B$ are due to ref. [11] having normalised $\Phi$ to unity during the RD era.

Finally, once we have solved for $I$, we find the GW signal by performing the integrals in eq. (2.2). The details of this integration can be found in appendix B, along with other details on our numerical techniques. We note that in this final integral we allow up to 5 percent error, as this is sufficient precision to compare with the gradual and rapid limits from previous work and make qualitative statements about the spectra and potential detectability.

It is useful for comparing back to phenomenological scenarios to derive a relationship between conformal time $\eta$ and temperature $T$ valid during the RD era [49]. From the Friedmann equations we have

$$\frac{aH}{a_{\rm eq}H_{\rm eq}}=\frac{a}{a_{\rm eq}}\sqrt{\frac{\rho}{2\rho_{\gamma,{\rm eq}}}}$$

(2.35)

where the subscript eq denotes evaluation at the late matter radiation equality time and $\rho_{\gamma}$ is the photon energy density. Near the reheating transition it is adequate to use eq. (2.12) which, during the RD era, gives $aH=(\eta-\eta_{R}/2)^{-1}$. Assuming entropy conservation $g_{s,\rm eq}a_{\rm eq}^{3}T_{\rm eq}^{3}=g_{s}a^{3}T^{3}$ where $g_{s}$ is the number of relativistic degrees of freedom of entropy. We reduce the ratio of energy densities by taking $\rho\propto gT^{4}$ where $g$ is the effective number of relativistic degrees of freedom of radiation. Taking $a_{\rm eq}H_{\rm eq}=k_{\rm eq}=0.01041\;{\rm Mpc^{-1}}$ [46, 50] gives

$$\frac{(\eta-\eta_{R}/2)^{-1}}{0.01041\;{\rm Mpc^{-1}}}=\frac{1}{\sqrt{2}}\left(\frac{g_{s,\rm eq}}{g_{s}}\right)^{1/3}\left(\frac{g}{g_{\rm eq}}\right)^{1/2}\frac{T}{T_{\rm eq}}.$$

(2.36)

We calculate $T_{\rm eq}=T_{0}(1+z_{\rm eq})=0.801\;{\rm eV}$ from the current CMB temperature $T_{0}=2.7255\;{\rm K}$ [51] and the redshift $z_{\rm eq}=3411$ [46]. Finally taking $g_{\rm eq}=3.38$ and $g_{s,{\rm eq}}=3.94$ [52] we obtain

$$\eta-\frac{\eta_{R}}{2}=58.2\;{\rm MeVpc}\left(\frac{g_{s}}{106.75}\right)^{1/3}\left(\frac{g}{106.75}\right)^{-1/2}T^{-1}.$$

(2.37)

We use this both to evaluate the reheating temperature, as well as to match our evolution of the Friedmann equations, eqs. (2.9) to (2.11), onto the standard RD era for $\eta\gg\eta_{R}$.

After the gravitational potential source decays away enough, the GW energy density parameter $\Omega_{\rm GW}$ comes to be constant during the RD era. For reference, we define the time at which $\Omega_{\rm GW}$ becomes constant to be $\eta_{c}$. The GWs will evolve further during the late MD era with the present value of the energy density parameter given by [12]

$$\Omega_{\rm GW}(\eta_{0},k)=0.39\left(\frac{g(\eta_{c})}{106.75}\right)^{-1/3}\Omega_{r,0}\Omega_{\rm GW}(\eta_{c},k)$$

(2.38)

where $\Omega_{r,0}$ is the present value of the radiation energy density parameter and $\eta_{0}$ is the current conformal time. This also takes into account the change in the effective relativistic degrees of freedom $g$ between when the GWs were sourced and now.

In this section, we will present the results of our improved calculation.

Our main results are presented in figures 3 and 4, which show our numerical calculations of the GW spectra induced from a scale invariant curvature spectrum for a variety of matter-radiation transition speeds. In figure 3, for each value of $\beta$ we we set $\eta_{\Gamma}$ such that the transition occurs at $\eta_{R}=1.164\times 10^{-9}\;{\rm Mpc}$, corresponding to a reheating temperature $T_{R}=100\;{\rm GeV}$. We also set $\Gamma_{\rm max}=500\Gamma_{R}$, where $\Gamma_{R}=1.590\times 10^{-14}\;{\rm GeV}$ is the constant decay rate that produces a transition at $\eta_{R}$, to ensure that for large $\beta$ the matter evaporates rapidly when the decay turns on. In contrast figure 4 demonstrates the same GW spectra for a different choice of the decay rate parameters. In this case we set $\Gamma_{\rm max}=7.948\times 10^{-12}\;{\rm GeV}$, the same as before, and fix $\eta_{\Gamma}=1.164\times 10^{-9}\;{\rm Mpc}$ so that the rapid limit is the same between the two plots. Now $\eta_{R}$ is no longer fixed and as $\beta$ decreases, the transition shifts to earlier times. The parameters of figures 3 and 4 are summarised in tables 1 and 2 respectively.

In addition to characterising the transition in terms of $\beta$ we can also compare the duration in regular time, $\tau$, directly to the Hubble rate at matter radiation equality, $H(\eta_{\rm eq})$. This gives us a more intuitive quantity to compare real models to, for example, in a gradual transition $\tau H\sim 1$ whereas a phase transition can complete with $\tau H\in(10^{-4},10^{0})$. We can compute $\tau$ from the conformal time

$$\tau=\int_{\eta_{i}}^{\eta_{f}}d\overline{\eta}\;a(\overline{\eta}).$$

(3.1)

We define the initial and final times $\eta_{i}$ and $\eta_{f}$ from the comoving matter energy density $\rho_{m}(\eta)(a(\eta)/a(\eta_{\rm MD}))^{3}$, where $\eta_{\rm MD}$ is the start of the MD era. The start of the transition $\eta_{i}$ is defined as when one percent of the matter has decayed and the end $\eta_{f}$ is when only a factor of $1/e$ remains.

Figures 3 and 4 are motivated by different cosmological model-building motivations: if a model has heavy matter whose decay is uncertain, then figure 4 may be more useful. On the other hand, if one is building a model by imposing a transition at a specific time, adjusting the type or amount of heavy matter, then figure 3 may be most useful. Both figures share the same qualitative features, however, demonstrating how the signal evolves between the slow and fast limits.

For comparison, the dashed, black curves in both figures show the rapid limit, constructed from the analytic approximation to the large scale (small $k/k_{\rm max}$) and resonant peak contributions to the spectrum, derived in the appendix of ref. [12]. Similarly, the solid, black curve denotes the gradual limit reconstructed from ref. [13] using their approximate fits to the scale factor, equation of state and gravitational potential. We see excellent agreement with our numerical code in the rapid limit with the nearly coincident blue and orange curves ($k=3000/\eta_{R}$ and $k=500/\eta_{R}$ respectively) perfectly matching the analytic approximation in the regions for which it is valid, both in the broad low $k/k_{\rm max}$ regime of the spectrum, as well as at the resonant peak. For the short scale part of the spectrum (large $k/k_{\rm max}$) away from the peak, they diverge from the approximate rapid limit because our result is a complete computation of the spectrum. This is consistent with the numerical results from ref. [12] which similarly diverge from the approximation in this region. We include both the orange $\beta=500/\eta_{R}$ curve and the blue $\beta=3000/\eta_{R}$ curve to demonstrate that by $\beta=500/\eta_{R}$ we have effectively converged to the rapid limit. The agreement is less precise in the gradual case, with the olive curve ($k=0$) matching the overall scaling of the dashed curve for large $k$, but differing at small $k$. We also see that our analysis does not have the oscillations at large $k$ values. Both of these will be discussed below.

A noteable difference between the curves plotted in figure 3 and figure 4 is that when $\eta_{\Gamma}$ is fixed, $\eta_{R}$ varies, and thus the time scale at which the transition occurs changes. Similarly, $k_{\rm max}$ is different for each curve and the spectrum is shifted to higher frequencies for earlier $\eta_{R}$. This dependence of $\eta_{R}$ and $k_{\rm max}$ on $\beta$ is summarised in table 2. However, since the spectrum is a function of $k/k_{\rm max}$ this is not immediately evident from figure 4.

The most important feature of these results is demonstrating the evolution of the resonant-like peak (the so-called poltergeist effect) as the matter-radiation transition time increases. First, we see that the peak does exist for finite $\beta$, and therefore it is not a numerical artefact of any of the step functions used in the analyses of ref. [12]. We see that as the duration of the transition increases, the entire spectrum becomes suppressed as would be expected. In fact we see that the red curve (corresponding to $\beta=50/\eta_{R}$) is already suppressed by roughly three orders of magnitude compared to the rapid limit. This corresponds to a transition duration of roughly $\tau H=8.64\times 10^{-2}$, demonstrating that to achieve the theoretical maximum strength signal, the transition needs to be faster than $\mathcal{O}(0.01)$ of the Hubble time.

In addition to this, the resonant peak broadens, which can be seen particularly in the $\beta=3000/\eta_{R}$ and $\beta=500/\eta_{R}$ lines, which are nearly coincident to the right of the resonant peak. The resonance can be understood as the amplification of GWs that are comoving with the sound waves produced in the plasma after the sudden transformation of matter into radiation. When the transition is spread out over a small time, the sound waves exist in a frequency band, amplifying a range of comoving GWs and broadening the peak. However, if the transition time is sufficiently long, then the overdensities dissipate instead of sourcing sound waves.

Next we discuss the gradual limit, and in particular, the differences between our analysis and ref. [13] that we mentioned above. For $\beta\lesssim 10/\eta_{R}$ (the brown curve and below) we see that for large enough $k/k_{\rm max}$ the spectra converges to the gradual limit. The significant suppression at large $k$ was previously noted in the literature [13], which has been understood as arising from the cross term in $I^{2}$ when $I$ is broken into MD and RD contributions,

	$\displaystyle I(u,v,k,\eta)$	$\displaystyle=k\int_{0}^{\eta_{R}}d\bar{\eta}\;\frac{a(\bar{\eta})}{a(\eta)}kG_{k}^{{\rm MD}}(\eta,\bar{\eta})f(u,v,k,\bar{\eta})+k\int_{\eta_{R}}^{\eta}d\bar{\eta}\;\frac{a(\bar{\eta})}{a(\eta)}kG^{\rm RD}_{k}(\eta,\bar{\eta})f(u,v,k,\bar{\eta})$
		$\displaystyle=I_{{\rm MD}}(u,v,k,\eta)+I_{\rm RD}(u,v,k,\eta).$		(3.2)

The GW spectrum itself can then be decomposed into a term from the RD era (sourced by $\overline{I_{\rm RD}^{2}}$), a term from the MD era (sourced by $\overline{I_{\rm MD}^{2}}$) and a cross term (sourced by $2\overline{I_{\rm MD}I_{\rm RD}}$), that is

$$\Omega_{\rm GW}=\Omega_{\rm RD}+\Omega_{\rm MD}+\Omega_{\rm cross}\ .$$

(3.3)

Since we numerically evolve all relevant quantities, and therefore don’t use piecewise expressions for the scale factor, Green’s function and gravitational potential, we do not need to decompose the signal into these components. However, to facilitate comparison to ref. [13], we have identified the corresponding contributions in our analysis (using $\beta=0/\eta_{R}$), and the results are shown in figure 5.

We see that for $k\gtrsim 0.03k_{\rm max}$ each contribution is comparable in magnitude with $\Omega_{\rm MD}\simeq\Omega_{\rm RD}\simeq-\Omega_{\rm cross}/2$ such that the resulting signal is many orders of magnitude smaller than any individual contribution. To understand why the cross term cancels so precisely, we note that in the large $k/k_{\rm max}$ regime, the gravitational potential undergoes significant decay before it oscillates with suppressed amplitude. Therefore, the contribution from the RD era $\Omega_{\rm RD}$ is nearly entirely sourced during the regime in which $\Phi$ is rapidly decreasing, before it even begins to oscillate. This is also the regime with the most significant contribution to the mixed term, as the matter density has not entirely disappeared yet.

This cancellation particularly demonstrates the importance of numerical precision in these calculations, and as can be seen in figures 3 and 4 our results in the large $k/k_{\rm max}$ regime are in agreement with version 3 of ref. [13]²²2We thank the authors of ref. [13] for updating their analysis after discussion..

Quite crucially though, once the power spectrum has converged to the gradual limit it does not oscillate with $k$ as it does in the results from ref. [13]. We observe that the oscillations in $k$ have a period of $2\pi/\eta_{R}$, indicating their origin is from the matching between the MD and RD eras. These oscillations arise both from discontinuities in the second derivatives of $a$ and $f$ as well as from the piecewise treatment of the Green’s function. This can be seen in figure 6, which shows that if we use our numerical solution for the scale factor in place of eq. (2.12) (but retain the piecewise approximation for the Green’s functions) the oscillations are suppressed, and then if we further replace the approximate Green’s function with our numerical solution they are absent. As such, our signals are free of the oscillatory artefact and highlight the true characteristic of the signal from a gradual transition.

As we mentioned, we also find that our numerical results differ considerably from the gradual limit for small $k/k_{\rm max}$, converging to the olive $\beta=0/\eta_{R}$ curve. This is due to the gradual limit employing a fit to $\Phi$ valid for $k\gtrsim 30/\eta_{\rm eq}$, where $\eta_{\rm eq}$ is the matter-radiation equality time of the transition in question (see eq. (3.10) and following discussion in ref. [13]). By numerically evolving $\Phi$ we do not face the same restriction and are able to account for oscillations in $\Phi$ during the MD era that become large for small $k$. This is evident in figure 5 where we see that $\Omega_{\rm RD}$ comes to dominate the spectrum for small $k$.

We also observe that for finite-duration transitions, the GW signal at large $k/k_{\rm max}$ converges to the slow transition regime faster than the GW signal at lower $k/k_{\rm max}$. This is also a consequence of the gravitational potential undergoing significant decay before it oscillates with suppressed amplitude at large $k/k_{\rm max}$, which as explained above leads to the precise cancellation between components. Physically, at these scales the matter overdensities have time to dissipate during the extended matter-radiation transition period and do not source the sound waves that resonantly enhance the GWs.

In figure 7 we show the GWs signals induced during transitions at $T_{R}=100\;{\rm GeV}$ of varying speeds by a realistic curvature spectrum as in eq. (2.3) with $A_{s}=2.1^{-9}$ and $n_{s}=0.97$ [46]. We contrast these with the sensitivities of upcoming GW experiments. We see that for a transition at $100\;{\rm GeV}$, LISA is the applicable detector and would only be sensitive to $\beta\gtrsim 100/\eta_{R}$ which correspond to transitions with duration $\tau H\gtrsim 4.490\times 10^{-2}$. Since the tilt of the curvature power spectrum is very close to one, as we move to different reheating temperatures, the peak frequency will change but the amplitude of the spectrum will remain fairly constant. Additionally changes in the effective degrees of freedom $g_{c}$ will also influence the amplitude of the spectrum, but this is only a small effect, scaling by a factor of $\mathcal{O}(1)$ at most. As such, only SKA and THEIA would be able to probe more gradual transitions with $\beta\gtrsim 50/\eta_{R}$ and $\tau H\gtrsim 8.641\times 10^{-2}$ occurring with reheating temperatures around $T_{R}\sim 1\;{\rm MeV}$. Other interferometers will perform similarly to LISA, being sensitive to transitions with duration $\tau H\gtrsim 4.490\times 10^{-2}$, with DECIGO probing reheating temperatures around $T_{R}\sim 30\;{\rm TeV}$ and the Einstein Telescope and Cosmic Explorer probing $T_{R}\sim 10^{6}\;{\rm GeV}$.

Finally, we contrast our approach with those of references [21, 24]. In these works, the finite width of the transition is included only via a normalisation factor $S$ in the gravitational potential $\Phi$, which describes the suppression of $\Phi$ before its oscillation. This suppression is determined by a numerical fit, parameterised by a Gaussian parameter $\sigma$ which describes the mass distribution of the black holes or Q-balls. Their analysis uses piecewise expressions for the scale factor, Hubble parameter, and Green’s function. In contrast, in our analysis we have replaced all three with smoothly varying expressions. Because of this, no further normalisation factor is needed in $\Phi$; its decline can be observed in the gradual line in figure 2. Because we consider a tanh decay profile to interpolate between fast and slow transitions, as opposed to a Gaussian distribution of decaying matter, it is difficult to compare our results directly; there is no direct equivalence between $\beta$ and $\sigma$ values. We note, however, that figure 6 of ref. [21] shows a resonance peak for all $\sigma$ values, whereas we clearly see the resonance disappear as the transition time increases. This is to be expected as they explicitly note that they only consider mass distributions sufficiently small that the MD and cross terms in eq. (3.3) are negligible. Therefore, this is the first work to explore the intermediate regime between the rapid and gradual transitions.

The GW signal from an early MD epoch has attracted a great deal of interest, particularly as new GW detectors are designed. Detailed calculations of finite-duration transitions between matter and radiation domination are necessary, particularly because a fast transition shows a resonant peak while a slow transition displays a very precise cancellation. In this work we have calculated the GW signal carefully for these intermediate transitions, avoiding all step-function approximations, whether in the scale factor or Green’s function. Our results show that the GW signal interpolates between the instantaneous limit and the gradual limit quickly enough that there is reason to be optimistic about the prospects of distinguishing between different sources of matter domination. Ultimately, one might be able to reconstruct three features of the scenario responsible for an equation of state: when the period of matter domination began, how long it lasted and how quickly it transitioned to radiation.

In this work, we have used eq. (2.8) to interpolate between the fast and slow transition regimes. Of course, physical scenarios, such as those with primordial black holes or Q-balls, will not have a tanh profile. We leave extending this analysis to these physically motivated scenarios for future work, as the calculation of the equivalent of eq. (2.8) is complicated due to gauge effects between the synchronous and Newtonian gauges.

MP was supported by an Australian Government Research Training Program (RTP) scholarship and a Monash Graduate Excellence scholarship (MGES). CB acknowledges support from the Australian Research Council via the Discovery project DP210101636.