Semi-analytic techniques for calculating bubble wall profiles

Abstract

We present semi-analytic techniques for finding bubble wall profiles during first order phase transitions with multiple scalar fields. Our method involves reducing the problem to an equation with a single field, finding an approximate analytic solution and perturbing around it. The perturbations can be written in a semi-analytic form. We assert that our technique lacks convergence problems and demonstrate the speed of convergence on an example potential.

1 Introduction

The decay of a false vacuum is a complex problem with numerous applications in cosmology [1,2,3,4,5,6] and is particularly important in the study of baryogenesis [7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26] (although there are mechanisms for producing the baryon asymmetry that do not require calculating the decay of the false vacuum [27,28,29,30,31,32,33]). Calculating tunneling rates is also an important problem in the study of vacuum stability [34,35,36,37,38,39,40].

While the physics of the tunneling process is qualitatively well understood [41], quantitatively it is a complicated problem that involves solving a set of highly nonlinear coupled differential equations usually requiring a numerical solution. The two techniques that are most commonly used to solving the tunneling problem are path deformation [42, 43] and minimizing the integral of the squared equations of motion for a set of parametrized functions [44, 45], although other methods also exist [46]. Exact solutions only exist for specialized cases [47, 48].

In this paper we offer a new approach by solving the tunneling problem semi-analytically. First, we give an analytic solution to an ansatz for an arbitrary potential. Since the ansatz is only an approximate solution, in the next step we derive a perturbative expansion that converges to the exact solution. For each term in the perturbative series we provide a semi-analytic solution. Since our starting potential is arbitrary, and because the convergence of the perturbative expansion is independent of the potential, our method is completely general.

To derive the ansatz we take advantage of the fact that the multi-field problem can be approximated by finding the solution to a single-field potential. The approximate tunneling solution can be found along the field direction that connects the true and false vacua. This single-field potential can be solved in terms of a single parameter. To improve the initial ansatz we compute correction functions to the ansatz in a manner analogous to Newton’s method of finding roots. The result is a perturbative series of corrections that are expected to converge quadratically.

The differential equations that define these corrections can be solved analytically in terms of eigenvalues of the mass matrix and a function of the initial ansatz. In doing so we use techniques that were recently employed to analytically solve number densities across a bubble wall [49]. Although the technique has elements in common with Newton’s method it does not share its trouble with null derivatives giving divergent corrections or oscillations around the solution. We also argue that the other problems with Newton’s method are generically not relevant to our method.

The layout of this paper is as follows. In Sect. 2 we give a brief overview of the false vacuum problem. In Sect. 3 we develop an ansatz form that approximately solves a general variety of multi-field potentials with a false vacuum, where the potential is specified by a single parameter. In Sect. 4 we derive the perturbative corrections to the ansatz forms and discuss the convergence. In Sect. 5 we use this method to solve a problem which can be directly compared with the literature. Concluding remarks are given in Sect. 6.

2 Fate of the false vacuum

Consider a potential of multiple scalar fields $V(\phi _i )$ with at least two minima. The trivial solution to the classical equations of motion is stationary extremizing the potential. This solution typically gives the field a non-zero vacuum expectation value and is responsible for giving standard model particles their mass via electroweak symmetry breaking. The other, less obvious solution is one where the fields continuously vary from one minimum to another. In this case, the false vacuum decays into the true vacuum via tunneling processes, and is termed the ‘bounce solution’ [41]. If this is achieved within a first order phase transition, regions of the new vacuum appear and expand as bubbles in space. In this paper we are interested in calculating the profile of the bubble, that is, the spacetime dependence of the bubble nucleation.

The spatial bubble profile is obtained by extremizing the Euclidean action

$$\begin{aligned} S_E = \int \mathrm d^d x\,\left[ \frac{1}{2}\left( \partial _\mu \varphi _i\right) ^2 + V(\varphi _i) \right] \end{aligned}$$

(1)

where $d=4$ for zero temperature tunneling and $d=3$ for finite temperature tunneling relevant to cosmological phase transitions. The nucleation rate per unit volume is

$$\begin{aligned} \Gamma = A(T) e^{- \frac{S_E}{T}} \end{aligned}$$

(2)

where A(T) is a temperature dependent prefactor proportional to the fluctuation determinant, T is the temperature and $S_E$ is the euclidean action for the bounce solution which satisfies the classical equations of motion. In the case of a spherically symmetric bubble the classical equations of motion are

$$\begin{aligned} \frac{\partial ^2 \varphi _i}{\partial \rho ^2} + \frac{2}{\rho } \frac{\partial \varphi _i}{\partial \rho } - \frac{\partial V}{\partial \varphi _i} = 0 \end{aligned}$$

(3)

and the bounce solution satisfies the conditions $\varphi _i(0)\approx v_i^\mathrm {true}$, $\varphi _i(\infty ) = v_i^\mathrm {false}$ and $\varphi _i'(0) = 0$.¹ Here $\rho$ is the ordinary 3D spherical coordinate, as we are considering finite temperature, and $v_i^\mathrm {true}$ and $v_i^\mathrm {false}$ are the vacuum expectation values of the field $\varphi _i$ in the true and false vacua, respectively. The equations of motion resemble the classical solution of a ball rolling in a landscape of shape-V with $\rho$ playing the role of time, but including a $\rho$-dependent friction term.

3 Approximate solution to the multi-field potential

3.1 Reducing to a single-field potential

The bounce solution can be approximated by the bounce solution of a single differential equation as follows [42, 43]. First make a shift of fields such that the false vacuum is at the origin in field space. The true vacuum is then at $v \hat{\phi }_1$ where $\hat{\phi }_1$ is a unit vector that points in the direction of the true vacuum. Then define a complete set of unit vectors orthogonal to $\hat{\phi }_1$ and rewrite the potential in the rotated basis $V(\varphi _1 , \varphi _2, \ldots ) \mapsto V(\phi _1, \phi _2, \ldots )$. Then consider the potential only in the $\hat{\phi }_1$ direction between the minima, $V(\phi _1, 0,\ldots )$. One can then solve the single equation of motion

$$\begin{aligned} \frac{\partial ^2 \phi _1}{\partial \rho ^2} + \frac{(d-1) }{\rho } \frac{\partial \phi _1}{\partial \rho } - \frac{\partial V(\phi _1,0,\ldots )}{\partial \phi _1} = 0 \end{aligned}$$

(4)

to derive an initial ansatz that approximately solves the full classical equations of motion. Let us therefore turn our attention to the most general renormalizable tree level potential with a single field,

$$\begin{aligned} V(\varphi ) = M^2 \varphi ^2 + b \varphi ^3 + \lambda \varphi ^4. \end{aligned}$$

(5)

An approximate expression for the effective action of a similar potential was found in Ref. [50]. The above potential can be rescaled $\phi = \varphi _\text {min} \varphi$ where $\varphi _\text {min}$ is the global minimum of the above potential. Then the rescaled potential has a global minimum at $\phi =1$. To ensure that the effective action is unaffected by this rescaling, we also make the replacement $\rho \mapsto \varphi _\text {min}\rho$. We paramatrize the rescaled potential as²

$$\begin{aligned} V(\phi ) = \frac{(4 \alpha -3 )}{2} E \phi ^2 +E \phi ^3 - \alpha E \phi ^4. \end{aligned}$$

(6)

Tunneling between two vacua requires the existence of a potential barrier or “bump” separating the two minima. As parametrized in Eq. (6), this type of barrier can only exist if $E<0$ and $\alpha \in (0.5,0.75)$.³ To illustrate this point, we present in Fig. 1 the potential in $\phi$ given in Eq. (6) for both the edge choices of $\alpha$ and the mean allowed choice, using several choices of E. One can see that, for $\alpha =\frac{1}{2}$, we have exactly the Mexican hat potential (albeit shifted to $\phi =0.5$) with degenerate minima, and for $\alpha = \frac{3}{4}$, there is no potential barrier between false and true minima.

Fig. 1

We present our rescaled scalar potential parametrized by E and $\alpha$ given in Eq. (6). Each panel displays $-E=1,2,5,10,20$. The top and bottom panels are the edge choices of $\alpha$, with $\alpha =\frac{1}{2}$ on top and $\alpha =\frac{3}{4}$ on the bottom. The middle panel has the mean value of $\alpha =\frac{5}{8}$

3.2 Developing ansatz solutions

In deriving an approximate solution to the potential in Eq. (6), we first note that the effective potential is proportional to E. Thus, one can factor |E| out of the equations of motion by further rescaling $\rho \mapsto \rho / \sqrt{|E|}$. Then the equations of motion only depend on $\alpha$.

Under the scaling we have introduced,

$$\begin{aligned} {\left\{ \begin{array}{ll} \varphi \mapsto \phi = \varphi _\text {min} \varphi \\ \rho \mapsto \dfrac{\varphi _\text {min}}{\sqrt{|E|}}\rho \end{array}\right. } \end{aligned}$$

(7)

the Euclidean action becomes

$$\begin{aligned} S_E = 4\pi \frac{\phi _m^3}{\sqrt{|E|}} \int \mathrm d\rho \, \rho ^2 \left[ \left( \frac{\partial \phi }{\partial \rho }\right) ^2 - \widetilde{V}(\phi )\right] \end{aligned}$$

(8)

where⁴ $\widetilde{V} \equiv V/|E|$, and we have integrated over the angular variables assuming isotropy. The integral in Eq. (8) must only depend on $\alpha$, as in

$$\begin{aligned} S_E = 4 \pi \frac{ \phi _m ^3}{\sqrt{|E|}} f(\alpha ). \end{aligned}$$

(9)

Meanwhile, we will approximate the rescaled field itself with the well-known “kink” solution [44, 45]

$$\begin{aligned} \phi \approx \frac{1}{2} \left( 1 - \tanh \left[ \frac{\rho - \delta (\alpha )}{L_w(\alpha )}\right] \right) \end{aligned}$$

(10)

parametrized by the offset $\delta$ and the bubble wall width $L_w$. Thus it remains to determine the $\alpha$-dependent functions $\delta (\alpha )$, $L_w(\alpha )$ from the kink solution, and $f(\alpha )$ in the Euclidean action.

We first evenly sample values of $\alpha$ within (0.5, 0.75), then numerically solve the full bubble profile using conventional techniques. Next, for each value of $\alpha$, we fit the kink solution given in Eq. (10) to the full solution, extracting $L_w$ and $\delta$. Lastly, we numerically integrate to find f in the Euclidean action. This results in a tabulation of values for $L_w$, $\delta$, and f, for each value of $\alpha$. Using the apparent $\alpha$ dependence and intuition from our parametrization of the potential, we find ansatz functional forms in terms of $\alpha$ for each of these parameters.

The offset $\delta$ should diverge at the boundaries $\alpha =0.5$ and $\alpha =0.75$, and is found to be quite small otherwise. It also appears to have approximate odd parity about the mean allowed value of $\alpha =0.625$. We modeled this with odd powers of non-removable poles at the boundaries of $\alpha$, along with an offset and a linear correction about the mean:

$$\begin{aligned} \delta (\alpha ) \approx \delta _0 + k\left( \alpha -\frac{5}{8}\right) + \sum _{n=1}^2 a_n \left[ \frac{\alpha -\frac{5}{8}}{\left( \alpha -\frac{1}{2}\right) \left( \alpha -\frac{3}{4}\right) }\right] ^{(2n-1)}.\nonumber \\ \end{aligned}$$

(11)

We then fit these parameters using the tabulated values.

The bubble wall width $L_w$ is dominated by two asymptotes. It diverges at $\alpha =\frac{3}{4}$, and become small as $\alpha \rightarrow \frac{1}{2}$. We used this form to model the asymptotic behavior,

$$\begin{aligned} L_w(\alpha ) \approx \ell _0 \left[ \left( \alpha -\frac{1}{2}\right) ^r + \frac{c}{\left|\alpha -\frac{3}{4}\right|^s}\right] . \end{aligned}$$

(12)

As before, the normalization $\ell _0$, the two exponents r and s, and the coefficient c are fit using the tabulated values from the full numerical calculation. Interestingly, we find almost exactly that $r=18$. (The full fitted parameters are given in Table 1.)

Table 1 The fitted values for the parameters that define the approximate ansatz functions for f which is used in the Euclidean action, the bubble wall width $L_w$, and the offset $\delta$ from the kink solution

The Euclidean action determined by $f(\alpha )$ diverges at $\alpha =\frac{1}{2}$ and is zero at $\alpha =\frac{3}{4}$. This is modeled by

$$\begin{aligned} f(\alpha ) =f_0 \frac{\left|\alpha -\frac{3}{4}\right|^{p}}{\left|\alpha -\frac{1}{2}\right|^{q}} \end{aligned}$$

(13)

where only a normalization parameter and exponents need to be fitted.

In Table 1 we present all the fitted values that go into these ansatz approximate forms. The numerical values computed for these functions as well as the resulting fits are given in Fig. 2. We did not estimate uncertainties in the full numerical calculations nor in the fitted parameters, though in principle this could be done. Thus we are not able to compute rigorous measures of the goodness of fits. As these fits are only used to form a base ansatz solution which then receives perturbative corrections, such an undertaking lies outside the scope of this work. We do, however, provide in Fig. 2 the residuals between the fitted curves and tabulated values.

Fig. 2

We present the fits to the offset $\delta$ (top panel) and the bubble wall width $L_w$ (middle panel) of the kink solution, and the integral f (bottom panel) appearing in the Euclidean action. The numerically computed values are presented along with the fitted curves, as well as the residuals

We note that |E| scales as $|b \phi ^3|$ so $S_E/T$ scales as $\frac{\phi _m}{T} \sqrt{\frac{\phi _m}{b}}$, where b is the cubic coupling of the unscaled field, as in Eq. (5). Also b controls the height of the barrier separating the two minima.

4 Perturbative solution

In the previous section, we developed fitted curves to estimate the parameters of the well known kink solution. In this section, we will take advantage of rescaling to compute convergent perturbative corrections. The process is largely analogous to Newton’s method for finding roots of functions. Here, we iteratively determine functional corrections to the ansatz form.

4.1 Perturbative corrections to the ansatz

We first note that along the trajectory in field space from the false vacuum to the true vacuum, the magnitude of any of the fields in $\phi = \{\phi _i (\rho )\}$ generically does not exceed the distance between the two minima. That is,

$$\begin{aligned} |\phi _i(\rho )|\lesssim |v_\text {true} - v_\text {false} |. \end{aligned}$$

(14)

If we rescale our fields as described in Sect. 3, the distance between the ansatz and the actual solution is bounded by 1, but is usually much smaller than 1. (This is illustrated with the concrete example presented in Sect. 5.) Let us call the ansatz to $\phi _i$, $A_i$ with correction $\varepsilon _i$, so that $\phi _i = A_i + \varepsilon _i$. Applying this to Eq. (3) yields

$$\begin{aligned}&\frac{\partial ^2 A_i}{\partial \rho ^2} + \frac{\partial ^2 \varepsilon _i}{\partial \rho ^2} + \frac{2}{\rho } \frac{\partial A_i}{\partial \rho } + \frac{2}{\rho } \frac{\partial \varepsilon _i}{\partial \rho } \nonumber \\&\quad = \left. \frac{\partial V(\phi )}{\partial \phi _i} \right| _{A} + \sum _j \left. \frac{\partial ^2 V(\phi )}{\partial \phi _i \partial \phi _j} \right| _{A} \!\!\!\!\varepsilon _j \end{aligned}$$

(15)

where $A\equiv \left\{ A_i\right\}$. We can then rearrange the above to separate terms that depend only on the ansatz forms from those involving the unknown correction functions, $\varepsilon _i$. This leaves us with a set of coupled inhomogeneous differential equations for $\varepsilon _i$:

$$\begin{aligned} \frac{\partial ^2 \varepsilon _i}{\partial \rho ^2} + \frac{2}{\rho }\frac{\partial \varepsilon _i}{\partial \rho } - \left. \frac{\partial ^2 V(\phi )}{\partial \phi _i \partial \phi _j} \right| _{A} \!\!\!\!\varepsilon _j = B_i(\rho ). \end{aligned}$$

(16)

Here the functions $B_i(\rho )$ are the inhomogeneous part of the differential equations for $\varepsilon _i$, and they are given by

$$\begin{aligned} B_i(\rho ) \equiv \left. \frac{\partial V(\phi )}{\partial \phi _i} \right| _{A} \!\!\!\!- \frac{\partial ^2 A_i}{\partial \rho ^2} - \frac{2}{\rho } \frac{\partial A_i}{\partial \rho }. \end{aligned}$$

(17)

One can see that the value of the functions $B_i$ represents how well the ansatz forms solve the equations of motion. This can be seen not only as the definition of $B_i$ are the equations of motion where the fields are taken to be $A_i$, but also because if $B_i$ were zero, then the differential equations for the corrections to the ansatz $\varepsilon _i$ would become homogeneous.

We can linearize and approximately solve these differential equations analytically by approximating the mass matrix by a series of step functions with a correction which we will use to form the convergent series of perturbations. Considering only the homogeneous case ($B_i=0$), Eq. (16) is solved by introducing $\varepsilon$ of the form

$$\begin{aligned} \varepsilon \sim \frac{ze^{\lambda \rho }}{\rho } \text { .} \end{aligned}$$

(18)

This will reduce the equation to an eigenvalue problem in the mass matrix. It is therefore convenient to define the homogeneous solutions with the notation

$$\begin{aligned} \varepsilon ^h_{ik} = \frac{z_{ik}e^{\lambda _k\rho }}{\rho }\end{aligned}$$

(19)

where the index i refers to field $\phi _i$, and the index k will run over the n eigenvalues of the mass matrix, and k is not summed over. Substituting $\varepsilon ^h_{ik}$ for $\varepsilon _i$ in Eq. (16) with $B_i=0$ yields

$$\begin{aligned} \sum _j M_{ij} z_{jk} = \lambda _k^2 z_{ik} \end{aligned}$$

(20)

where M is the mass matrix

$$\begin{aligned} M_{ij} = \left. \frac{\partial ^2{V}(\phi )}{\partial \phi _i\partial \phi _j} \right| _A \!\!\text { .} \end{aligned}$$

(21)

In this way, $z_{ik}$ is the ith element of the eigenvector of M that has eigenvalue $\lambda _k^2$. It is necessary, however, to use both positive and negative roots, $\pm \lambda _k$. Thus we must further introduce $\tilde{\lambda }_j$ and $\tilde{z}_{ij}$ with

$$\begin{aligned}&\tilde{\lambda }_1 = \lambda _1 \text {, } \tilde{\lambda }_2 = -\lambda _1 \text {, } \tilde{\lambda }_3 = \lambda _2 \text {, } \ldots \end{aligned}$$

(22)

$$\begin{aligned}&\tilde{z}_{i1} = z_{i1} \text {, } \tilde{z}_{i2} = z_{i1} \text {, } \tilde{z}_{i3} = z_{i2} \text {, } \tilde{z}_{i4} = z_{i2} \text {, } \ldots \end{aligned}$$

(23)

so that in $\tilde{\lambda }_j$ and $\tilde{z}_{ij}$, the index $j=1,2,\ldots ,2n$. Finally, we use the techniques described in [49] to arrive at the particular solution,

$$\begin{aligned} \varepsilon _i^\lessgtr = \sum _{j=1}^{2n}\sum _{k=1}^n \tilde{z}_{ij} \frac{e^{\tilde{\lambda }_j}\rho }{\rho }\left( \int _0^\rho t e^{-\tilde{\lambda }_j t} h_{jk} B_k^\lessgtr (t)\,\mathrm dt -\beta _j^\lessgtr \right) \text { .}\nonumber \\ \end{aligned}$$

(24)

In the above, the functions $B_k$ are defined in Eq. (17), and the constants $\beta ^\lessgtr _j$ are determined by boundary and matching conditions. The constants $h_jk$ that appear in the integrand are determined by the $2n^2$ constraint equations

$$\begin{aligned}&\sum _{j=1}^{2n} \tilde{z}_{ij} h_{jk} \tilde{\lambda }_j = \delta _{ik}, \end{aligned}$$

(25)

$$\begin{aligned}&\sum _{j=1}^{2n} \tilde{z}_{ij} h_{ik} = 0 \text { .} \end{aligned}$$

(26)

Each of the above equations are $n\times n$ matrix equations, giving $2n^2$ total constraints.

To account for the corrections to the mass matrix we relabel the solution we found $\varepsilon ^0 _i$, substitute into the differential equations $\varepsilon _i = \varepsilon ^0_i + \delta \varepsilon _i + \ldots$ and restore $\eta _{ij}(\rho )$ in the differential equations. Keeping only terms to first order we write

$$\begin{aligned}&\frac{\partial ^2\varepsilon _i}{\partial \rho ^2}+\frac{2}{\rho }\frac{\partial \varepsilon _i}{\partial \rho }-\left. \frac{\partial ^2V (\phi )}{\partial \phi _i\partial \phi _j} \right| _{A} \varepsilon _j = B_i(\rho ), \nonumber \\&\frac{\partial ^2\delta \varepsilon _i+\varepsilon ^0 _i}{\partial \rho ^2}+\frac{2}{\rho }\frac{\partial \delta \varepsilon _i+\varepsilon ^0_i}{\partial \rho }-\left( \delta \varepsilon _j+\varepsilon ^0 _j\right) \left( \bar{M}_{ij}+\eta _{ij}\right) = B_i (\rho ).\nonumber \\ \end{aligned}$$

(27)

Since $\varepsilon _i ^0$ solve the initial differential equations in terms of $\bar{M}_{ij}$ by definition we can make an immediate simplification. Keeping only terms up to first order in our expansion we then write

$$\begin{aligned} \frac{\partial ^2\delta \varepsilon _i}{\partial \rho ^2}+\frac{2}{\rho }\frac{\partial \delta \varepsilon _i}{\partial \rho }-\delta \varepsilon _j \bar{M}_{ij}-\eta _{ij} \varepsilon ^0_j= & {} 0,\nonumber \\ \frac{\partial ^2\delta \varepsilon _i}{\partial \rho ^2}+\frac{2}{\rho }\frac{\partial \delta \varepsilon _i}{\partial \rho }-\delta \varepsilon _j \bar{M}_{ij}= & {} \eta _{ij} \varepsilon ^0_j. \end{aligned}$$

(28)

This once again is a set of coupled linear differential equations which has the same form as the original set of differential equations. So the solution is the same as before with $\eta _{ij} \varepsilon _j ^0$ replacing $B_i$. One then finds the series of $\delta \varepsilon _k$ up to a desired tolerance to find each value of $\varepsilon _i$. This process continues until the equations of motion are satisfied up to the desired tolerance.

4.2 Observations on convergence

Newton’s method is well known to have four major issues. In our analogous form, these issues would be:

• 1. If the initial ansatz function is too far from the true function, convergence will be slow.

• 2. Oscillating solutions where $\varepsilon ^{(n)}(\rho ) \approx - \varepsilon ^{(n+1)}(\rho )$.

• 3. Divergent corrections that arise in Newton’s method if the function’s derivative becomes undefined or zero. The equivalent issue will be discussed in detail below.

• 4. Being in the wrong basin of attraction and converging to the wrong function.

We will demonstrate that our method as applied here does not suffer from these problems with the exception of issue 4, where in principle a local minimum could be closer to the initial ansatz than the closest bounce-like extrema. This, however, is a limitation of all other known algorithms for finding bubble wall profiles. Meanwhile, we have already demonstrated in Sect. 4.1, our that our algorithm is free from the first problem as the guess of the initial ansatz ensures that the error functions are generically bounded by 1 (but should be much less than 1). For the remaining two issues, a little more care is needed.

Let us examine the issue of oscillating solutions. Let the updated function be

$$\begin{aligned} A_i^{(n)} = A_i^{(0)} + \sum _{k=1}^{n-1}\varepsilon _i^{(k)}, \end{aligned}$$

(29)

where $A_i^{(0)}$ is the initial ansatz and $\varepsilon _i^{(k)}$ are the correction functions. Suppose that, for field $\phi _i$, the successive correction functions begin oscillating at iteration n, so that $\varepsilon _i^{(n)} = -\varepsilon _i^{(n+1)}$. But this means that the equations of motion for $A_i^{(n+2)} = A_i^{(n)} + \varepsilon _i^{(n)} + \varepsilon _i^{(n+1)}$ can be written before Taylor expanding as

$$\begin{aligned}&\frac{\partial ^2 [A_i^{(n)} + \varepsilon _i^{(n)} +\varepsilon _i^{(n+1)}]}{\partial \rho ^2} + \frac{2}{\rho }\frac{\partial [A_i^{(n)} + \varepsilon _i^{(n)} + \varepsilon _i^{(n+1)}]}{\partial \rho }\nonumber \\&\qquad + \left. \frac{\partial V(\phi )}{\partial \phi _i }\right| _ {\left\{ A_k^{(n)} + \varepsilon _k^{(n)} + \varepsilon _k^{(n+1)}\right\} }\nonumber \\&\quad = \frac{\partial ^2 A_i^{(n)}}{\partial \rho ^2 } + \frac{2}{\rho }\frac{\partial A_i^{(n)}}{\partial \rho } + \left. \frac{\partial V}{\partial \phi _i} \right| _{\left\{ A_1^{(n+2)},\ldots ,A_i^{(n)}, A_{i_1}^{(n+2)},\ldots \right\} } = 0.\nonumber \\ \end{aligned}$$

(30)

Thus, in the case of a single field, an oscillating solution means that corrected field at the order where oscillation begins has solved the equations of motion exactly. In the multi-field case, as the fields $\phi _{j\ne i}$ converge without oscillating corrections, the changes to the derivative of the potential energy will diminish and thus will resemble the single-field case. In the case that more than one field has begun to receive oscillating corrections, this could prevent a rapid convergence but does not necessarily preclude it as the equations for the fields are still coupled.

In Newton’s method of finding roots, a major issue is when the derivative of the function becomes zero or undefined. The closest analogy to our method is the case where the mass matrix $\left. \frac{\partial ^2 V(\phi )}{\partial \phi _i \phi _j} \right| _{A^{(n)}}$ becomes zero or singular. In fact this is not an issue, as we can demonstrate. In the case that the mass matrix is zero, the differential equations become

$$\begin{aligned} \frac{\partial ^2 \varepsilon _i}{\partial \rho ^2} + \frac{2}{\rho } \frac{\partial \varepsilon _i}{\partial \rho } = B_i(\rho ). \end{aligned}$$

(31)

This is easily solved as

$$\begin{aligned} \varepsilon _i = \beta _0 + \frac{\beta _{-1}}{\rho } +\int ^\rho _0 \frac{\mathrm{{d}}y}{y^2} \int ^y_0 x^2 B(x)\,\mathrm dx \end{aligned}$$

(32)

where $\beta _0$ and $\beta _{-1}$ are both zero if the mass is zero everywhere. The case to consider is when the mass matrix is singular. This in fact does arise quite typically at some spatial points, but this is not a problem because the matrix inverse is not needed, and zero eigenvalues can be treated easily by using singular value decomposition or strategic placement of the step functions. Also, this issue is avoided if the full inhomogeneous differential equation is directly solved numerically.

5 Comparison with a solved example

We apply our method with the sample potential given in [43]

$$\begin{aligned} V(x,y) = (x^2 + y^2) \left[ 1.8(x-1)^2 + 0.2(y-1)^2 -\delta \right] .\nonumber \\ \end{aligned}$$

(33)

For $\delta = 0.4$ the potential deforms quite dramatically from the initial Ansatz so the convergence will be slower than for a typical case. We make a rotation in field basis $(x,y) \mapsto (u,v)$ such that u traces a straight line path from the origin to the global minimum and v is of course orthogonal to u. Our one dimensional potential is then given writing the potential in the rotated basis and setting v to zero. We then rescale such that the minimum is at $u=1$ and then we divide by |E| to get

$$\begin{aligned} \frac{V(u,0)}{|E|} = 0.36 u^2 - u^3 + 0.57u^4. \end{aligned}$$

(34)

We then use our analytic formulas to write the ansatz and make the appropriate rescalings to $u(\rho )$ and $\rho$ such that the ansatz is the solution to the original 1D potential. In the (x, y) basis the ansatz is

$$\begin{aligned} x(\rho )= & {} 1.046 \left( 1-\tanh \left[ \frac{\rho -0.437}{1}\right] \right) ,\end{aligned}$$

(35)

$$\begin{aligned} y(\rho )= & {} 1.663 \left( 1-\tanh \left[ \frac{\rho -0.437}{1}\right] \right) . \end{aligned}$$

(36)

Note that the wall width is only equal to 1 due to the rescaling. We have to sanitize our initial ansatz to set the derivative to zero as $\rho \mapsto 0$ or the correction diverges due to the $\phi ^\prime /t$ term in the differential equations. To achieve this we subtract from our initial ansatz

$$\begin{aligned} \delta x(\rho )= & {} \left. \frac{\partial x}{\partial \rho } \right| _{\rho = 0} \exp \left[ - \left. \frac{\partial x}{\partial \rho } \right| _{\rho = 0} \rho \right] , \end{aligned}$$

(37)

$$\begin{aligned} \delta y(\rho )= & {} \left. \frac{\partial y}{\partial \rho } \right| _{\rho = 0} \exp \left[ - \left. \frac{\partial y}{\partial \rho } \right| _{\rho = 0} \rho \right] . \end{aligned}$$

(38)

If one uses a small amount of step functions to approximate the spacetime dependent mass matrix one can find that the corrections $\delta \varepsilon _i$ are slowly converging. In particular it is useful to have step functions for regions where $m_{12}(\rho ) = 0$ and $m^2 _i <0$ as the functional form of the solutions changes in these regions. In the former the differential equations decouple for a region, for the latter, some of the exponents $\alpha _i$ are imaginary (but the $\varepsilon _i (\rho )$ remains real).

In Fig. 3 we show each iteration of the trajectory in the $(x(\rho ), y(\rho ))$ field space, along with the numerical trajectory as derived in [43]. The algorithm essentially converges after 3 perturbations. In Fig. 4 we show the x and y fields starting with the base ansatz forms $x^{(0)}$ and $y^{0}$, and then including the first three perturbative corrections, denoted by

$$\begin{aligned} \phi ^{(n)}(\rho ) = \phi ^{(n-1)}(\rho ) + \varepsilon _\phi ^{(n)}(\rho ),\quad \text {with}\; \phi = x, y. \end{aligned}$$

(39)

Figure 4 also includes the error functions $\varepsilon _\phi ^{(n)}(\rho )$ to illustrate the overall and diminishing magnitude of corrections to the fields in successive perturbations.

Fig. 3

The tunneling trajectory in the space of the x and y fields at the level of the base ansatz (solid straight line), and including the first three iterative corrections (solid curved lines). Also included is the numerical result, independent of our semi-analytic method, from Ref. [43] (dashed curve)

Fig. 4

The base ansatz form, first three iterations of corrections and the full numerical result of the x and y fields are presented in the top left and top right panels, respectively. The first three iterative corrections to the ansatz form of the x and y fields are given in the bottom left and bottom right panels, respectively

In Fig. 5 we show the error function to the ansatz for the x field, $B_x(\rho )$, which arises from the inhomogeneous part of Eq. (16). The error function is given for the bare ansatz solution of $x(\rho )$ and the first three perturbative corrections. We point out that the magnitude of the error is reduced by roughly a factor of 5–10 from each perturbative correction, and that the error function for $x^{(3)}(\rho )$ has reduced in magnitude by a factor of 300 compared to that of $x^{(0)}(\rho )$.

Fig. 5

The error function, B, of the ansatz solution $x^0$ when applied to the field equations, and the same for the ansatz solution including the first three perturbative corrections denoted by $x^{(n)} = x^{(n-1)} + \varepsilon ^{(n)}$

6 Conclusion

In this work we presented a new method to calculate the bubble profile in a bounce solution for a multi-field potential with a false vacuum. The method uses fitted functions to estimate the parameters of the single-field kink solution which is used as an ansatz form. It then applies this form to the full multi-field potential, which receive perturbative correction functions that are reduced to elementary numerical integrals. We have argued that the perturbative series of corrections should converge quadratically, and is immune to the issues of the analogous Newton’s method. This method is shown to be effective in solving a toy model with two scalar fields.