Decay of the time derivative of solutions of the Klein-Gordon equation in decelerating expanding space-times Suppose that we have a model of a universe* given by a flat FLRW metric.* In short, the model universe has $n\in\mathbb N$ dimensions, is homogeneous, isotropic and its expansion is governed solely by the scale factor $a:]0,\infty[\to]0,\infty[$. In the limiting case (the "big bang", $t=0$), we would have $a(0)=0$.
The Klein-Gordon equation (also known as the wave equation) in this model universe is given by (see arXiv:1801.08944, formula (14))
\begin{equation}\tag{*}\label{*}\ddot\phi(t, x)+\frac{n\dot a(t)}{a(t)}\dot\phi(t,x)-\frac1{a^2(t)}\Delta\phi(t,x) =0,\end{equation}
where $\Delta=\Delta_x$ is the Laplace operator and the solution is a (for instance smooth) function $$\phi:]0,\infty[\times\hat M\to\mathbb R,$$ where $\hat M$ is the space form of the model universe. (I will only look at $\hat M=(\mathbb R/\mathbb Z)^3$, and $\hat M =\mathbb R^3$.)

I want to prove suitable decay of $\partial_t\phi$ for any solutions to \eqref{*}. In arXiv:1801.08944, appendix A, they give a very informal heuristic for why they expect that, when $a(t)=t^p$ for some $p\in]0,1[$ (which means decelerating expansion), then we have $$\lVert\partial_t\phi(t,\cdot)\rVert_{L^\infty(\hat M)}\lesssim a(t)^{-\frac{n+1}2}.$$
My goal will be to prove this. $\phi$ will always be assumed to

*

*in the Schwartz space $S(\mathbb R^n)$ if $\hat M=\mathbb R^n$;

*smooth if $\hat M=(\mathbb R/\mathbb Z)^n$.

My question: How can I do this?
My attempt: See my answer below.

* For my definition of the FLRW metric, see for instance this question of mine (only dimension $n=3$) or formula (12) of arXiv:1801.08944 (any dimension $n\in\mathbb N$).
 A: I will only treat the special case $n=3$ and $p=\frac 23$ (the Einstein-De-Sitter universe). I encourage the reader to try to generalize this to arbitrary $n\in\mathbb N$ and $p\in]0,1[$.

By the above regularity assumptions, it is valid to Fourier expand,* i.e. to write $$\phi(x, t) = \int_{\mathbb R^3} c_{\mathbf k}(t) \exp(2\pi i\langle \mathbf k, x\rangle)\,\mathrm d\mathbf k$$
if $\hat M =\mathbb R^3$ and $$\phi(x, t) = \sum_{\mathbf k\in\mathbb Z^3} c_{\mathbf k}(t) \exp(2\pi i\langle\mathbf k, x\rangle)\,\mathrm d\mathbf k,$$
if $\hat M=(\mathbb R/\mathbb Z)^3$, where, in both cases, $$c_{\mathbf k}(t) = \int_{\hat M} \phi(x, t)\exp(-2\pi i \langle\mathbf k,x\rangle)\,\mathrm d{\mathbf k}.$$

The Fourier modes $c_{\mathbf k}$ satisfy the ODE $$\ddot c_{\mathbf k}(t)+\frac2t\dot c_{\mathbf k}(t)+k^2 t^{-4/3} c_{\mathbf k}(t)=0.$$
Following appendix A of arXiv:1801.08944, by defining

*

*$\tau = \tau(t)=3 t^{1/3}$;

*$\tilde c_{\mathbf k}(\tau(t))= c_{\mathbf k}(t)$;

*$\tilde d_{\mathbf k}(\tau) = \left(\frac{\tau}3\right)^2 \tilde c_{\mathbf k}(\tau)$;

one obtains the new ODE ($k$ denotes the absolute value of $\mathbf k\in\hat M$; Note also that the complex conjugate of $c_{\mathbf k}$ is given by $c_{-\mathbf k}$ since $\phi$ is real) $$\tilde d_{\mathbf k}''(\tau)+\left(k^2-\frac2{\tau^2}\right)\tilde d_{\mathbf k}(\tau)=0.$$
We also have $$\tilde d_{\mathbf k}(\tau(t)) = \frac{\tau(t)^2}9 \tilde c_{\mathbf k}(\tau(t)) = t^{2/3} c_{\mathbf k}(t).$$

The general solution to this is (TODO: Say how to solve it, the idea is to introduce a new function $\tilde e_{\mathbf k}(\tau) =\tau^{-1/2} \tilde d_{\mathbf k}(\tau)$), for $k\neq 0$, given by $$\tilde d_{\mathbf k}(\tau) = C_1 \sqrt{\tau} J_{3/2}(k\tau) + C_2\sqrt{\tau} Y_{3/2}(k\tau),$$
where $J, Y$ are Bessel functions of the first and second kind, respectively.** The constants $C_1, C_2$ also depend on $\mathbf k$. Note that $Y$ has a non-removable singularity at $0$, which is why this solution is not well-defined at $k=0$.
$C_1, C_2$ can be determined by setting $\tau = 1$ and $\tau = 2$. Let $\vec c =\vec c(\mathbf k)= \begin{pmatrix} C_1\\C_2\end{pmatrix}$ and $\vec d =\vec d(\mathbf k)=\begin{pmatrix} \tilde d_{\mathbf k}(1)\\\tilde d_{\mathbf k}(2)\end{pmatrix}$, then we have $$\vec c = A^{-1} \vec d,$$ where $$A=A(k) = \begin{pmatrix} J_{3/2}(k) & Y_{3/2}(k)\\\sqrt{2} J_{3/2}(2k) & \sqrt{2} Y_{3/2}(2 k)\end{pmatrix}.$$
Note that $A(k)$ is only invertible for $k$ in the complement of a discrete subset of $\mathbb R$ (TODO: Prove this), but one can take a limit since the Fourier modes are continuous, so this is not a big problem.
We can now write for the $k\neq 0$ not in the discrete set from before, $$\tilde d_{\mathbf k}(\tau) = \sqrt{\tau} \left\langle A(k)^{-1} \vec d(\mathbf k), \begin{pmatrix} J_{3/2}(k\tau)\\ Y_{3/2}(k\tau)\end{pmatrix}\right\rangle.$$

Remark. Setting for instance $\vec d(\mathbf k) = (1,1)^\top$, one gets $$\tilde d_{\mathbf k}(\tau)=\frac{-2 \left(k^2 \tau +1\right) \sin (k-k \tau )+\left(2 k^2 \tau +1\right) \sin (2 k-k \tau )+k (\tau -2) \cos (k (\tau -2))-2 k (\tau -1) \cos (k-k \tau )}{\tau \left(2 k^2 \sin (k)+\sin (k)-k \cos (k)\right)}$$
which equals for fixed $\tau$ and $k\to 0$ $$\frac{\frac{6+\tau^3}3 k^3 + O(k^5)}{\frac{7\tau k^3}3+O(k^4)},$$ which is $$\frac{6+\tau^3}{7\tau}+O(k^2).$$ In particular, this fits very nicely with the fact that the general solution to (*) with $k=0$ is $$C_1 \tau^2 + C_2/\tau.$$

The rest is done in my answer about decay of solutions to the wave equation $\ddot\phi(t, x)+\frac{n p}{t}\dot\phi(t,x)-t^{-2p}\Delta\phi(t,x)=0$.

Remark. Using properties of the Bessel functions, it seems that my argument can be expanded to all $p\in]0,1[$. In particular, it would also apply to different physical theories, including for example 11-dimensional supergravity, for which $p=\frac{1+\sqrt{21}}{10}$ according to arXiv:1809.04724, formulas (12), (13).

* Cf. https://en.wikipedia.org/wiki/Fourier_inversion_theorem and https://en.wikipedia.org/wiki/Convergence_of_Fourier_series and the references contained therein.
** See https://en.wikipedia.org/wiki/Bessel_function.
