of their energy density spectrum (at sub-horizon scales) $\Omega_{\text{GW}}$, defined as the GW energy density $\rho_{\text{GW}}$ per unit logarithmic comoving wavenumber interval, normalized to the critical density $\rho_{\text{crit}}=3m_p^2H^2$ [@Caprini:2018mtu], $$\begin{aligned}
\Omega_{\text{GW}}(\tau, k) &\equiv\dfrac{1}{\rho_{\text{crit}}}\dfrac{d\rho_{\text{GW}}(\tau,k)}{d\ln k} = \frac{k^2}{12a^2(\tau)H^2(\tau)}\Delta_h^2(\tau,k)\,,  \label{energyspectrum}\end{aligned}$$ It is customary to factorize the tensor power spectrum at arbitrary times as a function of the primordial inflationary spectrum $\Delta_{h,\text{inf}}^2(k)$ \[c.f. Eq. (\[inf\])\] by means of a transfer function [@Boyle:2005se] $$\begin{aligned}
\Delta_h^2(\tau,k)\equiv T_h(\tau,k) \Delta_{h,\text{inf}}^2(k)\,,~~~~ T_h(\tau,k) \equiv {1\over 2}\left(a_k\over a(\tau)\right)^2\,,\label{transferfunc}\end{aligned}$$ which characterizes the expansion history between the moment of horizon re-entry $\tau = \tau_k$ of a given mode $k$, defined as $a_kH_k \equiv k$ with $a_k \equiv a(\tau_k)$, $H_k\equiv H(\tau_k)$, and a later moment $\tau > \tau_k$. For the power spectrum today we will use the notation $T_h(k) \equiv T_h(\tau_0,k)$. Note that the factor $1\over2$ in Eq. (\[transferfunc\]) is simply due to averaging over harmonic oscillations of the modes deep inside the horizon.

If we assume that immediately after inflation, the Universe became radiation dominated (RD) with equation of state $w = 1/3$, the resulting present-day GW energy density spectrum is (quasi-)scale invariant, for the frequency range corresponding to the modes crossing the Hubble radius during RD. Setting $n_t = 0$ and averaging 