Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Mathematica Cosmic inflation using mathematica

  1. Apr 27, 2017 #1
    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.
    Image.jpg

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

    rho.jpg
     
  2. jcsd
  3. May 2, 2017 #2
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Cosmic inflation using mathematica
  1. Mathematica to use it (Replies: 2)

Loading...