Numerically solving scalar field model of dark energy

Question

Not sure if I should be asking this in maths or physics section so please excuse me if asked in the wrong section. I am trying to numerically integrate the equations obtained in the exponential potential scalar field model of dark energy. The set of equations I have obtained is:

\begin{align} x &= tH_{0} \\ y_{1} &= \frac{a}{a_{0}} \\ y_{2} &= \frac{\psi}{\psi _{0}}\\ y_{3} &= \frac{dy_{2}}{dx} \end{align} which on differentiating w.r.t x I get:
\begin{align} {y_{1}}\prime &= y_{1}\left(\frac{\Omega_{m}}{y_{1}^3} + \frac{y_{3}^2}{6.0} + \frac{V(y_{2})}{2.H_{0}^2}\right) \\ y_{2}\prime &= y_{3} \\ y_{3}\prime &= -3\frac{y_{1}\prime}{y_{1}}y_{3} - \frac{1}{H_{0}^2}\frac{\partial v(y2)}{\partial y_{2}} \end{align} My doubt here is I am using RK4 method for solving this and when I write the increment should I increment the independent variable $x$ too?Like when I write formulae for $k_1,\,k_2,\,k_3,\,k_4$ etc I am incrementing $y_1,\,y_2,\,y_3$. As far as I know it is not needed as these ODE are autonomous equations yet when I try to plot them I am not getting the expected result.

Another doubt is I am plotting it w.r.t scale factor( $a=0$ to $a=1$) and when I plot $y_1$ which is $y_{1} = \frac{a}{a_{0}}$ shouldn't I be getting a linear curve or am I missing on some major concept?

Answer 1

Let's look at what $y_1 = \frac{a}{a_0}$ means.

Since $a = a_0 e^N$, where $N$ is the e-fold number, $\frac{a}{a_0} = e^N$ and so $y_1 = e^N$.

When you plot $y_1$ against $N$, you do indeed get the green curve shown in the plot in your question (using a calculator to plug in values of $N$ to get $e^N$ gives a rough approximation).

I think the next step would be for you to plot the Hubble parameter $H$ against redshift $z$ (you can find $H$ by using $H = \frac{\dot{a}}{a}$). Depending on the values you are using for the slope of the potential (generally called $\lambda$ in the literature) and the initial values of the potential $V_0$ and the field itself ($\phi$ and $\dot{\phi}$)*, you should be able to find a value of $H_0$ consistent with observational data i.e. $H_0 \approx 70$ kms$^{-1}$Mpc$^{-1}$.

Put $\Omega_m$ and $\Omega_\phi$ on a separate plot too. If you plot these against redshift you should again be able to find a set of model parameters where $\Omega_m \approx 0.3$ and $\Omega_\phi \approx 0.7$, again in line with current observations.

Also, if you want to make life a little easier for yourself, use from scipy import integrate and then integrate.odeint. This is a nice function in Python and saves you from having to solve the equations numerically.

---

*I was coupling the field to matter, but found that values of $\lambda = 0.23$, $V_0 = 0.76 \times 10 ^{-120}$, $\phi_0= 0.7$ and $\dot{\phi_0}= 0$ worked nicely. I was matching the Planck data, so this gave $H_0 = 67.8$ kms$^{-1}$Mpc$^{-1}$.