Message Boards Message Boards

Use an algorithm to obtain numerical solutions for the EW vacuum 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 online. In the Appendix, "A Numerical computation of the bounce solution", it outlines an algorithm involving a shooting method to solve a boundary value problem.

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

Mp = 2.435*10^18;
α = 1.4*10^-5;
β = 6.3*10^-8;
γ = -0.013;

The appendix outlines the algorithm, where Eq.(2.23) is referring to the differential equation (and boundary conditions) stated above.

This is my attempt thus far:

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 = 2.3;
n = 0;
k = 0;
\[Epsilon] = 10^-10;
\[CapitalDelta] = 1;
\[CapitalLambda] = 1;
Bvalues = {};
absdiffvalues = {};
diffvalues = {};
\[Delta]values = {};
eqn0 = B + B^3/8 (γ + (α/2)*Log[B] + α*(Log[B])^2 + β*(Log[B])^3 + β*(Log[B])) x^2;



    While[\[CapitalLambda] > \[Epsilon] && n <= 500,
  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 = First[
    NDSolve[{y''[x] + 3 y'[x]/x == DU, y[zero] == initialcon1, 
      y'[zero] == initialcon2}, y[x], {x, zero, inf}, 
     Method -> {"ExplicitRungeKutta"}]];
  x2sol =  x^2*y[x] /. sol;
  \[CapitalDelta] = (x2sol /. x -> inf) - (x2sol /. x -> ref);
  \[CapitalLambda] = Abs[\[CapitalDelta]];
  AppendTo[diffvalues, \[CapitalDelta]];
  AppendTo[absdiffvalues, \[CapitalLambda]];
  If[n >= 3,
   B0 = Part[Bvalues, n];
   B1 = Part[Bvalues, n - 1];
   B2 = Part[Bvalues, n - 2];
   If[! (B0 > B1 > B2 || B2 > B1 > B0),
    k++;
    \[Delta] = \[Delta]/10;
    ];
   ];
  If[\[CapitalDelta] > \[Epsilon],
   B += \[Delta];
   ];
  If[\[CapitalDelta] < -\[Epsilon],
   B -= \[Delta];
   ];
  AppendTo[Bvalues, B];
  n++;
  If[Divisible[n, 500],
   Print[n];
   Print[\[CapitalDelta]];
   Print[B];
   Print[k];
   Print[\[Delta]];
   Print[initialcon1];
   Print[initialcon2];
   ];
  ];

Unfortunately it seems to get to a value of delta = -0.000129013 but then seems to not get any better no matter how many iterations it runs through. In the paper it seems as if it only takes around 100 iterations in this case to reach a suitable solution. Is there something I am missing, i.e. to search for appropriate initial guess, changing the numerical method?

Any help on this matter would be greatly appreciated.

POSTED BY: Benjamin Leather
Answer
3 months ago

I'm unsure but I think you should check that you've mixed arbitrary precision numbers with fixed precision ones (ie, 3 and 3.0).

If your running an algorithm that computes extra digits of precision per (loop/pass), you can't afford not to be "picky" about which your using when as far as your NDSolve it looks to be a simple diff'eq, but I have no time to do that work.

POSTED BY: John Hendrickson
Answer
3 months ago

Group Abstract Group Abstract