# Numerically solving scalar field model of dark energy

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? • What language are you using for the numerical integration? (Or are you doing this by hand??) Because it's relatively easy to solve these equations using the scipy integrate function in Python- I did exactly this for my Master's dissertation. – astronat Jul 7 '17 at 0:00
• @astronat I am using python and I wrote the entire algorithm for solving. Problem is when I try to plot matter density and scalar density the plots dosent make any sense the eos does. Can you help me with this? Should I edit the question and include the entire code used by me? – tsv Jul 7 '17 at 3:28
• Can you add a picture of plot of the matter and scalare field density parameters? – astronat Jul 7 '17 at 7:27
• @astronat Have added the plot in the question. – tsv Jul 7 '17 at 7:34
• It sounds to me that you may have an error in your RK4 coding. I suggest to try solving a 2nd order ODE that can be solved analytically and check that you have no coding errors. – Lewis Miller Jul 7 '17 at 12:37

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}$.

• Thank you. This did clear things up though I am still stuck on using scipy(throws up some error-am using python for the first time), but did solve it using the numerical method. – tsv Jul 7 '17 at 11:25
• No probs. If you like, take a look at the appendix in my dissertation where I've got a listing of the code I wrote to do this (the relevant bit is on page 43). – astronat Jul 7 '17 at 11:33