s}}(\tau - \tau_*)\,,~~~~~~
a_{\rm RD}H_{\rm RD}\tilde{\tau}_r(\tau) \equiv 1+a_{\rm RD}H_{\rm RD}(\tau - \tau_{\rm RD})\,. \label{tautilde}\end{aligned}$$ From now on, and without loss of generality, we fix $a_* = 1$ and $\tau_* = 0$ at the end of inflation (onset of SD).

During the stiff epoch the GW equation of motion reads, using , $$\begin{aligned}
h_k^{\prime\prime}+\frac{2\alpha_{\rm s}}{\tilde\tau_s}h_k^{\prime}+k^2h_k=0\,, \label{eomconstantw}\end{aligned}$$ where $'$ denotes derivatives with respect to $\tilde\tau_s$ (since $d\tilde{\tau}_s=d\tau$, we do not distinguish between the derivatives with respect to $\tilde{\tau}_s$ or $\tau$). Requiring that the tensor mode function must match the inflationary spectrum $h_k(\tau) = h_k^{\text{inf}}$ and $h_k^\prime(\tau) = 0$ in the superhorizon limit $k\tau \ll 1$, the solution to during the stiff period when $\tau_{*}\leq\tau\leq\tau_{\text{RD}}$ is $$\begin{aligned}
h_k^{\text{(stiff)}}(\tau)=\Gamma\left(\alpha_{\rm s}+\frac{1}{2}\right)\left(\frac{2}{\alpha_{\rm s}y}\right)^{\alpha_{\rm s}-{1\over2}}J_{\alpha_{\rm s}-{1\over2}}(\alpha_{\rm s}y)~h_k^{\text{inf}}\,, ~~~~~~y \equiv {k\over aH} = {k\tilde\tau_{s}(\tau)\over \alpha_{\rm s}}\,,  \label{stiffsol}\end{aligned}$$ where $J_\nu(x)$ is the Bessel function of the first kind. Using the small argument limit of the Bessel function, $J_{\nu}(x) \simeq x^\nu / (2^\nu \Gamma(\nu+1))$ for $x \ll 1$, we obtain $h_k^{\text{(stiff)}}(\tau) \simeq h_k^{\rm inf}$ when $k \ll  aH$, as it should. The superscript $^{(\text{stiff})}$ indicates that the solution applies in the stiff epoch. We will use analogous notations in what follows.

Using , the GW equation of motion during the 