In cosmic inflation, we have the equation of motion for the inflaton given by,(adsbygoogle = window.adsbygoogle || []).push({});

$$\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)

**Physics Forums - The Fusion of Science and Community**

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Mathematica Cosmic inflation using mathematica

Have something to add?

Draft saved
Draft deleted

Loading...

Similar Threads - Cosmic inflation using | Date |
---|---|

Using Maxima to plot error in Fourier series | Nov 5, 2017 |

Matlab Solving equation with integration using MATLAB | Aug 4, 2017 |

What Software do you use to do your math manipulations? | May 21, 2017 |

ANSYS APDL Inflated Balloon Volume | Aug 21, 2011 |

**Physics Forums - The Fusion of Science and Community**