Multifield Polygonal Bounces

We propose a new approach for computing tunneling rates in quantum or thermal field theory with multiple scalar fields. It is based on exact analytical solutions of piecewise linear potentials with many segments that describes any given potential to arbitrary precision. The method is first developed for the single field case in 3 and 4 space-time dimensions and demonstrated on examples of classical potentials as well as the calculation of quantum fluctuations. A systematic expansion of the potential beyond the linear order is considered, taking into account higher order corrections, which paves the way for multiple scalar fields. We thereby provide a fast semi-analytical tool for evaluating the bounce action for theories with an extended scalar sector.


I. INTRODUCTION
Stability of the vacuum and phase transitions in the early universe are subjects of deep interest to particle physics and cosmology. The non-perturbative problem of the tunneling among two vacua was developed in seminal works [1][2][3][4] for single scalar field theories. The problem of evaluating the lifetime of such metastable states was solved by computing the action of a semiclassical instanton solution, called the bounce, interpolating between the two minima. The form of the bounce was proven to have O(D) invariance under general conditions [5] for D > 2 in flat spacetime.
While finding the bounce in four dimensions is needed to assess the stability of the vacuum, an analogous calculation becomes important at finite temperature. The bounce action in D = 3 dimensions sets the probability of bubble nucleation [6] and controls the quality of the contingent phase transitions. Moreover, the shape of the field solution, e.g. the size and thickness of the bubble is directly related to the power spectrum of gravitational waves [7,8] (see [9] for a recent analysis).
Computing the bounce action involves solving a nonlinear second order differential equation with a friction term dependent on D. Finding an analytical solution in a closed form is in general impossible for an arbitrary potential. However, an approximation can be found in the thin-wall regime [2] and examples of exactly soluble potentials include the binomial, logarithmic [10] and quartic one [11]. In most occasions the calculation of the bounce is thus performed numerically. For renormalizable single field potentials, one can use rescaling to define a single parametric problem and solve it by the usual shooting method [12,13]. Moreover, it is possible to derive an absolute lower bound on the bounce action [14][15][16][17] and to * victor.guada@ijs.si † alessio.maiezza@irb.hr ‡ miha.nemevsek@ijs.si provide estimates based on a tunneling potential [18] as well as machine learning techniques [19]. A remarkably simple example of a soluble bounce is the linear potential. This is the basis for our discussion that builds on the work of Duncan and Jensen (DJ) [20] in which two linear segments are combined into a triangular potential barrier. The shooting is transformed into an algebraic problem that is solved analytically in D = 4. This approximation was studied in [21] for single field and [22] for multi-field potentials. Another combination of two segments, one with a linear and other with a quartic potential was considered in [23]. The analytical continuation of the triangular solution from Euclidean to Minkowski space was developed in [24].
Finding the Euclidean action becomes harder when an arbitrary number of fields is considered. As shown recently [25], the bounce still keeps the O(D) invariance. Nevertheless, finding the path in field space and computing the bounce with multi-field potentials is significantly more challenging. The main difficulty with the usual shooting approach is finding the fine-tuned initial field value in the multidimensional field space, especially close to the thin wall limit, and integrating the system of coupled differential field equations.
In this work we propose a new approach to obtain the bounce solution that is based on a generalization of the DJ calculation, namely gluing an arbitrary number of linear segments into a polygonal potential and solving the resulting system for any D and any number of fields. By increasing the number of segments one can approximate any potential that admits a bounce solution, with an arbitrary precision and thereby obtain the relevant action.
The polygonal method enables one to work out bounces within non-analytic potentials, even when the usual approaches may have issues with stability.
In §II we review the basics of vacuum tunneling, introduce the polygonal method and construct the single field bounce solution. We discuss how this approach is employed in §III, where the relative convergence is evaluated on selected problems and contact is made with the existing tools. In §IV we show how these bounce solutions are used in the calculation of the decay rate pre-factor from one loop quantum fluctuations. In §V we extend the method beyond the linear approximation and pave the way for §VI, where the multi-field case is developed. We conclude with an outlook in §VII and leave details to appendices: dimensions other than D = 3, 4 are covered in A, the two segment calculation is expanded in B and further details on root finding can be found in C.

A. Bounce redux
Let us recall the basic features of vacuum transitions in field theory. We consider a single real scalar field ϕ in D dimensions, subject to an arbitrary potential V (ϕ) with non-degenerate minima, shown on the left of FIG. 1.
The probability of tunneling from one ground state to another is proportional to the Euclidean action S D . We assume the D dimensional solution to be O(D) symmetric [5] for any number of fields [25] where ρ 2 = t 2 + x 2 i is the Euclidean radius that sets the size of the bubble.
The bounce is an instanton solution of the Euler-Lagrange equation that interpolates between the minima of V and therefore obeys the appropriate boundary con- where d i V is the derivative of V with respect to ϕ i . The field starts at ϕ i0 with zero velocity and rolls down to a stop in the false vacuumφ iN at ρ = ∞. The usual shooting procedure involves numerically integrating the bounce Eq. (2) and varying ϕ i0 until the boundary conditions are met. In this procedure, care should be taken when numerically evaluating ρ → 0 or ∞. Conversely, notions of zero and infinity are not relevant for polygonal bounces below.

B. Polygonal bounces
In this work we introduce the polygonal bounce (PB) by generalizing the approach of [20]. Instead of the generic potential with two minima, let V (ϕ) be approximated by a polygonal piecewise linear approximation, as shown in FIG. 1.
Let us first develop the idea for the single field case, dropping the field index i and introducing the segment index for the field valuesφ s , s = 1, . . . , N , such that the two minima reside atφ 1,N . The values of the potential areṼ s = V (φ s ) and the linear segments are For linear V s , the exact solution of (2) on the section s is with D > 2. Two dimensions require minor modifications derived in A. Because we are dealing with a finite number of segments, the solution either a) starts from ϕ 0 at ρ = 0 withφ 1 = 0, which gives b) or waits atφ 1 until ρ = R 0 , which translates into Regardless of the initial condition, the field in the final section ϕ N −1 stops in the second minimumφ N at some final radius R N −1 such that where a 0 = a N = 0, because the first derivatives are zero in the minima. Thus there is no issue with the ρ → 0 limit: in case a) the singularity of the friction term is regulated by b 1 = 0, while in the case b) there is no singularity to start with and R 0 is non-zero. Similarly, the role of ρ → ∞ is taken over by the final radius R N −1 that is finite and numerically under control, see C for details.
The total Euclidean action of the bounce is then a sum of linear parts with the integrated kinetic and potential pieces which is valid for both instances, a) and b), with the understanding that R −1 = 0 and in case a) R 0 = 0. To determine the action above the field segments ϕ s (ρ) of the bounce need to be computed. To this end, a segmentation of {φ s } is set up, such that given the V (φ s ), the a s parameters are fixed by (3). We shall return to the choice of segmentation procedure in §III A below. What remains to be calculated are the v s , b s and the unknown radii R s , s = 0, . . . , N − 1.
We now demonstrate that solving the PB is a single variable problem, i.e. once the initial radius is known, the entire solution is determined. The free parameters are fixed by matching conditions required to glue neighbouring linear bounces into a single smooth solution, as in FIG. 1. There are three conditions, two for the field value ϕ s (R s ) =φ s+1 = ϕ s+1 (R s ) to match onto the initial segmentation at R s and another one for the derivativė These three conditions per segment precisely determine the unknown v s , b s and R s . Therefore, one can increase the number of sections at will without introducing additional free parameters.
The two-segment N = 3 problem can be solved analytically in some instances, as shown in B. With more segmentation points, one can transform the system into a single variable problem that can be solved numerically. For some particular D, further simplifications are possible.
Let us derive the recursion relations for R s (v s , b s ) such that they can be computed numerically. We first derive v s and b s by subtracting (12) from (13) and using (14) v The individual radii can be solved directly from (12) with δ s =φ s+1 − v s . Resulting Eq. (17) is a fewnomial with simple closed form solutions The radii corresponding to D = 2, 6, 8 can be found in Eqs (A7)-(A9) of A. This concludes the analytical setup of the PB construction. a. Derrick's theorem for piecewise actions. A well known result due to Derrick [47] is the relation between the integrated kinetic and potential parts in (10) and (11).
We will use this theorem to find the PB solution and to test the goodness of the approximation, so let us recall its essential point. For the action to remain minimal upon rescaling the argument of the solution to ϕ(ρ/λ), the following identity has to hold S (λ) For piecewise actions, such as the PB under consideration, the above identity is modified because (9) becomes a sum of finite integration intervals. While rescaling ρ → ρ/λ has no effect on integration limits in the continuous limit, rescaling the finite intervals R s → R s /λ in (9) introduces a manifest λ dependence. As a result, and similarly for V is quickly recovered. At the same time, one can use the continuous version of (21) with the input potential (not the polygonal approximation) to verify the goodness of the polygonal solution. This is shown on the right side of FIG. 11 in C, where about a permille level is achieved with N = 400 segments.

III. EVALUATING POLYGONAL BOUNCES
A. Implementation a. Overview. Let us turn to the implementation of the PB method. In the work of [20], the bounce equations were cast into an algebraic system and solved in a closed form. The approach followed here instead is to recursively compute the bounce parameters and solve a single boundary condition equation.
The boundary equation is obtained by combining (7) with (15) and setting s = N − 1, which leads to valid for all D. Because the R s are already solved for, the final condition for v N −1 holds automatically. Alternatively, one can use the relation in (21) with the polygonal potential, and look for the solution of In order to solve the boundary equation, either (25) or (26), one has to find the initial radius R s from which the subsequent v s , b s , R s are computed recursively until the boundary condition is satisfied. This is the algebraic analog of the shooting method used to solve (2) directly. Adding more segmentation points improves the accuracy of the approximation, but does not exponentially increase the computational burden, timing scales linearly with N .
b. Segmentation. To set up the polygonal potential approximation, one chooses a set of field values {φ s } that interpolate between the positions between which the tunneling happens, as exemplified in FIG. 1. Throughout this work we assume the original potential V (ϕ) to be non-pathological in the sense that it admits at least one bounce solution between these two values 1 .
To describe an arbitrary potential, enough segments should be taken to capture all the non-linearities with desired precision. In addition, the action converges faster if the segmentation is tailored to a specific potential, i.e. if the density of points increases close to the extrema. This geometrical insight is a particular feature of the polygonal approach and allows for intuitive understanding of the problem prior to the actual calculation of the bounce.
For a sufficiently large N the specific choice of coverage is not relevant, the naïve uniform distribution reproduces any reasonable potential when N → ∞ and converges smoothly to the final value. In this limit, the resolution of ∆φ s is small enough such that ϕ 0 always falls aboveφ 1 and only case a) persists. This is to be expected because such limit is equivalent to the original problem in (2) where R 0 → 0 and only ϕ 0 matters.
c. Computing the initial bounce radius. With a given segmentation at hand one has to find the initial radius R in that solves the boundary equation. Actually, the task can be simplified by a priori isolating the field segment on which the solution exists.
One can see from the right panel of FIG. 1 that the list of Euclidean radii {R s }, must be real, positive and growing (the true minimum is on the left by convention). On the other hand, Eq. (25) contains a number of nested roots and becomes progressively non-linear as N grows and generically admits complex solutions for the radii.
Let us demonstrate that the final radius R N −1 becomes imaginary as R in is varied across the true solution. This can be understood by noticing that the discriminant (20)  Furthermore, note that in both cases a) and b) one only needs to solve for R in , from which the initial field value ϕ 0 can be determined. In case b) this is merely the position of the minimum ϕ 0 =φ 1 , while in case a), it is obtained from R in and (12) From here one can infer the interval for R in ∈ [0, R max in ] by setting ϕ 0 to the lower and upper boundary of the segment in (27). The way to find the segment with the solution a priori is therefore to evaluate the final radius from these two limiting R in and checking whether it becomes imaginary, as illustrated in FIG. 2.
Once the segment containing the solution has been found, one can proceed to solve the polygonal bounce by does not exist, since the field cannot wait at the true minimum. Instead, the choice of the exit point, i.e.φ 1 must be deep enough for the field, starting from ϕ 0 , to roll down to the false minimum.
Schematic overview of finding the PB. The segment with the solution (in this example s = 2 and Rin = R2) can be found by evaluating the PB on the boundaries of R min 2 = 0 and R max 2 and checking that the imaginary part of the final radius RN−1 becomes non-zero. Finally, the solution of R2 is found such that the scaling parameter λ → 1.
solving either (25) or (26). Another approach is to take advantage of the fact that the bounce solution depends solely on R in . This is a dimensional parameter, which can therefore be rescaled by the optimal amount computed from (26), which essentially aims to minimize the action. For example, one may begin with R max in , compute the corresponding λ, which in general will be different from 1, and proceed by iteration from R in = λR max in . This procedure converges in a few iterations to a permille level. Alternatively, one can solve (26) with standard root finding algorithms.
By increasing the number of segments, the initial radius (e.g. R 0 in case b)) decreases until R in = 0, when the domain of the solution disappears and one has to switch to the next segment. This agrees with (2), as does the fact that the final radius R N −1 grows steadily to infinity when N → ∞, see FIG. 12.

B. Examples, convergence and comparisons
a. Linearly displaced quadratic potential is the benchmark potential to test the PB method. It is defined as in the work of Coleman [2] and shown on the left panel of FIG. 1. For convenient numerical evaluation, we set λ = 0.25, v = 1; other points in parameter space can be obtained by rescaling [13]. For such choice of parameters, varying ε from 0.01 to 0.08 covers all the regions of interest, starting from thin wall regime of small ε, going to well separated minima until the second minimum disappears. We now apply the PB method to the potential in (28), employing the homogeneous segmentation for simplicity. The first results are the ϕ 0 and R 0 that attempt to solve (25). The solution for R 0 varies with N , therefore we show the behavior of R  The smaller ε is, the closer one goes towards the thin wall regime, where the field needs to wait close to the minimum. This means R 0 remains sizeable for higher values of N and one needs to introduce many segments for R 0 to reach zero, as clear from FIG. 3. On the other hand, the transition from b) to a) happens faster when ε increases. Finally, when ε is large enough, the transition eventually disappears and we are left with case a) right from the start at N = 3.
The number of dimensions also has an impact on the transition from b) to a), as seen in FIG. 3. Keeping ε fixed, the transition in D = 4 occurs for higher N with respect to D = 3. This is expected because the damping term in (2) is proportional to D and thus becomes more important in higher dimensions.
The final step after obtaining R 0 or ϕ 0 is to compute the main object of interest: the Euclidean action S D in (9)   In the limit of ε 0 one ends up in the thin wall regime, and therefore N = 3 has to produce the correct result of [20], in agreement with the inset of FIG. 4. With increasing ε, the potential in (28) will eventually lose the second minimum. For any potential close to this threshold, the resolution of the homogeneous segmentation has to be precise enough to describe the local maximum, otherwise the solution cannot exist a priori. This is precisely what happens in FIG. 4 for ε = 0.08, the N = 4 segmentation is too rough to possess an intermediate maximum. In general, the approximation worsens for 4 ≤ N < O(10), which is an artefact of the assumed uniform segmentation. Conversely, for higher N , the action starts to converge rapidly and the rate is faster in case b) for smaller , where the shooting method instead becomes increasingly unstable.
The initial approximations with small N s, shown in the inset, are already quite close to the end result, and are valid at about 10% level. It is clear that the N = 3 segmentation always underestimates the action and this simple approximation becomes progressively better as ε decreases. On the other hand, as N increases, the method starts to overestimate the bounce and converges to the final result from above. Even for moderate N = 10 the accuracy of the estimation is below 10% and goes below the permille level when N = 200. The convergence is slightly faster for N = 3, moreover the rate of convergence can be improved by choosing an appropriate segmentation.
To compare the PB method to existing methods, we show the results of other approaches in FIG. 4. The other three calculations are the usual shooting method of Eq. (2) and the out-of-the-box results from Cosmo-Transitions [29] and AnyBubble [33] packages. Note that in these examples all the methods agree within a few permille level.
b. Bi-quartic potential is another example of a simple but non-trivial exact solution [23] that glues two quar- tic field functions. Since the bi-quartic bounce is computed analytically, it can serve as an additional test of the polygonal method.
The bi-quartic potential is parameterized by ε ϕ 4 representing the gap in the potential difference, which varies from the thin to thick wall regime, as shown in FIG. 5. The presence of the cusp creates issues for standard approaches based on the shooting method, due to the nondifferentiable potential. On the other hand, the polygonal method turns out to be quite robust and the solution can always be found. Nevertheless, for smooth convergence it is convenient to employ a bi-uniform segmentation on both sides of the cusp.
The resulting bounce action is shown on the right of FIG. 5 and goes below percent level accuracy with O(100) field segments. The dashed lines also provide the comparison to CosmoTransitions 3 .
c. Other potentials and additional minima. We also tested the PB approach on the potential in (28), corrected with the logarithmic field dependence. With such a deformed potential, the calculation proceeds exactly as before and the polygonal approximation works as expected.
The method was successfully applied on examples with further intermediate local minima. Such a situation may arise when more fields are involved and a particular path in field space is chosen. The prototype of N = 5 with two triangles and a single additional minimum can be considered as the minimal setup illustrating such situations. As long as the bounce exists, i.e. if the intermediate minimum is not too deep, (25) gives a consistent real solution.

IV. QUANTUM FLUCTUATIONS WITH POLYGONAL BOUNCES
The simplicity of the semi-classical polygonal solution can be exploited also for computing the quantum corrections, i.e. the prefactor of the decay rate, originally 3 We were unable to recover the value of the action by using Cos-moTransitions out-of-the-box. Instead, we extrapolated the field bounce solution to manually compute the action shown on FIG. 5. Still, this procedure failed for lower values of ε ϕ 4 closer to the thin wall regime. Moreover, computing the bounce using AnyBubble was not possible for any value of ε ϕ 4 . derived in [3]. A number of studies on the prefactor proposed different numerical methods in D = 3 [36,37] and D = 4 [38,39], and recent progress has been made on precision calculations in presence of gauge interactions [42], scale invariant instantons and extended gauge theories [43]. On the other hand, not many explicit analytical results on the prefactor are available with a notable exception of the thin wall limit [40].
The total decay rate at one loop is where S 4 is the semi-classical action computed from the bounce solution ϕ(ρ) and det is the determinant of the fluctuation operator, i.e. the product of its eigenvalues, with zeros removed. Finally, δ 4 is the perturbative one loop counterterm of the action that absorbs the renormalization infinities.
In computing the determinant, we follow the work of Dunne [39], where the fluctuation operator O, i.e. the second variation of the action, is decomposed in a multipole expansion due to the O(4) symmetry where V (ρ) is the rewritten form of (28) with the removed asymptotic value of 1, such that for the fluctuations around the true vacuum V (φ 1 ) = 0.
Instead of computing all the eigenfunctions ψ l (ρ) of O l and summing the corresponding eigenvalues, it is convenient to use the Gel'fand-Yaglom theorem [45] that relates the ratio of determinants to the value of the ratio of eigenfunctions evaluated at infinite Euclidean time The calculation of the pre-factor splits in two parts: the low l region up to an arbitrary l ≤ L O(10) and the high l region, going to infinity. In the low l part the ratio of determinants R l is computed by solving the partial differential equation for R l , because the solutions ψ free l for the fluctuations around the true vacuum are known Bessel functions ψ free l = I l+1 (ρ)/ρ. The bounce solution ϕ(ρ) determines the shape of V (ρ), and in the low l regime, the contribution to the decay rate is finite and proportional to the sum of the log of all the ratios of determinants: On the other hand when l > L 1 one can solve for the R l using the WKB approximation [39], which is regularized with the proper counter terms in δ 4 (with optional higher orders [46] for faster convergence). This high-l part of the rate, i.e. − ln Γ hi is with the three relevant integrals given by These integrals are straightforward to compute analytically and the total prefactor contribution is the sum of the low and high l pieces. ρ The crucial component in computing R l is of course the semi-classical bounce solution. In FIG. 6 we show the resulting R l using the precise numerical shooting solutions and the PB approximation with the minimal N = 3 and the more precise N = 50.
The values at infinity R l (∞) agree with the expectation of a single negative eigenvalue for l = 0, four-fold degenerate zero for l = 1 and the rest of l ≥ 2 being positive. This is true for the precise shooting procedure, however the N = 3 PB bounce produces a number of negative eigenvalues, while for N = 50 the correct spectrum is recovered. This happens because the semi-classical solution is not approximating the exact potential with sufficient precision, the proof for one negative and multiple zero eigenvalues [44] (and the entire calculation of the fluctuations) relies on the fact that the semi-classical action is extremised. Nevertheless, blithely summing the absolute values of R l (∞) gives a rather precise (and very simple) estimate of the decay rate prefactor, as seen in TAB. IV.
The crude N = 3 approximation fails when α 1, however it works well in the thin wall limit when α → 1 and all of the approaches coincide, as shown on the right panel of FIG. 6. In this section we develop a general procedure of including non-linear corrections to the PB. This is done by setting up a systematic procedure based on the Taylor expansion of the potential and then building the new bounce solution perturbatively on the PB ansatz.
Higher order corrections describe non-linear features that are not there in the leading approximation, for example around the extrema of V where the linear part of the potential vanishes. Although the PB solution is formally exact when N → ∞, the nonlinear corrections may enhance the convergence of the action, depending on the type of the potential and the order to which we correct.
a. Generalities. Consider the complete bounce solution expanded around the PB: ϕ = ϕ P B + ξ, such that the correction to the potential is evaluated on the PB background and the bounce equation becomes where α is an arbitrary linear part. The bounce correction ξ is then given by Evaluating the above integral I for an arbitrary δdV and computing the unknown parameters of ξ is involved and basically equivalent to the numerical integration of (2). However, a systematic expansion of the potential and linearization simplify this approach considerably. b. Perturbation. On a given segment, the potential can be expanded in Taylor series aroundφ s where the constants ∂V s , ∂ 2 V s , . . . are determined by matching the values and (higher) derivatives of V . When N increases, the segmentation becomes arbitrarily dense and thus the terms beyond the linear one in (43) become progressively negligible.
To illustrate this point, we expand V to second order where dṼ s stands for the derivative of the original potential evaluated atφ s . This is the additional information required from the original potential in order to get to the next-to-leading order. The α s coefficients are thereby fixed and the inclusion of the quadratic correction improves the fit of the potential near the extrema, as seen from FIG. 7. Moreover, with a large N , one has α s a s as clear from (44), which is consistent with the assumption of perturbativity.
With this approximation of the potential, the nonhomogeneous part of the correction is which can be evaluated for D = 3, 4 where the arbitrary integration constants ρ 0 , ρ 1 were chosen to simplify the expression for I s without loss of generality because they can be absorbed in ν s , β s . The remaining task is to compute the unknown coefficients ν s , β s and the new matching radii by requiring the solution to be continuous and differentiable as in the PB case. Given that ϕ P B and its matching radii are already close to the actual solution, the new radii have to be close to the previous ones Following the same procedure as in the PB construction above, we set up the modified initial, final and matching conditions for the correction ξ. These conditions are then perturbatively linearized in r s to get the recursion relations for the parameters and similarly a linear equation for the radius correction at each segment is Following the same logic as in the PB case above, we compute the initial radius correction r in by solving the linear equation that satisfies the final matching condition. Being a linear equation, this additional step does not require significant computing time but improves the accuracy of the action and speeds up convergence.
c. Improved action. To understand the effect of second order corrections, we reconsider the usual displaced quartic potential and show the improved action in FIG. 8. The correction significantly improves the approximation of the action by nearly an order of magnitude improvement for any given N and ε. In other words, to achieve the same level of accuracy one needs to consider half as many segments.
Because the polygonal bounce perturbation requires only to solve a linear equation, the computational cost of computing the bounce solution with a given accuracy is reduced significantly. Moreover, the final result of the bounce field configuration is again given in the form of segmented analytical functions, which allows for further manipulation.

VI. MULTI-FIELD POLYGONAL BOUNCES
Computing the false vacuum decay rate with multiple scalar fields faces a number of technical difficulties. These are related to the fact that the Euclidean action is not a minimum but a saddle point. In terms of the bounce solution, one has to look for the (fine-tuned) initial condition in the higher dimensional field space and then integrate the coupled system of differential equations, usually numerically.
Existing approaches to this problem [26][27][28][29][30][31][32][33][34][35] address these challenges in various ways. In general, solutions where the shooting and path deformation are decoupled exhibit oscillatory (and therefore slower) path convergence, multifield shooting face non-linear scaling with the number fields, and most approaches have difficulties with thin wall regimes and provide purely numerical output of the bounce field configuration, as well as the Euclidean action.
The PB solution overcomes a number of these shortcomings and provides a framework with the following features.
a) The multifield PB field solution remains as simple as in the single field case in (4). It is therefore fast to evaluate numerically and is retained upon iteration. The final result has a closed analytical form, which allows for further manipulation.
b) The solution is built iteratively, where a single iteration takes into account the curvature in field space by explicitly solving the ρ dependence and simultaneously deforms the path. This eliminates the oscillatory behavior and the solution converges quickly, within O(1) iterations, see FIG. 10.
c) The method works very well in the thin wall limit, which is usually problematic due to severe finetuning. This feature is directly inherited from the single field case and is due to the fact that we are solving for the Euclidean time ρ variable and not in ϕ space. Of course, the method works equally well (see again FIG. 10) in the thick wall regime; moreover it is applicable to cuspy and unstable potentials, as well as paths with multiple minima. d) Finding the path in field space boils down to a coupled system of ordinary linear equations that scales linearly with the number of fields and number of segments. The procedure converges very close to the final path even with a few -O(1) segments.
One can switch to more segments in the final step only to ensure sufficient precision in the longitudinal direction, depending on the desired precision of the action.
e) It works for any space-time dimensions D > 2 (with D = 2 in the Appendix (A)), in particular it is simple to consider D = 3, 4, which are most relevant for physical applications.
A. Constructing multi-field polygonal bounces Let us describe the generalization of the PB approach to an arbitrary number of scalar fields. The starting point is a single field PB solutionφ is where i is the field index i = 1, . . . , n f and s = 1, . . . , N is the segment point. The ansatz is obtained from a selection of initial points in the multi-field spaceφ is , for instance by segmenting a straight line connecting the two minima, as in the left panel of FIG. 9, and computing the corresponding longitudinal PB, seen on the right panel of FIG. 9.
We then consider an expansion around the initial estimate, such that ϕ is (ρ) =φ is + ζ is . This produces a set of coupled bounce equations for each field direction The idea here is to look for a solution of the field expansion ζ which is of the polygonal type where a is corresponds to the leading constant expansion of the gradient of the potential around some deformed path, defined byφ is +ζ is . This is the main difference in contrast to the single field case: the position in field space is not fixed a priori and one has to allow for the segmentation to move in field space. The gradient parameters a is can be linearized in terms of the displacementζ js with a symmetric average It is crucial that the gradient in (56) is expanded beyond the constant leading order up to O(ζ) that includes the second derivative of the potential. This is needed to properly describe curved paths in field space.
Matching. To fix the remaining parameters of the ζ solution in (54), the field has to match onto the deformed path. We choose to match toζ at the fixed radii R s , computed from the initial longitudinal polygonal ansatz. This can be done for all the R s , except for the initial R i0 and final ones R iN −1 , which are free parameters for each field direction i.
The field values of the ansatzφ is are continuous from one section to another, while the derivatives may not be. The matching of derivatives at R s then gives the recursion relation for b is and field continuity, together with (57) provides the re- Initital/final conditions. In case a) the initial endpoint is free to move, however the solution starts at ρ = R i0 = 0 with a vanishing derivative, therefore In case b) the initial endpoint does not move and we Here we expanded the initial and final radii R i0 = R 0 (1 + r i0 ) and R iN −1 = R N −1 (1 + r iN −1 ) to leading order in r i0,N −1 , in order to maintain a linear system. As for the derivatives,φ whereφ i1 = 0 andφ i1 = 8ā i1 follows from (53). In summary we have the following conditions The final task is to solve this linear system. The initial conditions are solved in terms of which determines ζ i1 that has to be fixed toζ i2 at R 1 . The recursion relations (57) and (58) then provide the polygonal ansatz for ζ is , to be fixed ontoζ is+1 This continues until the final segment where the endpoint does not move anymoreζ iN = 0, in agreement with (61).
The final equation to be solved is then theζ iN −1 condition in (63). By construction, (54) keeps the same polygonal form in ρ, therefore it is simple to iterate and converges once the path in field space does not change anymore, i.e.ζ is 0.

B. Examples and path convergence
Let us consider a simple two field potential that has multiple solutions for spontaneous symmetry breaking vevs ϕ i = v i . The metastable minima are in general of different depths with V (v 1 ) = V (v 2 ), which allows for the local false vacuum to decay into the global minimum by traversing the field space along the bounce solution.
To illustrate the multi-field PB method, we choose two exemplary points in the parameter space to cover both non-trivial cases: a) and b). Specifically, we take µ Again, the convergence of the action can be improved by taking into account also the ρ dependence of the PB ansatz, similar to the single field extension defined above. It is also possible to solve the multifield bounce equation by solving for ζ dynamically and gluing the corresponding Bessel functions. This is a somewhat tedious task that requires local field rotations and is beyond the scope of the current work, but a similar semi-numerical approach was done in D = 3 by [32].
Finally, the path converges to the final one without oscillations, in contrast to [29] where the ρ dependence of transverse field directions was dropped, effectively neglecting the kinetic term. Since we use an explicit solution in (54), the dynamical term of the curved path is taken into account. This happens also in [32], where the field construction is slightly more involved, requiring local rotations and evaluation of Bessel functions.

VII. CONCLUSIONS AND OUTLOOK
An efficient and fast approach for calculating the false vacuum tunneling rate is developed for arbitrary potentials with any number of fields up to the desired precision. The method is based on the simple, well-known exact solution [20] that is extended to any number of segments, space-time dimensions and number of scalar fields.
Usually, the simple single field problem of finding the bounce is solved by shooting -numerically integrating the bounce equation and looking for the correct initial condition. Here instead, the differential equations are solved exactly and are glued into a single continuously differentiable field. The boundary conditions can be solved exactly and the field solution is computed recursively. The remaining initial/final conditions are highly non-linear but can be solved by iterative use of Derrick's theorem or numerical root finding.
In contrast to numerical integration, the PB solution is given by segmented polynomials. This allowed for simple analytical manipulation, such as including corrections of higher orders in the potential expansion, quantum or thermal fluctuations, expanding to more fields and it ultimately reduces the computational cost. Because the one field solution depends on a single dimensional parameter, which is the initial radius defined on some initial segment, the fine-tuning of initial conditions is avoided. This is advantageous especially in the thin wall regime, where the usual shooting procedure struggles.
The method was applied to a number of single field examples, from the simplest displaced quartic potential to more involved cases, such as the bi-quartic potential. The resulting bounce action converges quickly with N O(10) and reaches a permille level precision as seen in FIG. 4, where the comparison with existing tools is made. The semiclassical bounce solution was also employed in the calculation of the one-loop quantum corrections, i.e. the prefactor of the decay rate.
The simplest polygonal potential can serve as an ansatz to be perturbatively deformed in order to describe the remaining non-linearities. These are generically important close to the extrema and their inclusion improves the convergence of the bounce action, as seen from FIG. 8. The ability of perturbative expansion allows for the generalization to the multi-field case. The main challenge with respect to the single field case is finding the path in field space. The PB approach solves it by starting from an initial polygonal ansatz that is iteratively deformed by solving the bounce equations at the leading order. Path deformation is solved by a linear system and converges very quickly without oscillations such that the action is recovered to arbitrary precision within a few iterations.
In summary, we find that the PB method is a robust, precise and reliable way of computing the semi-classical tunneling rate for any given potential. This approach describes the false vacuum decay in flat space time, however the solution can also be used in curved space-time within a small gravitational field approximation [48,49]. The PB solution and its extension can thus provide a tool with an analytical insight in characterizing stable vacua of theories with multiples scalar fields [50][51][52][53][54][55], describing bubble nucleation and the quality of potential first order phase transitions as well as the related spectrum of gravitational waves.

ACKNOWLEDGMENTS
We would like to thank Borut Bajc for valuable discussions and suggestions. The work of VG was supported by the Slovenian Research Agency's young researcher program under the grant No. PR-07582. AM was partially supported by the H2020 CSA Twinning project No. 692194, "RBI-T-WINNING". MN was supported by the Slovenian Research Agency under the research core funding No. P1-0035 and in part by the research grant J1-8137. MN acknowledges the support of the COST actions CA15108 -"Connecting insights in fundamental physics" and CA16201 -"Unraveling new physics at the LHC through the precision frontier". and the recursion relations in (15) are modified by applying the replacement of (A2) to v s . The radii in two dimensions are solved by where W (z) is the product log function that returns the solution of w to the equation z = w e w for a given z.
The polygonal bounce setup for D = 6, 8 closely follows the procedure outlined in II above, apart from the solution of the radii fewnomial in (17). Indeed, the two closed form solutions for D = 6, 8 are Single-field. The simplicity of having only three points allows for some further progress. In particular, the shooting in ϕ 0 for case a) can be carried out analytically for any D.
a) Here, b 1 = 0 and R 1 is easy to solve from (12), while R 2 follows from (25) The recursion for v s in (15) and the final condition for v 2 (7) for N = 3 give For case a) to be consistent, the final solution should obey ϕ 0 >φ 1 .
b) Plugging the initial/final conditions of (6) and (7) into (25), (12) and (15) gives This system can be reduced to a single non-linear equation that can be solved numerically for any D. However, D = 4 is special, here a simple closed form solution can be obtained. The above equations can be rewritten as Expressing R 2 0 = R 1 (R 1 − ∆ 2 ) and R 2 2 = R 1 (R 1 + ∆ 3 ) and plugging R 0,2 into (B10) gives Multi-fields. The minimal multi-field case with N = 3 can be carried out analytically up to a single n 2 f linear system. The initial conditions in (61) and (62) with recursion relations (57) and (58) give  (25). Right: The continuous version of Derrick's theorem (21) with T computed with the PB and V from the input potential in (28). The normalized quantity acts as a test of convergence and goodness of approximation.
This leaves us with three equations for r i0 , r i2 andζ i2 Inserting r i0 from (B14) into (B15) gives a linear system forζ i2 that can be solved using the explicit form of a i1,2 (ζ i2 ) given in (55). Onceζ i2 is given, r i2 follows from (B16), which concludes the calculation of ζ.
Remarkably, this simple estimate already gives a rather good approximation for the path in field space, the main inaccuracy in the bounce action is due to the poor estimate of the ρ dependence.
Appendix C: Real radii and root finding a. Real radii. The radii solutions in Eqs. (A7)-(A9), as well as those in (18), (20) above, allow for a number of branches. The ones chosen above are such that the resulting R s are real and positive. Moreover, the slope of the potential a s has to be appropriately factorized in the expressions above in order to maintain the reality of R s during the transition through the maximum of V when a s flips the sign. This choice of signs also ensures that the radii of segments below the initial ϕ 0 automatically remain 0, i.e. R s = 0 forφ s < ϕ 0 .
b. Root finding. The starting point for root finding is to determine the real domains of the initial parameters ϕ 0 and R 0 for a) and b) cases, respectively. This defines the region of parameter space where a consistent solution can be searched for. To illustrate this point, we show the behavior of the final radius with respect to R 0 and ϕ 0 in FIG. 11. It is curious that the solution to the matching equation in (25) lies precisely on the edge of the real domain. In order to implement the root searching numerically, one has to define a starting estimate for R 0 or ϕ 0 . It turns out that for case a) the more stable option is to choose the initial estimate for ϕ 0 close to the false vacuum ϕ 0 φ 1 , while in the case b) the N = 3 result gives a fairly reliable starting point. Moreover, the behavior of case b) root finding convergence is in general more stable with respect to case a).
The behavior of ϕ 0 that solves the polygonal bounce in case a), is shown on the left of FIG. 12, where the field is normalized to the position of the false minimum inφ 1 . Notice that as ε decreases, the solution gets closer toφ 1 and eventually crosses over to case b). The smaller N approximation typically underestimates the final value and oscillates towards the limiting value, which is an artefact of the segmentation.
Note also that for ε = 0.05(0.04), the solution for case a) does not exist until N 10(70) when the segmentation becomes refined enough for the method to work and which is precisely when R 0 becomes non-zero in FIG. 3. Another particularity related to the segmentation happens with ε = 0.07 in D = 4 where we start in case a) for N = 3, 4, switch to case b) and return back to a) at N = 8.
The right panel of FIG. 12 shows the extent of the non-trivial part of the bounce field solution in the ρ dimension, i.e. the final radius R N −1 , normalized to the N = 3 approximation. Above this radius, the bounce solution remains constant as in FIG. 1. As we expect to get back to (2) in the continuous limit, the R N −1 should go to infinity when N increases, which is evident from the right panel of FIG. 12. As discussed above, the R N −1 is a finite and numerically well defined quantity that regulates the infinity of ρ. In particular, the extent to which the final radius grows is surprisingly small. Even for a large number of points N ∼ 400 where the bounce action is already quite precise, the final radius is merely about 50% larger than the initial estimate from N = 3.