equation as $$\label{eq:3:heq}
h_{ij}''(\tau,\vec{x})-\nabla^2h_{ij}(\tau,\vec{x})
=
16\pi Ga(\tau)T_{ij}^{\mathrm{TT}}(\tau,\vec{x}).$$ where prime denote the differential with respect to the conformal time. Assuming that the typical physical scale is much smaller than the horizon scale, we neglected the terms including $a''$. Here, $T_{ij}^{\mathrm{TT}}(\tau,\vec{x})$ is computed by applying the projection operator in the momentum space as $$\begin{split}
 &T_{ij}^{TT}(\tau,\vec{k})=\Lambda_{ij,kl}(\hat{k})T_{ij}(\tau,\vec{k})
  =\Lambda_{ij,kl}\left\{\partial_k\phi\partial_l\phi\right\}(\tau,\vec{k}),\\
  &\Lambda_{ij,kl}(\hat{k}) 
  = P_{ik}(\hat{k})P_{jl}(\hat{k})  -\frac{1}{2}P_{ij}(\hat{k})P_{kl}(\hat{k}),\\
  &P_{ij}(\hat{k}) = \delta_{ij}-\hat{k}_{i}\hat{k}_{j},
\end{split}$$ where $\left\{\partial_k\phi\partial_{l}\phi\right\}(\tau,\vec{k})$ is the Fourier transform of $\partial_{k}\phi(\tau,\vec{x})\partial_{l}\phi(\tau,\vec{x})$.

We assume that the source term is nonzero only for the time interval $ {\tau}_{i}\leq{\tau}\leq{\tau}_{f}$. In this case we can solve (\[eq:3:heq\]) by using the Green’s function and obtain $$\label{eq:3:h}
{h_{ij}}({\tau},\vec{ k})
= { A_{ij}}(\vec{ k}){ \sin\left[k\left({ \tau}-{ \tau}_{f}\right)\right]}
 +{ B_{ij}}(\vec{ k}){ \cos\left[k\left({ \tau}-{ \tau}_{f}\right)\right]}
 \hspace{5ex}{ \mathrm{for}}\hspace{1ex}{\tau}\geq{\tau}_{ f},$$ where $$\begin{split}
&{ A_{ij}}(\vec{ k})
 ={\frac{-16\pi G}{k}\int^{{\tau}_{f}}_{{ \tau}_{i}}d{\tau'}
    \cos\left[k\left({\tau}_{f}-{\tau}'\right)\right]a({ \tau}')
        T_{ij}^{TT}}({ \tau}',\vec{ k}),\\
&{ B_{ij}}(\vec{ k})
 ={\frac{16\pi G}{k}\int^{{\tau}_{f}}_{{ \tau}_{i}}d{\tau'}
    \sin\left[k\left({\tau}_{f}-{\tau}'\right)\right]a({ \tau}')
        T_{ij}^{TT}}({ \tau}',\vec{ k}).
\end{split}$$ Here, $h_{ij}(\tau,\vec{k})$ is the Fourier transform of $h_{ij}(\tau,\vec{x})$.

The energy density of the gravitational waves is given by $$\rho_{\mathrm{gw}}(\tau)
=
\frac{1}{32\pi Ga^4(\tau)}
\braket{
h_{ij}'(\tau,\vec{x}) h_{ij}'(\tau,\vec{x})
}_V
=
\frac{1}{32\pi Ga^4(\tau)V}
\int\frac{d^3k}{\left(2\pi\right)^3}
h_{ij}'(\tau,\vec{k})h_{ij}^{\ast'}(\tau,\vec{k}),$$ where $\langle \cdots\rangle_V$ denotes the average over the spatial volume $V$. Taking the time average over a period of the oscillations, we obtain $\rho_{\mathrm{gw}}(\tau)$ as $$\begin{split}
