# Mathematica Cosmic inflation using mathematica

1. Apr 27, 2017

### Whitehole

In cosmic inflation, we have the equation of motion for the inflaton given by,

$$\ddot\phi + 3H(1 + Q)\dot\phi + V_\phi = 0$$

the Friedman equation is given by

$$H^2 = \frac{1}{3 M_p^2}(\frac{1}{2} \dot\phi^2 + V + \rho_r)$$

where $M_p$ is the reduced Planck mass. The differential equation for the radiation density is given by

$$\dot\rho_r + 4H\rho_r = 3 Q H \dot\phi^2$$

Now for a given $V$, I want to solve these differential equations using NDSolve, let $V = \lambda (\phi^2 - \sigma^2)^2$

Let $\phi[t] = y[t]$, $\rho_r[t] = \rho[t]$, $Q=100$

Mp = 2.4353*10^18; (* Reduced Planck mass = 2.4353*10^18 GeV *)
m = 1.8*10^13; (* Inflaton mass = 1.8*10^13 GeV *)
Rm = (1.8*10^13 )/(2.4353*10^18); (* Rescaled inflaton mass *)
$\sigma$ = 2.24*10^19; (*Energy scale of symmetry breaking = 2.24*10^19 GeV*)
$R\sigma$ = (2.24*10^19)/(2.4353*10^18); (*Rescaled energy scale of symmetry breaking*)
$\lambda$ = ((1.8*10^13)/(2.4353*10^18))^4;
tf = 10^9;
sol[a_, b_] := NDSolve[{dy'[t] + 303 H[t] dy[t] + 4 \lambda (y[t]^2 - R\sigma^2) y[t] == 0, H[t] == Sqrt[(0.5 dy[t]^2 + \lambda (y[t]^2 - R\sigma^2)^2 + \rho[t])/3], \rho'[t] + 4 H[t] \rho[t] == 300 H[t] dy[t]^2, y'[t] == dy[t], y[0] == -a, dy[0] == b, \rho[0] == Rm^4}, {y, dy, H, \rho}, {t, 0, tf}]
sol[1, 1/100000]

This will generate interpolating functions. I need to determine the proper initial conditions $y[0] = -a$ and $tf$ ($dy[0] = b$, but it can be set to $10^{-6}$) for a fixed $Q$, in this case suppose I choose $Q=100$, so that I can get the correct values for the quantities.The initial conditions and $tf$ should be adjusted so that,

($N = \int Hdt$ where I represented $N$ by N1x in the code)
H[t] == Sqrt[(0.5 dy[t]^2 + \[Lambda] (y[t]^2 - R\[Sigma]^2)^2 + \[Rho][t])/3];
ndsol[a_, b_] := NDSolve[{D[N1x[t], t] == Evaluate[H[t] /. First@sol[a, b]], N1x[0] == 0}, N1x, {t, 0, tf}]

This plot should reach N1x = 60 after the plot/calculation,

Manipulate[Plot[Evaluate[N1x[t] /. First@ndsol[a, b]], {t, 0, tf}], {{a, 1}, 0, 100, Appearance -> "Labeled"}, {{b, 1/100000}, 0, 0.00009, Appearance -> "Labeled"}]

and the plot of $T$ vs. $t$ and $H$ vs $t$ generates a plot like the image below. ($~T = \Big(\frac{3 \rho_r}{10 \pi^2}\Big)^\frac{1}{4}$)
Manipulate[Plot[{Evaluate[H[t] /. First@sol[a, b]], Evaluate[(3 \[Rho][t]/10 Pi^2)^(1/4) /. First@sol[a, b]]}, {t, 0, tf}, PlotRange -> Automatic], {{a, 1}, 0, 50, Appearance -> "Labeled"}, {{b, 1/100000}, 0, 0.00009, Appearance -> "Labeled"}]

My code generates an error "For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions" for some $y[0] = -a$ and $tf$, so I'm having a hard time determining those two parameters. Actually, I don't know what tf to actually use, I'm just guessing. Is there a way to tell which $tf$ I should use so that the only thing I need to vary is $y[0]$?
A general rule for determining $y[0]$ is that it shouldn't be too far from the minimum of the potential $V$ but also not too close.

An example of my plot where $T \rightarrow yellow$ and $H \rightarrow blue$ (you can't even see it)

2. May 2, 2017

### PF_Help_Bot

Thanks for the thread! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post? The more details the better.