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

Matthew Pearce  Lauren Pearce  Graham White  Csaba Balázs

###### Abstract

An epoch of matter domination in the early universe can enhance the primordial
stochastic gravitational wave signal, potentially making it detectable to
upcoming gravitational wave experiments. However, the resulting gravitational
wave signal is quite sensitive to the end of the early matter-dominated epoch.
If matter domination ends gradually, a cancellation results in an extremely
suppressed signal, while in the limit of an instantaneous transition, there is
a resonant-like enhancement. The end of the matter dominated epoch cannot be
instantaneous, however, and previous analyses have used a Gaussian smoothing
technique to account for this, and consider only a limited regime around the
fast transition limit. In this work, we present a study of the enhanced
gravitational wave signal from early matter domination without making either
approximation and show how the signal smoothly evolves from the strongly
suppressed to strongly enhanced regimes.

##  1 Introduction

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] and non-topological solitons [22, 23]. Currently,
there are no constraints on the existence or duration of such an epoch,
provided that it ended before nucleosynthesis [24].

Since an early MD epoch can arise in such a varied array of models [25, 26,
27, 28, 23, 29, 30, 31, 32, 33, 22, 34, 35, 36, 37, 38, 21] (for a review, see
ref. [39]), 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, 23], 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 A.

##  2 Modelling end of early matter domination

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, 40, 41, 42]

$\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. [43]. 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.111We 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 actually occur at time $\eta_{R}$ when the decay becomes
efficient compared to the Hubble rate $\Gamma\sim\mathcal{H}$. As such
$\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 [44]

$\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.

Figure 1: The time dependence of the scale factor in a gradual transition. The
numerical solution to eq. (2.11) in black is contrasted to the instantaneous
approximation eq. (2.12) in dashed red and dotted blue for the two choices of
$\eta_{R}$ as discussed in the text. The errors shown are the fractional
errors in the each fit compared to the numerical computation, given by
$\frac{|a-a_{\rm f}|}{a}$

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. In the the conformal Newtonian gauge,
neglecting anisotropic stress of radiation, the perturbations for a given
wavenumber $k$ evolve as [44]

$\displaystyle\delta_{m}^{\prime}$
$\displaystyle=-\theta_{m}+3\Phi^{\prime}-a\Gamma(\eta)\Phi$ (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)$
(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)

with 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 [45]

$\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.

Figure 2: Time dependence of the gravitational potential $\Phi$ for wavenumber
$k=100/\eta_{R}$. The evolution for a rapid transition with
$\beta=3000/\eta_{R}$ in blue is contrasted with that of a gradual transition
with $\beta=0/\eta_{R}$ in dashed orange. The error is that of the analytic
approximation during the RD era defined as $|\Phi-\Phi_{A}|/\Phi_{\rm amp}$
with $\Phi_{\rm amp}$ as defined in eq. (2.29). We note that $\eta_{\rm osc}$
differs in the rapid and graduate cases, explaining why the error appears at
different times.

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 A, 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 [46]. 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}}$ [43, 47] 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}$ [48] and the redshift $z_{\rm
eq}=3411$ [43]. Finally taking $g_{\rm eq}=3.38$ and $g_{s,{\rm eq}}=3.94$
[49] 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.

##  3 Results

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.

$\beta$ |  $\eta_{\Gamma}$ (Mpc) |  $\Gamma_{\rm max}$ (GeV) | $\tau H$
---|---|---|---
$3000/\eta_{R}$ | $1.157\times 10^{-9}$ | $7.948\times 10^{-12}$ | $2.685\times 10^{-3}$
$500/\eta_{R}$ | $1.158\times 10^{9}$ | $7.948\times 10^{-12}$ | $9.926\times 10^{-3}$
$100/\eta_{R}$ | $1.171\times 10^{-9}$ | $7.948\times 10^{-12}$ | $4.490\times 10^{-2}$
$50/\eta_{R}$ | $1.192\times 10^{-9}$ | $7.948\times 10^{-12}$ | $8.641\times 10^{-2}$
$25/\eta_{R}$ | $1.242\times 10^{-9}$ | $7.948\times 10^{-12}$ | $1.623\times 10^{-1}$
$10/\eta_{R}$ | $1.420\times 10^{-9}$ | $7.948\times 10^{-12}$ | $3.423\times 10^{-1}$
$5/\eta_{R}$ | $1.750\times 10^{-9}$ | $7.948\times 10^{-12}$ | $5.291\times 10^{-1}$
$1/\eta_{R}$ | $4.612\times 10^{-9}$ | $7.948\times 10^{-12}$ | $8.062\times 10^{-1}$
$0$ | - | $3.179\times 10^{-14}$ | $8.373\times 10^{-1}$
Table 1: Parameters associated with figure 3, where we fixed $\eta_{R}=1.164\times 10^{-9}\;{\rm Mpc}$, which corresponds to $T_{R}=100\;{\rm GeV}$ and $k_{\rm max}=3.867\times 10^{11}\;{\rm Mpc^{-1}}$. When $\beta=0/\eta_{R}$ this uniquely determines $\Gamma_{\rm max}$. In the opposite limit, to have a rapid transition, it is necessary to have $\Gamma_{\rm max}\gg H$ and thus we use larger values for nonzero $\beta$. $\beta$ |  $\eta_{R}$ (Mpc) |  $T_{R}$ (GeV) | $k_{\rm max}\;({\rm Mpc^{-1}})$ | $\tau H$
---|---|---|---|---
$3000/\eta_{\Gamma}$ | $1.170\times 10^{-9}$ | $99.42$ | $3.844\times 10^{11}$ | $2.663\times 10^{-3}$
$500/\eta_{\Gamma}$ | $1.169\times 10^{-9}$ | $99.52$ | $3.848\times 10^{11}$ | $9.877\times 10^{-3}$
$100/\eta_{\Gamma}$ | $1.157\times 10^{-9}$ | $100.6$ | $3.890\times 10^{11}$ | $4.516\times 10^{-2}$
$50/\eta_{\Gamma}$ | $1.136\times 10^{-9}$ | $102.5$ | $3.962\times 10^{11}$ | $8.837\times 10^{-2}$
$25/\eta_{\Gamma}$ | $1.087\times 10^{-9}$ | $107.1$ | $4.141\times 10^{11}$ | $1.722\times 10^{-1}$
$10/\eta_{\Gamma}$ | $9.237\times 10^{-10}$ | $126.0$ | $4.873\times 10^{11}$ | $4.003\times 10^{-1}$
$5/\eta_{\Gamma}$ | $6.637\times 10^{-10}$ | $175.3$ | $6.793\times 10^{11}$ | $6.617\times 10^{-1}$
$1/\eta_{\Gamma}$ | $1.387\times 10^{-10}$ | $838.9$ | $3.215\times 10^{12}$ | $8.300\times 10^{-1}$
$0$ | $7.400\times 10^{-11}$ | $1572$ | $6.048\times 10^{12}$ | $8.259\times 10^{-1}$
Table 2: Parameters associated with figure 4, where we fixed
$\eta_{\Gamma}=1.164\times 10^{-9}\;{\rm Mpc}$ and $\Gamma_{\rm
max}=7.948\times 10^{-12}\;{\rm GeV}$. Figure 3: GW power spectrum induced
from a scale invariant power spectrum ($n_{s}=1$). we fixed $\eta_{\Gamma}$ in
eq. (2.8) to produce a transition occurring at $T_{R}=100{\rm GeV}$. The
curves are labelled by the value of $\beta$. We compare our results to the
rapid limit from ref. [12] in solid black and the gradual limit from ref. [13]
in dashed black. Figure 4: GW power spectrum induced from a scale invariant
power spectrum ($n_{s}=1$). We fixed $\eta_{\Gamma}=1.164\times 10^{-9}\;{\rm
Mpc}$ in eq. (2.8) so that the most rapid transition occurs at $T_{R}=100{\rm
GeV}$. The curves are labelled by the value of $\beta$. We compare our results
to the rapid limit from ref. [12] in solid black and the gradual limit from
ref. [13] in dashed black.

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, we have included the rapid limit in both figures, constructed
from the fit in the appendix of ref. [12], as the solid black curve.
Similarly, the dashed 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 blue curve ($k=3000/\eta_{R}$) perfectly
matching the analytic approximation, both in the broad low $k/k_{\rm max}$
regime of the spectrum, as well as at the resonant peak. In addition, we show
in orange the $\beta=500/\eta_{R}$ curve, which is very slightly suppressed in
comparison to the rapid limit, indicating that the $\beta=3000/\eta_{R}$ curve
has truly converged to the 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.

Figure 5: The GW signal $\Omega_{\rm GW}$ for a gradual transition with
$\beta=0/\eta_{R}$ decomposed into contributions from the MD and RD eras and
the cross term.

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]222We 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.

Figure 6: GW power spectrum comparing the gradual limit for different levels
of approximation. The solid black curve shows the limit reproduced from ref.
[13] with an instantaneous approximation for $a$ and the Green’s function. The
dash-dotted blue curve show our numerical result if we were to use the same
piecewise approximation for the Green’s function as in ref. [13]. Lastly the
dashed orange curve shows our fully numerical results.

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$ [43]. 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}$.

Figure 7: The GW signal induced from a transition at $T_{R}=100{\rm GeV}$ with
a realistic curvature power spectrum as in eq. (2.3). The signals are labelled
by the value of $\beta$. We see only the $\beta=500/\eta_{R}$ and
$\beta=100/\eta_{R}$ spectra would be detectable at LISA. We compare to
sensitivity curves of future GW experiments for LISA observing for four years
[6], the Einstein Telescope observing for one year [50, 51, 46], the Cosmic
Explorer [52], DECIGO with three units observing for one year [50, 53], THEIA
observing for twenty years [54, 55] and SKA [56].

Finally, we contrast our approach with those of references [21, 23]. 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.

##  4 Concluding remarks

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.

## References

  * [1] LIGO Scientific, VIRGO collaboration, _Characterization of the LIGO detectors during their sixth science run_ , _Class. Quant. Grav._ 32 (2015) 115012 [1410.7764].
  * [2] LIGO Scientific, Virgo collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102 [1602.03837].
  * [3] LIGO Scientific, Virgo collaboration, _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101 [1710.05832].
  * [4] LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT collaboration, _Multi-messenger Observations of a Binary Neutron Star Merger_ , _Astrophys. J. Lett._ 848 (2017) L12 [1710.05833].
  * [5] A. Kosowsky, M.S. Turner and R. Watkins, _Gravitational waves from first order cosmological phase transitions_ , _Phys. Rev. Lett._ 69 (1992) 2026.
  * [6] C. Caprini et al., _Detecting gravitational waves from cosmological phase transitions with LISA: an update_ , _JCAP_ 03 (2020) 024 [1910.13125].
  * [7] A. Vilenkin, _Gravitational Field of Vacuum Domain Walls and Strings_ , _Phys. Rev. D_ 23 (1981) 852.
  * [8] D.I. Dunsky, A. Ghoshal, H. Murayama, Y. Sakakihara and G. White, _GUTs, hybrid topological defects, and gravitational waves_ , _Phys. Rev. D_ 106 (2022) 075030 [2111.08750].
  * [9] H. Assadullahi and D. Wands, _Gravitational waves from an early matter era_ , _Phys. Rev. D_ 79 (2009) 083511 [0901.0989].
  * [10] L. Alabidi, K. Kohri, M. Sasaki and Y. Sendouda, _Observable induced gravitational waves from an early matter phase_ , _JCAP_ 05 (2013) 033 [1303.4519].
  * [11] K. Kohri and T. Terada, _Semianalytic calculation of gravitational wave spectrum nonlinearly induced from primordial curvature perturbations_ , _Phys. Rev. D_ 97 (2018) 123532 [1804.08577].
  * [12] K. Inomata, K. Kohri, T. Nakama and T. Terada, _Enhancement of Gravitational Waves Induced by Scalar Perturbations due to a Sudden Transition from an Early Matter Era to the Radiation Era_ , _Phys. Rev. D_ 100 (2019) 043532 [1904.12879].
  * [13] K. Inomata, K. Kohri, T. Nakama and T. Terada, _Gravitational Waves Induced by Scalar Perturbations during a Gradual Transition from an Early Matter Era to the Radiation Era_ , _JCAP_ 10 (2019) 071 [1904.12878].
  * [14] G.D. Coughlan, W. Fischler, E.W. Kolb, S. Raby and G.G. Ross, _Cosmological Problems for the Polonyi Potential_ , _Phys. Lett. B_ 131 (1983) 59.
  * [15] B. de Carlos, J.A. Casas, F. Quevedo and E. Roulet, _Model independent properties and cosmological implications of the dilaton and moduli sectors of 4-d strings_ , _Phys. Lett. B_ 318 (1993) 447 [hep-ph/9308325].
  * [16] T. Banks, D.B. Kaplan and A.E. Nelson, _Cosmological implications of dynamical supersymmetry breaking_ , _Phys. Rev. D_ 49 (1994) 779 [hep-ph/9308292].
  * [17] T. Banks, M. Berkooz and P.J. Steinhardt, _The Cosmological moduli problem, supersymmetry breaking, and stability in postinflationary cosmology_ , _Phys. Rev. D_ 52 (1995) 705 [hep-th/9501053].
  * [18] M. Pospelov, A. Ritz and M.B. Voloshin, _Secluded WIMP Dark Matter_ , _Phys. Lett. B_ 662 (2008) 53 [0711.4866].
  * [19] N. Arkani-Hamed, D.P. Finkbeiner, T.R. Slatyer and N. Weiner, _A Theory of Dark Matter_ , _Phys. Rev. D_ 79 (2009) 015014 [0810.0713].
  * [20] A. Berlin, D. Hooper and G. Krnjaic, _PeV-Scale Dark Matter as a Thermal Relic of a Decoupled Sector_ , _Phys. Lett. B_ 760 (2016) 106 [1602.08490].
  * [21] K. Inomata, M. Kawasaki, K. Mukaida, T. Terada and T.T. Yanagida, _Gravitational Wave Production right after a Primordial Black Hole Evaporation_ , _Phys. Rev. D_ 101 (2020) 123533 [2003.10455].
  * [22] G. White, L. Pearce, D. Vagie and A. Kusenko, _Detectable Gravitational Wave Signals from Affleck-Dine Baryogenesis_ , _Phys. Rev. Lett._ 127 (2021) 181601 [2105.11655].
  * [23] S. Kasuya, M. Kawasaki and K. Murai, _Enhancement of second-order gravitational waves at Q-ball decay_ , _JCAP_ 05 (2023) 053 [2212.13370].
  * [24] M.A. Amin, M.P. Hertzberg, D.I. Kaiser and J. Karouby, _Nonperturbative Dynamics Of Reheating After Inflation: A Review_ , _Int. J. Mod. Phys. D_ 24 (2014) 1530003 [1410.3808].
  * [25] K. Harigaya, K. Inomata and T. Terada, _Gravitational wave production from axion rotations right after a transition to kination_ , _Phys. Rev. D_ 108 (2023) L081303 [2305.14242].
  * [26] D. Borah, S. Jyoti Das, R. Roshan and R. Samanta, _Imprint of PBH domination on gravitational waves generated by cosmic strings_ , _Phys. Rev. D_ 108 (2023) 023531 [2304.11844].
  * [27] T.C. Gehrman, B. Shams Es Haghi, K. Sinha and T. Xu, _The primordial black holes that disappeared: connections to dark matter and MHz-GHz gravitational Waves_ , _JCAP_ 10 (2023) 001 [2304.09194].
  * [28] G. Domènech and M. Sasaki, _Gravitational wave hints black hole remnants as dark matter_ , _Class. Quant. Grav._ 40 (2023) 177001 [2303.07661].
  * [29] N. Bhaumik, A. Ghoshal, R.K. Jain and M. Lewicki, _Distinct signatures of spinning PBH domination and evaporation: doubly peaked gravitational waves, dark relics and CMB complementarity_ , _JHEP_ 05 (2023) 169 [2212.00775].
  * [30] C. Chen, A. Ota, H.-Y. Zhu and Y. Zhu, _Missing one-loop contributions in secondary gravitational waves_ , _Phys. Rev. D_ 107 (2023) 083518 [2210.17176].
  * [31] J.L. Rey Idler, _Primordial black holes from inflation and their gravitational wave signals_ , Ph.D. thesis, U. Autonoma, Madrid (main), 2022.
  * [32] N. Bhaumik, A. Ghoshal and M. Lewicki, _Doubly peaked induced stochastic gravitational wave background: testing baryogenesis from primordial black holes_ , _JHEP_ 07 (2022) 130 [2205.06260].
  * [33] J. Kozaczuk, T. Lin and E. Villarama, _Signals of primordial black holes at gravitational wave interferometers_ , _Phys. Rev. D_ 105 (2022) 123023 [2108.12475].
  * [34] M.R. Haque, D. Maity, T. Paul and L. Sriramkumar, _Decoding the phases of early and late time reheating through imprints on primordial gravitational waves_ , _Phys. Rev. D_ 104 (2021) 063513 [2105.09242].
  * [35] G. Domènech, V. Takhistov and M. Sasaki, _Exploring evaporating primordial black holes with gravitational waves_ , _Phys. Lett. B_ 823 (2021) 136722 [2105.06816].
  * [36] I. Dalianis and C. Kouvaris, _Gravitational waves from density perturbations in an early matter domination era_ , _JCAP_ 07 (2021) 046 [2012.09255].
  * [37] G. Domènech, C. Lin and M. Sasaki, _Gravitational wave constraints on the primordial black hole dominated early universe_ , _JCAP_ 04 (2021) 062 [2012.08151].
  * [38] N. Bhaumik and R.K. Jain, _Small scale induced gravitational waves from primordial black holes, a stringent lower mass bound, and the imprints of an early matter to radiation transition_ , _Phys. Rev. D_ 104 (2021) 023531 [2009.10424].
  * [39] G. Domènech, _Scalar Induced Gravitational Waves Review_ , _Universe_ 7 (2021) 398 [2109.01398].
  * [40] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada and T.T. Yanagida, _Inflationary primordial black holes for the LIGO gravitational wave events and pulsar timing array experiments_ , _Phys. Rev. D_ 95 (2017) 123510 [1611.06130].
  * [41] K.N. Ananda, C. Clarkson and D. Wands, _The Cosmological gravitational wave background from primordial density perturbations_ , _Phys. Rev. D_ 75 (2007) 123518 [gr-qc/0612013].
  * [42] D. Baumann, P.J. Steinhardt, K. Takahashi and K. Ichiki, _Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations_ , _Phys. Rev. D_ 76 (2007) 084019 [hep-th/0703290].
  * [43] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [44] V. Poulin, P.D. Serpico and J. Lesgourgues, _A fresh look at linear cosmological constraints on a decaying dark matter component_ , _JCAP_ 08 (2016) 036 [1606.02073].
  * [45] D.S. Gorbunov and V.A. Rubakov, _Introduction to the Theory of the Early Universe_ , World Scientific Publishing Company (2011), 10.1142/7873.
  * [46] K. Inomata and T. Nakama, _Gravitational waves induced by scalar perturbations as probes of the small-scale primordial spectrum_ , _Phys. Rev. D_ 99 (2019) 043511 [1812.00674].
  * [47] Planck collaboration, _Planck 2015 results. XIII. Cosmological parameters_ , _Astron. Astrophys._ 594 (2016) A13 [1502.01589].
  * [48] D.J. Fixsen, _The Temperature of the Cosmic Microwave Background_ , _ApJ_ 707 (2009) 916 [0911.1955].
  * [49] L. Husdal, _On Effective Degrees of Freedom in the Early Universe_ , _Galaxies_ 4 (2016) 78 [1609.04979].
  * [50] C.J. Moore, R.H. Cole and C.P.L. Berry, _Gravitational-wave sensitivity curves_ , _Class. Quant. Grav._ 32 (2015) 015014 [1408.0740].
  * [51] M. Maggiore et al., _Science Case for the Einstein Telescope_ , _JCAP_ 03 (2020) 050 [1912.02622].
  * [52] D. Reitze et al., _Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO_ , _Bull. Am. Astron. Soc._ 51 (2019) 035 [1907.04833].
  * [53] S. Kawamura et al., _Current status of space gravitational wave antenna DECIGO and B-DECIGO_ , _PTEP_ 2021 (2021) 05A105 [2006.13545].
  * [54] The Theia Collaboration, C. Boehm, A. Krone-Martins, A. Amorim, G. Anglada-Escude, A. Brandeker et al., _Theia: Faint objects in motion or the new astrometry frontier_ , 1707.01348.
  * [55] J. Garcia-Bellido, H. Murayama and G. White, _Exploring the early Universe with Gaia and Theia_ , _JCAP_ 12 (2021) 023 [2104.04778].
  * [56] G. Janssen et al., _Gravitational wave astronomy with the SKA_ , _PoS_ AASKA14 (2015) 037 [1501.00127].
  * [57] I. Bogaert, _Iteration-free computation of gauss–legendre quadrature nodes and weights_ , _SIAM Journal on Scientific Computing_ 36 (2014) A1008 [https://doi.org/10.1137/140954969].

##  Appendix A Technical Details of Code

There are several systems of differential equations that need to be solved
numerically in the process of calculating the GW spectra. These are the
Friedmann equations, eqs. (2.9) to (2.11), the GW Green’s function, eq. (2.5),
and the evolution of the perturbations, eqs. (2.17) to (2.21). For this
purpose we use Odeint from the Boost C++ library. Specifically we use the
runge_kutta_dopri5 solver with an absolute tolerance of $10^{-200}$ and a
relative tolerance of $10^{-11}$. The absolute tolerance of $10^{-200}$ is
necessary to prevent the solver stalling once $\rho_{m}$ decays to values
close to the smallest representable machine number.

To perform integrals numerically we used Gauss-Legendre quadratures which we
generate using the algorithm described in ref. [57]. For the integral over
$\bar{\eta}$ in eq. (2.4) we note that once $\Phi$ begins to oscillate it does
so with period $2\pi\sqrt{3}/k$ and we ensure that the integrand is sampled
densely enough that there are approximately 50 points per period.

For the integral over $u$ and $v$ in eq. (2.2) we first transform to the
variables $s=u-v$ and $t=u+v-1$ which gives [11]

$\overline{\mathcal{P}_{h}(\eta,k)}=2\int_{0}^{\infty}dt\int_{-1}^{1}ds\left(\frac{t(2+t)(s^{2}-1)}{(1-s+t)(1+s+t)}\right)^{2}\overline{I^{2}(u,v,k,\eta)}\mathcal{P}_{\zeta}(uk)\mathcal{P}_{\zeta}(vk)$
(A.1)

where the remaining factors of $u$ and $v$ can be expressed in terms of $s$
and $t$ as $u=(t+s+1)/2$ and $v=(t-s+1)/2$. This form has several benefits.
Firstly, the limits of the inner integral over $s$ are no longer dependent on
the outer integration variable. Additionally the integrand is symmetric in $s$
which we will exploit to only have to perform the integral from 0 to 1.
Lastly, $I$ has a pole at $t=\sqrt{3}-1$, which is responsible for the
resonant peak in the rapid transition spectra.

We found that the integrand was only very mildly dependent on $s$ and so we
could perform the integral over $s$ from 0 to 1 accurately with only 5
quadrature points. Increasing the number of quadrature points in $s$ to 10
only impacted the result at less than the $\mathcal{O}(0.01\%)$ level in both
the rapid and gradual cases. For the integral over $t$ we divide the
integration region into $0\leq t<\sqrt{3}-1$, $\sqrt{3}-1<t\leq 2(\sqrt{3}-1)$
and $t>2(\sqrt{3}-1)$. To ensure we account for the important contribution
from the pole we sample the regions above and below the pole with 40
quadrature points each. The remainder of the integral was computed with 20
more quadrature points. We found that this distribution of quadratures could
have as much as 5 percent error for some values of $k$ for a rapid transition
when compared to using 200 points either side of the pole and 100 in the
remainder. The error in the gradual case was negligible due to $I$ becoming
approximately independent of $s$ and $t$. Using this larger number of
quadrature points would significantly increase the run-time and a 5 percent
error is adequate to compare with the gradual and rapid limits from previous
work and make qualitative statements about the spectra and potential
detectability.