Real-time dynamics of false vacuum decay

We investigate false vacuum decay of a relativistic scalar field initialized in the metastable minimum of an asymmetric double-well potential. The transition to the true ground state is a well-defined initial-value problem in real time, which can be formulated in nonequilibrium quantum field theory on a closed time path. We employ the non-perturbative framework of the two-particle irreducible (2PI) quantum effective action at next-to-leading order in a large-N expansion. We also compare to classical-statistical field theory simulations on a lattice in the high-temperature regime. By this, we demonstrate that the real-time decay rates are comparable to those obtained from the conventional Euclidean (bounce) approach. In general, we find that the decay rates are time dependent. For a more comprehensive description of the dynamics, we extract a time-dependent effective potential, which becomes convex during the nonequilibrium transition process. By solving the quantum evolution equations for the one- and two-point correlation functions for vacuum initial conditions, we demonstrate that quantum corrections can lead to transitions that are not captured by classical-statistical approximations.

In general, the decay of a metastable vacuum consists of a first-order phase transition where the order parameter changes from a metastable phase to a stable phase.The transition may occur through nucleation of bubbles caused by quantum or statistical fluctuations.Bubbles below a critical size collapse due to the overwhelming surface energy cost, while larger bubbles grow rapidly and eventually fill space to complete the phase transition.
The conventional approach to computing the bubble nucleation rate in field theory involves solving the equations of motion in Euclidean signature to find the saddle point of the path integral (the bounce).While this semiclassical method is expected to work well for strongly first-order phase transitions, when the jump in the order-parameter field at the phase transition is sufficiently large, it has limitations for weaker transitions and breaks down for strongly coupled quantum systems (see, e.g., [21][22][23][24]).Furthermore, this method is limited to describing bubble nucleation and does not capture other processes like spinodal decomposition and resonant nucleation, which can also cause false vacuum decay [25,26].
This work aims to go beyond semiclassical methods for the computation of the decay rate.In principle, false vacuum decay is a well-defined initial-value problem in real time, which can be formulated in nonequilibrium quantum field theory on a closed time path [27].The problem involves the computation of the time dependence of the order-parameter field from its initial value in the metastable phase to its final value in the stable phase.
Here, we consider a relativistic scalar quantum field theory with quartic self-interaction in 3+1 space-time dimensions.We compute the time-dependence of the field using the two-particle irreducible (2PI) effective action at next-to-leading order (NLO) based on a nonperturbative large-N expansion [28][29][30].In the high-temperature regime, we also compute the real-time dynamics from an approximation based on classical-statistical field theory simulations on a lattice [31][32][33].
This allows us to establish a "dictionary" between semiclassical characterizations of false vacuum decay and nonequilibrium field theoretical descriptions in terms of time-dependent order-parameter fields and correlation functions.Importantly, to formulate a quantitative comparison of the decay rate, we properly take into account fluctuations using an effective potential, building upon prior work in this direction [20,[34][35][36][37].
This work aims to identify relevant and measurable information on the decay that can be extracted and constructed from correlation functions, which is particularly practical for situations in which semiclassical concepts, such as the identification and counting of bubbles, become ambiguous.This is especially relevant given the increased experimental possibilities of tabletop experiments.These systems offer a very high degree of controllability and a genuinely quantum-mechanical nature, making them unique platforms for analogue simulations of quantum field dynamics [38,39], including the decay of a false vacuum, see, e.g., experimental proposals for false vacuum decay using ultracold atoms [40][41][42][43][44][45][46], quantum spin chains [47,48], and a recent experiment using a ferromagnetic superfluid [49].
This paper is organized as follows.After introducing the false vacuum decay problem in Sec.II, in Sec.III, we summarize relevant aspects of the Euclidean framework to describe phase transitions.Sec.IV considers the dynamics for a high-temperature situation, where classical-statistical fluctuations are expected to dominate over quantum fluctuations.In Sec.V, we establish a link between Euclidean and real-time decay rates in a system at high temperatures.Section VI addresses the problem of real-time dynamics in the quantum regime at low temperatures, where we discuss the importance of genuine quantum corrections to the real-time dynamics using the 2PI effective action approach.Finally, we conclude and give an outlook in Sec.VII.Several appendices provide details on quantum field theory out of equilibrium and the derivation of the time evolution equations.

II. THE PROBLEM OF FALSE VACUUM DECAY IN FIELD THEORY
The system we consider is a (3+1)-dimensional scalar field theory with classical action 1 for a scalar field φ(x 0 , x) depending on time x 0 ≡ t and spatial coordinates x, whose potential is given by The parameters m 2 , λ, h are real and positive.The potential with h = 0 exhibits two degenerate minima at constant values of φ given by φ ± (h = 0) = ± 6m 2 /λ separated by a potential barrier.We consider h ̸ = 0 to lift this degeneracy, such that a double-well potential arises with two minima whose free-energy difference is nonzero, The local minimum φ + denotes the "false vacuum" in a first-order phase transition to the global minimum, or "true vacuum," φ − .False vacuum decay is an initial-value problem, starting from a nonequilibrium initial condition (the metastable state).If one only considers classical dynamics, the time evolution of the field is governed by the classical equation of motion, i.e., 1 Throughout this work we use units where FIG. 1. Sketch of the classical potential V (φ) in (2).
if spatial homogeneity is assumed.If the initial field value is close to the local minimum φ + , then the field may not reach the global minimum without additional kinetic energy to overcome the potential barrier.In quantum field theory, however, the corresponding state is always expected to decay towards the state with the lowest available free energy.In general, the description of the transition requires the consideration of quantum and/or statistical fluctuations.Therefore, both the initial value φ + , as well as the dynamical evolution Eq. ( 3), and thus the final field value φ − , will receive corrections due to quantum-statistical fluctuations.At very low temperatures or vacuum initial conditions, the transition is dominated by quantum fluctuations, where tunneling plays an essential role.
In quantum-statistical field theory, the initial pure or mixed state of the system can be described by a density operator ϱ(t 0 ) at a given time t 0 .The goal is to understand the temporal evolution of the system from its initial state at t 0 to its subsequent states at later times.The expectation value of any observable O(φ) and its dynamical evolution can be derived from the functional integral [50] whose definition involves the Schwinger-Keldysh closed time path [27], see Appendix A for more details.For instance, for O(φ) = φ the path integral (4) would encode the evolution equation for the field expectation value ϕ(t) = ⟨φ⟩(t) which, for a spatially homogeneous system, now involves the quantum-statistical fluctuations in the initial state and the corrections to the classical evolution Eq. ( 3).The time-dependent expectation value ϕ(t) evolves during the transition from its initial value ϕ + = ϕ(t 0 ) to its final value ϕ − = ϕ(t) at sufficiently later times t > t 0 .In general, the functional integral ( 4) is too complicated to be solved numerically without further approximations.Because of the rapidly oscillating integrand ∼e iS , standard numerical techniques based on impor-tance sampling cannot be applied, and one has to resort to suitable approximate descriptions, which we turn to in the following sections.

III. CONVENTIONAL EUCLIDEAN APPROACH TO FALSE VACUUM DECAY
This section summarizes the most important calculation steps within the conventional Euclidean approach to false vacuum decay.This is to prepare for the comparison with the real-time analysis in the following sections.The central observable here is the decay rate per unit volume.The latter consists of the probability of nucleation of bubbles of the new phase per unit volume per unit time.

A. Bounce
Let us consider a theory defined by the classical action (1).The Euclidean action is obtained by analytically continuing time to imaginary values by performing the Wick rotation τ = it.It is written as S E = −iS and is a functional of the Euclidean field φ E (τ, x): (5) In the Euclidean action, the sign preceding the potential term is reversed compared to the Minkowskian action.Due to this "inverted potential," the Euclidean action allows for nontrivial classical solutions connecting the two minima.One of these solutions is crucial in describing false vacuum decay, as described in the following.
To compute the decay rate per unit volume, one can expand the Euclidean action around the saddle-point solution φ B (τ, x), which satisfies the classical Euclidean field equation of motion together with the boundary conditions where β = 1/T is the inverse temperature.
The field configuration φ B , often called bounce or instanton, is associated with the critical bubble within the metastable phase, the interior of which consists of the stable phase.We rewrite the stationarity condition (7) by assuming four-dimensional rotational symmetry of the solution in the limit of vanishing temperature.Thus, the bounce equation in terms of the hyperradial coordinate ρ 2 = x 2 + y 2 + z 2 + τ 2 simplifies to the one-dimensional differential equation where prime denotes d/dφ E .This equation describes the motion of a particle in the inverted potential and subject to a "friction force" 3 ρ dφ E dρ .On the other hand, in the high-temperature limit, the leading contribution may be obtained from the dimensionally reduced theory by performing the integration over the Euclidean time coordinate, i.e., Here, S 3 [φ E ] denotes the rescaled action of a threedimensional field theory for a τ -independent field.In this case, the critical bubble configuration has a threedimensional rotational symmetry.It can be parametrized by φ B (r), with r 2 = x 2 + y 2 + z 2 , the solution of the bounce equation similar to the previous case: In both the high-T and low-T cases, the bounce solutions satisfy the boundary conditions and for the respective radial coordinate.The corresponding explicit expression for the high-T Euclidean action in polar coordinates is given by

B. Nucleation rate
The knowledge of the bounce action allows one to estimate the decay rate.We now consider fluctuations around the non-trivial saddle point, described by ( 8) and (10), respectively.The actual calculation of the one-loop fluctuation determinant involves one negative eigenmode, several zero modes and infinitely many positive modes.Restricting to the limits of low-and high-temperature cases, the nucleation rate per unit volume for the former case can be written as [15] The prime on the determinant indicates that the zero eigenvalues, i.e., the translational modes, have been integrated out (but the negative eigenmode is kept) [20].In this case, due to the symmetry of the solution, there are four zero modes, each bringing a factor of ( S E 2π ) 1 2 and resulting in the factor ( S E 2π ) 2 in (14).Note that the decay rate per unit volume is exponentially suppressed by S E .
In the high-T case, the decay rate per unit volume is expected to be given by [18] − indicates the single negative eigenvalue.The different exponent 3/2 in the prefactor compared to the low-T case is due to the three translational modes in the high-temperature regime.
In typical cases, the action appearing in the exponents in (14) or (15) is significant, such that the exponential is very small.Therefore, the order of magnitude of Γ/V 3 is determined predominantly by the instanton action and the temperature.In this work, we concentrate on this exponential part and disregard the prefactor, which may be expected to be of the order of unity.
In order to make practical use of the above formulae, it is necessary to calculate the bounce.However, relying solely on a saddle point approximation around the bounce obtained from the classical potential V may not provide an accurate description.The impact of finite temperature can be significant because thermal fluctuations modify the potential, leading to a finitetemperature effective potential.In other cases, the classical action may not exhibit two distinct minima, but quantum corrections introduce an additional local minimum in the effective potential [51,52].Therefore, incorporating the quantum-statistical fluctuations into a modified effective potential can offer a more efficient description [20, 35-37, 53, 54].
For the underlying real-time problem of false vacuum decay, due to its nonequilibrium nature, it can also be necessary to consider the effective potential as time dependent.As a consequence, the bounce and the decay rate per unit volume can, in principle, vary with time.Therefore, a single decay rate may not be enough to fully capture the nonequilibrium process's complexity.We consider in the next section the nonequilibrium realtime problem in a high-temperature scenario to address these questions.

IV. REAL-TIME FALSE VACUUM DECAY IN THE CLASSICAL-STATISTICAL REGIME
In this section, we analyze the evolution in time of a field theory characterized by the classical action (1).Initially, the field is trapped in the false minimum of the potential (2).We assume that the system is prepared in an initial state at a high temperature, represented by the density operator ϱ(t 0 ).This choice allows us to approximate the dynamics through classical-statistical field theory evolution [55].
The idea behind the classical-statistical field description is that the field is evolved according to the classical real-time equations of motion.The quantum expectation values are then obtained by averaging over the statistical ensemble with the initial correlation functions determined by the quantum theory.

A. Dynamics of the field
We introduce dimensionless quantities such that the parameters in the potential can be reduced to only one dimensionless parameter by rescaling as The dimensionless parameter h indicates how much the potential is tilted.We prepare the system's initial state with a spatially homogeneous field average close to the potential's local minimum.More precisely, we set for the field at initial time ϕ(t = 0) = ϕ 0 , φ(t = 0) = 0.The initial fluctuations are encoded via the initial two-point correlation function.These fluctuations change the effective potential of ( 2) and, as a result, the field is typically displaced from the local minimum of the effective potential at initial time.
The subsequent time evolution of the system can be divided into several phases.It typically takes some time for the transition to the global minimum of the effective potential to begin.During the early-time dynamics of the field average or one-point function, the system oscillates around the local minimum of the effective potential.As the energy is converted into fluctuations with nonzero momentum, these oscillations gradually decrease in amplitude.During this period, a resonating momentum band gets excited and at later times, the oscillations decay.In the second, much longer phase, the system evolves towards a prethermal state with a field value close to the local minimum.Prethermal means that the system has many equilibrium properties but is not yet thermalized [56], as the field is still trapped in the false minimum.In particular, the field obeys the equipartition of energy between momentum modes, as expected for a classical system.This prethermalization has been verified in our simulations, where we observed that the average kinetic energy of each mode is T eff /2.This allows us to define an effective temperature T eff of the system as where the brackets indicate the ensemble average and the average over volume is implied.In this second phase, the system effectively loses the details of the initial conditions.At the same time, the dynamics depends on only a few parameters: the couplings of the effective potential, which differ from the bare potential due to fluctuations, and the total energy density.Consequently, the specific initialization of the fluctuations in the system becomes less important at later times.
For h ̸ = 0, the field eventually transitions from the local to the global minimum.The timescale at which the transition takes place depends on the value of h.In particular, for large values of h, the transition can occur before establishing prethermalization.In contrast, for small values of h, the transition occurs when the system is already nearly thermal.The strength of the fluctuations or, equivalently, the effective temperature is another parameter determining the transition time.
Numerical results are obtained from simulations on a cubic lattice with periodic boundary conditions.The number of lattice sites per spatial dimension and the spacing are denoted as N x and a x , respectively, such that the linear lattice size is L = N x a x .We set a t = 0.01a x , N x = 64, h = 0.055.The initial conditions for the field were chosen as where ϕ 0 ≡ ⟨φ(0, x)⟩ = 1.99.We choose the initial distribution as f p (0) = 2.225 ω p (0), and ω p (0) = p 2 + m 2 init is the initial frequency and c p satisfies [55].For these values, the time evolution of ϕ(t), the equaltime statistical connected two-point function, defined as 2 , and the effective temperature T eff (t) is shown in Fig. 2. The black solid lines show the correlation functions obtained from averages over n = 200 runs.The gray dashed lines indicate a subset of single realization runs, which are not averaged, for comparison.One observes a short phase of decaying oscillations of ϕ(t) at very early times.Next, in the more prolonged phase before the transition, the field value decreases almost linearly, together with the temperature.We emphasize that this stage of the dynamics depends on the choice of the initial conditions.The white noise initial spectrum, considered here, results in energy flow in momentum space towards the infrared region.This leads to the observed slow growth of F (t, t), which, in turn, decreases the value of ϕ(t).Simultaneously, the temperature T eff time slightly decreases.Eventually, the transition to the true minimum begins, resulting in the faster decrease of ϕ(t).
The decay of the false vacuum can proceed via bubble nucleation in position space.The formation of bubbles can be observed in the evolution of individual classicalstatistical simulation runs.To illustrate this, snapshots of the evolution for a two-dimensional section of the classical field are shown in Fig. 3 using a single run without averaging.At the beginning (top left), the background is filled with the metastable phase (lighter shade) in the presence of fluctuations (darker shade).Finally, in the last portion of the trajectory, the field oscillates around the true ground state and starts thermalizing (bottom right).In the snapshots of the real-space configurations, a bubble of the new phase exceeding the critical size appears (top right) and starts expanding (bottom left).In the following subsection, we focus on how to extract the decay rate using the notion of an effective potential.

B. Time-dependent effective potential
To calculate the exponential factor that governs the decay rate per unit volume Γ/V 3 in (15), one has to compute the effective temperature for the system T eff and the shape of the bounce φ B , which in turn depends on the (effective parameters of the) potential V (φ).As described in the following, we consider a method for ex-tracting the effective potential directly from the lattice simulations.In a regime with slow field evolution, such as realized in the prethermal regime described above, kinetic energy contributions should be small, i.e., φ ≪ mϕ and φ ≪ m 2 ϕ.We may extract an effective free energy in terms of the macroscopic field ϕ.
More precisely, up to the time when the first bubble is nucleated, we consider a free energy or effective potential V eff (ϕ) of the same form as the bare one, but with effective couplings: It involves the effective mass m eff and effective coupling λ eff .We also allow for a constant source term ∼J linear in the field to probe the effective potential at different values of ϕ.Specifically, for each external source J, we denote the value the field settles at a given time t in ϕ = ϕ(J) = ⟨φ(J)⟩.This field value corresponds to the stationary point of the effective potential, where its derivative with respect to ϕ vanishes, i.e., ∂V eff (ϕ)/∂ϕ = 0. Inversion reveals a one-to-one correspondence between each value of ϕ and the corresponding linear term, denoted as J = J(ϕ), explicitly given by Our analysis for extracting the effective potential is similar to the one in [34], where the time-evolution projection of the microscopic equation of motion onto the zero momentum sector is considered.In particular, they consider the volume-averaged field value in a single run.As a consequence, the field is subject to a stochastic force in addition to the deterministic one.Our analysis, instead, is based on tracking the ensemble average field rather than a single run averaged over the volume.As a result, the stochastic component vanishes, leaving only the deterministic parts.Also, the second time derivative of the field becomes significantly smaller because the oscillations average out within the ensemble.The linear term as a function of the field expectation value is illustrated in Fig. 4 using circles for various times, before and after the transition.
Using (21), we can thus estimate the effective couplings m 2 eff and λ eff from our numerical simulations.The effective parameters are determined by minimizing the sum of the squares of the residuals.The solid lines depicted in Fig. 4 represent the fitted values of the derivative of the effective potential.We observe that the effective couplings exhibit a slow evolution over time.
While the extracted data points in Fig. 4 appear to be well described by (20) in the outer convex regions of the potential, this is not the case in the inner nonconvex part of the polynomial ansatz.In the inner region, indicated by the dashed lines in Fig. 4, the convexity of the effective potential implies a nonanalytic form beyond (20).Accordingly, the data points extracted from the simulations deviate from the fit in the inner region.The FIG. 4. Linear source J as a function of the field expectation value ϕ at various times t.The solid lines represent the fitted results obtained in the prethermal state, i.e., before decay occurs, using (21).The data points in the inner region of the potential deviate from the fit, as indicated by the dashed lines.This observation reveals a convex shape for the inner part of the potential.
coexistence of two hysteresis phases suggests the presence of this inner region and is typically associated with the notion of "Maxwell's construction."

V. COMPARISON OF THE REAL-TIME DECAY WITH THE EUCLIDEAN APPROACH
This section compares the real-time evolution of false vacuum decay explored in the previous section and the Euclidean approach discussed in Sec.III.We focus on the decay rate per unit volume, the fundamental quantity for our analysis.In the Euclidean approach, the rate per unit volume for the specific high-temperature scenario can be approximated as Γ/V 3 ∼ e −S3/T eff , where S 3 represents the three-dimensional bounce action, and T eff is the effective temperature.
Before analyzing the numerical data, it is essential to understand how the decay rate is affected by the timescales of the system: the prethermalization timescale (the timescale for the prethermalization processes in the metastable phase to complete) and the transition timescale (the timescale for the nucleation of a bubble to occur).Our central question focuses on the possibility of capturing the true real-time dynamics through a single rate and how to extract it from the real-time data.In the case of prethermalization times significantly longer than the transition timescale, the decay rate will be independent of time.However, for practical, computational or experimental observations, it is essential to choose parameter values that allow us to observe the transition at a reasonable time.As a result of this choice, the two timescales may not be so different, leading to some time dependence of the decay rate.

A. Extracting the decay rate from the effective potential
To extract the relevant timescales for computing the rates, we identify the onset of the decay process at approximately t 0 ≈ 800.This is evident from the data points deviating from the polynomial fit, as shown in Fig. 4. Our primary focus is calculating the rate per unit volume around this time t 0 , including a small time interval before and after.Note that at much later times the field reaches the true minimum, and the decay process is already complete, rendering the concept of a decay rate irrelevant.Conversely, the decay rate is expected to be extremely small for much earlier times before the decay starts.We will investigate this point in the following subsections.
Following Sec.III, we determine the profiles of the critical bubble φ B (r) for various times.Rarely can (10) be solved analytically.Thus, we numerically compute the associated bounce using the so-called shooting method [15,57].The solutions are shown in Fig. 5, where a clear time dependence of the bounce profile is observed.For each time, we calculate the corresponding bounce action S 3 (t) by inserting the profile into (13), and the timedependent effective temperature T eff (t) is determined using (18).The resulting ratio S 3 (t)/T eff (t) is displayed in the inset of Fig. 5. Notably, this ratio decreases with time and only converges to a time-independent value at significantly later times.
The corresponding decay rate, defined as (Γ/V 3 ) 1/4 , has the units of inverse time.The result is illustrated in Fig. 8 as the gray solid line.This decay rate is depicted for the aforementioned time range and will be used for later comparison.We note again that this result does not consider any corrections to the rate that come from the preexponential factor.

B. Extracting decay rates from real-time dynamics
We have developed two methods to compute the decay rate from our real-time simulations.

Survival probability
The first method we employ to find the decay rate involves monitoring the formation of bubbles as they nucleate in specific regions of space.To identify the nucleation of bubbles, we conduct a series of classicalstatistical simulations on cubic lattices of various volumes V 3 = 128 3 , 256 3 , 512 3 .We maintain identical parameters in each simulation and initiate the system in the metastable state.We then observe the elapsed time until the system grows a bubble of the stable phase, denoting it the nucleation time t.To ensure the statistical significance of our results, we perform a total of n = 125 simulations for each volume.To resolve a critical bubble, we consider a small lattice cube and compute the average field value inside it for each point on the lattice.The size of the cubes is chosen according to the expected radius of the bubbles, estimated to be around r ∼ 10 (see Fig. 5).The primary purpose is to prevent the doublecounting of the same bubble during the nucleation process.While the specific length of the cubic boxes is not critical, it is chosen to ensure that the results are independent of this scale.The system is considered to have decayed when one (or more) of these cubic boxes exhibit an average field value greater than the value ϕ max ≈ 0.
We create a time distribution based on collected data of bubble nucleation events and then calculate the cumulative distribution, also known as survival probability.The survival probability is defined as p(t) ≡ n surv (t)/n, where n surv (t) represents the number of simulations that have not nucleated an expanding bubble up to a given time t.
By utilizing the cumulative distribution instead of the decay time distribution, as done in [58], we gain an advantage in studying the distribution's tails with fewer realizations.This approach involves summing up the decay times until a given time, resulting in a smoother curve less influenced by fluctuations in individual decay times.In Fig. 6, we present the rescaled survival probability p(t) V /V3 on a log-linear scale as a function of time for various volumes V 3 .The rescaling of the distributions is obtained using a reference volume V = 128 3 .The inset of the figure displays the same functions before the rescaling.Notably, noticeable variations can observed across different volumes.In larger volumes, bubble nucleation occurs earlier, and time distributions are narrower.This observation aligns with the intuition that larger systems are more likely to nucleate bubbles sooner.In the hypothetical scenario of an infinitely large lattice volume, the distribution would resemble a step function akin to a theta function.
Importantly, Fig. 6 demonstrates that, after the rescaling, all curves lie approximately on top of each other.This is consistent with the prediction that the decay rate per unit volume must be volume independent.
Using a range of grid sizes, we could observe a broader range of times compared to what would have been possible with a single grid.This combined distribution allows for a collective analysis.Since a single decay rate is not feasible for the considered choice of the parameters, we instead consider a time-dependent Γ(t) and extract it using The outcome of this analysis, i.e., the corresponding decay rate (Γ(t)/V 3 ) 1/4 , is presented in Fig. 8, where the decay rates for V 3 = 128 3 , 256 3 , 512 3 are shown in the same colors as in Fig. 6.The standard deviation estimates the error.

Field expectation value
The second method for obtaining the decay rate involves studying the evolution of the one-point function during the transition from one minimum to the other.Fig. 7 illustrates the behavior of the one-point function ϕ(t) for various volumes V 3 .It displays the expected volume dependence for small volumes and volume independence for larger volumes, specifically when V 3 ≥ 400 3 .The gray region indicates the transition time range for the large volume limit.
We define the time-dependent ratio of expectation val-500 750 1000 1250 1500 1750 Time: t 7. One-point function ϕ(t) as a function of time t for different lattice volumes V3, ranging from 64 3 to 600 3 .In the large volume limit, the one-point function becomes volume independent.
ues as In this context, the one-point function ϕ and the minima ϕ + and ϕ − are all time-dependent quantities.This fact arises from the system settling into the local minimum of the potential ϕ + after an initial transient early time regime.Over time, the potential itself evolves, gradually shifting the minimum ϕ + towards smaller values, as observed in the previous section.Likewise, the global minimum ϕ − changes over time as all fluctuations must settle and thermalize after the transition.To determine the values of ϕ + and ϕ − , we employ polynomial fitting of the one-point function.We fit over a selected time range during which the field settles in the respective minimum of the potential, forming nearly a flat line in the field.By extrapolating these fits to the entire time range after the early time regime, not relevant to our analysis, we accurately determine the values of ϕ + and ϕ − .This approach ensures that p ϕ (t) takes the value of 1 in the false vacuum and 0 in the true vacuum, as expected for a probability function.The probability p ϕ (t) can be expressed as (for a detailed derivation, see [20,59]) where, in the most general form, the function I(t) is given by In this expression, t 0 indicates the starting time of the transition.The rate is multiplied by the volume of the bubble of radius r(t, t ′ ) formed at time t ′ .The radius at time t of a bubble nucleated at time t ′ with velocity v w is 700 800 900 1000 Time: FIG. 8. Decay rate (Γ(t)/V3) 1/4 (from the survival probability), (Γ ϕ (t)/V3) 1/4 (from the evolution of the one-point function) and (Γ(t)/V3) 1/4 (from the bounce) as a function of time.
approximately given by r(t, t ′ ) ≈ v w (t−t ′ ), neglecting the bubble size at nucleation.Note that the previous analysis focuses exclusively on bubble nucleation and expansion while ignoring bubble collisions.
Finally, from (25), one can compute the decay rate per unit volume as where I (4) is the fourth derivative of the function I(t), and we set v w = 1 for simplicity.Remarkably, this result is independent of the particular choice for t 0 .The resulting curve for the corresponding decay rate (Γ ϕ (t)/V 3 ) 1/4 is shown in Fig. 8 as the black solid curve.We display the decay rate over a period of time starting when the rate starts to be non-negligible and ending when about half of the volume shifts to the new phase, around p ϕ ≈ 0.5.
Before we conclude, we can formulate simple estimates to understand the order of magnitude of the rate.At this stage, we neglect the time dependence of Γ ϕ for simplicity.Let us consider the transition time denoted by ∆t ≡ t f − t 0 , i.e., the time required to transition from the start (at t 0 ) to completion (at t f ).In that case, the general expression (25) simplifies to and therefore, for t = t f , I(t f ) ≈ 1, leading to which gives us the transition time estimate in terms of the rate per unit volume as In the case that the decay proceeds by a single bubble nucleating and expanding, I simplifies to which gives Therefore, the inverse timescale of the transition, i.e., the rate Γ ϕ = ∆t −1 , can be estimated as Γ ϕ ≈ p ′ ϕ (t)/p ϕ (t).

C. Comparison
Now that we have all the necessary elements, we can compare real-time decay rates and those computed using the instanton method.
Fig. 8 summarizes the decay rates obtained using the three methods.
Using the survival probability method, we can combine the results obtained using different volumes.Remarkably, within the overlapping regimes, the corresponding rate (Γ(t)/V 3 ) 1/4 is volume independent and nicely forms a single line.
Furthermore, the figure exhibits an intriguing trend, common to all methods, showing that the decay rates increase exponentially in time.Both real-time methods, the survival probability one and the one-point function method, show good quantitative agreement.In contrast, the instanton method exhibits a slight shift from the other two curves, possibly attributed to the neglected pre exponential factors.Including this prefactor was beyond the scope of our analysis.Nevertheless, the prefactor is expected to be less sensitive to the time-dependent shape of the potential compared to the exponential term.This is indeed reflected in the fact that both the instanton and the real-time methods predict a similar time evolution of the decay rate, and the quantitative agreement of the decay rate is quite satisfactory.
We verified that when using the thin wall approximation of the bounce, the agreement between the two methods becomes worse.This is reflected in both the magnitude of the decay rate as well as its time dependence.

VI. TOWARDS THE QUANTUM REGIME OF REAL-TIME FALSE VACUUM DECAY
In the previous section, we focused on investigating a situation where the decay of the false vacuum takes place at high temperatures.This scenario allows for the application of classical-statistical simulations to describe the dynamics.In this context, the decay was observed to occur through bubble nucleation, and we analyzed realtime observables.
In this section, our aim is to formulate the problem of false vacuum decay directly in quantum field theory.
The description is based on the nonequilibrium 2PI effective action on a closed time path.It describes the time evolution in terms of dynamical equations for one-and two-point correlation functions of the field.

A. Quantum effective action
The 2PI effective action, denoted by Γ[ϕ, G], is formulated in terms of a one-point function, ϕ(x), and a two-point function, G(x, y).It can be expressed as follows [55]: where S C [ϕ] is the classical action integrated over the Schwinger-Keldysh contour, G −1 0 (ϕ) is the inverse classical propagator, and Γ 2 [ϕ, G] is the sum of all 2PI graphs with full propagator lines.
The quantum equations of motion can be derived by taking the functional derivative of the effective action Γ[ϕ, G] with respect to the macroscopic field ϕ(x) and the propagators G(x, y).The stationary conditions are then obtained by setting these functional derivatives to zero, resulting in a system of coupled equations.
The equations can be expressed by decomposing the two-point function into its statistical component F (x, y) and spectral component ρ(x, y) as2 [60] and are given by field and the fluctuations and is given by The self energies Σ ρ , Σ F , and Σ (0) are defined as follows: where the dependence of the self energies on ϕ and G is implied.The right-hand sides of the evolution equations involve memory integrals, which take into account the whole time history of the evolution.
To solve the set of coupled equations ( 34) and ( 35), the self-energies must be truncated.A powerful nonperturbative truncation scheme is provided by the 2PI large-N expansion [28,29].Already the NLO of this expansion is known to be able to describe various phenomena observed for nonequilibrium scalar fields, such as instabilities [61,62] and thermalization [63,64].Even though the convergence of this expansion relies on a large value of N , the method has been successfully applied to models with a moderate number of field components and even for single-component fields.
In the following subsection, we present some numerical results for the dynamics of the false vacuum decay within the NLO large-N expansion, focusing on a single field component.The explicit form of the equations of motion is presented in Appendix B, along with the corresponding expressions for the self-energies within the chosen approximation scheme.To fully determine the system's time evolution described by the equations of motion, the initial conditions must be specified, as outlined Appendix B [(B11) and (B12)].

B. 2PI NLO evolution for vacuum initial conditions
In this section, we consider general nonequilibrium initial conditions where-unlike before-the field does not reach a prethermalized state before the transition.Therefore, it is not possible to estimate the decay rate per unit volume using the same assumptions as in the previous analysis.
In such circumstances, we focus on the dynamics of field correlation functions.By studying these correlation functions, we can track the decay of the field from the false minimum without referring to notions such as bubbles, which are well defined only in the semiclassical limit.
A powerful advantage of the 2PI effective action approach is that it allows a clear comparison between the quantum dynamics and the classical-statistical approximation.In order to detect both qualitative and quantitative indications of quantum effects in the dynamics, we run the simulations twice with identical initial conditions within the NLO large-N approximation.First, by simulating the entire quantum dynamics and, second, only the statistical-classical limit of the dynamics (for further explanation, see Appendix C).The results are shown in Fig. 9, where the dynamics of the one-and twopoint functions is plotted.The full quantum dynamics (classical-statistical limit of the dynamics) is indicated by dashed lines (solid lines).As can be observed, the field's decay in the full quantum dynamics occurs faster than in the corresponding classical-statistical limit.This observation suggests that the quantum part in the dynamics, although challenging to isolate due to the complexity of the equations, speeds up the transition.
Next, we explore a variation in the initial displacement of the one-point function from the false minimum.We employ m 2 = −4, λ = 10, h = 0.01, and V 3 = (0.5 × 64) 3 .We keep these parameters fixed and vary the initial value of the macroscopic field ϕ 0 .The numerical results are depicted in Fig. 10: the full quantum evolution (the classical-statistical approximation) is indicated by dashed lines (solid lines).We start with an initial condition for the field such that its amplitude generates enough fluctuations that lead to the decay of the field.In this case, the field crosses the barrier almost immediately for very large initial field oscillation amplitudes.Then, we progressively reduce the field's initial amplitude to have fewer fluctuations generated by the oscillations.The field experiences parametric resonance followed by damped oscillations and remains in the metastable state for a longer period of time.
In the case of the full NLO evolution, the field oscillates and, eventually, depending on the initial conditions, transitions.In contrast, in the NLO classical-statistical limit, there is a threshold value (|ϕ 0 | > 0.5) beyond which the system does not decay at all in our simulations but continues oscillating in the well. 3Therefore, in this case, we can say that quantum effects aid in the decay process.
Based on these observations, we can conclude that even if the system does not possess sufficient energy to surpass the potential barrier in the presence of classical-statistical fluctuations, quantum fluctuations drive the transition to the other side.A similar behavior was observed also in the case of initial conditions with a small displacement of the one-point function from the minimum but, instead, stronger fluctuations encoded in the two-point function.

C. Assessing NLO N = 1 approximation accuracy
Having compared the classical-statistical and the quantum dynamics of the false vacuum decay within the NLO approximation, we now want to check the accuracy of this approximation for N = 1.To do so, we simulate the full classical-statistical dynamics using the same parameters as in Fig. 9 (m 2 = −4, λ = 30, √ λϕ 0 /m = 0.3, f p = 0).The lattice volume was chosen large enough such that the results are volume independent.
The results from the classical-statistical simulations are illustrated in Fig. 11 by solid lines.The overall behavior of the dynamics of the one-point function is qualitatively the same as in the previous cases, with the NLO classical-statistical usually showing a faster evolution than the full classical statistical approach for this choice of parameters.The two-point function also shows different behavior: in full classical-statistical simulations, the graph shows a clearer peak during the field decay, which is not observed in NLO classical-statistical simulations.
One reason for this discrepancy of the large-N expansion for N = 1 is that higher-order contributions such as NNLO are no longer suppressed.We investigate the discrepancy between the two methods further in Appendix D. There we present a comparison between classical-statistical and NLO classical-statistical time evolution of the one-point function for a different choice of initial conditions for the fluctuations in terms of box initial conditions.

VII. CONCLUSIONS AND OUTLOOK
In our study, we presented a real-time approach to investigate the phenomenon of false vacuum decay in quantum field theory.We focused on the dynamics of a relativistic scalar field in 3 + 1 dimensions, initialized near a local minimum of a tilted double-well potential.
First, we considered a scenario where strong fluctuations are present.We employed the classical-statistical approximation to simulate the dynamics of the system in this regime.By this, we demonstrated that the realtime decay rates are comparable to those obtained from the bounce solution of the conventional imaginary-time approach.In general, we found that the decay rates are time dependent.Correspondingly, we extracted a timedependent effective potential, which becomes convex during the nonequilibrium transition process.
Second, we extended our analysis to the case of strongly coupled systems with low occupancy, where quantum effects are expected to play an important role.We used the nonequilibrium quantum evolution equations derived from the irreducible effective two-particle action truncated at NLO in 1/N to study the dynamics of the field and two-point functions.We demonstrated that quantum corrections can lead to transitions that are not captured by classical-statistical approximations.
Our results help bridging the nonequilibrium quantum field theory and the Euclidean approach to false vacuum decay.In contrast to the latter, the real-time method allows one to address the initial-value problem with the different stages of the subsequent transition.False vacuum decay is often not captured in terms of a (constant) decay rate but may be more efficiently characterized in terms of a time-dependent effective potential, which describes the transition from the metastable phase to the stable minimum.The notion of a time-dependent effective potential for a quantum description of false vacuum decay should be accessible in controlled experiments with ultracold quantum gases by measuring higher-order correlation functions [38,39].allows us to rearrange the fields in terms of these linear combinations (A2) All time-ordered connected correlation functions can be obtained by taking functional derivatives of the Schwinger functional W [J, R] defined as Z[J, R] = e iW [J,R] , with respect to the source terms and then setting the sources to zero.Specifically, the one-point function, or macroscopic field, reads and the connected two-point function, or propagator, Appendix B: Evolution equations The equations presented in Sec.VI are derived from first principles and capture the whole quantum evolution.However, due to their complexity, they generally cannot be solved without resorting to approximations.Therefore, we introduce a systematic nonperturbative approximation scheme known as the 1/N expansion.This appendix presents the evolution equations derived from the 2PI effective action under the 1/N expansion, truncated at NLO. See, e.g, [66] for more details.
To illustrate different approximation schemes in subsequent sections, we will consider an example of a quantum field theory involving a real, N -component scalar field φ a (x) with an O(N )-symmetric classical action corresponding to a potential as in (2) for a single component field.We specifically concentrate on the case of N = 1 from now on.The explicit form of the self energies in (34) are given by + I F (x, y)ρ(x, y) + P ρ (x, y)F (x, y) + P F (x, y)ρ(x, y) .

(B4)
In the previous equations, we denote the statistical and spectral components of the summation functions as where we set t 0 = 0, and with the field dependent resummation functions given by B8) and the field dependence is contained in the functions H ρ (x, y) ≡ −ϕ(x)ρ(x, y)ϕ(y) .
For practical computations, we assume spatial isotropy and homogeneity.Therefore, the integrals in momentum space will only depend on a single radial momentum.
Initial conditions must be prescribed to solve (34) and (35).The initial conditions for the macroscopic field can be written as Gaussian initial conditions for the statistical two-point function considered in this work correspond to (B12) In Sec.VI, we have chosen the initial quasiparticle occupation in vacuum, i.e., f p = 0.The initial conditions for the spectral function are determined by its antisymmetry, ρ(t, t, |p|) = 0, and the equal-time commutation relations of the underlying bosonic quantum theory.The source code used for the 2PI simulations is accessible in the Zenodo repository under [67].The numerical time evolution is computed using a spatial grid of N x = 64 lattice points and a lattice spacing of a x = 0.5, unless stated differently.Moreover, the time step is chosen to be a t = 0.01a x , guaranteeing energy conservation at the level of a few percent.

Appendix C: Classicality condition
This appendix provides additional details on the classical-statistical field theory framework.In this approach, the field is treated classically, and the dynamics involves sampling over initial conditions evolving each realization in real time according to the equation of motion of the classical field.Observables are obtained by averaging over classical trajectories.
A generating functional analogous to the one discussed in (A1) can be formulated for computing correlation functions within classical-statistical field theory.This generating functional enables us to calculate various statistical quantities and correlations.For a more comprehensive and detailed explanation of the subject, see [55].
The exact time evolution equations for ρ cl and F cl share the same structure as in (34), featuring self-energies Σ ρ cl , Σ F cl .In particular, we are interested in seeing how the previously discussed full 2PI equations of motion differ from their classical-statistical approximation.It is important to consider how the equations governing the evolution of the system are affected in this limit and the effects on the dynamics, which is explored in Sec.VI.One finds that the classical self-energies are obtained from the respective quantum ones as [55] Σ The classical-statistical self-energies are derived from their quantum counterparts by neglecting terms involving the product of two spectral functions ρ compared to the product of two statistical propagators F .Importantly, this condition does not require thermal equilibrium and can be applied to systems undergoing out-ofequilibrium dynamics, irrespective of the magnitude of the macroscopic field [33].
The comparison between quantum and classicalstatistical dynamics is made more transparent by rescaling the macroscopic field and statistical correlation function.Specifically, we have the following rescaling ϕ(x) → √ λϕ(x) , F (x, y) → λF (x, y) , (C2) whereas the spectral function ρ(x, y) remains unchanged.Similarly, the statistical self-energy rescales as Σ F (x, y) → λΣ F (x, y).Because of this rescaling, the classical self-energies also show an important property: they allow us to remove self-interaction influence from the dynamics.As a result, the coupling constant only affects the initial conditions.However, quantum corrections break this property and introduce terms proportional to ρ 2 , which become more significant as the coupling strength increases.

Appendix D: Starting with box initial conditions
This appendix investigates the qualitative and quantitative differences between the full classical-statistical and NLO classical-statistical dynamics in the same spirit of Sec.VI C.This time, we make a different choice of initial conditions.In order to skip the parametric resonance stage, it is convenient to start with a fixed nonvanishing background field ϕ 0 and initial occupation numbers of the form These initial conditions correspond to a box in momentum space up to the momentum scale Q = 1, with amplitude n 0 .We set ϕ 0 = 1.65/ √ λ, which is close to the false minimum of the potential for h = 0.9.We use V 3 = (32 × 1.5) 3 for the grid.In the left panel of Fig. 12, we compare the extracted transition timescale ∆t (that is, the time when the field becomes negative) resulting from the full classical-statistical simulations averaging over 100 runs (black triangles) and the simulations of the 2PI 1/N NLO classical-statistical (red dots).For high initial occupation numbers, the transition time of the full classical-statistical simulations is slower than the one of the NLO classical-statistical simulations.This is due to higher order contribution to the dynamics (NNLO, 0.8 0.9 1.0 1.1 Initial occupancy: λn 0 etc.) as discussed in the main text.Going towards lower initial occupation numbers, the transition timescales of the NLO classical-statistical simulations become slower than the full classical-statistical case, showing around λn 0 ≈ 0.9 an asymptotic value where the transition timescale diverges.To get more insight into the decay dynamics, we show in the right panel of Fig. 12 the time evolution of the field for three selected initial conditions from the left panel.The solid (dashed) lines indicate the full classical-statistical (NLO classical-statistical) field evolution.We observe that the full classical-statistical and NLO classical-statistical do not match well.In particular, the NLO classical-statistical dynamics is characterized by a sharp decay and shows oscillations around the true vacuum after the transition.The full classicalstatistical dynamics starts decaying earlier, and the transition takes a longer time to occur.This discrepancy requires further clarification.

FIG. 2 .
FIG.2.Time evolution of the one-point function ϕ(t), the equal-time statistical propagator F (t, t), and the effective temperature T eff (t).The black solid lines represent the correlation functions obtained from averages over n = 200 runs, while the gray dashed lines show a subset of single realization runs.

FIG. 3 .
FIG.3.Snapshots of the evolution for the classical field φ(t, x) from a two-dimensional slice using a single realization run: (top left) metastable phase, (top right) bubble nucleation, (bottom left) bubble expansion, and (bottom right) stable phase.

FIG. 5 . 4 .
FIG.5.Profiles of the bounce φ B (r) as a function of the radius r are shown for various selected times t, as presented in Fig.4.The inset displays the corresponding ratio between the bounce action S3 and the effective temperature T eff , obtained using(18).

FIG. 9 .
FIG.9.Transition from a false vacuum computed using 2PI NLO (dashed lines) and its approximation to the classicalstatistical dynamics (solid lines).Left: One-point function ϕ(t).Right: Equal-time statistical two-point function F (t, t).The simulations were realized for different initial conditions for the parameter controlling the asymmetry of the potential h.We employ m 2 = −4, λ = 30, √ λ m ϕ0 = 0.3.

FIG. 11 .
FIG.11.Transition from a false vacuum in full classical-statistical simulations.Left: Evolution of the one-point function ϕ(t) for different linear couplings h.Right: Equal-time statistical two-point function F (t, t) for different linear couplings h.We employ the same parameters as in Fig.9, and V3 = (1 × 256)3 , for various h as shown in the legend.Each curve is an ensemble of 20 runs.

9 FIG. 12 .
FIG. 12.Comparison of the full classical-statistical with NLO classical-statistical simulations for box initial conditions (D1).Left: Transition timescale ∆t in terms of the initial occupancy n0 for NLO classical-statistical (red dots) and full classicalstatistical simulations (black triangles) obtained averaging over 100 runs.The NLO classical-statistical simulations indicate around λn0 ≈ 0.9 an asymptotic value where the transition timescale diverges.Right: Evolution of the one-point function ϕ(t) for three selected initial conditions shown in the left panel.The dashed lines indicate the full classical-statistical and NLO classical-statistical are indicated by solid lines.We employ h = 0.9, √ λ m ϕ0 = 1.65.