BubbleDet: A Python package to compute functional determinants for bubble nucleation

We present a Python package, BubbleDet, for computing one-loop functional determinants around spherically symmetric background fields. This gives the next-to-leading order correction to both the vacuum decay rate, at zero temperature, and to the bubble nucleation rate in first-order phase transitions at finite temperature. For predictions of gravitational wave signals from cosmological phase transitions, this is expected to remove one of the leading sources of theoretical uncertainty. BubbleDet is applicable to arbitrary scalar potentials and in any dimension up to seven. It has methods for fluctuations of scalar fields, including Goldstone bosons, and for gauge fields, but is limited to cases where the determinant factorises into a product of separate determinants, one for each field degree of freedom. To our knowledge, BubbleDet is the first package dedicated to calculating functional determinants in spherically symmetric background


Introduction
Functional determinants are ubiquitous in field theory, encoding the effects of quantum or statistical fluctuations.Starting from the path integral, they arise whenever one makes a semiclassical or saddlepoint approximation, and hence appear in a wide range of physical phenomena.For example, functional determinants play a central role in vacuum decay [5][6][7], thermal bubble nucleation [8][9][10], the study of solitons [11], and baryon number violation through anomalies [12,13].
The semiclassical or saddlepoint approximation to a path integral is the infinitedimensional generalisation of Laplace's method.For a set of fields ϕ i (x) with action S[ϕ i ], this takes the form Dϕ i e −S [ϕi] ≈ A e −B . ( The largest contribution to the result, B, is equal to the action evaluated at a saddlepoint.The prefactor A comes from fluctuations around the saddle, which at leading order arise quadratically in the action.Effectively the path integral becomes an infinite product of Gaussian integrals that, in turn, results in an infinite product of eigenvalues-a functional determinant [5,8].
The computation of the saddlepoint action B has received much attention.For vacuum decay, many different algorithms have been proposed [14][15][16][17][18][19], and there are at least five software packages dedicated to this computation [1,[20][21][22][23].On the other hand, to our knowledge there are no existing public codes capable of computing the functional determinant.
The computation of A is the focus of this article.Typically, it has been assumed that A is of negligible numerical importance in comparison with e −B , and hence one may be content with an estimate for A based on dimensional analysis.This assumption is often motivated by the exponential form of equation (1), and the familiar rapidity of exponential growth.
Despite appearances, the functional determinant A generically takes an exponential form in field theory.While there is undoubtedly some arbitrariness in the distinction between exponential and prefactor, below we argue that field-theoretic functional determinants naturally give exponentials.
• The prefactor is the sum of one-loop vacuum Feynman diagrams, including disconnected diagrams.This is equal to the exponential of the sum of connected vacuum diagrams, A = e connected .The divergences of the connected diagrams cancel against the one-loop counterterms in the action, leaving a finite remainder.
• The prefactor is the integral over the phase space of fluctuations around a saddlepoint, and hence it is the exponential of their entropy.For a field theory, these fluctuations are an infinite number of harmonic oscillators.The entropy of a single oscillator is S i = O(1), and the total entropy is extensive: A = e i Si .
• Consider a simple example, such as a constant background field.In this case B = x V 0 and A = e − x V1 , where V 0 and V 1 are the tree-level potential and one-loop effective potential respectively.
As a consequence of this exponentiation, after scaling out dimensions, the typical magnitude of log A relative to B is the same as for any other one-loop correction to the tree-level. 1 Furthermore, for phase transitions the tree-level action is typically fine-tuned, thus making the relative impact of A even greater.It is however important to stress that perturbation theory may still work since higher loop corrections are suppressed by additional powers of couplings, and it is but the tuned tree-level action at fault.
To date there have been several calculations of functional determinants in various field-theoretic models.Analytic results have been attained in one spatial dimension [24], in the thin-wall limit [25][26][27][28][29][30], and in a scaleless potential [6,7].Outside these simplifying cases, computations have been carried out numerically [31][32][33][34][35][36][37][38].Calculations of functional determinants are universally involved, and consequently they have yet to achieve widespread usage, either for phenomenological applications or for comparisons to new theoretical methods.Motivated by this, we introduce BubbleDet, a Python package that automatically calculates bosonic functional determinants in spherically symmetric scalar field backgrounds.

Quick start
The simplest way to install BubbleDet is to use a Python package manager, such as the Package Installer for Python (PIP) or Conda.To install with PIP, run the following in a Linux or Unix (including macOS) terminal or in a Windows command prompt 2   $ pip install BubbleDet Alternatively, to install with Conda, use $ conda install -c conda -forge BubbleDet Once installed, BubbleDet can be imported just as any other Python package.
To use BubbleDet to compute functional determinants, one must pass a precomputed bounce solution.Since BubbleDet is written in Python, it is straightforward to obtain and pass the bounce solution from CosmoTransitions [1].This is however optional, and BubbleDet can be used together with any bounce solver.
A number of examples are included with the package, the simplest of which, first example.py,demonstrates the computation of the functional determinant for the vacuum decay of a pure scalar field with a quartic potential.Further examples demonstrate how to use the package in a variety of settings, including for a symmetrybreaking vacuum transition in the Abelian Higgs model and for a high-temperature phase transition in a Yukawa model.All the examples can be viewed in the online documentation, at https://bubbledet.readthedocs.io/,where one can also find comprehensive documentation of all the package functions.
3 Theoretical background

Overview of the problem
A metastable or false vacuum state can decay through bubble nucleation.This is a semiclassical process, dominated by a critical bubble or bounce, which leads to a local escape from the false vacuum.In quantum field theory, the study of this process was initiated by Coleman and Callan in the late 70s [5,14], building on earlier work in classical statistical field theory [8,39].
More recently, the predicted metastability of the electroweak vacuum state has necessitated the development of quantitatively reliable predictions for the decay rate [6,7].This has been bolstered by the possibility of observing gravitational waves from a cosmological phase transition in the early universe [40], for which accurate predictions of the bubble nucleation rate are required to determine the peak frequency and amplitude of the signal [41,42].
Consider a single real scalar field in d-dimensions with Euclidean action 2 Use pip3 in place of pip on systems where pip refers to the Python 2 instance.where we have kept all indices lowered to emphasise the Euclidean signature of the metric.Let us assume that V (ϕ) has one (metastable) minimum at ϕ = ϕ F and a deeper one at ϕ = ϕ T , such as in Figure 1.The metastable, or false, vacuum will then decay to the stable, or true, vacuum with a rate per unit volume given at zero temperature by [5,7,14,[43][44][45]] Here det ′ signifies that the translational zero modes are omitted from the determinant, and ϕ b denotes the bounce solution, also referred to as the critical bubble.We assume that the bounce only depends on the radial coordinate r ≡ √ x µ x µ and satisfies [14] δS subject to the boundary conditions, Here we have introduced the shorthand ∂ ≡ ∂ ∂r .
These boundary conditions are preserved by the functions included in the functional determinant, as they can be considered additive fluctuations about the background.The fluctuations are therefore regular at the origin and go to zero at infinity.
Note that equation (3) only incorporates 1-loop corrections-via the determinant-and higher loop corrections are omitted in this work.

Multi-field models
In more complicated models, the full functional determinant runs over all the fields.In principle these fields can mix, either through the mass matrix, or through derivative terms.However, in the present work we will assume that the functional determinants can be diagonalised in field space.We also assume that that there is only one background field ϕ(r) that is coordinate dependent.
As a concrete example, consider first a scalar theory with a global U(1) symmetry and with Euclidean Lagrangian If we now expand around a radially-symmetric background ϕ(r) that extremises the action, Φ = 1 √ 2 (ϕ(r) + H(x) + iG(x)), we find where the dots denote terms higher order than quadratic in fluctuations.In this case, the quadratic part of the Lagrangian, and consequently the functional determinant, factorises into Higgs and Goldstone pieces.
The generalization of equation ( 3) to this model is then where V X and J X are volume and Jacobian factors arising due to zero modes.The former is the volume of the space of zero modes, and the latter is the Jacobian for the transformation to collective coordinates [43,46].They are discussed further in Section 3.3 and in Appendix D. For the Higgs, we have that . Note that we need not take the modulus of the Goldstone determinant, as only the Higgs determinant contains a negative mode.
Extending this model, consider next the inclusion of an n-component scalar field, coupled as where the index a runs over 1, 2, . . ., n.Assuming the χ field does not take a background expectation value, the Lagrangian expanded to quadratic order reads There are no zero modes for the χ field, as it does not break any symmetry.It also does not mix with the Higgs or Goldstone fluctuations at quadratic order.So, the full nucleation rate reads In all the examples above, the total one-loop contribution to the decay rate is a product of independent functional determinants, each of which takes the form where n is the number of such field degrees of freedom, and J and V are the Jacobian and volume factors for the zero modes.Since we are interested in the rate per unit volume, we always remove the volume associated with space-time translations.The W factor in equation ( 13) denotes a field-dependent mass squared, for example For the above models, a list of W functions determines the Lagrangian at quadratic order.However, in models with off-diagonal terms at quadratic order, W should be promoted to a matrix in field space.
The inclusion of vector fields, such as in the electroweak theory, inevitably leads to off-diagonal terms in field space, through mixing with Goldstone fields [7,38,47,48].Accounting for such mixing terms goes beyond the scope of this first version of BubbleDet.However, one can approximate the full vector-field one-loop determinant by dropping the off-diagonal terms.This can be expected to capture the correct order of magnitude of the functional determinant.Vector fields are discussed further in Appendix B. Nevertheless, while we do not consider full multi-field determinants in this work, several such computations exist in the literature [6,7,38,[49][50][51][52].
Finite temperature transitions At high temperatures it is possible for fields to borrow energy from the thermal bath and escape from a metastable state.In addition, if the nucleating field evolves parametrically slower than the inverse temperature, we can describe its dynamics classically in real time.As such, J.S. Langer's framework of classical nucleation theory is applicable [8,39,53,54], and the rate factorises into dynamical and statistical parts For thermal nucleation in d + 1 spacetime dimensions, the statistical factor A stat coincides with the vacuum decay rate in d dimensions given by equation ( 3), and thus can be computed directly with BubbleDet.Albeit with one caveat: Thermal corrections from nonzero Matsubara modes should be included in the tree-level potential when computing A stat [10,53,55].
The dynamical factor A dyn contains dissipative effects and, unlike A stat , it is not expected to exponentiate.In Langer's framework, which assumes Langevin dynamics, the dynamical factor is equal to the real-time growth rate of a critical bubble divided by 2π .This can be expressed as [39,56,57] in terms of the negative eigenvalue of the functional determinant λ − , and the Langevin damping coefficient η.The negative eigenmode is the lowest eigenvalue of the Higgs fluctuation operator O H (ϕ b ), and corresponds to isotropic growth or shrinking of the bubble.This identification receives corrections at higher orders [54].The computation of λ − can be carried out using BubbleDet, but that of η requires additional real-time input.Setting η = 0 yields the approximation of Ref. [9], though in this limit the saddlepoint approximation is expected to break down [56].

The Gelfand-Yaglom theorem
To find the rate in equation ( 3) we need to evaluate the functional determinant det −∇ 2 + W (r) .For a constant scalar field ϕ one finds the usual effective potential [58,59], however, closed analytical expressions are in general not available for a spatially varying field.Instead numerical techniques are required.As an initial step it is useful to exploit the spherical symmetry and to expand all eigenfunctions in spherical harmonics: where is the degeneracy factor for the orbital quantum number l. Dependence on the magnetic orbital quantum number (normally denoted m) is trivial; it is completely accounted for by the degeneracy factor.The spherical Laplacian is To compute the determinant for given l we use the Gelfand-Yaglom theorem, which in our case states that [60][61][62][63] det where the ψ l b,F (r) satisfy the differential equations with the boundary condition ψ l b,F (r) ∼ r l as r → 0. Note that these equations for ψ l b,F are initial value problems, whereas the corresponding eigenfunctions satisfy boundary value problems.Since W (∞) is a constant, the differential equation for ψ l F (r) can be solved analytically.

One-loop correction to the action
As discussed in Section 3.2, the problem of calculating the rate in equation ( 3) is reduced to solving the differential equations for each l.Given the bounce, ϕ b (r), these equations can be solved numerically.There are however a few complications.First, the determinant vanishes if there are zero modes, so these have to be removed.Second, in practice we can only solve equation (21) for a finite number of l's.And third, the product-or equivalently a sum in the exponent-in equation ( 17) is generically ultraviolet divergent.Let us deal with these problems in turn.
Zero modes If we have a pure scalar theory, all zero modes occur either in the l = 0 or in the l = 1 determinant.The procedure to remove zero modes is equivalent for the two cases so we focus on the latter, and refer to Appendix D.3 for the l = 0 case.For a single scalar, equation (22) with l = 1 gives which has the solution ψ 1 b (r) ∝ ∂ϕ b (r).Note that the determinant indeed vanishes since lim r→∞ ∂ϕ b (r) = 0.
Formally one can remove these zero modes-which arise because the bounce breaks the translation symmetry-by using collective coordinates [43,46].This means that we re-express eigenfunctions that generate the symmetry as a coordinate shift for all other eigenfunctions.Then, since everything is transitionally invariant, the integration over all possible translations gives the d-dimensional volume V; we also have to include a Jacobian factor, J , since we changed variables.Our job now is to find the value of the determinant once zero modes have been removed.
To do this we follow [31,34,36,64] and deform the equation to which effectively shifts all eigenvalues by k 2 .The determinant without zero modes is then reproduced by the following limit After taking the k → 0 limit one finds [36] S[ϕ b ] 2π Here ϕ ∞ is defined from the asymptotic behaviour of ϕ b (r) as r → ∞ via In addition, from equation ( 4) we find that For potentials in which the nucleating, or Higgs, field is massless in the metastable phase, W (∞) = 0, the calculation in [36,64] needs to be modified.In this case the asymptotic behaviour of the bounce follows from the m F → 0 + limit of equation (27), here assuming d > 2. 3 While the derivation differs from the massive case, the final result for the determinant agrees with equation (26).It is worked out in Appendix D.2.
Analytic solution for large l For asymptotically large l, the computation of the determinant simplifies.Physically we can think of l as the classical orbital momentum l ∼ ⃗ p × ⃗ r; this means that the source at the origin-the critical bubble background-becomes less significant when l ≳ mR, where R is the bubble radius and m is the particle mass.So for large l we can use a WKB approximation to solve equation (21) analytically [35].To do this we use R.E.Langer's method [65] and define Ψ(x) = r d/2−1 ψ(r) together with a change of variables to x = log r.Equation ( 21) is then equivalent to For large l we can solve this equation in powers of l −1 .We leave the details to Appendix E and merely quote the leading order solution [35] log where we have defined In practice we use the Gelfand-Yaglom theorem to solve equation ( 21) numerically for l up to some l max .We then use the WKB approximation to solve equation (21) analytically for all remaining l's.

Divergent sums and renormalization
From equation (30) we see that for large l The degeneracy factor scales as l d−2 , so in general there is a divergent sum for all d ≥ 2. Take for example d = 4.In this case we also need O l −3 terms: Because the degeneracy factor is deg(4; l) = (l + 1) 2 , both terms diverge once we sum over l.To handle this we use dimensional regularization and set d = 4 − 2ϵ.The only sum with a pole is of the form In the MS-scheme we should multiply this expression by exp(γ)µ 2 4π ϵ and add counterterms-these can directly be read off from the (constant-field) 4 Coleman-Weinberg potential [58] : which fixes V ct (ϕ).Here C = 1 for Higgs and Goldstone determinants and C = d − 1 for vector determinants.After integrating over the volume we find Note that this contribution enters the rate as e −Sct [ϕ] .
Putting everything together and using gives where (C, a) = (1, 0) for scalars, and (C, a) = (d − 1, 1/(d − 1)) for vectors.Similar terms appear in all even-numbered dimensions, albeit they multiply different terms in the WKB approximation.To renormalize our theory in d dimensions we need all terms up to l −d+1 in the WKB approximation.The full arrows show a strong connection between components which is structurally essential, while the dashed arrows show a weaker connection, such as a particular instance of usage.This first example can be found at https://bubbledet.readthedocs.io/.

The BubbleDet code
The overall structure of the BubbleDet package is outlined in Figure 2. The package defines three Python classes: • BubbleConfig describes the background field.To initialise an object of this class requires a scalar potential V (ϕ), the false vacuum ϕ F , a bubble profile ϕ b (r), and the dimension d. • ParticleConfig describes a fluctuating particle.It is initialised by the function W (ϕ), the spin of the particle s, the number of its internal flavour or colour degrees of freedom n, and a flag denoting the type of zero modes present.• BubbleDeterminant computes the functional determinant.It is initialised by one BubbleConfig instance and a list of ParticleConfig instances.
For more details, see the documentation and examples.
Once an instance of the class BubbleDeterminant is initialised, the most important class method is which computes a functional determinant in its totality, regularised in MS, and with zero mode factors included where necessary.In this convention, the result is an additive correction to the action.Here the index i runs over the ParticleConfig list.For scalar fields the factor dof(d, 0, n i ) = n i , and for vector fields dof(d, 1, n i ) = (d − 1)n i (see Appendix B).Note that while Jacobian factors are included in all cases where zero modes exist, we do not include the (infinite) volume factor for the zero modes corresponding to translations.As a consequence, exp(−findDeterminant()) has units of mass d for the Higgs determinant.
For classically conformal models the size of the bounce needs to be integrated over to obtain the rate.In BubbleDet this integration does not change the mass-dimension of the rate.See the discussion in Appendix D.4.
For thermal transitions, one should pass the keyword thermal=True.In this case, the dynamical prefactor term − log A dyn is added to the output, given in terms of the negative eigenvalue and assuming zero damping coefficient.
The BubbleDeterminant class also contains a number of ancillary methods for computing different specific parts of the prefactor, such as the negative eigenvalue. Figure 2 shows a component diagram of the structure of the package, and how it can be used in conjunction with CosmoTransitions to compute the vacuum-decay rate for our introductory example.

Application of the Gelfand-Yaglom method
From the Gelfand-Yaglom theorem, the determinant for a given orbital quantum number l is given by the ratio ψ l b (r)/ψ l F (r) at r = ∞.While this ratio is finite, both ψ l b (r) and ψ l F (r) grow exponentially at large r, which can lead to overflow for floating-point numbers.To avoid this, one can work directly with the ratio,6 This function solves the following initial value problem with T l (0) = 1, ∂T l (0) = 0, and where we have defined in terms of modified Bessel functions I.We integrate the first step of the initial value problem, from 0 to δr, by using a Taylor expansion around the origin, and making use of the initial conditions; see Appendix C.1.We then proceed by using the fourth-order Runge-Kutta method to evolve from δr to some r max , the largest radius at which the bounce profile is given.Errors on T l (∞) arise from two main sources: discretisation errors due to non-zero radial steps δr > 0, and the error due to non-infinite r max < ∞.The former error scales as (δr) 4 as long as the bubble profile is known to this accuracy. 7The error due to the finite maximum radius is smaller than any power of 1/r max , for m F > 0, and for sufficiently large r max , and we estimate it from the value of the derivative at the final step, ∂T l (r max ).The total error on T l (∞) is simply estimated as the larger of these two sources of error.
For the massless case m F = 0, the solution converges relatively slowly as r → ∞ so additional methods are adopted to accelerate convergence; see Appendix C.2.

Extrapolating the sum over orbital quantum number
After extracting the l = 0 and 1 modes, which need separate consideration, the logarithm of the complete functional determinant involves a sum over l from 2 to ∞.The large l divergences of this sum cancel against the one-loop counterterms of equation (34).For a single scalar field, the renormalised sum is thus Figure 3: Convergence of the WKB approximation at large l.Shown is the relative difference between the complete log T l , calculated using the Gelfand-Yaglom theorem, and a number of WKB approximations to it.Higher order WKB approximations are used to accelerate the convergence of sums over l.The specific model is a real scalar with potential given by Eq. ( 50), at α = 0.5 and in d = 3.The data for this plot is produced by the example wkb.py.
To avoid divergences in intermediate computations, and to speed up convergence, we add and subtract a WKB approximation of T l : The WKB factors log T (WKB) l are individually finite.Including them in the summand cancels the large l dependence of log T l up to o(1/l), which ensures that the sum converges.In practice, this sum is truncated at some large value l max ≫ 1, with the residual scaling as an inverse power of l max , and the term labeled ≈ 0 in equation ( 42) is dropped.
To accelerate the convergence of this sum we utilise sequence acceleration, adopting two different methods: the epsilon algorithm [67] (which implements an iterated Shanks transformation [68]), and a fit extrapolation, based on fitting a polynomial in 1/l.For the latter, the order of the polynomial is chosen to minimise χ 2 per degree of freedom.Typically the epsilon extrapolation performs better than the fit extrapolation when l max is relatively small, but performs worse for larger l max as numerical errors build up.We choose the final extrapolated result as that with smaller estimated error.
The terms log T (WKB) l are computed within the WKB approximation, the expressions for which are collected in Appendix E. Higher order terms of the WKB expansion are suppressed by higher powers of 1/l, so that one need only compute a finite number of WKB terms to attain a given rate of convergence in 1/l.By carrying out this approach to higher and higher orders, one can achieve much faster convergence of the sum, as demonstrated in Figure 3.
Note that in higher dimensions the sum over l becomes increasingly ultraviolet divergent, so higher orders of the WKB expansion are needed to renormalise the sum and to accelerate convergence.This is because the degeneracy factor grows faster in higher dimensions, as deg(d; l) ∼ ld−2 where l = l + d 2 − 1, while the expansion for log T l takes the same form in all dimensions, where R is the radius of the critical bubble.When subtracting terms with the WKB approximation we therefore need to know (at least) the first ⌊ d 2 ⌋ terms to cancel all divergences.
An important numerical consideration is that one must compute each log T l to higher accuracy in higher dimensions, as the necessary cancellations are more delicate.For example in d = 4 and choosing l max ∼ 25, we need to know the numerically obtained log T l to a relative accuracy of l max /l 3 max ∼ 10 −3 .For d = 6 we also need the next term, so we need to know log T l to l max /l 5 max ∼ 10 −6 .In BubbleDet we utilise the WKB approximation up to and including log T (9) l ; see Appendix E. We do so for all dimensions, thereby ensuring a relatively rapidly converging sum in lower dimensions.Conversely, to achieve a fixed accuracy requires larger l max in higher dimensions.

Finding the asymptotic behavior
Here, we will discuss the implemented numerical method for finding log ϕ ∞ , which is defined by Once we have determined log ϕ ∞ , it is straightforward to deal with zero modes [36,69].
See for example equation ( 26).An estimate for log ϕ ∞ can be found by fitting the numerical bounce to equation (44) for some sufficiently large radius.
Still, this must be done carefully since a numerically obtained profile will not exactly follow the asymptotic behaviour in equation (44).Large deviations may happen for example when a shooting algorithm stops as the bounce tail crosses the metastable vacuum.As the very tail end of a numerical bounce has little effect on the corresponding action, its accuracy is often not a high priority.
The need for caution is illustrated in Figure 4: The orange line is obtained by directly estimating log ϕ ∞ at different radii.The bubble is produced with CosmoTransitions in a model in three dimensions with scalar potential, Note how the estimate begins to diverge at large radial distances, r > 10.Also, there are oscillations, visible on the right pane, which follow from the use of cubic splines in CosmoTransitions.The peaks of the oscillations correspond to the best estimates as they are the nodes of the splines.
To evade these problems at the tail-end, we have constructed an algorithm that is robust even when the tails are inaccurate.In the following, we focus on the massive case, m F > 0. An example result from the algorithm is plotted in blue in Figure 4, along with the direct estimates.
The algorithm takes two user-given parameters: tail, which determines how close to the end of the profile to fit, and log phi inf tol, which gives an estimate of the relative uncertainty on the tail of the profile.
The main algorithm then finds log ϕ ∞ by extrapolating from the chosen tail; see Figure 12 and the surrounding discussion.If there are no such suitable tail points, if Figure 4: Estimates for log ϕ ∞ by directly fitting the asymptotic tail behavior, Equation (44), to points on the numerical critical bubble solution for the potential in Equation (45), and the result from BubbleDet with uncertainties.The right panel zooms in on the region of the left pane in which the estimate is nearly constant.
for example the bounce profile is too short, the code resorts to a fail-safe algorithm, which instead performs the fit at r = √ r 1/2 r max , where r max is the largest radius.The error is then estimated by comparison to fits at radial distances (r 2 1/2 r max ) 1/3 and (r 1/2 r 2 max ) 1/3 .In case the parameter log phi inf tol does not correctly reflect errors in the bubble profile, the more sophisticated main method is run a second time but with log phi inf tol modified to give at most the error of the fail-safe method.Together with the first run of the main method and the fail-safe method, this gives three estimates for log ϕ ∞ , of which that with the least error is returned.
We refer to Appendix C.3 for a more detailed description of the fitting algorithm, and for a discussion of the massless case, m F = 0.

Computing the negative eigenvalue
For thermal nucleation, the decay rate factorises into a product of statistical and dynamical parts, as is shown in equation (15).While the statistical part can be computed through application of the Gelfand-Yaglom method, this is not so for the dynamical part.In the latter, the negative eigenvalue of the Higgs operator O H (ϕ b ) appears in combination with the real-time damping rate.For computation of the negative eigenvalue, BubbleDet provides the function findNegativeEigenvalue().
The negative eigenvalue, λ − < 0, is defined by the following eigenvalue problem, where ∂f (r) → 0 as r → 0 + and f → 0 as r → ∞.Note, that here we have used the information that the negative eigenmode is spherically symmetric, l = 0.
The eigenvalue problem can be approximated as a discrete matrix equation, Here, the matrix M is a discretization of the linear differential operator in Eq. (46), such that all of the rows (indices i) correspond to individual locations of r, and N is the number of discrete points in the numerical critical bubble.The code obtains a numerical estimate of the negative eigenvalue by passing M to scipy.sparse.linalg.eigs.This estimate is then improved by an extrapolation: The implemented derivatives, discussed below, are accurate to such an order that the errors in λ − decrease as N −4 .Obtaining the estimates with one half and one third of the points allows for an extrapolation to N → ∞ and also an error estimation.
The boundary condition at infinite radius, r → ∞, cannot be implemented accurately in the discretized version.To estimate the size of the corresponding error, in the code there are two alternative boundary conditions implemented at the maximal numerical radius, r max : ∂f (r max ) = 0 and f (r max ) = 0 . ( In the limit r max → ∞, these both reduce to the correct boundary condition for the negative eigenmode.Deviations are exponentially small for a large enough maximal radius, ∝ exp(−2 m 2 F − λ − r max ).Values with errors are computed for each of the two boundary conditions, extrapolated to N → ∞.The final result and error are given as the midpoint and extent of the combined uncertainty range.
We refer to Appendix C.4 for more details.

Requirements on the input bounce profile
Before computing the functional determinant, one must solve for the bounce.While BubbleDet does not provide this functionality, there are a number of other packages which are built for this purpose: CosmoTransitions in Python [1], AnyBubble and FindBounce in Mathematica [20,23], and BubbleProfiler and SimpleBounce in C++ [21,22].As a word of caution, for computation of the functional determinant, the bounce profile has to be known rather accurately, and out to rather large r.Unlike the bounce action, the functional determinant depends strongly on the large r asymptotics of the bounce, and further fluctuations with high orbital quantum numbers probe the short scale structure.With CosmoTransitions we have found that using tolerance parameters xtol=1e-9 and phitol=1e-9 typically yields an accurate enough bounce profile for computing the functional determinant reliably.Note that this is significantly more accurate than CosmoTransitions's default values.

Comparisons to literature
We have carried out extensive tests of the BubbleDet code against the literature, comparing against all the explicit results we could find.We have generally found agreement within quoted errors.
A number of authors have numerically computed the one-loop vacuum decay rate for a single real scalar field in a generic quartic tree-level potential [34,38,64], Here the subscript in V 4 refers to the highest power of ϕ present in the potential.By using the result in Appendix A.1 the action can be brought to the form where α = λm 2 η −2 and β = m 6−d η −2 .This dimensionless form is quite useful since the determinant is-up to a − d 2 log β factor-only a function of α.Note that a term linear in ϕ can be removed by a shift ϕ → ϕ + const.
Refs. [34,64] have carried out the computation in this potential for d = 4.We find agreement to better than 1% with the tabulated results of Ref. [34] for the full d V (ϕ) of parameters studied, α ∈ [0.01, 0.99].With Ref. [64], we find good agreement with their figures 5 and 6, though we disagree on their figure 7.
We have also compared the result from BubbleDet to exact analytical results for classically scaleless potentials; for example the four-dimensional potential V (ϕ) = − 1 4 λϕ 4 [70,71].Previous results exist in four dimensions [6,7,30], and as an additional crosscheck we have also derived analytical results for three and six dimensions; to our knowledge these results don't exist in the literature, but the derivation is identical to the four-dimensional case and won't be repeated here 8 .A summary of the results can be found in table 1.
Figure 5 compares the results from BubbleDet to the exact results for scaleless potentials given in table 1.There we have introduced S 1 to denote the one-loop correction to the effective action, i.e. minus the logarithm of the functional determinant prefactor.In all cases we find agreement at better than 1%, excepting a very narrow range where S 1 goes through zero in d = 4, and in most cases the agreement is better than 0.1%.For this model, exact results are also available for the decomposition of the determinant into orbital quantum numbers l, and for fixed l we find even better relative agreement, as shown in the example file unbounded.py.
In d = 3, Ref. [38] numerically computed the Higgs, Goldstone and vector field determinants for the same potential, given by equation (49), as well as for which, after using the result in Appendix A.2, becomes For the V 4 potential in d = 3 we find agreement with the fit functions of Ref. [38] to the 1% level, and for the V 6 potential we find agreement to the 1% to 5% level, except for where S 1 goes through zero, near α ≈ 0.4.This agreement matches the expected accuracy of the fits.

Thin-wall limit
In the thin-wall limit there are a number of analytic results for the one-loop vacuum decay rate, in d = 2 [27,28], d = 3 [26,30] and d = 4 [25,30], all of which use the potential of equation ( 50) up to trivial scalings and shifts.We also used a method similar to that in [30] to derive new analytical results for d = 5, 6 and 7.
Figure 5: The Higgs determinant, for given R, compared between BubbleDet and the analytical results for the scaleless potentials given in Table 1.The left-hand plot shows the relative error as a function of R. Note that larger errors in d = 4 arise because S 1 (R) ≈ 0; this can be seen clearly in the right-hand plot where the relative difference is instead normalized by S 1 (R = 1).All lines use the renormalization-scale µ = R −1 and λ = 0.1.The analytical bounces, given in table 1, are supplied to BubbleDet, which then calculates the one-loop action using default settings.
With our potential conventions, and in the MS renormalisation scheme with the MS scale equal to the mass, these results read In Figure 6 we show how our numerical results approach the analytic thin-wall results in the approach to the thin-wall limit.Fitting a cubic function of α to the data yields extrapolations for α → 1 which agree with the expected analytical values to better than 1% accuracy for dimensions d = 2, 3, . . ., 6.For d = 7, the increased sensitivity to shorter distance fluctuations requires including higher orbital quantum numbers (see Section 4.2).This in turn prevents us from reaching 1 − α ≲ 0.05 while keeping numerical errors under control, and hence significantly worsens the accuracy of the α → 1 extrapolation in d = 7, which agrees to only about 10% with the analytic result.
Note that the thin-wall limit is a particularly difficult regime numerically, as to resolve the scale hierarchy m F R ≫ 1 requires very high orders in the sum over orbital quantum number l, which leads to accumulation of numerical errors.As a consequence, we expect BubbleDet to perform better, and to achieve higher accuracy away from the thin-wall limit, where analytic results are lacking.
For d = 1, there is one Euclidean time dimension and zero spatial dimensions, and hence we are treating tunnelling in quantum mechanics.At α = 1, the potential of equation ( 50) corresponds to the double-well potential, which admits kink or soliton solutions [24].These are quantum mechanical instantons.The functional determinant in a kink background gives the one-loop correction to the splitting between the even Figure 6: The approach to the thin-wall limit in d = 2, 3, . . .7. Unfilled circles show the analytic results valid in the limit α → 1. Crosses show data computed using BubbleDet, and the lines show cubic fits to this data.We also show the nearest equivalent for d = 1, the one-loop correction to the energy eigenvalue splitting.The data for this plot can be generated using our example scripts thinwall.pyand kink.py.
and odd energy eigenstates [69,72].Evaluating this analytically, one finds With default settings, BubbleDet reproduces this one-loop result to a relative accuracy of 0.001%.This is shown in Figure 6.Similar agreement is found for the sine-Gordon theory, V s-G (ϕ) = 1 − cos(ϕ), which also admits kink solutions in d = 1 and for which S 1 = 1 2 log π − 2 log 2.

Derivative expansion
For fluctuations of some scalar field χ, much heavier than the background Higgs field ϕ, a derivative expansion of the functional determinant may be possible.This is an expansion in powers of the ratio of masses squared, m 2 ϕ /m 2 χ .Consider, for example, so that m 2 χ = g 2 ϕ 2 .For sufficiently large g 2 , the functional determinant for χ should be well approximated by a derivative expansion.At leading order (LO) and next-toleading order (NLO), the result takes the form where V 1 and Z 1 are the contribution to the effective potential and field normalisation for ϕ due to fluctuations of the χ field.Expressions for Z 1 in arbitrary dimensions can be constructed from the results of Ref. [66].
In Figure 7 we plot the relative difference between the full functional determinant and its LO and NLO approximations within the derivative expansion, for dimensions d > 2. For d = 2 the derivative expansion in this model is infrared divergent (as r → ∞) at NLO, hence we do not include it. 9For d = 3 the derivative expansion is Figure 7: Deviations from the derivative expansion for the determinant of the χ field; see equation (56).Here the value of m 2 ϕ /m 2 χ is given in the true vacuum.Crosses connected by dashed lines show the LO derivative expansion, while unfilled circles connected by full lines include the NLO order corrections.For small m 2 ϕ /m 2 χ , convergence appears to arrest at about 10 −5 , especially in higher dimensions.This is because one must reach very large values of the angular quantum numbers l, and in so doing numerical errors accumulate (here on default tolerance settings).infrared divergent at next-to-next-to-leading order (NNLO), and in fact there exists a term between NLO and NNLO which is invisible to the derivative expansion; see for example Ref. [53].This explains why the slope of the NLO line in Figure 7 is approximately 3/2 and not 2.In general as the coupling g 2 is increased, agreement with the derivative expansion improves, reaching 0.001% at NLO in dimensions where the derivative expansion is under control.This provides a nontrivial test of BubbleDet, which works as expected.

Speed tests
For speed comparisons we compare the time it takes to calculate the determinant relative to the time it takes for CosmoTransitions to obtain the bounce solution.
From Figure 8 we see that, on average, it only takes twice as long to find the bubble determinant as the bounce.Note that our use of a sequence acceleration, combined with using more terms in the WKB approximation, is crucial.

Applications
As argued in the introduction, functional determinants in quantum field theory exponentiate, and can thereby yield important corrections to nucleation or decay rates.Here we illustrate this, applying BubbleDet to a range of scenarios.50) with d = 3. Computations were carried out on a laptop with a 3.2 GHz Apple M1 processor.Each data point represents the average of 300 runs.The required time to calculate the bounce action with CosmoTransitions is shown in orange, and the required time for the determinant is shown in blue.The latter increases in the approach to the thin-wall limit due to the necessity of reaching larger orbital quantum numbers, l max ≫ m F R, where R is the bubble radius.

Thermal bubble nucleation and dimensional reduction
following Euclidean action (in d = 4) The notation is standard and follows Ref. [53].The thermal evolution of the longwavelength modes of the field Φ are described by a dimensionally-reduced effective field theory (in d = 3), which is where to leading order ϕ is equal to the zero Matsubara mode of Φ divided by √ T , and the effective parameters are By shifting ϕ → ϕ + const and scaling, this model can be made to agree with the conventions of our V 4 .However, here we work directly with the dimensionful parameters of the effective field theory.As given in equation ( 15), the thermal nucleation rate is the product of two terms.The statistical part A stat is the vacuum decay rate for the 3-dimensional EFT.For  57).Here we show the variation with temperature for one physical parameter point.The data for this plot can be generated using our example script thermal.py.
the dynamical part A dyn we neglect damping.The total decay rate is then The fluctuation determinant runs over the degrees of freedom of the EFT, i.e. the 3-dimensional scalar field ϕ.The fermion contributes to the nucleation rate through the temperature dependence of the effective parameters.Note that there is no factor of 1/T in the exponent, as factors of T have been absorbed into the effective parameters and field.
In Figure (9) we show the nucleation rate for an example thermal history of this model, with parameters {s, m 2 , g, λ, m Ψ , y} = {0, −1, 0.3, 0.1, −0.2, 0.3}.At this parameter point the critical temperature is T c ≈ 8.42|m|, and the one-loop correction to the jump in the order parameter ∆⟨ϕ⟩ c is about 8% [73].However, the relative importance of loop corrections grows as the system supercools.Figure (9) shows that the functional determinant can yield sizeable corrections to the rate, and that these eventually become nonperturbative in the approach to spinodal decomposition, at T sp ≈ 8.19|m|.

Gauge symmetry breaking
Let us consider the Abelian Higgs theory in d = 4, with Lagrangian given by equation (81).To simplify the analysis [48], we include a (ϕ * ϕ) 3 operator in the potential to generate a tree-level barrier between symmetric and broken phases; see equation (51).
The full nucleation rate is then given by (see Appendix B) We remind the reader that we have dropped off-diagonal gauge-Goldstone mixing.
Figure 10 shows the tree-level and one-loop corrections to the nucleation rate for the parameter point {m 2 , λ, c 6 } = {1, −0.15, 0.01}, varying the gauge coupling.as the contribution to the 1-loop effective action from the gauge boson.As can be seen, this dominates over the scalar contributions for m H /m W ≲ 1.The data for this plot is for fixed scalar parameters m 2 , λ and c 6 , and varying gauge coupling g 2 .It can be generated using our example script symmetry breaking.py.
For mass ratios below about m H /m W ≲ 2/3, or gauge couplings above g 2 ≈ 1.3, the functional determinant of the gauge field grows larger than the tree-level bounce action.For such large gauge couplings, the gauge determinant is exponentially enhanced by the ratio of gauge to scalar masses, or equivalently by g 2 /λ ≫ 1.In this case the broken perturbative expansion can be repaired by integrating out the gauge field before solving the bounce equation [74].

Analogue false vacuum decay in d = 1 + 1
In studies of analogue false vacuum decay, a nucleating scalar field with an effective relativistic energy-momentum relation can arise as an angular degree of freedom in cold atom setups. 10In this case a trigonometric potential arises [76][77][78] The authors of Ref. [78] investigated the effects of renormalisation in classical stochastic lattice simulations of vaccuum decay, studying parameter points λ = 2, v2 ∈ [2, 10], and we presume m 2 = 1.They suggested a definition for a lattice effective potential aiming to capture "the renormalization effects present in the fluctuation determinant", though the latter was not directly computed.Using the lattice effective potential led to a roughly e 30 increase of the nucleation rate relative to tree-level, which was almost independent of v over the range studied.
For the parameter points studied in Ref. [78], we have used BubbleDet to further investigate their proposal.In support, we indeed find that the functional determinant leads to an increase in the nucleation rate which varies little with v over the range studied.However, we find the magnitude of the increase to be significantly smaller, being only a factor of ∼ e 7 .This calculation can be found in our example script analogue.py.A full comparison would require matching regularisation schemes.However, recent work has cast doubt on the continuum limit of these classical stochastic lattice simulations [79,80].

Effects of potential shape
To better understand the effects of potential shape on the functional determinant, it is useful to work with the dimensionless potentials introduced above, which have been put into a common form.The potentials are linear in α, and α → 0 corresponds to the thick-wall limit, while α → 1 corresponds to the thin-wall limit.In addition, the determinants only depend nontrivially on α, and not β.
In addition to the polynomial V 4 and V 6 potentials discussed above, we introduce the following logarithmic potential This potential arises in classically scale-invariant models, and its thermal evolution typically exhibits strongly supercooled phase transitions [81,82].Equivalently (see Appendix A.3) we can consider the dimensionless potential where α ∈ (0, 1) and the factor in-front of the action is We also rewrite the trigonometric potential from equation ( 63) in a dimensionless form by using the result in Appendix A.4: where (1 − α) = 1/λ 2 and the prefactor is All the potentials we consider are collected in Appendix A, together with their scaled dimensionless forms.In Figures 11a, 11b, 11c, and 11d we show the magnitude of the functional determinants for each potential relative to the corresponding treelevel bounce actions, focusing on d = 4.Note that we have scaled out β, as well as any group-theoretic factors from symmetry breaking.We see that the behavior, and size, of the functional determinants varies greatly depending on the potential.This is no surprise, as the bounce solution, and the corresponding functional determinants, explore the global properties of the potential, not just the minima.
Notably, for the V 4 and V cos potentials, the ratio S 1 /S 0 , for the Higgs determinant, tends towards a constant in the thin-wall limit (α → 1), while it grows as (1 − α) −1 for the V 6 and V log potentials.This is because in the thin-wall limit the leading behavior is determined by the (free-energy) pressure difference between the false and true vacua: ).The quartic and cosine potentials are special because V ′′ (ϕ), and thus the pressure, coincide in the two phases at the critical temperature, due to a symmetry between the two phases.So, similarly to the tree-level action, log det O H behaves as R d−1 .For the other potentials, and particles, this is not the case.And in general all determinants scale as R d ; the coefficient can even be found from the one-loop effective potential.

Conclusions
BubbleDet enables computing functional determinants for a wide range of applications.This is particularly important for the quantitative reliability of studies of cosmological phase transitions.The common approximation A ≈ T 4 for the nucleation prefactor (see equation ( 1)) is the origin of one of the largest sources of theoretical uncertainty in current predictions of the gravitational wave signal, as revealed through residual renormalisation scale dependence [41,42].BubbleDet marks a step in overcoming this theoretical uncertainty, in preparation for planned gravitational wave observatories such as LISA [83].
The code has been thoroughly tested, and reproduces a large number of results found in the literature.It is robust, and has been shown to yield sub-percent errors Figure 11: Size of Higgs, Goldstone, and vector determinants for various potentials.Figure 11a corresponds to the potential in equation (50); figure 11b corresponds to the potential in equation ( 52); Figure 11c corresponds to the potential in equation (65); and Figure 11d corresponds to the potential in equation (66).All figures use the dimensionless gaugecoupling g = 0.5.In addition, the unit of all plots is in β −1 ; so increasing β from 1 to 2 scales all curves by a factor of 1 2 .Note that the scales are different for each plot.
from the thick to the thin-wall limits and for a wide range of potentials.There is rarely need to tweak meta-parameters.All the results shown in the present paper were computed using default method and tolerance options. 11ubbleDet can be installed just as any other Python library, and is available in the PyPi and conda-forge repositories.It is easy to interface with CosmoTransitions, so that existing scripts using this library can be straightforwardly upgraded to include functional determinants.Computing a functional determinant typically takes less than about a second on a standard laptop, fast enough to do so for parameter scans of models beyond the Standard Model.
The most important extensions for BubbleDet are to allow for multiple background fields, and for couplings between fluctuating field degrees of freedom, i.e. off-diagonal quadratic terms in the action.These extensions are planned for a second version of BubbleDet, and would allow one to tackle, for example, multi-scalar extensions of the Higgs sector.

A Potentials used in the paper
To compare the effects of different potential functions, we shift our fields and coordinates such that every tree-level action is of the form where α ∈ (0, 1) and α → 0 + corresponds to the thick-wall limit; and α → 1 − to the thin-wall one.We further ensure that the scaled tree-level potentials are linear in α In this form the Higgs determinant is independent of β, up to a − 1 2 log β contribution from each zero mode.The tree-level bounce action therefore dominates if β ≫ 1, and perturbation theory works well.

A.1 A ϕ 4 potential
Consider the potential We always have the freedom to re-scale coordinates and to shift our fields.Meaning that we can eliminate three parameters in favour of two dimensionless ratios.For example, in the literature it is common to use which means that the action and potential become which, after scaling x → m −1 x, ϕ → m|λ| −1/2 ϕ, becomes

A.3 A logarithmic potential
Consider the potential We perform the following re-definitions After which we find the dimensionless potential where Here W is the Lambert W function and β is the overall factor multiplying the action.

A.4 A trigonometric potential
Take the potential To rewrite this potential in our preferred form we use to find where (1 − α) = 1/λ 2 and the overall factor is

B Vector fields
As a minimal example, consider the gauged version of the U(1) model of equation 6, the Abelian Higgs model, where We choose the class of Fermi gauges, with gauge-fixing parameter ξ and ghost field c12 .
If we again assume a radially-symmetric background ϕ(r), and expand in fluctuations about this, we find For this model, the quadratic part of the action is not diagonal in this field basis, as one can see from the terms on the last line.
The first version of BubbleDet is not able to accommodate such off-diagonal terms.However, counting physical degrees of freedom, one expects roughly that the result should be expressible in terms of d−1 massive vector degrees of freedom, one Goldstone degree of freedom, and one scalar. 13This counting is apparent in the one-loop Landaugauge effective potential V 1 (ϕ), for which the off-diagonal derivative terms are zero, where we have defined From this, one expects For the one-loop contribution to the rate, the off-diagonal derivative terms are not small, though we expect that neglecting them should nevertheless give a reasonable indication of the order of magnitude.
For models with larger gauge groups, or more scalars, one can likewise neglect any mixing.The final result is then always in a form similar to equation (85).

C Algorithm details C.1 Initial value problem
For the initial value problem in the application of the Gelfand-Yaglom method, one must be mindful to avoid 1/r singularities at the origin.We integrate the first step of the initial value problem, from r = 0 to δr, using so that After the first step, there are no further coordinate singularities, and we update with the fourth order Runga-Kutta algorithm.

C.2 Gelfand-Yaglom method for massless bounces
Bounces that behave as ϕ b ∼ r −(d−2) for large r require special treatment.In this case the Gelfand-Yaglom solution converges slowly, and to get a decent approximation the grid must be pushed to ever-larger r.   2) .We then solve equation ( 21) numerically from r = 0 to r = r max , and use the solution at r = r max to solve equation ( 21) analytically from r = r max to r = ∞.The asymptotic behaviour of ∆W , defined in equation (31), depends on the form of the potential.To handle general cases we perform a fit to ∆W ≈ W ∞ r −a∞ for large r, for which equation ( 21) can be solved analytically.This procedure significantly improves convergence even with a small r max .We similarly improve the WKB approximation by performing various integrals analytically from r = r max to r = ∞.

C.3 Algorithm for fitting ϕ ∞
The basic idea in the main method for the massive case is to find a quantity for which the direct estimates of log ϕ ∞ behave linearly.This way, one can use the linearity to extrapolate to infinite radius.The algorithm uses the following quantity, Note that the possible numerical inaccuracies of the bounce solution at large radii are packed very near zero in terms of ∆V quad /V quad .This follows from the approximately logarithmic relation with the radius, for some c = O(1) depending on the form of the potential.See Figure 12 for an example.
The algorithm chooses four points (shown as gray squares in Figure 12), which are then used to extrapolate to ∆V quad V quad = 0 .
This corresponds to taking the r → ∞ limit.The four points are chosen according to the two user-given parameters, with default values tail=0.007and log phi inf tol=0.001.These values were tested to be stable across different models.
In an idealised case without any numerical inaccuracies, the four points are chosen to have values tail, 2 tail, 3 tail and 4 tail, as illustrated in Figure 12.However, to account for numerical inaccuracies in the profile and potential, the algorithm performs consistency tests, which may lead to choosing different points.
Starting from r = r max and working inwards towards r = 0, the algorithm chooses the subset of points that pass the following criteria: • Both |ϕ − ϕ F | and ∆V quad /V quad must be larger than the next point away from the centre of the profile.• The floating point errors due to the subtraction in equation ( 90) are very small in comparison to tail.
If any of these conditions fail, the algorithm discards that and all points with larger r.Then, it begins again from the next point, working inwards.Once a point is deemed eligible, the algorithm chooses it as the first chosen point ϕ 1 if ∆V quad /V quad > tail.The rest of the points are chosen if they are above the tail value set by the previous chosen point ϕ i : In addition, the second chosen point ϕ 2 has to have a small enough estimated error due to the proximity of the end of the profile.One can estimate this as This would be exact if the tail satisfied the linearised equation of motion exactly with the boundary condition that ϕ(r max ) = ϕ F .For the second chosen point, we ensure that this is less than log phi inf tol, as a relative error.Overall, this procedure improves the stability of the extrapolation to numerical inaccuracies in the profile or potential.
The algorithm then extrapolates to infinite radius, performing both a quadratic and a linear fit, and using equation (94) to normalise the χ 2 .The final error for the main method is chosen to be the maximum of two estimates: the difference between linear and quadratic fits, and the square root of the covariance in the linear fit.
Note that when improving the accuracy of the numerical critical bubble while keeping tail and log phi inf tol fixed, the error for the resulting log ϕ ∞ does not converge to zero, but rather to a constant.To shrink the errors further requires shrinking the tail and log phi inf tol parameters, in addition to improving the bubble.However, these errors in log ϕ ∞ are typically overshadowed by other sources of errors in the computation of the full determinant.
The algorithm for the massless case is quite different.The power-like behaviour of the asymptotic field makes a direct fit in r more feasible than the massive case.However, errors from the tail-end of a bubble profile only decrease as a power of (r max − r) and not exponentially.The algorithm performs a fit to ϕ(r) Here, the constant ∆ϕ is the other solution to the linearised equation of motion around the metastable phase.It does not adhere to the boundary condition and hence picks up deviations from the correct asymptotic behavior.We estimate the error based on two sources: The square root of the covariance in the fit, ∆ cov ϕ ∞ , and also the error ∆ const ϕ ∞ from the constant ∆ϕ ̸ = 0 being nonzero.The total error estimate for log ϕ ∞ is then This can become infinite or complex, for example when the tail crosses the initial phase at r max .In these cases, the fit is discarded, and the result of the fail-safe algorithm, described in Section 4, is used instead.More generally, the result with the smaller error is returned.

C.4 Algorithm for finding the negative eigenvalue
Here we discuss further details of the implemented discretization.The first and second derivatives in the differential operator are discretized in M so that their values would be exact if f i was a fourth-order polynomial around the evaluation point.This leads to the discretization error for the negative eigenvalue decreasing as N −4 , when increasing the number of points, N .Hence, the derivative and the second derivative require information from five points of j around the evaluation point, i. (Compare with five free parameters in a fourth-order polynomial.)The points are chosen as j ∈ {i − 2, . . ., i + 2}.Thus, the ith row of the discretized matrix looks like M ij = (0, . . ., 0, The chosen implementation of the derivatives is not straightforwardly possible near the boundaries i = 0 and i = N − 1, as the index j would go out of bounds.Near the center of the bubble profile, the matrix M is modified if the profile would continue beyond i = 0 as which implements the boundary condition ∂f = 0 at the origin.At the end of the profile, the modification of the matrix depends on the chosen boundary condition.It corresponds to the continuation of the profile beyond f N −1 as either where the former corresponds to the zero derivative and the latter to the zero value.Finally, we want to discuss one subtle point in the discretisation at r = 0: The second term in the differential operator in equation (46) appears singular.This is however not a problematic singularity.Due to the condition that ∂f (r) → 0 as r → 0 + , one can show that lim Thus, for the first row of the matrix M , corresponding to r = 0, one can use the discretized version of ∂ 2 instead of that of the apparently singular operator D Volume, Jacobians, and removing zero modes

D.1 Breaking translational invariance
Zero modes arise if the bounce breaks a symmetry of the full theory.In particular, consider small deviations about the bounce: where ξ a (x) are a complete set of eigenvectors that we normalize as d d xξ a ξ b = 2πδ ab , and the integration measure is Π a dc a .Let us start with translational symmetries.Since the bounce isn't invariant under coordinate shifts x µ → x µ + a µ , we expect d zero modes; which we denote by ξ µ (x).The idea with collective coordinates is that we can equivalently express these zero modes as a coordinate shift of the bounce14 : As such we can identify ξ µ (x) = N −1 ∇ µ ϕ b (r) and c µ = a µ N from equation (102).To fix the normalisation N we use which implies that N 2 = S[ϕ b ]/(2π).The Jacobian for the transformation of the integration measure for all d modes is then In effect all zero modes are absorbed by the bounce, and the integral over the zero modes yields where V is the spacetime volume.
Translation zero modes arise in the Higgs determinant at l = 1 in the sum over angular quantum number-as ∇ µ ϕ b (r) transforms as a vector.The full result for integration over modes at l = 1 is given in equation (26).

D.2 Translations for massless Higgs potentials
The derivation of equation ( 26), for translational zero modes assumes that the nucleating scalar is massive in the false vacuum, m 2 F = V ′′ (ϕ F ) > 0. As discussed around equation (24), this was done by deforming the l = 1 differential equation into: This shifts all eigenvalues by k 2 , and as such we can remove the final zero mode by using However, this procedure does not have a smooth m F → 0 limit.To see this, let us follow Ref. [64] and express ψ 1,k b (R) via where is the normalized zero mode.In the massive case, we can then solve for ψ 1,k b (R) by using that ψ F +k 2 and ψ 1 b (r) ∼ e −rmF at large r.This needs modification for m F = 0, as the behaviour of the solutions changes.
To proceed it is useful to add k 2 also for the false-vacuum solution and to directly where U 1,k is the consequent k 2 deformation of equation (40).The solution at k 2 = 0 is leading to T 1 (∞) = 0 as expected.Now, similarly to Ref. [64] our strategy is to integrate equation (109) with a suitable factor.To get something useful we choose The r d+1 factor is crucial for the equation to solve: It is designed to cancel the d+1 r term in U 1,k , when integrating by parts.We can now move derivatives from T 1,k to T 1 to arrive at Taking R much larger than any intrinsic scale from ϕ b , so that the V ′′ (ϕ b ) term can be dropped from equation (109), we find the large r asymptotics, for some constants a, b.Using this we can solve for T 1,k (∞) to leading order in k 2 : All in all we find the same expression as in the massive case.As an example, take d = 4 and V (ϕ) = − 1 4 λϕ 4 ; the Fubini-Lipatov instanton is then [70,71] where R > 0 is an arbitrary parameter.With this "bounce" we identify and thus which is in agreement with [6,7].

D.3 Internal symmetries
We next consider Goldstone bosons.In this case our theory is invariant under the action of some internal symmetry group G; we write the linearised action of this group on our scalars as where a runs over the group generators, 1, . . ., dim(G), and i, j run over the indices of the representation in which ϕ transforms.The bounce will generically only be invariant under a subset of these generators T ∈ H ⊂ G, while the other, broken ones, satisfy where we have introduced barred indices a to run over the broken generators T a ∈ G/H.In general then we expect dim(G/H) zero modes.To remove these zero-modes we again express the linear shift in terms of a basis of functions The broken generators rotate one solution, ϕ i b , to another solution.So we can identify ξ i a (x) = N −1 G T a ij ϕ j b (r) and c a = N G ε a .The normalization factor for the Goldstones is and the corresponding Jacobian is J G = N n G G where n G = dim(G/H) is the number of Goldstone bosons.By absorbing these zero modes as an arbitrary rotation of ϕ i b we are left with an integration over the broken subgroup: where V G = vol(G/H) is the volume of the broken subgroup.So for example, a U(1) group that is broken completely gives V G = 2π.Note that V G depends on the group manifold, and not just the algebra, so it is sensitive to the global structure of the gauge group, such as division by a discrete group; see for example Refs.[91,92].These Goldstone zero modes arise at l = 0 in the sum over angular quantum number.To compute the full l = 0 result, we can follow the method of Ref. [64], just as for translation zero modes; see equation (26).The result is

D.4 Dilatations
In models with classical scale invariance there is a zero mode arising at l = 0 from the breaking of scale transformations or dilatations, ϕ(x) → s d−2 2 ϕ(xs), s > 0. For an infinitesimal transformation we take s = 1 + ε with |ε| ≪ 1.In this case the zero mode is ξ s (r) = d−2 2 ϕ b (r) + r∂ϕ b (r), and the Jacobian factor for the transformation to the collective coordinate is Formally J s is infinite-the zero mode is not normalisable-but the full determinant will be finite once we remove the zero mode [7].
Like the integration over translations, the integration over scale factors V s = ∞ 0 ds is divergent in a standard saddle-point evaluation of the path integral.As a consequence, we omit the factor of V s in BubbleDet.However, unlike translation invariance, scale invariance often does not survive quantisation, being broken by loop corrections.One may therefore obtain a finite result by integrating over the scale factor after having computed the functional determinant for the other modes [7].As the instanton itself breaks scale invariance, practically this requires computing the functional determinant for each value of the scale factor and integrating the result.
Proceeding as before, we deduce, The result in equation ( 129) is finite, and for the V (ϕ) = − 1 4 λϕ 4 potential in d = 4 we find in agreement with [7].
Finally, we note that in the literature the dilatation collective coordinate is sometimes chosen differently.For the V (ϕ) = − 1 4 λϕ 4 potential in d = 4, it is often chosen as the scale R in the Fubini-Lipatov instanton, With R as the collective coordinate, the dilatation integration measure is dR, while with our choice of collective coordinate, s, it is ds = dR/R.The total integrated results must agree, so the integrand with s as the collective coordinate is R times the integrand with R as the collective coordinate.

D.5 Special conformal transformations
In addition to dilatations, we can also have zero modes associated with special conformal transformations: For d = 4 we can see from the analytic "bounce" that ψ sc ∝ ∂ϕ b ; so no new zero modes appear.This also holds for conformal bounces in d = 3 and d = 6.

Figure 1 :
Figure 1: Schematic example potential showing a (metastable) false vacuum at ϕ F and a (stable) true vacuum at ϕ T .

Figure 2 :
Figure 2: Component diagram of the BubbleDet package, here shown when BubbleDet is used together with CosmoTransitions.The full arrows show a strong connection between components which is structurally essential, while the dashed arrows show a weaker connection, such as a particular instance of usage.This first example can be found at https://bubbledet.readthedocs.io/.
l) log T l − log T l) log T l − log T (WKB) l ≈0 .

192π 3 λ 2 6Table 1 :
log R + 7 log λ − 18 5 log µR − 16.1573 Results for scale-invariant models with unbounded potentials (λ > 0) in d = 3, 4 and 6.Here S 0 and S 1 denote the tree-level and one-loop effective action respectively, i.e. S 0 = B and S 1 = − log A in equation (1).The numerical factors in S 1 , correspond to various combinations of Riemann zeta functions and Euler-Mascheroni constants.The full rate per unit volume is given by Γ = d(log R)e −S[R] .

For a simpleFigure 8 :
Figure 8: Time required to calculate the Higgs determinant for the potential given in equation (50) with d = 3. Computations were carried out on a laptop with a 3.2 GHz Apple M1 processor.Each data point represents the average of 300 runs.The required time to calculate the bounce action with CosmoTransitions is shown in orange, and the required time for the determinant is shown in blue.The latter increases in the approach to the thin-wall limit due to the necessity of reaching larger orbital quantum numbers, l max ≫ m F R, where R is the bubble radius.

Figure 9 :
Figure 9: Thermal prefactor corrections for a Yukawa model; see equation.(57).Here we show the variation with temperature for one physical parameter point.The data for this plot can be generated using our example script thermal.py.

Figure 12 :
Figure12:The linear behavior of the log ϕ ∞ estimates.The result is obtained from a linear extrapolation from the fit points, which are chosen according to the tail parameter.The failure of estimates happens near zero.Error estimates are plotted but not visible.