Non-minimal coupling ξ - numerical bounce solutions (shooting method, false vacuum decay)

I have the potential below: $$V(\phi)=-\frac14 a^2(3b-1)\phi^2+\frac12 a(b-1)\phi^3+\frac14 \phi^4 +a^4c$$ This potential has 2 minima, the false vacuum $$\phi_f=0$$ which tunnels to the global minimum, the so-called true vacuum, at $$\phi_t=a$$. The parameter $$c$$ at the potential is responsible for the character of our initial false vacuum ($$c>0$$ for transitions from de Sitter (dS) space and $$c=0$$ for transtions from Minkowski space). Here we will take $$a=1$$ and $$b=1/10$$.

I want to solve numerically the two-equations problem:

$$\ddot\phi+3\frac{\dot\rho}{\rho}\dot\phi-\xi\phi R=\frac{dV}{d\phi}$$ where $$\xi$$ is the non-minimal coupling and the dots are $$d/d\tau$$, and

$$\ddot\rho=\frac{\kappa\rho}{3(1-\kappa\xi\phi^2)}\bigg(-\dot\phi^2-V+3\xi\bigg(\dot\phi^2+\ddot\phi\phi+\dot\phi\phi\frac{\dot\rho}{\rho}\bigg)\bigg)$$

with boundary conditions: $$\dot{\phi}(0)=\dot{\phi}(\tau_{max})=0,\\ \rho(0)=0,\\ \rho(\tau_{max})=0\quad\text{(for dS false vacuum)},\\ \rho(\tau_{max})=\rho_{max}\ne 0\quad\text{(for Minkowski false vacuum)}$$

Also, one more important boundary condition, not refered in the paper but it is obvious at the numericals, is the following: $$\phi(\infty)=\phi(\tau_{max})=\phi_f=0$$

Here $$R$$ is the curvature scalar and it is defined as: $$R=-6\Bigg(\frac{\ddot\rho\rho+\dot\rho^2-1}{\rho^2}\Bigg)=\frac{\kappa\rho}{(1-\kappa\xi\phi^2)}\bigg(\dot\phi^2+4V-6\xi\bigg(\dot\phi^2+\ddot\phi\phi+3\dot\phi\phi\frac{\dot\rho}{\rho}\bigg)\bigg)$$

I tried to solve this problem by adjusting to my case, lets say $$\xi=0.1$$ an algorithm made by jdp in my previous question about the curved case problem ($$\xi=0$$). But, no light in the horizon, my free-15days trial has been expired, maybe i have to renew it, my cracked edition does not support the manipulate command. I post the following for your help below:

and for the bubble radius $$\rho$$:

• Also, my bad coding skills for a dS transition ($$c=0.05$$) with e and e0 taking place of jdp's a and a0 and x is $$\tau$$ and $$\kappa=8\pi G=8\pi M_{Pl}^2$$. Ηere I took $$\kappa=0.03$$ from another paper, maybe a Manipulate command could fix this parameter:
u[a_][b_][c_][phi_[x_]]:=-1/4a^2(3b-1)phi[x]^2+1/2a(b-1)phi[x]^3+1/4phi[x]^4+a^4c;
R=k/(1-k*xi*phi[x]^2)*(phi'[x]^2+4u-6xi(phi'[x]^2+phi[x]phi''[x]+3phi'[x]*phi[x]*r'[x]/r[x]));

solve[phi_, r_, x_, xmin_, xmax_, e_, e0_, fv_, k_,xi_,R_, u_] :=
Module[{due, eq},

due = D[u, phi[x]] /. phi[x] -> e;
eq["phi"] = phi''[x] + 3 r'[x]/r[x] phi'[x] -xi*phi[x]*R- D[u, phi[x]] == 0;
eq["r"] = r''[x] + k*r[x]/(3(1-k*xi*phi[x]^2)) (phi'[x]^2 + u-3xi(phi'[x]^2+phi''[x]phi[x]+phi'[x]phi[x]r'[x]/r[x])) == 0;
eq["ic"] = {phi[xmin] == e + 1/8 due xmin^2,
phi'[xmin] == 1/4 due xmin, r[xmin] == xmin, r'[xmin] == 1};
ParametricNDSolve[
Flatten[{eq["phi"], eq["r"], eq["ic"]}], {phi, r}, {x, xmin,
xmax}, {e}]
];

{xmin,xmax,e0,fv,k,xi,a,b,c} = {.01,15,.95,0,.03,.1,1,.1,0.05};
pnds = solve[phi, r, x, xmin, xmax, e, e0,fv,k,xi,R,
u[a][b][c][phi[x]]]

e = (e /. FindRoot[phi[e][xmax] == fv /. pnds, {e, e0}])

Plot[Evaluate[{phi[e][x], r[e][x]/10} /. pnds], {x, xmin, xmax},
PlotLegends -> {"phi[x]", "r[x]/10"}]



I think I covered everything. Sorry for the great extent of my question. Please help me if you can. HEEEEEEELP!!!

• From the picture above it looks like xmin=0, xmax=35. Is it correct? Mar 9 at 3:09
• yeap, but to avoid singularity your start from a very small lets say xmin=0.0001 Mar 9 at 12:32

We don't need parameter $$\kappa$$ since we can express scalar curvature as follows R = -6/r[x]^2 (r''[x] r[x] + r'[x]^2 - 1). With this definition numerical model generates some solution

{xmin, xmax, e0, fv, k, xi, a, b, c} = {.01, 15, .95, 0, .03, .1,
1, .1, 0.05};
u = -1/4 a^2 (3 b - 1) phi[x]^2 + 1/2 a (b - 1) phi[x]^3 +
1/4 phi[x]^4 + a^4 c;
R = -6/r[x]^2 (r''[x] r[x] + r'[x]^2 -
1)(*k/(1-k*xi*phi[x]^2)*(phi'[x]^2+4u-6xi(phi'[x]^2+phi[x]phi''[\
x]+3phi'[x]*phi[x]*r'[x]/r[x]))*);

due = D[u, phi[x]] /. phi[x] -> e;
eq1 = phi''[
x] + (3 r'[x]/r[x] phi'[x] - xi*phi[x]*R - D[u, phi[x]]) == 0;
eq2 = r''[x] +
k*r[x]/(3 (1 - k*xi*phi[x]^2)) (phi'[x]^2 + u -
3 xi (phi'[x]^2 + phi''[x] phi[x] +
phi'[x] phi[x] r'[x]/r[x])) == 0;
ic = {phi[xmin] == e + 1/8 due xmin^2, phi'[xmin] == 1/4 due xmin,
r[xmin] == xmin, r'[xmin] == 1};

sol = ParametricNDSolve[
Flatten[{eq1, eq2, ic}], {phi, r}, {x, xmin, xmax}, {e}];

sol1 = FindRoot[phi[e][xmax] == fv /. sol, {e, e0}]


Visualization

Plot[Evaluate[{phi[e][x], r[e][x]/10} /. sol1 /. sol], {x, xmin,
xmax}, PlotLegends -> {"phi[x]", "r[x]/10"},
PlotLabel -> Row[{"e = ", e /. sol1}]]


Update 1. Actually we can compute solution with k=1 to compare with paper. This solution is different from above. It is not clear how physicists can describe this discrepancies.

 {xmin, xmax, e0, fv, k, xi, a, b, c, k} = {.01, 15, .95,
0, .03, .1, 1, .1, 0.05, 1};
u = -1/4 a^2 (3 b - 1) phi[x]^2 + 1/2 a (b - 1) phi[x]^3 +
1/4 phi[x]^4 + a^4 c;
R = k/(1 - k*xi*phi[x]^2)*(phi'[x]^2 + 4 u -
6 xi (phi'[x]^2 + phi[x] phi''[x] + 3 phi'[x]*phi[x]*r'[x]/r[x]));

due = D[u, phi[x]] /. phi[x] -> e;
eq1 = phi''[
x] + (3 r'[x]/r[x] phi'[x] - xi*phi[x]*R - D[u, phi[x]]) == 0;
eq2 = r''[x] +
k*r[x]/(3 (1 - k*xi*phi[x]^2)) (phi'[x]^2 + u -
3 xi (phi'[x]^2 + phi''[x] phi[x] +
phi'[x] phi[x] r'[x]/r[x])) == 0;
ic = {phi[xmin] == e + 1/8 due xmin^2, phi'[xmin] == 1/4 due xmin,
r[xmin] == xmin, r'[xmin] == 1};

sol =
ParametricNDSolve[
Flatten[{eq1, eq2, ic}], {phi, r}, {x, xmin, xmax}, {e}];

sol1 = FindRoot[phi[e][xmax] == fv /. sol, {e, e0}]

(*Out[]= {e -> 0.949034}*)

Plot[
Evaluate[{phi[e][x], r[e][x]/10} /. sol1 /. sol], {x, xmin, xmax},
PlotLegends -> {"phi[x]", "r[x]/10"},
PlotLabel -> Row[{"e = ", e /. sol1}]]


• I used the second definition of R because of this saying in the paper In our definition of R second power of r appears in the denominator so there is a possibility of divergence here. Therefore, it is much more convenient for numerical calculations to express the scalar curvature using the Friedman equations as: (go to eq.2.10 to check it). So i assumed the 2nd def is more suitable for my numerical. But the solution looks fine, i'll check it and i comment on it Mar 9 at 12:35
• Thanks a lot for your answer!! I'll try to do some manipulate staff too like @jdp Mar 9 at 12:41
• And we've got a possibility for a second ρ=0 in the dS transition, (fig 5 left panel). I don't have access to mathematica for the moment, when I'll have got I check your solution, xmin and xmax, the shape of ρ in the case where c=0.05 must be like fig.5 left panel around x=25. Also, I'll do some manipulation like jdp did in the first question Mar 9 at 15:32
• @JenniferDerleth See Update 1 to my answer. Mar 10 at 3:03
• yeah, this looks like a de Sitter to Minkowski space transition, the shape of ρ is fine, thanks, it seems good. Today i'll have the time to check it too – Mar 10 at 14:29