# Numerical solution to Mukhanov-Sasaki equation

I am trying to figure out how to solve the Mukhanov-Sasaki equation to compute the power spectrum of an inflation potential that exhibits an ultra slow-roll phase that gives rise to an enhancement in the power spectrum, suitable for primordial black hole formation. In terms of the efold variable, the MS equation is

$$\frac{d^{2}\zeta_{k}}{dN^{2}}+\left( 3-\epsilon_{1}+\epsilon_{2} \right)\frac{d\zeta_{k}}{dN}+\left( \frac{k}{aH} \right)^{2}\zeta_{k}=0$$ with the Bunch-Davies vacuum as the initial condition

$$u_{k}=\lim_{k>>aH}\frac{e^{-ik\tau}}{\sqrt{2k}}$$ I have computed the Hubble flow parameters but now I am unsure how to solve this ODE for all the modes of interest. I know that I need to consider them deep inside the horizon and evolve them to horizon exit where they freeze over, however, I am unsure how to proceed.

• What specifically are you unsure about: how to initialize the equation for each $k$ in the Bunch-Davies limit? How to solve it once you've done that? – bapowell Nov 26 '18 at 15:46
• Yes, that's what I want to know. I know that I have to apply the BD vacuum to all the modes some efolds before horizon exit and let it evolve sufficiently long until it reaches a constant value. – Judas503 Nov 26 '18 at 16:35

You need to solve the MS equation separately for each $$k$$, initialized far enough in the past so that each mode is approximately in the Bunch-Davies limit. I've found that choosing $$N_i$$ so that $$k = \mathcal{O}(100)a(N_i)H(N_i)$$ is sufficient for this. Note that there is a different $$N_i$$ for each $$k$$. Here, $$a(N) = a_0 \exp(N_0-N)$$, where you need to pick a value for $$a_0$$. For example, you could set $$a_0$$ such that a particular scale, like the quadrupole, is at horizon crossing at $$N_0$$, i.e. $$k_{\ell = 2} = a_0 H_0$$.
Then, you just evolve each mode forward in time until inflation ends. The largest-scale modes will have frozen-out long before the end of inflation, but this way you can compute $$P(k)$$ by evaluating each individual mode $$u_k$$ on an equal-time slice. Alternatively, you could evolve each mode until they are well-outside the horizon, say, when $$k < aH/100$$ or something. The modes freeze-out quite soon after horizon exit so applying a generous cut-off like this should work well.
• Thanks. I have one last question. In the paper I'm studying, $k=0.05Mpc^{-1}$ has been chosen to leave the horizon at N=0. So if I want to evolve sufficiently far enough, then N will take negative values (I have only solved the inflation background for N 0 to 56). – Judas503 Nov 26 '18 at 18:27
• So what about smaller momentum modes, e.g. $k = 0.002 \, h{\rm Mpc}^{-1}$? These cross the horizon at $N<0$ in your model as well. – bapowell Nov 26 '18 at 19:09