False vacuum bounce solution (curved space) shooting method problems

Question

Merry Christmas, I think I can still say it. I am back this time for my problem in a simpler case, without coupling. I have a problematic code now and I want your help. I want to solve numerically the below coupled equations:

$\tilde{\phi}''+\frac{3\tilde{\rho}'}{\tilde{\rho}}\tilde{\phi}'=\frac{dU}{d\phi}$

and

$\tilde{\rho}''=-\frac{8\pi}{3}μ^2\tilde\rho\Bigg(\tilde{\phi}'^2+\tilde{U}\Bigg)$

the primes are $d/dx$ and μ the mass scale $μ=Μ/Μ_{planck}$

with boundary conditions for $\tilde{\phi}$ to be: $\lim_{x\rightarrow \infty}\tilde{\phi}(\tilde{x})=\tilde{\phi}_{fv}\sim 1,\qquad \frac{d\tilde{\phi}}{d\tilde{x}}\Bigg|_{\tilde{x}=0}=0$

and for $\tilde{\rho}:$

$\tilde{\rho}(0)= 0,\qquad \tilde{\rho}'(0)= 1$

with the following potential

$ \tilde{U}(\tilde{\phi})=\frac18 (\tilde{\phi}^2-1)^2+\frac{\tilde{\epsilon}}{2}(\tilde{\phi}-1)$

I'm following the method of this paper https://arxiv.org/pdf/2205.03140v2.pdf (if you have questions about my coding go to the APPENDIX A.6 and check the method's description in the flat and curved space cases, it is quick 3 pages ), where they use Taylor expansion to transform the BCP to a numerically solvable ICP. I want to solve for small $\epsilon$ lets say 0.11.

Here it is. My amateur, super problematic code. I post it as detailed as I can in order to be clear to you.

ro = r[x]; fi = f[x]; (*ρ(x) and φ(x)*)
rod = D[ro, x];
rodd = D[ro, {x, 2}];
fid = D[fi, x];
fidd = D[fi, {x, 2}];
eps = 0.11 (*ε*); xmin = 0.01; xmax = 12; fv = 0.93; 
f0 = -1.05;(*α in the paper, I put it equal to true vacuum value ~ -1, 
because I want the graph to start from there*)
eq1 = fidd + (3 rod/ro) fid - (1/2) fi (fi^2 - 1) - eps/2 == 0;
eq2 = rodd + (8 Pi/ 3) mrat^2 ro (fid^2 + (1/8) (fi^2 - 1)^2 + (eps/2) (fi - 1)) == 0;
(*initial and boundary conditions*)
incon1 = f[xmin] == f0 + (1/8) ((1/2) f0 (f0^2 - 1) + eps/2) xmin^2;
incon2 = f'[xmin] == (1/4) ((1/2) f0 (f0^2 - 1) + eps/2) xmin;
incon3 = r[xmin] == xmin;
incon4 = r'[xmin] == 1;
bcon1 = f[xmax] == fv;
bcon2 = f'[xmin] == 0;
bcon3 = r[0] == 0;
bcon4 = r'[0] == 1;
solution = NDSolveValue[{eq1, eq2, bcon1, bcon2, bcon3, bcon4}, {fi, ro}, {x,xmin,xmax},
Method -> {"Shooting", "StartingInitialConditions" -> {incon1, incon2, incon3, incon4}}];
fsol = Part[Evaluate[f[x] /. solution], 1];
Plot[fsol, {x, xmin, xmax}, AxesLabel -> {x, φ}];
(* Manipulate Command for "f0" and "xmin-xmax"??*)

• I think I need Manipulate command (for f0 and xmin-xmax?), but I can't use it in my Mathematica because it is cracked. When I used it shows an error.

• For a solution I expect an instanton, which a good guess is $\phi(x)\sim \tanh(smthing)$. Check the graph ($\xi=0$ our case and $\eta$ is our $x$)

• When I run the code, stiff errors arise. I observed that in a number of shooting method questions here too, but I have no idea how to overcome these stiff problems.

• I have already checked the answer in the flat case question here with a little different boundary conditions:How to program efficient undershoot/overshoot , I can solve the flat space problem too. Also, I can solve the curved space problem with another method, I have a code solving it via minimization of the energy functional, but I have to make a specific hypothesis for the function ρ(x) in the metric. So I want to solve it in a more general way via the equations of motion posted here.

I THINK I COVERED EVERYTHING. GIVE ME YOUR LIGHTS PLEASE !!!! ITS IMPORTANT!! HAVE A NICE HOLIDAY TOO!!.

Any ideas how to fix the boundary condition phi[xmax]=fv???

I'm leaving the code in its initial form in the case something its going wrong with the mass scale- eps change.

Answer 1

I deleted most of my earlier answer after reading the paper, and deleted or folded earlier comments into this answer.

The authors point out that the problem isn't a Cauchy initial value problem (IVP), but a boundary value problem (BVP), so they're using a shooting method. For flat space, they keep the phi'[0] boundary condition, but replace the phi[infinity] = fv condition (which makes the original problem a BVP instead of an IVP) with phi[0] = a, 'a' being a parameter to be varied. You now have a Cauchy IVP.

Due to the singularity at the origin, they need to shift their initial values over to xmin, which they do via a Taylor expansion. This gives you a new set of IVP's at xmin. You can now solve the resulting problem as a function of parameter 'a', then vary 'a' to match the condition at phi[infinity] = fv to desired precision.

For the curved space case, they follow a similar approach, but also shift the r variable over to xmin using a similar Taylor expansion, which they don't show the derivation for, just the result. This is how they get r[xmin] = xmin, etc. So their "shooting method" is partially implicit in the derivation they performed, except for the dependence on a parameter 'a'. I.e. they've restructured the BVP as a new IVP, and shifted the IVP equations over to xmin instead of 0 to avoid the singularity.

One thing that perhaps wasn't stressed enough was that this shift to xmin needs to be applied to both the boundary and initial value conditions. Three of the four are the same - only the condition on phi at infinity is swapped out to convert from a BVP to an IVP. The shooting method will attempt to test its intermediate results against the supplied boundary conditions, and if they remain singular, the process won't converge.

Once these steps are performed, you now have a Cauchy problem which can be solved using the equations and the initial values they derived at xmin, using ParametricNDSolve without any shooting method selected, and parameterized by 'a'. The result can be passed to FindRoot with an initial guess for 'a', and solved to match the boundary condition at "infinity" for phi, i.e. match to fv.

You could also try using the shooting method option for NDSolve instead, using the shifted initial and boundary conditions. The shooting tutorial mentions that FindRoot is used internally to seek convergence. There is an advantage to using ParametricNDSolve/FindRoot separately however; if you run into problems, it's easier to see what's going wrong. As long as ParametricNDSolve runs successfully, it's possible to use Plot/Table/etc. to obtain information about what your solution looks like even if FindRoot fails, as will be seen below.

After running an earlier version of this code, it was found that no solution was possible with the specified parameters. After further discussion/investigation with the OP, we found that the variable 'eps' was being used in two different ways. It controls the shape of the potential function the OP uses, but the authors of the referenced appendix in the link use it as the ratio of the scalar field mass to the Planck mass. I've changed the name of this parameter to massRatio in the code below.

It also became apparent that it would be useful to be able to try out different equations for the potential function. To simply this process, and avoid having to re-code equations by hand, I modified the earlier code to allow defining potential functions externally and passing them into the body of the solver. The necessary derivatives are computed internally and used to generate the appropriate equations for each potential.

This is the original potential function:

u["baseline"][eps_][phi_[x_]] := 
  1/8 (phi[x]^2 - 1)^2 + eps/2 (phi[x] - 1);

The [phi[x]] term at the end allow various derivatives to be performed. The [eps_] is for the constant parameter used in the equation; additional parameters may be added here if needed. The u["baseline"] is a way to label trial potentials. Any short descriptive string can serve as a label. It's also possible to omit the label and just use u0, u1, etc.

Here's the main function:

solve[phi_, r_, x_, xmin_, xmax_, a_, a0_, fv_, mrat_, u_] :=
  Module[{dua, eq},
   (* Compute the derivative of the potential wrt phi evaluated at \
phi=a *)
   dua = D[u, phi[x]] /. phi[x] -> a;
   eq["phi"] = phi''[x] + 3 r'[x]/r[x] phi'[x] - D[u, phi[x]] == 0; 
   eq["r"] = r''[x] + 8 Pi/3 mrat^2 r[x] (phi'[x]^2 + u) == 0;
   eq["ic"] = {phi[xmin] == a + 1/8 dua xmin^2, 
     phi'[xmin] == 1/4 dua xmin, r[xmin] == xmin, r'[xmin] == 1};
   ParametricNDSolve[
    Flatten[{eq["phi"], eq["r"], eq["ic"]}], {phi, r}, {x, xmin, 
     xmax}, {a}]
   ];

Call it with

 {xmin,xmax,a0,fv,mrat,eps} = {.01,12,-1.05,.93,.25,.11};
    pnds = solve[phi, r, x, xmin, xmax, a, a0,fv,mrat,
u["baseline"][eps][phi[x]]]

Then solve for the optimal 'a' value:

a = (a /. FindRoot[phi[a][xmax] == fv /. pnds, {a, a0}])

and plot the result (I scaled r down to better fit both curves on the plot and added the Epilog to the plot at 'fv' to serve as a reference marker.):

Plot[Evaluate[{phi[a][x], r[a][x]/10} /. pnds], {x, xmin, xmax}, 
         PlotLegends -> {"phi[x]", "r[x]/10"},Epilog ->  Line[{{0,.93},{xmax,.93}}]]

There was also a request for a way to run this code with Manipulate for various inputs. In order to get acceptable output whether FindRoot converges or not, I wrapped the call with Enclose/ConfirmQuiet, which will throw a Failure box to the surrounding enclose if messages are generated.

  Manipulate[
 Module[{pnds, phi, r, x, a}, 
  pnds = solve[phi, r, x, xmin, xmax, a, a0, fv, mrat, 
    u["baseline"][eps][phi[x]]];
  {Plot[Evaluate[{phi[a0][x], r[a0][x]/10} /. pnds], {x, xmin, xmax}, 
    PlotLegends -> {"phi[x]", "r[x]/10"}, 
    Epilog -> Line[{{0, fv}, {xmax, fv}}]], 
   Enclose[ConfirmQuiet[
     a = (a /. FindRoot[
         phi[a][xmax] == fv /. pnds, {a, a0}])]]}], {{xmin, .01}, .01,
   1}, {{xmax, 12}, 1, 31}, {{a0, -1.05}, -1.1, 0}, {{fv, 0.93}, 
  0.093, .93},
 {{eps, .11}, -.2, .2},
 {{mrat, .25}, .0001, 1.}]

Here's a sample output:

[![enter image description here][1]][1]

It's also possible to verify that the function does in fact give valid solutions. Write the equations/boundary/initial conditions as lhs == 0, choose a value for the parameter a, generate interpolating functions, and plug these into the lhs of each equation. The equations should give zero for all x in the range, and the initial/boundary conditions at their endpoints. Here's code to test the result:

Manipulate[
 Module[{pnds, phi, r, x, a, eq1, eq2, ic, ifn, subst, sol},
  pnds = 
   solve[phi, r, x, xmin, xmax, a, a0, fv, mrat, 
    u["baseline"][eps][phi[x]]];
  eq1 = phi''[
     x] + (3 r'[x]/r[x]) phi'[x] - (1/2) phi[x] (phi[x]^2 - 1) - eps/2;
  eq2 = r''[
     x] + (8 Pi/3) eps^2 r[
      x] (phi'[x]^2 + (1/8) (phi[x]^2 - 1)^2 + (eps/2) (phi[x] - 1));
  ic = {phi[xmin] - (a + (1/8) ((1/2) a (a^2 - 1) + eps/2) xmin^2), 
    phi'[xmin] - ((1/4) ((1/2) a (a^2 - 1) + eps/2) xmin), 
    r[xmin] - xmin, r'[xmin] - 1,
    {phi[xmax] - fv}};
  ifn = {phi[a0], r[a0]} /. pnds;
  subst = Thread[{phi, r} -> ifn];
  sol = {{eq1, eq2}, ic} /. a -> a0 /. subst;
  Column[{Plot[sol[[1, 1]], {x, xmin, xmax}], 
    Plot[sol[[1, 2]], {x, xmin, xmax}], sol[[2]]}]
  ]
 , {{xmin, .01}, .01, 1}, {{xmax, 12}, 1, 31}, {{a0, -1.05}, -1.1, 
  0}, {{fv, 0.93}, 0.093, .93},
 {{eps, .11}, -.2, .2},
 {{mrat, .25}, .0001, 1.}]

I'm seeing numerical round off errors on the order of 5x10^-6 overall, so the solutions appear valid.