Numerical solution of the Mukhanov-Sasaki equation

[Rafid Mahbub]

Hi,

I am trying to figure out how to solve the Mukhanov equation numerically in Mathematica, but have some problems dealing with it. In terms of the number of efolds, the Fourier modes satisfy the following ODE in terms of the Hubble slow roll parameters:

$$ \frac{d^{2}u_{k}}{dN^{2}}+(1-\epsilon_{H})\frac{du_{k}}{dN}+\left[ \frac{k^{2}}{\mathcal{H}}+(1+\epsilon_{H}-\eta_{H})(\eta_{H}-2)-\frac{d}{dN}(\epsilon_{H}-\eta_{H}) \right]u_{k}=0 $$
the solutions of which give the power spectrum $$ \mathcal{P}_{\mathcal{R}}=\frac{k^{3}}{2\pi^{2}}|\frac{u_{k}}{z}|^{2}_{k<<\mathcal{H}} $$
Now the Hubble slow roll parameters depend on the solution of the inflaton's equation which I know how to solve. I am a bit concerned in how to deal with $k$ in the Mukhanov equation and then in the evaluation of the power spectrum. For this problem, the usual Bunch-Davies vacuum is assumed in the asymptotic past-
$$ u_{k}\rightarrow \frac{e^{-ik\tau}}{\sqrt{2k}} $$

 

[bapowell]

Suppose you want the largest-scale modes ($k$ corresponding to the CMB quadrupole) to be exiting the horizon at $N=60$. To initialize these largest-scale modes, go back in time until you are "close" to the BD limit, say, when $k = 100aH$. The corresponding value of $N$ should be your $N_i$ in the mode equation for that $k$.