- #1
Whitehole
- 132
- 4
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)
$$\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)