Message Boards Message Boards

An algorithm to obtain appropriate numerical solutions for the bounce

GROUPS:

I am looking to replicate an algorithm explained in the paper "Impact of new physics on the EW vacuum stability in a curved spacetime background" by E. Bentivegna, V. Branchina, F. Continoa and D. Zappal` to numerically solve the boundary-value problem encountered when trying to find the bounce solution in false vacuum decay. It can be found here: https://arxiv.org/abs/1708.01138

In the Appendix, "A Numerical computation of the bounce solution", it outlines an algorithm involving a shooting method to solve the following boundary problem:

$y^{\prime\prime}(x) + \frac{3}{x}y^{\prime}(x) = \frac{dU}{dy};\\ y(\infty) = 0,\ y^{\prime}(0) = 0, \\U(y) = \frac{1}{4}y^{4}(\gamma +\alpha\ln^{2}y + \beta\ln^{4}y),$

where ${}^{\prime}$ indicates differentiation with respect to $x$.

The value of the constants in the equation are as follows:

Mp = 2.435*10^18;
\[Alpha] = 1.4*10^-5;
\[Beta] = 6.3*10^-8;
\[Gamma] = -0.013;
\[Lambda]6 = 0;

The appendix outlines the algorithm as follows, where Eq.(2.23) is referring to the differential equation (and boundary conditions) stated above. Appendix Part 1 Appendix Part 2 Appendix Part 3

I am attempting to recreate this algorithm in Mathematica but my attempts thus far of simply trying numerically integrate the differential equation using a shooting method with boundary conditions specified in the Appendix has just yielded errors owing to finding complex infinities.

Any help with this matter would be greatly appreciated.

Attachments:
POSTED BY: Benjamin Leather
Answer
1 month ago

Hello,

I don't no if my code is correct!. Plot is smilar to your attachment.

α = 1.4*10^-5;
β = 6.3*10^-8;
γ = -0.013;
U = Re[D[1/4*y[x]^4*(γ + α*Log[y[x]]^2 + β*Log[y[x]]^4), y[x]]]

sol = With[{ϵ = 1/10000}, NDSolve[{y''[x] + 3/x*y'[x] == U, y[5000] == 0, y'[ϵ] == 0}, y, {x, ϵ, 5000}, 
Method -> {"BoundaryValues" -> {"Shooting", "StartingInitialConditions" -> {y[ϵ] == 2, y'[ϵ] == 0}}}]]

Plot[y[x] /. sol, {x, 1/10000, 200}, PlotRange -> All]

enter image description here

Regards,MI

POSTED BY: Mariusz Iwaniuk
Answer
1 month ago

Thank you for your response. Unfortunately, this is not the answer I am looking for with regards to recreating the algorithm stated in the paper and because of this the figures for the profile of the bounce solution are not the same.

POSTED BY: Benjamin Leather
Answer
1 month ago

Here is my updated code to see if this helps people. It still doesn't yield the required results but the algorithm follows the same procedure as in the paper.

zero = 10^-10;
inf = 10^9;
ref = inf - 10^-3;
\[Delta] = 0.01;
DU = Re[D[
    1/4 y[x]^4*(\[Gamma] + \[Alpha]*(Log[y[x]])^2 + \[Beta]*(Log[
           y[x]])^4), y[x]]];
eqn = y''[x] + 3 y'[x]/x == DU;
B = 0.070;
n = 0;
k = 1;
\[Epsilon] = 10^-10;
\[CapitalDelta] = 1;

\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\) = 1;
Bvalues = {};
absdiffvalues = {};

While[
\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\) > \[Epsilon],
  eqn0 = B + 
    B^3/8*(\[Gamma] + (\[Alpha]/2)*
        Log[B] + \[Alpha]*(Log[B])^2 + \[Beta]*(Log[
           B])^3 + \[Beta]*(Log[B])^4) x^2;
  initialcon1 = eqn0 /. x -> zero;
  initialcon2 = D[eqn0, x] /. x -> zero;
  sol = NDSolve[{y''[x] + 3 y'[x]/x == DU, y[inf] == 0, 
     y'[zero] == 0}, y[x], {x, zero, inf}, 
    Method -> {"BoundaryValues" -> {"Shooting", 
        "StartingInitialConditions" -> {y[zero] == initialcon1, 
          y'[zero] == initialcon2}}}];
  x2sol = x^2*y[x] /. (First@sol);
  \[CapitalDelta] = (x2sol /. x -> inf) - (x2sol /. x -> ref);

\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\) = 
   Abs[\[CapitalDelta]];
  AppendTo[absdiffvalues, 
\!\(\*OverscriptBox[\(\[CapitalDelta]\), \(~\)]\)];
  If[n >= 3,
   B0 = Part[Bvalues, n];
   B1 = Part[Bvalues, n - 1];
   B2 = Part[Bvalues, n - 2]
   ];
  If[\[CapitalDelta] > \[Epsilon],
   If[B1 <= B0 <= B2  , 
     k++;
     B += \[Delta]/(k*10);
     ];,
   B += \[Delta]/(k*10);
   ];
  If[\[CapitalDelta] < -\[Epsilon],
   If[B2 <= B0 <= B1  ,
     k++;
     B -= \[Delta]/(k*10);
     ];,
   B -= \[Delta]/(k*10);
   ];
  AppendTo[Bvalues, B];
  n++;
  If[Divisible[n, 500],
   Print[n];
   Print[\[CapitalDelta]];
   Print[B];
   Print[k];
   Print[initialcon1];
   Print[initialcon2];
   ];
  ];

Any help with this problem would still be greatly appreciated.

POSTED BY: Benjamin Leather
Answer
1 month ago

Group Abstract Group Abstract