Semi-analytic techniques for calculating bubble wall profiles

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.

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 a e-mail: sujeet.akula@coepp.org.au b e-mail: csaba.balazs@monash.edu c e-mail: graham.white@monash.edu 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.

Fate of the false vacuum
Consider a potential of multiple scalar fields V (φ 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 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 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 and the bounce solution satisfies the conditions 1 Here ρ is the ordinary 3D spherical coordinate, as we are considering finite temperature, and v true i and v false i are the vacuum expectation values of the field ϕ 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 ρ playing the role of time, but including a ρ-dependent friction term.
3 Approximate solution to the multi-field potential

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φ 1 whereφ 1 is a unit vector that points in the direction of the true vacuum. Then define a complete set of unit vectors orthogonal toφ 1 and rewrite the potential in the rotated basis V (ϕ 1 , ϕ 2 , . . .) → V (φ 1 , φ 2 , . . .). Then consider the potential only in theφ 1 direction between the minima, V (φ 1 , 0, . . .). One can then solve the single equation of motion 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, An approximate expression for the effective action of a similar potential was found in Ref. [50]. The above potential can be rescaled φ = ϕ min ϕ where ϕ min is the global minimum of the above potential. Then the rescaled potential has a global minimum at φ = 1. To ensure that the effective action is unaffected by this rescaling, we also make the replacement ρ → ϕ min ρ. We paramatrize the rescaled potential as 2 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 α ∈ (0.5, 0.75). 3 To illustrate this point, we present in Fig. 1 the potential in φ given in Eq. (6) for both the edge choices of α and the mean allowed choice, using several choices of E. One can see that, for α = 1 2 , we have exactly the Mexican hat potential (albeit shifted to φ = 0.5) with degenerate minima, and for α = 3 4 , there is no potential barrier between false and true minima.

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 ρ → ρ/ √ |E|. Then the equations of motion only depend on α.
Under the scaling we have introduced, the Euclidean action becomes where 4 V ≡ V /|E|, and we have integrated over the angular variables assuming isotropy. The integral in Eq. (8) must only depend on α, as in Meanwhile, we will approximate the rescaled field itself with the well-known "kink" solution [44,45] φ parametrized by the offset δ and the bubble wall width L w .
Thus it remains to determine the α-dependent functions δ(α), L w (α) from the kink solution, and f (α) in the Euclidean action.
We first evenly sample values of α within (0.5, 0.75), then numerically solve the full bubble profile using conventional techniques. Next, for each value of α, we fit the kink solution given in Eq. (10) to the full solution, extracting L w and δ. Lastly, we numerically integrate to find f in the Euclidean action. This results in a tabulation of values for L w , δ, and f , for each value of α. Using the apparent α dependence and intuition from our parametrization of the potential, we find ansatz functional forms in terms of α for each of these parameters.
The offset δ should diverge at the boundaries α = 0.5 and α = 0.75, and is found to be quite small otherwise. It also appears to have approximate odd parity about the mean allowed value of α = 0.625. We modeled this with odd powers of non-removable poles at the boundaries of α, along with an offset and a linear correction about the mean: We then fit these parameters using the tabulated values. The bubble wall width L w is dominated by two asymptotes. It diverges at α = 3 4 , and become small as α → 1 2 . We used this form to model the asymptotic behavior, As before, the normalization 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.) The Euclidean action determined by f (α) diverges at α = 1 2 and is zero at α = 3 4 . This is modeled by 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.
We note that |E| scales as |bφ 3 | so S E /T scales as 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.

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 con- Here, we iteratively determine functional corrections to the ansatz form.

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 φ = {φ i (ρ)} generically does not exceed the distance between the two minima. That is, 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 φ i , A i with correction ε i , so that φ i = A i + ε i . Applying this to Eq. (3) yields where A ≡ {A i }. We can then rearrange the above to separate terms that depend only on the ansatz forms from those involving the unknown correction functions, ε i . This leaves us with a set of coupled inhomogeneous differential equations for ε i : Here the functions B i (ρ) are the inhomogeneous part of the differential equations for ε i , and they are given by 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 ε 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 ε of the form 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 where the index i refers to field φ i , and the index k will run over the n eigenvalues of the mass matrix, and k is not summed over. Substituting ε h ik for ε i in Eq. (16) with B i = 0 yields j M i j z jk = λ 2 k z ik (20) where M is the mass matrix In this way, z ik is the ith element of the eigenvector of M that has eigenvalue λ 2 k . It is necessary, however, to use both positive and negative roots, ±λ k . Thus we must further introducẽ λ j andz i j with so that inλ j andz i j , the index j = 1, 2, . . . , 2n. Finally, we use the techniques described in [49] to arrive at the particular solution, In the above, the functions B k are defined in Eq. (17), and the constants β ≶ j are determined by boundary and matching conditions. The constants h j k that appear in the integrand are determined by the 2n 2 constraint equations 2n j=1z i j h ik = 0 .
Each of the above equations are n×n matrix equations, giving 2n 2 total constraints. To account for the corrections to the mass matrix we relabel the solution we found ε 0 i , substitute into the differential equations ε i = ε 0 i + δε i + . . . and restore η i j (ρ) in the differential equations. Keeping only terms to first order we write Since ε 0 i solve the initial differential equations in terms of M i j by definition we can make an immediate simplification. Keeping only terms up to first order in our expansion we then write 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 η i j ε 0 j replacing B i . One then finds the series of δε k up to a desired tolerance to find each value of ε i . This process continues until the equations of motion are satisfied up to the desired tolerance.

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 ε (n) (ρ) ≈ −ε (n+1) (ρ).

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 where A (0) i is the initial ansatz and ε (k) i are the correction functions. Suppose that, for field φ i , the successive correction functions begin oscillating at iteration n, so that ε . But this means that the equations of motion for A can be written before Taylor expanding as 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 φ j =i converge without oscillating corrections, the changes to the derivative of the potential energy will diminish and thus will resemble the singlefield 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 ∂ 2 V (φ) ∂φ i φ j 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 This is easily solved as where β 0 and β −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.

Comparison with a solved example
We apply our method with the sample potential given in [43] V (x, y) = (x 2 + y 2 ) 1.
For δ = 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) → (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 We then use our analytic formulas to write the ansatz and make the appropriate rescalings to u(ρ) and ρ such that the ansatz is the solution to the original 1D potential. In the (x, y) basis the ansatz is 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 ρ → 0 or the correction diverges due to the φ /t term in the differential equations. To achieve this we subtract from our initial ansatz If one uses a small amount of step functions to approximate the spacetime dependent mass matrix one can find that the corrections δε i are slowly converging. In particular it is useful to have step functions for regions where m 12 (ρ) = 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 α i are imaginary (but the ε i (ρ) remains real).
In Fig. 3 we show each iteration of the trajectory in the (x(ρ), y(ρ)) 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 φ (n) (ρ) = φ (n−1) (ρ) + ε (n) φ (ρ), with φ = x, y. Figure 4 also includes the error functions ε (n) φ (ρ) to illustrate the overall and diminishing magnitude of corrections to the fields in successive perturbations.
In Fig. 5 we show the error function to the ansatz for the x field, B x (ρ), which arises from the inhomogeneous part of Eq. (16). The error function is given for the bare ansatz solution of x(ρ) 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) (ρ) has reduced in magnitude by a factor of 300 compared to that of x (0) (ρ).

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.