Bubble nucleation and quantum initial conditions in classical statistical simulations

Classical-statistical lattice simulations provide a useful approximation to out-of-equilibrium quantum field theory, but only for systems exhibiting large occupation numbers, and only for phenomena that are not intrinsically quantum mechanical in nature. In certain special circumstances, it can be appropriate to initialize such real-time simulations with quantum-like zero-point fluctuations. We will revisit these points, and investigate reports that quantum bubble nucleation rates in 1+1 dimensions can be computed through the classical evolution of such a quantum-like initial condition. We find that although intriguing, the reported numerical agreement between classical-statistical simulations and the quantum nucleation rate in 1+1 dimensions is a coincidence, which is not specific to this choice of initialisation, is parameter and lattice cut-off dependent and disappears as the number of space-dimensions increases from 1+1 to 2+1


The classical-statistical approximation
The classical-statistical approximation (CS) to real-time quantum field dynamics consists in replacing the evolution of the quantum operators (such asφ(x, t)) by classical dynamics of an ensemble of random initial conditions. The ensemble is taken to reproduce the initial correlators of the quantum system, and each random member of the ensemble is evolved by means of the classical equations of motion. The expectation values of observables are then computed as averages over the ensemble.
The CS approximation is reliable only when the occupation numbers (particle numbers) n k of the fields are large, n k 1 (se for instance [2]). For a typical massive scalar field, the field equation reads dφ (x, t).
(1.1) with V nl denoting non-linear self-interactions. When non-linearities are small, individual momentum modes behave as harmonic oscillators, and we may write for the field and canonical momentum operatorŝ Then the occupation number of a field momentum mode is given by where by standard convention, we have taken |f k | 2 (0) = 1/ω k , ω 2 k = k 2 + m 2 . The zeropoint fluctuations (corresponding to the zero point energy of a harmonic oscillator) is the "1/2", while the excitations of the system above the vacuum are the n k . The classical limit corresponds to n k 1 2 . This argument relies on the particle numbers n k which is a quasi-particle concept, valid at weak coupling. The argument may be generalised and made more precise in the context of the Keldysh formalism and Kadanoff-Baym equations for the real-time correlation functions (see for instance [3], and a related discussion in [2]). It is convenient to introduce the "statistical" and "spectral" propagators The real-time evolution of correlators may be expressed through diagram expansions in terms of F and ρ, both for quantum and classical field theory [3,4] (see also [5] for an explicit application and comparison). In the classical approximation, certain diagrams turn out to be absent 1 so that whenever the quantum theory self-energy contains a combination of the form 2 in the classical theory the same diagram has only Σ classical F 2 (x, y).
(1.8) 1 In terms of the Keldysh field basis of φ cl and φq, for instance in λφ 4 -theory, the 3-φq vertex is absent, and any diagram involving this vertex. 2 For different theories, diagrams and self-energies, prefactors may vary. For instance, in λφ 4 -theory, the sunset diagram produces 3F 2 − ρ 2 /4 in the self-energy component for ρ, while F 2 − 3ρ 2 /4 appears in the self-energy component for F [6].
Hence, the classical approximation is good, whenever ρ 2 can in fact be neglected, F 2 ρ 2 . For weak coupling and in the quasi-particle picture, F (n k + 1/2)/ω k , ρ 1, in which case the criterion for classicality again amounts to n k 1. A third fully non-perturbative derivation of the classical-statistical approximation follows directly from the Keldysh contour path integral [7] 3 . In short, whereas the quantum result comes about from averaging over the field variables at all times on the Keldysh contour (all "paths"), the CS approximation amounts to only averaging over the field variables at the initial time, corresponding to the ensemble of initial conditions. The authors of [7] proceed to show that the phenomenon of tunneling in quantum mechanics may be computed from the complete path integral or the Schrödinger equation yielding the same result, while the CS approximation fails to correctly reproduce the tunneling rate. Similarly, the CS approximation fails to describe the famous quantum violation of Bell or Leggett-Garg inequalities [8], even in the free-field limit. In fact, the CS approximation may only be used to compute certain "classical" observables. 4 One obvious example is that the fundamental commutator [φ(x), π(y)] − = iδ d (x−y), which is inherently "quantum", vanishes in the CS approximation. 5 The failure of the CS approximation to reproduce the tunneling rate in quantum mechanics (as given by solving the Schrödinger equation) is a standard (see [9] for a recent analysis in the present context).

Classical-Statistical simulations and the "half "
There is nothing to prevent us from performing CS computations from any initial condition, provided we are able to somehow generate the configurations making up the initial ensemble. One example is the classical thermal equilibrium-like state, parametrized by some temperature T , n k + 1/2 = T /ω k , which up to corrections from non-linear interactions is a fixed point of the dynamics. But evolving some generic initial ensemble amounts to classical field theory from a non-equilibrium initial state, not necessarily with any connection to a quantum system.
As we have seen, only for large occupation numbers (or large F ) can a CS computation be expected to yield a good approximation to the quantum result, and only for appropriate observables. Also, the approximation can in many cases only be expected to hold for a finite time, since the late time asymptotic state is the classical equilibrium 6 rather than the quantum equilibrium.
Most physical phenomena are dominated by some characteristic momentum range, and the spectrum of momentum modes split up into regions with large ( 1), moderate ( 1) and small (< 1) occupation numbers. As long as the modes relevant for the phenomenon of interest are highly occupied, the expectation is that classical dynamics will give a reliable 3 The primary focus of [7] is the subsequent evaluation of the path integral in the Picard-Lefschetz formalism, but the relation to the CS approximation is independent of that further application. 4 The ones involving the φ cl -field and not the φq-field in the Keldysh basis. 5 Note that in actual classical field theory, one may define objects with similar properties, such as the Poisson bracket of canonical variables. But although much of the concrete numerical computation is the same, conceptually classical field theory and the CS approximation to quantum fields are distinct. 6 Which is badly defined in the continuum limit, but has meaning on a finite lattice with a momentum cutoff.
result when applied to all modes. Typical examples include (near-to-)equilibrium systems at high temperature (see for instance [10,11]), large objects such as topological defects or sphalerons (for instance [12,13]), as well as high-occupancy phenomena such as resonances and instabilities (for instance [14][15][16][17][18][19][20][21]). For a few very specific cases, a very special initial condition has been employed dubbed "the quantum half" [14,22,23]. The prescription is to represent an initial quantum vacuum state with n k + 1 2 = 1 2 by an ensemble of classical initial conditions, and evolve the system classically from there. In most cases, this is very problematic, since n k 1 is certainly not satisfied, and the energy density of the initial state is cut-off dependent and divergent.
Another issue is that while the true quantum dynamics ensures that the zero-point fluctuations stay put in each mode [24], allowing only the exchange of the n k between modes, the classical dynamics does not distinguish between the n k and the 1/2 excitations, and will allow all to be exchanged. Extracting energy from the zero-point fluctuations in this way is an unphysical effect, which is negligible if n k + 1/2 is anyway large, but may be very important when n k + 1/2 1/2.
However, one property can make it reasonable to describe the quantum dynamics of a quantum-like "half" initial condition by the CS approximation: For non-interacting fields, the operator field equations are linear, as described above allowing us to expand the Heisenberg field operators as independent time-independent harmonic oscillators. To compute numerically we only need to solve for f k (t), while the a k , a † k are time-independent operators containing the information about the initial state. Since the evolution is linear, it makes no difference whether we evolve from the initial condition f k (0) = 1/ √ ω k and multiply by n k + 1/2 at the end, or whether we classically evolve an ensemble of initial conditions φ k (0) with the property that φ k (0)φ † k (0) = (n k + 1 2 )/ω k . This is the CS approximation, and so for a non-interacting field, the approximation to the evolution is exact, irrespective of n k 7 .
This means that for systems, where for some reason the occupation numbers grow large while still in the linear regime (for small coupling, say), we are allowed to initialise the classical system in the quantum-vacuum like state n k = 1/2, and evolve the system using classical equations of motion throughout; at early times because the system is linear, at late times because the system has large occupation numbers. We only require that occupation numbers grow large before self-interactions become important.
To summarize: The ideal prescription to simulate a phenomenon arising from a quantum vacuum initial condition is to 1) start off with 1/2 in all modes, 2) evolve them all with the (quantum, but equivalently classical) linear equations until non-linear self-interactions become important, 3) discard all the modes that have not by then acquired large occupation numbers, and only 4) continue the now classical evolution of the highly occupied modes. Various levels of adherence to these rules can be argued for on a case-by-case level. Examples, where this applies include: • The primordial perturbations responsible for cosmological structure formation. These are the zero-point fluctuations of a weakly coupled scalar field, that grow because the accelerated expansion of space introduces non-adiabatic evolution of the modes [25]. This is one instance of the phenomenon known as "squeezing" of an initial vacuum state. Observations show that non-Gaussianities are minute, and so the entire earlytime evolution from vacuum fluctuations (1/2) to non-vacuum (n k 1) may be simulated using (almost linear) classical evolution.
• Resonant preheating after inflation arises when at the end of inflation, the oscillating inflaton mean field is in resonance with certain field modes (whether of another field or the inflaton itself). This resonance amplifies these modes from an initial vacuum state to large occupation numbers [15,16] (for a recent review, see [26]). Since the selfinteraction is usually quite small (λ 10 −12 or smaller for many inflation models), occupation numbers can grow very large before non-linearities become important. And so the CS approximation is valid all the way from the quantum vacuum initial state.
• Tachyonic preheating (or spinodal decomposition) occurs in hybrid inflation-type models, where a negative curvature of the potential V triggers an instability of certain modes k 2 + V < 0 [17]. These modes grow exponentially, until self-interactions become important. If the self-interaction is small, the classical evolution again holds from the initial quantum vacuum state (when the evolution is linear), and also in the subsequent non-linear regime, because occupation numbers are by then 1 (see for instance [22,23]).
• Certain plasma instabilites in gauge theories can also be described as unstable modes, at weak coupling [18,19]. As these acquire large occupation numbers, the CS approximation can be applied also in the context of the approach to thermal equilibrium in heavy-ion collisions.
A final point worth mentioning is that the classical regime with large occupation numbers does not imply that one particular classical realization (one member of the ensemble) is singled out. All observables must be computed as statistical expectation values over the whole classical ensemble of configurations, which is then expected to reproduce well the expectation values over the wave function (or density matrix) of the quantum system.

Classical simulations of vacuum decay
A quantum system at zero temperature in a local potential minimum (a "false" vacuum) may decay into a state in the global minimum (the "true" vacuum) through quantum mechanical tunneling. In the Euclidean formulation of quantum field theory the transition is described by an instanton [27], and from a path integral point of view, the transition is mediated by non-classical paths, paths that do not satisfy the classical equations of motion. The transition rate is straightforwardly computed in quantum mechanics, but is substantially harder to extract in quantum field theory.
In [1], an approximate agreement was reported between the instanton computation of the transition rate in 1+1 space-time dimensions, and the CS evolution of a vacuum ("half") initial state in the unstable vacuum. The authors were surprised and intrigued by their result, since tunneling is precisely the type of very quantum processes, where one would expect the CS approximation to fail. Indeed in quantum mechanics (field theory in 0+1 dimensions), the CS approximation does fail to reproduce the quantum tunneling rate [7].
Classical simulations of bubble nucleation are natural in the context of a finite-temperature phase transition, where the initial state is described by the finite temperature distribution of occupation numbers above the unstable vacuum. Then the transition is a classical effect whereby there is some (Boltzmann) probability that the ambient thermal fluctuations manage to spontaneously form a true-vacuum bubble, large enough to make it over the potential barrier and expand to eventually fill the whole of space (we will return to this point in more detail below). It follows that the finite-temperature bubble nucleation rate can in principle be computed by classically evolving all field configurations starting in the local potential minimum, and then averaging them over the initial Boltzmann distribution, schematically Γ Finite T = P Boltzmann,T [configuration] × transition rate of the configuration. (1.10) The result of [1] would suggest that the quantum tunneling rate follows from the same set of classical trajectories, but averaged over the quantum vacuum-like initial distribution [configuration] × transition rate of the configuration.
(1.11) This is a surprising result, and warrants further scrutiny. In particular, since classical evolution conserves energy, it would imply that quantum tunneling is simply the classical evolution of the subset of the initial condition ensemble, that have enough energy to nucleate a bubble.
In [9], the numerical computations of [1] were reproduced, although crucially it was pointed out that to get the reported agreement between numerical and instanton results, a "fudge" factor had to be introduced. The agreement occurs for 1/2 which amounts to rescaling the zero-point fluctuations from n k = 1 2 to 1 8 . The authors of [9] then carried out similar simulations of different initial conditions, and different models as well as an analysis related to cut-off dependence and renormalisation. The conclusion remained, that only when rescaling the amplitude of the initial conditions by an arbitrary factor < 1 is it possible to find the approximate agreement reported in [1].
We will expand further on that analysis, and show that the reported agreement is indeed a coincidence to do with the choice of the parameters of the model, the lattice cutoff and the fudge factor, and that it is not specific to the "half" initial condition. We will also generalise the simulations to 2+1 dimensions, and show that there is no agreement there. We will see that there are some essential differences between nucleation in 1+1 and higher dimensions.

Tunneling and Bubble Nucleation
Consider a potential V with two non-degenerate minima, with a barrier in-between. If the system is initially in the local minimum with highest energy, a transition may occur whereby the system moves to the global minimum with lowest energy ("false vacuum decay"). Energetically, it is very expensive for the field to move across the barrier in all of space simultaneously. Instead, one local region of space (a bubble) is created with the field in the global vacuum inside, in the local minimum outside, and with the field continuously interpolating between the two on the boundary (the wall).

Classical Bubble Nucleation
Classical bubble nucleation is the process by which random classical fluctuations (for instance in equilibrium at a temperature T ) by chance organise themselves into such a bubble. This happens all the time, but most bubbles are so small, that they collapse again. The energy criterion controlling the process is the balance between the energy cost of creating the bubble wall, interpolating between vacua, and the energy gain from the inside of the bubble having a lower potential energy than when the bubble is not there. In the simplest approximation one may write where σ is the surface tension, the energy associated with the interpolating field wall, and ∆V is the difference in potential at the two minima V global − V local (which is negative). In 1+1 dimensions, the volume is the distance between walls, 2R, while the surface is just a factor of 2 (2 walls), In order for a transition to happen, a random fluctuation has to occur that creates a pair of walls. Once these walls are established, there is no further energy cost in increasing the size of the bubble. The total energy is linearly decreasing with increasing R. We define the critical energy and the critical radius to be In 2+1 dimension, things are qualitatively different. Now which is maximised to give the saddle point solution In most cases, a random fluctuation does not acquire this critical radius, and the transition does not complete. The bubble shrinks again. But occasionally, a critical-size bubble is generated, which then continues to grow. In 3+1 dimensions, we have so that Throughout, we have assumed that the bubble is spherical, since this maximises the volume/area. There will be subleading contributions from many other near-spherical configurations.
In thermal equilibrium, the bubble nucleation rate is then proportional to the Boltzmann probability of a large enough random fluctuation Dividing by the volume V (not to be confused with the potential) and t normalises the rate to unit volume and time, respectively. A more detailed numerical analysis along the lines of (1.10) allows the direct computation of this quantity [28][29][30].
In a non-thermal environment, for instance a state with some non-thermal occupation numbers n k , the probability of creating such bubbles will depend on the state. As for the thermal equilibrium state, it may require that random multi-wavelength fluctuations manage to organise themselves into a large enough bubble configuration. But one could also imagine a state with only long-wavelength fluctuations (say, of size R crit ), in which critical-size bubbles are ubiquitous.
There is also the possibility that the state (whether thermal or not), simply has an energy density (much) larger than the height of the potential barrier. Then the system hardly notices, that the minima are separated, and will not need to minimise energy into a spherical bubble to perform the transition. In this case, transitions are common and fast. If the energy density is larger than |∆V |, one may also get transitions back again.
Finally, there is the possibility that the entire physical volume has too little energy to even make a single critical bubble. This is only a practical issue in a numerical simulation of a finite volume, and hence finite total energy. Then a transition will never happen, if the dynamics are classical and energy conserving.

Quantum Bubble Nucleation
Quantum tunneling is most apparent in situations where a barrier separates two local minima of the potential, and the energy of the state is smaller than the height of the barrier. In quantum mechanics (field theory in 0+1 dimensions), starting in one minimum, one may straightforwardly solve for the wavefunction of the system, giving a non-zero probability of finding the particle inside, and on the other side of the barrier. In time, there is an ever increasing probability for the particle to be measured in the other minimum. In the case when the second minimum is in fact the global minimum, we speak of vacuum decay.
In field theory, the analogous process can also be interpreted in terms of Euclidean instanton paths, famously in [27]. This instanton is a 4-D spherically symmetric saddle point of the Euclidean action. One may again write down To a good approximation, the rate of tunneling may then be written as but keeping in mind that this is the saddle point action rather than an energy, and that no temperature is involved.

The wall tension σ
Whereas ∆V is simply the difference between potential minima, computing the wall tension σ in the general case requires knowledge of the wall profile. For classical nucleation, an approximation is found by solving the (spherically symmetric) equation of motion for a static field profile interpolating between the two minima: For d = 1, this looks like time evolution in the potential −V , and is usually solved by numerical means (shooting) [31]. Then one may compute the wall tension as (2.12) In the limit when the wall is much thinner than the size of the bubble, the term (d − 1)/r may be neglected. Then it is not necessary to know the detailed shape of the wall, as one may rewrite (2.12) into which is easily computed, at least numerically. For the 4-dimensional instanton we must first rotate to Euclidean space, the saddle point equation in d + 1 dimensions becomes which in 4-dimensional spherical coordinates is equivalent to Eq. (2.11), in one dimension higher. Hence for a thin wall, the calculation of the wall tension proceeds in exactly the same way. This does not directly imply a relation between the tunneling rates, since E crit and S E are very different objects.

A convenient toy model potential
Following [1] we will focus on a specific potential, defined by the action It is parameterized by three quantities, λ, φ 0 and V 0 . For λ > 1 the periodic potential has global and local minima at φ = 2nπφ 0 and φ = (2n + 1)πφ 0 , respectively, with integer n. The potential is chosen to have V (φ local ) = 0 and ∆V = −2V 0 , and we define the masses The height of the potential barrier separating the two minima is given by (2.18) We will follow [1] and set λ = 1.2. The potential is therefore parametrized by m f and φ 0 . In this parametrization φ 0 fixes the location of the local vacuum but also influences the relative height of the potential barrier. We show in Fig. 1 the potential for example sets of parameters. We will compute the bubble nucleation rate primarily as a function of φ 0 , and from the potential alone, we expect the rate to decrease with increasing φ 0 .

Numerical implementation
We discretize the action on a space-time lattice, and solve the classical equation of motion, (2.20) A symplectic integrator scheme is used to ensure energy conservation for long simulation times. The lattice has periodic boundary conditions and the number of lattice sites per dimension and the spacing are denoted as N x and a, giving the linear lattice size L = N x a. We recast the lattice action in lattice units, whereby all dimensionfull quantities appear in dimensionless versions by means of the lattice spacing as am f , a d+1 V 0 , a 2 k 2 and so on. Consequently, the dispersion relation on the lattice is determined by the discretized Laplacian and given by where for each spatial dimension i, k i = n i 2π Nx for n i = −N x /2 + 1, ..., N x /2. The quantity am f then defines the lattice cut-off, since if the maximum momentum is aΛ π then the cut-off in physical units is Λ/m f = π am f . As am f decreases, the cut-off increases. We will in the following only explicitly write out powers of a when needed.
The quantum-like initial conditions are Gaussian distributed field fluctuations defined by These vacuum fluctuations are added to a homogeneous field placed initially at the local minimum φ(x) = πφ 0 . The "fudge factor" was introduced by [9] to parametrically fit tunneling rates to instanton results. = 1 is the physical value that mimics a quantum vacuum state, whereas other values have no obvious physical interpretation. As we will see, and consistent with [9], the apparent agreement between CS results and the instanton rate arises for 0.5. In preparation of the later discussion, it is instructive to generate a single initial condition φ(x), and simply compute the distribution of local field values. Fig. 2 shows a histogram superimposed on the potential. For a fudge factor of = 0.5, we see that the entire field configuration is inside the false vacuum initially. But for = 1, already at the initial time, the field is on the other side of the potential barrier in some small parts of space.
Following [1,9], as the simulation proceeds, we monitor the observable cos(φ/φ 0 ) , where . refers to the ensemble average, to define whether a configuration has transitioned to one of the neighboring global minima. For homogeneous configurations at the local/global minima this observable takes the value −1 or +1. A configuration is then said to have transitioned if where ∆ t=0 is the standard deviation of the same observable cos(φ/φ 0 ) at the initial time. Given an ensemble of N configurations, we define N surv (t) to be the number of these configurations that by a given time t have not yet transitioned. We then perform a fit to the form 24) where N 0 refers to the starting point of the fit. Then Γ is the bubble nucleation rate. Typically it takes some time before the configurations begin to transition. The fit was therefore performed from a time when N 0 was 60 percent of the total configurations in the simulation. An example is shown in Fig. 3.
-12 - We first consider the system in 1+1 dimensions, exactly as in [1,9]. We compute the nucleation rate for different values of φ 0 and for different values of the fudge factor, shown in Fig. 4. The number of configurations is N = 100, which was sufficient to get convincing results. In all our figures, the (bootstrap) statistical error bars are roughly the same size as the plotting symbols. The instanton estimate of the quantum tunneling rate is obtained via the expression [9] Where S B refers to bounce action computed with the tool CosmoTransitions [31]. Making allowance for possible small differences in fitting procedure and numerical implementation, this reproduces the results of [1] and [9], which may be summarized as follows: In 1+1 dimensions, CS simulations of a quantum vacuum-like initial ensemble produces a bubble nucleation rate of a similar order of magnitude as the quantum instanton result, at least for φ 0 ≤ 1.25. However, this agreement is achieved from a quantum-like initial condition not with occupation numbers of 1/2, but instead 1/8, hence a fudge factor of 1/2 [9]. In fact, tuning the fudge factor down from 1, one may achieve different levels of agreement with the instanton result at different values of φ 0 . The "half" initial condition (fudge factor 1) overestimates the quantum nucleation rate by several orders of magnitude, and has a weak dependence on the shape of the potential, φ 0 . In Fig. 5 we consider another choice of initial condition, namely the classical equilibrium mentioned above While the quantum vacuum has constant occupation number for all modes, in the classical equilibrium they are suppressed in the UV. The energy density is however still divergent as the cut-off am f goes to zero. We perform the same simulation procedure as previously, but now for different values of the parameter T . We see that just as we did for the fudge factor , we may also tune T to a semi-quantitative agreement with the instanton nucleation rate, in this case T = 0.1m f . Since the initial conditions correspond to a divergent energy density in the continuum, it is also prudent to test the robustness of our results to changing the cut-off, in our parametrization the quantity am f . The result of this procedure is shown in Fig. 6 for = 0.5. We see that while giving a weaker effect than varying , changing the cut-off is an alternative way of tuning the rate to match the instanton rate. As might be expected, smaller am f corresponding to larger cut-off and more energy in the system leads to a larger nucleation rate. We may consider introducing a mass counterterm (or more generally, renormalise the potential [32]) to counter the effect of the divergent initial condition. But because the zero point fluctuations do not stay put in classical simulations, this is difficult to achieve (see for instance [24]), and does not in itself solve the problem of a divergent energy being available for tunneling. The particular potential considered here is also not readily renormalisable. Figure 6: The dependence of the nucleation rate on the lattice cut-off in 1+1 dimensions, for quantum-= 0.5 initial conditions and comparing to the instanton rate.

Looking for bubbles and energy considerations in 1+1 dimensions
The nucleation rate for = 1 has a weak dependence on φ 0 , and similarly for = 0.5 for small φ 0 . We can begin to understand this at least qualitatively by considering the energy density of the configurations.
The top left panel of Fig. 7 shows the time-evolution of the observable cos(φ/φ 0 ) for individual configurations for = 1 at φ 0 = 1.2. The bottom left panel is one field configuration in space at different times, labelled by the value of cos(φ/φ 0 ) at that time. We see that all configurations transition through the threshold value −0.6 almost immediately, and that the field configurations have many nuclei and bubbles. This is an example of an initial condition with an energy density ρ larger than the potential barrier V max . There is no need for the configuration to randomly organise itself into a critical bubble for the transition to take place. In contrast, the right-hand panels of Fig. 7 show a simulation at = 0.5 and φ 0 = 1.2. Here, the transitions happen as an exponential decay. Also, field configurations evolve over time from a single, initially small, bubble (light blue) to a larger and larger bubble.
To make this more explicit, we compute the total energy and average energy density of the configurations. Fig. 8 shows the dependence of the average energy density on φ 0 and . The grey shaded region is where the average energy density is smaller than the potential barrier, while above, the energy density is larger than the barrier. Roughly speaking, one would expect the rate of nucleation to be exponentially suppressed only in the grey region, as a critical bubble needs to emerge through a stochastic process. And one would expect the rate to be unsuppressed everywhere else. Of course, individual field configurations are inhomogeneous, multiple nuclei complicate the picture, and some out-of-equilibrium initial states may have special properties enhancing nucleation. And so the boundaries of the grey region should be considered fuzzy.
We see that for = 1, we only enter the grey region far beyond the range of the figure. And that for = 0.5, we enter the region around φ 0 = 1.05, corresponding roughly to where the exponential dependence on φ 0 kicks in in Fig. 4.
Since the energy density is approximately ∝ 2 and the potential barrier ∝ φ 2 0 , the criterion for entering the grey region ρ = V max amounts to φ 0 ∝ . The proportionality constant in the case depicted here happens to be 2.05, and so the rate for the = 1 initial condition of physical relevance only becomes exponentially suppressed around φ 0 = 2.1, where the instanton rate is Γ L φ 2 0 V 0 = 4.7 × 10 −9 . In Fig. 9 we have extended the range of Fig. 4 to include the exponentially suppressed region for = 1. We see again that the CS approximation overestimates the nucleation rate by several orders of magnitude.
We can attempt to compute the wall tension in 1+1 dimensions from further developing the naive model of Sec. 2.1. We note that if the bubble is really in the global minimum inside the bubble (cos(φ/φ 0 ) = 1) and in the local minimum outside (cos(φ/φ 0 ) = −1), Figure 8: The dependence of energy density on φ 0 and . We also see a small dependence on the cut-off am f .
where 2R is the wall separation, and R hence the radius of the bubble. We now compute numerically the energy of bubbles, where we by hand force the interior and exterior to be We fit the parameters a and b for each critical bubble and relate to σ and ∆V via In this way, an estimate for σ can be obtained by extrapolating E Bubble to cos φ φ 0 = −1. Fig. 10 shows an example fit to a critical bubble obtained with φ 0 = 0.5, = 0.17. Fig. 11 shows the corresponding energy of the critical bubble, 2σ, for different values of φ 0 . This then is the minimal energy required for a configuration to classically transition to the global minimum. If the volume, cut-off, and is such that the energy is smaller than this value, the evolution of this interacting scalar field initially in an out-of-equilibrium state, will eventually drive the system to the classical equilibrium state in the local minimum.  Since the occupation numbers of the modes play an important role in the CS approximation, we also compute these through

Energy depletion and thermalization in 1+1 dimensions
In Fig. 12 we show the occupation numbers for a set of modes in time 8 . We see that initially, the occupation numbers (3.6) are indeed 2 /2, and as the nucleation is triggered (within a time of a few in mass units), they increase as potential energy is converted into excitations. The energy is mostly deposited in IR modes. As discussed in the preceding section, we can engineer an initial configuration with total energy less than the E crit which will never transition. An example of this is shown in Fig. 13. For very long time, we expect the particle numbers to slowly reorganise themselves into a classical thermal spectrum. Clearly, in 1+1 dimensions, this is an extremely long time, longer than we are able to simulate. This also implies that the nucleation rates that we have computed so far indeed arise from the quantum-like initial state. It is not such, that the initial state first thermalises to the equilibrium, after which the nucleation takes place. We will return to this point when considering 2+1 dimensional simulations.
4 Generalising to 2+1 dimensions. We now perform 2+1 dimensional simulations completely analogously to the 1+1 dimensional case. We discretize the 2+1 dimensional action on a quadratic lattice of size N 2 x . The scale is still set by the mass am f , and the other dimensionless combinations are now a 3 V 0 and a 1/2 φ 0 . The instanton prediction is obtained by generalizing the 1+1 dimensional case, and is given by where the bounce action S B is again obtained from CosmoTransitions [31].
In Fig. 14 we show the transition rates for different values of φ 0 for lattice simulations with fudge factors = 1, = 0.8, = 0.5, respectively, and again comparing to the instanton prediction. It is clear that the CS rate even when applying a fudge factor vastly overestimates the quantum tunneling rate. Quantum bubble nucleation and false vacuum decay cannot be modelled in this way. Simulations for large values of φ 0 and = 0.5 did not transition to a global minimum at all.   15 shows the nucleation rate for different values of the cut-off am f for = 0.5. As in 1+1 dimensions a higher cutoff (low am f ) results in higher rates, but the dependence is much stronger in 2+1 dimensions than it was in 1+1 dimensions. In 2 spatial dimensions, the number of UV modes grows much faster as the cut-off increases, and the initial energy density then also increases faster.

Looking for bubbles and energy considerations in 2+1 dimensions
We will again take a closer look at the configurations close to the transition. Fig. 16 shows two sets of simulations with high and low transition rates, respectively. The simulations shown in the left-hand panels have = 1.0 and φ 0 = 1.2 while the right-hand panels correspond to = 0.5 and φ 0 = 1.2.
We observe a clear difference in that on the left-hand side, transitions happen almost immediately, and the field configurations display many small bubbles nucleating close to each other. As in 1+1 dimensions, this is a case of the energy density being larger than the potential barrier. On the right-hand side, we have an exponential decay, with just a single bubble nucleating in the entire volume.
In a similar way as in 1+1 dimensions, we compute the initial energy density and compare it to the potential barrier. This is shown in Fig. 17 where again the grey area corresponds to parameter combinations where the energy density is smaller than the barrier, and a bubble must be created for the transition to take place.
We can again estimate the wall tension by computing the energy of 2-dimensional bubbles, but by hand fixing the inside and outside to the local and global minimum values. Fig. 18 shows the energy as a function of the radius of a growing bubble. As we argued in Sec. 2.1, the critical bubble in 2+1 and higher dimensions have a non-zero R crit , in contrast to the 1+1 dimensional case. From the quadratic fit in the figure we can tentatively estimate the critical bubble size to be Rm f = 15 − 20.

Energy depletion and thermalization in 2+1 dimensions
Finally, we will consider the evolution of the occupation numbers of the fields, also in 2+1 dimensions. Fig. 19 shows the occupation numbers as defined in  In 1+1 dimensions, we saw that the initial quantum-like distribution is essentially unchanged up until the transition takes place. But in 2+1 dimensions, even before the transition happens the dynamics have begun redistributing the energy to approach the thermal equilibrium state. For φ 0 = 1.2 this process does not have time to complete, but for φ 0 = 1.4, the transition rate is so small, that the system thermalises, reaching an asymptotic state. This would not happen in the true quantum system. It seems that in 1+1 dimensions, the time scales are such that classical nucleation is always much faster than thermalization. While in 2+1 dimensions, kinetic equilibration is often well underway by the time the transition happens. This ordering of time scales is dependent on the potential (the strength of self-interactions), the initial condition ( , say) and the cut-off (am f ).

Conclusions
Motivated by the intriguing possibility proposed in [1], that classical-statistical simulations could have something to say about quantum vacuum decay, we have investigated such simulations, both in 1+1 and 2+1 dimensions. The conclusion is disappointing, although perhaps not wholly unexpected.
As also demonstrated in [9], the reported approximate agreement between the instanton calculation and the CS simulations is there, but as we have seen it arises through arbitrarily adjusting the parameters of the initial conditions (the amplitude , cut-off am f ), and is also not specific to the quantum-like state with equal occupation numbers in all modes (thermal initial conditions work just as well, when tuning T ). In fact, the actual, = 1, "half" initial condition intended to be mimic the zero-point fluctuations of the false vacuum produces a nucleation rate several orders of magnitude larger than the instanton nucleation rate, also in the range of φ 0 , where the energy density is smaller than the barrier. In addition, obtaining even approximate agreement between CS simulations and the instanton result is specific to 1+1 dimensions. In 2+1 dimensions, the CS simulations consistently overestimate the nucleation rate by many orders of magnitude. We also attempted simulations along the same lines for the physically relevant case of 3+1 dimensions, but the nucleation rate is then far below our numerical reach, and advanced Monte-Carlo techniques are likely required to compute also the classical rate [28][29][30].
As mentioned in the introduction, the CS-approximation may be derived directly as a limit of the full real-time path integral [7]. It is only a good approximation for interacting quantum evolution for large occupation numbers, and even then only when computing "classical" observables. Quantum vacuum decay is both inherently quantum and by construction has an initial condition with occupation numbers 1. Such initial conditions can only reliably be simulated in the CS approximation for very small coupling, when the evolution equations are (approximately) linear. But we have seen that even the proposed "half" initial condition probes the non-linear regions of the potential considered here (Fig.  2).
We conclude that computing quantum tunneling rates in field theory beyond [27] remains a difficult task, which cannot be simulated using classical dynamics of an ensemble of configurations. It likely requires non-perturbative numerical methods at the level of the path integral, known to be challenging for real-time systems out of equilibrium (although see [6,24,33,34] and [7]). Fortunately, in almost all cases, phase transitions involve nonvacuum initial states, for which the quantum rate is insignificant compared to the classical nucleation rate. And classical nucleation rates may be computed using CS simulations or stochastic evolution in effective theories [28,29].