Analog vacuum decay from vacuum initial conditions

Ultracold atomic gases can undergo phase transitions that mimic relativistic vacuum decay, allowing us to empirically test early-Universe physics in tabletop experiments. We investigate the physics of these analog systems, going beyond previous analyses of the classical equations of motion to study quantum fluctuations in the cold-atom false vacuum. We show that the fluctuation spectrum of this vacuum state agrees with the usual relativistic result in the regime where the classical analogy holds, providing further evidence for the suitability of these systems for studying vacuum decay. Using a suite of semiclassical lattice simulations, we simulate bubble nucleation from this analog vacuum state in a 1D homonuclear potassium-41 mixture, finding qualitative agreement with instanton predictions. We identify realistic parameters for this system that will allow us to study vacuum decay with current experimental capabilities, including a prescription for efficiently scanning over decay rates, and show that this setup will probe the quantum (rather than thermal) decay regime at temperatures $T\lesssim10\,\mathrm{nK}$. Our results help lay the groundwork for using upcoming cold-atom experiments as a new probe of nonperturbative early-Universe physics.

Since the pioneering early work of Coleman and collaborators [1][2][3], false vacuum decay (FVD) has primarily been studied using instanton methods, in which one obtains a semiclassical approximation of the decay rate by solving the equations of motion in imaginary time.These methods are made tractable by imposing O(d+1) symmetry on the resulting Euclidean 'bounce' solutions which describe the bubble nucleation event (with d the number of spatial dimensions).However, this symmetry assumption is broken on dynamical and/or inhomogeneous spacetimes that are relevant to cosmology, and precludes us from studying interesting and observationally important issues such as correlations between multiple bubbles [27,28].Furthermore, additional assumptions are required to interpret the instanton in real time; specifically, it is assumed that a critical bubble 'appears' at some instant in time.This prevents any study of the precursors of such an event in terms of the real-time dynamics of the field.
Recently, a promising new method for addressing these questions has emerged: the use of ultracold atomic Bose gases as quantum simulators of relativistic bubble nucleation [29][30][31][32][33][34][35][36][37][38][39][40].These systems exhibit coherent quantum behavior on scales that can be directly imaged in the laboratory, and can be manipulated into mimicking the dynamics of a Klein-Gordon field in a potential with true and false vacua.Cold-atom experiments have already been successfully used to study discontinuous phase transitions in quantum fields [41][42][43][44][45][46], including nonrelativistic thermal vacuum decay [47].Atomic simulators of relativistic FVD are now under active development by several groups, offering the prospect of studying vacuum decay in real time and in a controlled and reproducible manner, with the promise of new insights that complement those from long-established Euclidean techniques.These insights could have a transformative impact on our understanding of the early Universe, potentially helping to answer some of the most fundamental questions in cosmology, such as why there is more matter than antimatter [11][12][13], and whether our observable Universe is embedded in a larger 'multiverse' [6][7][8][9][10].
Previous analyses of these analogs have focused on their classical equations of motion, showing that these are equivalent to the Klein-Gordon equation for a relativistic field in the appropriate limit.Here we go further by calculating the spectrum of quantum vacuum fluctuations in the analog false vacuum state.This fluctuation spectrum is a crucial input for lattice simulations of the cold-atom system, in which the fluctuations are represented as classical stochastic variables in order to obtain a semiclassical approximation of the decay process.These simulations are our main theoretical tool for guiding the development of the analog experiments, and ultimately for helping us interpret the experimental data.
After describing our proposed analog system in Sec.II, we show in Sec.III that the false-vacuum fluctuation spectrum matches that of a Klein-Gordon field on scales where the classical analogy holds.This result was not guaranteed by the existing classical analogy, and thus provides further evidence for the suitability of this system as a relativistic analog.After an exhaustive search of the cold-atom literature, we identify a homonuclear potassium-41 mixture as the most promising experimental setup, and in Sec.IV we present a realistic set of parameters for a 1D realization of this system.This includes a protocol for scanning over parameters that allows us to vary the decay rate while keeping all other scales in the effective relativistic theory fixed.In Sec.V we then carry out a suite of semiclassical lattice simulations of this system, using our results for the fluctuation spectrum to generate realistic vacuum initial conditions.We verify that the field undergoes exponential decay as expected, and that the decay rate scales exponentially with the amplitude of the initial fluctuations, in qualitative agreement with the instanton prediction.Finally, in Sec.VI we explore the impact of finite temperatures on the decay rate, and argue that current experimental technologies can probe the regime of quantum rather than thermal decays.We summarize our results in Sec.VII, and discuss avenues for further development of this work.

II. THE ANALOG FALSE VACUUM
In this section we review the essential details of the analog FVD system we are interested in, as first proposed by Fialko et al. [30], and subsequently studied in Refs.[31][32][33][34][35][36][37].This system consists of a two-component Bose-Einstein condensate (BEC), with each atomic species described by a complex bosonic field ( The operators ψ † i (x) and ψi (x) create and annihilate atoms of species i in the position eigenstate |x⟩, respectively.Their amplitudes therefore determine the local number density of each species, ni (x) = ψ † i ψi , while their phases φi (x) encode coherent wavelike behavior and interference effects.The dynamics of these fields are described by the Hamiltonian which consists of a nonrelativistic kinetic term for each species, as well as a quartic self-interaction of strength g > 0 due to repulsive s-wave contact interactions between atoms.This interaction sets the characteristic energy scale of the BEC, E = gn, where n = ⟨n⟩ is the mean number density.The integral in Eq. ( 2) is over a finite spatial volume V that is either one-or two-dimensional, with the BEC confined tightly along the remaining dimensions, rendering them nondynamical.
We have specialized here to the case where both species have equal masses (m 1 = m 2 = m), equal intraspecies scattering (g 11 = g 22 = g), and zero interspecies scattering (g 12 = g 21 = 0).These conditions can be realized in practice by letting our two species be two different hyperfine states of the same atomic isotope, and applying an external magnetic field at the zero-crossing of a Feshbach resonance in the interspecies channel g 12 [31,48].Another possibility is to trap a single atomic species in a double-well potential; the atoms in each of the two wells then act as the two species, and only scatter with other atoms in the same well [49,50].
The Hamiltonian (2) excludes the usual external potential term that describes the trapping of the atoms along the extended direction(s).Our proposed experiment uses a 'box trap' which effectively approximates an infinitewell potential [51,52], so that the given Hamiltonian is accurate inside the trap.This is desirable for simulating relativistic physics as it maintains translation invariance in the interior region, with a near-homogeneous density profile.The density rapidly tapers to zero at the walls of the trap on a characteristic scale called the healing length, For the experimental parameters we consider here, this scale is smaller than the size of the BEC by a factor of 500 (see Table I).We therefore treat the field as homogeneous with periodic boundaries throughout this paper, as in most previous studies of this system [30][31][32][34][35][36][37][38].(This setup is also a reasonable approximation to a 1D ring trap, as used in e.g.Ref. [53].) Extending our results below to include the box trap and corresponding boundary conditions requires a calculation of the full spectrum of inhomogeneous eigenmodes, which has yet to be carried out for this system.We will present this calculation and its impact on bubble nucleation in an upcoming companion paper.
The two condensed species are coupled via a linear interaction term in the Hamiltonian, which allows atoms of species 1 to convert into species 2 (and vice-versa) at a rate ν that undergoes rapid modulation at some angular frequency ω, where ϵ ≪ 1 and λ = O(1) are dimensionless constants.
In the setup with two hyperfine states, this coupling is introduced by applying a modulated radio-frequency (rf) field; in the double-well case, ν instead represents the tunneling rate between the two wells.We integrate out the fast oscillation to obtain an effective Hamiltonian Ĥeff Potential for the analog relativistic field φ, as given by Eq. (10).There are stable 'true vacuum' (TV) states at every even integer value of φ/(πφ0).For λ > 1 there are also metastable 'false vacuum' (FV) states for every odd integer value.
that is valid on timescales much longer than ω −1 [54].At linear order in ϵ, we find This time-averaged picture fails to capture the presence of Floquet instabilities induced in modes whose natural frequencies are close to the driving frequency ω [32,34].
One expects that setting ω sufficiently large (i.e., making the wavelengths of the unstable modes sufficiently short) will cause these instabilities to be quenched by damping effects on small scales; however, the exact nature of this process is still an open question.The relevance of the effective Hamiltonian (6) for quantum simulation comes from considering the field which is proportional to the relative phase between the two species. 3On scales much larger than the healing length, the classical equation of motion for this degree of freedom is identical to that of a relativistic scalar field,  (20) for the relative modes of the analog system.Right panel : Fluctuation power spectrum for the effective relativistic field φ.Both quantities interpolate between being Klein-Gordon-like (13) in the IR (ξk ≪ 1) and nonrelativistic (21) in the UV (ξk ≫ 1).The vertical dashed line in each panel indicates the crossover between these two regimes.Both use our fiducial parameters, given in Table I.
where we identify the 'speed of light' as Note that in reality this is the sound speed of phonons in the BEC, which is roughly eleven orders of magnitude smaller than the speed of light in vacuum.However, as we see below, it plays exactly the same role as the speed of light in the effective relativistic theory that emerges on large scales.The potential appearing in Eq. ( 8) is ) As shown in Fig. 2, this contains a series of true vacua at φ tv /φ 0 = 2jπ, j ∈ Z, and for λ > 1, a series of false vacua at φ fv /φ 0 = (2j + 1)π.These correspond to the two atomic species being in phase and in antiphase, respectively; the linear coupling means that there is an additional energy density of order ϵgn 2 associated with being in antiphase, while the modulation generates an effective potential barrier that makes this state metastable.Increasing the amplitude of the modulation via λ creates a deeper potential barrier, and increases the mass of fluctuations in the false vacuum,

III. QUANTUM FLUCTUATIONS IN THE FALSE VACUUM
We have reviewed the known result that, on scales much larger than the healing length, an atomic Bose-Bose mixture can reproduce the classical equation of motion of a Klein-Gordon field (8) with a false vacuum potential (10).However, vacuum decay is inherently quantummechanical, so it is important to test whether these systems are also analogous at the quantum level.Here we perform this test by calculating the power spectrum of fluctuations in the false vacuum state |Ω fv ⟩, where φk are the Fourier modes4 of the effective relativistic field (7).Below we find that, on scales much larger than the healing length (ξk ≪ 1), this spectrum asymptotically matches that of the corresponding Klein-Gordon field, with corrections suppressed by powers of (ξk) 2 and ϵ.To derive this result, we adopt the standard mean-field approximation [55] in which each atomic field consists of small quantum fluctuations around a highly-occupied classical condensate wavefunction, The factor (−1) here reflects the fact that the two species are in antiphase in the false-vacuum state.We expand around a homogeneous mean-field wavefunction, whose phase evolves at a rate set by the chemical potential, µ = (1 + ϵ)gn.To study the dynamics of the fluctuations, it is convenient to remove this time evolution with a canonical transformation ψi → e iµt/ℏ ψi .This modifies the Hamiltonian to Expanding this new Hamiltonian to quadratic order in the fluctuations, we find that it can be written as with K 0 a constant energy offset associated with the meanfield solution, and separate terms K± governing the total and relative fluctuation modes, with the normalization chosen such that the modes obey canonical bosonic commutation relations.The field we are interested in is defined solely in terms of the relative modes, and at linear order in the fluctuations is given by We can therefore ignore the dynamics of the total modes for now, given that they are decoupled in the linear regime.
(We return to them in Sec.VI, as they play a significant role in the presence of thermal noise.)To calculate the power spectrum (12), we must determine the eigenstates of the relative Hamiltonian K− and identify |Ω fv ⟩ as the lowest-lying of these states. 6We can do this by writing the Hamiltonian in diagonalized form, so that each normal mode, described by the ladder operators âk , â † k , acts as an independent harmonic oscillator.The false vacuum |Ω fv ⟩ is then identified as the state annihilated by âk for all wavenumbers k.In Appendix A we identify the appropriate Bogoliubov transformation relating the normal modes to the relative atomic field modes ψ− k , ψ− † k .The energy associated with excitations of the normal modes is given by which, on scales much larger than the healing length (ξk ≪ 1), reduces to the dispersion relation ( 13) of a Klein-Gordon field of the same false vacuum mass (11) we found in our classical analysis of the equations of motion.We can directly evaluate the power spectrum ( 12) by writing the Fourier modes φk in terms of the normal modes âk and using standard ladder operator identities.In the same IR limit as before we find Eq.( 13), which is exactly what we expect for the corresponding Klein-Gordon field.
We already know from our classical understanding of the system that the relativistic analogy breaks down on scales much smaller than the healing length (ξk ≫ 1).In this limit, we recover a white-noise fluctuation spectrum and the usual nonrelativistic dispersion relation, The former represents an excess of power at small scales compared to the Klein-Gordon spectrum ( 13), due to nonrelativistic, high-momentum excitations of individual atoms.The interpolation between this regime and the Klein-Gordon-like results on large scales is shown in Fig. 3.

IV. EXPERIMENTAL PARAMETERS
Our results for the false vacuum power spectrum are a general feature of the modulated Bose-Bose mixture system described in Sec.II, regardless of any particular experimental realization.In this section, we describe a concrete set of experimental parameters (summarized in Table I) that is achievable with current cold-atom experiments, and which will allow us to probe the physics of relativistic vacuum decay.
As highlighted in Sec.II, among the key requirements for our system are that both atomic species have equal masses (m 1 = m 2 ), equal intraspecies scattering lengths (a 11 = a 22 ), and negligible interspecies scattering (a 12 = 0). 7It is easy to select equal masses by using two hyperfine states of the same atomic isotope (i.e., a homonuclear mixture).However, the conditions on the scattering lengths are more difficult to arrange.It is possible to set a 12 to zero by applying an external magnetic 7 These 3D scattering lengths a ij determine the corresponding 1D interaction strength, g ij = 2ℏω ⊥ a ij , where ω ⊥ is the frequency of the transverse harmonic potential, Vtrap = is satisfied at the zero-crossing of a 12 .The 41 K resonance specified above is therefore the optimal candidate system for simulating relativistic vacuum decay.
The main technical challenge with this setup is that the resonance has a width of only 155.8 mG [59], necessitating a very high level of magnetic field stability in order to stay at the zero-crossing of a 12 , as illustrated in Fig. 4. Nonetheless, this level of stability is achievable with current experimental technologies.In particular, Borkowski et al. [67] have recently demonstrated magnetic field stability at the level of ∼ 2 ppm in a cold-atom experiment.For our proposed system this corresponds to |a 12 | ≤ 0.53 a 0 (where a 0 = 5.292 × 10 −11 m is the Bohr radius).This is less than 1% of the mean intrastate scattering length a = 60.24 a 0 , which should be sufficient precision for our purposes.
Given the 3D scattering properties of the two atomic species, the behavior of the effective 1D system is set by the number of condensed atoms, the size of the trap along the elongated and transverse directions, and the strength and modulation of the applied radio-frequency field.We have explored this parameter space with the goal of maximizing the natural condensate energy scale gn relative to the thermal energies k B T that can be achieved in current experiments, as this will allow us to investigate the regime of quantum (rather than thermal) decays.At the same time, we have ensured that this energy scale is not so high that transverse modes of energy ℏω ⊥ are excited, where ω ⊥ is the frequency of the harmonic trapping potential in the transverse directions.(We plan to test this explicitly in future work with 3D simulations that resolve the transverse directions.) In order to facilitate comparisons with instanton predictions (which are challenging to calibrate at any single point in parameter space), it is useful to vary the system parameters to scan over a broad range of bubble nucleation rates.The instanton decay rate per unit volume in this model scales as where n ≡ ξ d n is the dimensionless condensate number density (i.e., the number of atoms in a region of size equal to the healing length).In d = 1 dimensions the dependence on ϵ vanishes, and the decay rate is thus primarily controlled by n.This parameter also sets the size of fluctuations in the field relative to the characteristic value φ 0 , We find that it is possible to vary n while keeping the energy scale gn (and therefore all other dimensionless parameters of the system) fixed, by simultaneously increasing the number of atoms of each species N and decreasing the transverse trapping frequency ω ⊥ .This allows us to perform a controlled test of how the bubble nucleation rate scales with the amplitude of the initial fluctuations.
Our proposed parameters are summarized in Table I.We vary n by a factor of 5, which is sufficient to see a significant variation in the decay rate.As we show in Sec.VI below, the energy scale gn here is large enough that the quantum-decay regime is readily accessible to current or near-future experiments.

V. LATTICE SIMULATIONS
Part of the value of our results on the vacuum power spectrum in Sec.III is that they can be used as an input for semiclassical lattice simulations of the cold-atom system.These simulations are a powerful tool for exploring the real-time dynamics of bubble nucleation, and are a crucial ingredient for developing and interpreting analog FVD experiments.The key idea is to encode the nonclassical nature of the problem in the initial conditions of the simulation, by drawing an ensemble of random field realizations that sample vacuum fluctuations around the homogeneous false vacuum state [68].These realizations are then evolved forward by numerically integrating the classical equations of motion.This approach is widely used in the context of atomic physics and quantum optics (where it is referred to as the 'truncated Wigner approximation' [69][70][71]), and also underpins cosmological lattice simulations of inflation and preheating [72][73][74][75][76][77][78][79][80][81][82][83] as well as vacuum decay [68,84,85].
It is common for lattice simulations of cold-atom systems to initialize the fluctuations using a white-noise power spectrum (21) [30,31,33,34], particularly in situations where the processes of interest are insensitive to the precise form of this spectrum.Bubble nucleation, however, is extremely sensitive to the statistics of the initial fluctuations, as different initial states can decay at exponentially different rates.(For example, we see from Eq. ( 22) that there is an exponential sensitivity on n.)The vacuum fluctuation spectrum derived above is therefore a crucial ingredient for realistic simulations of analog vacuum decay.
In this section we use a suite of lattice simulations to study bubble nucleation from vacuum initial conditions in the 1D cold-atom system described in Sec.IV.We extract decay rates for different values of the fluctuationamplitude parameter n, and verify that the rates depend exponentially on this parameter, in agreement with the scaling found in the instanton approach.We perform the same test with white-noise initial conditions, and find decay rates that are globally larger than in the vacuum case.This confirms that vacuum decay in semiclassical lattice simulations is indeed sensitive to the statistics of the initial fluctuations, and that for the cold-atom system these must be correctly specified using Bogoliubov theory, as we have done here.We additionally investigate the conservation of the Noether charges of the effective Klein-Gordon theory in our simulations of the cold-atom system, as these are a useful diagnostic for the faithfulness of the relativistic analogy.

A. Code setup
We use a Fourier pseudospectral code with an eighthorder symplectic time-stepping algorithm [86] (see Appendix B for details), and work in units where the atomic mass m, healing length ξ, and sound speed c are set to unity (which is equivalent to also setting ℏ = gn = 1).Our simulations work at the level of the time-dependent Hamiltonian (4), resolving the modulation of the interspecies coupling so that we can test for the emergence of the effective time-averaged dynamics.
We simulate a system with the experimental parameters specified in Table I.In code units, this setup is realized by .Left panel : Survival probability for the false vacuum state as a function of time, as estimated using ensembles of 1024 simulations for each curve.We scan over the dimensionless number density n = ξn to probe a broad range of decay rates.The gray shaded region (0.01 ≤ Pr(survive) ≤ 0.5) is used to fit an exponential decay rate Γ for each curve (shown as dashed lines).Right panel : Dimensionless decay rate per unit volume as a function of n, computed for both vacuum and white-noise initial conditions.Both curves are well-described by a linear fit, as expected from instanton calculations.The white-noise case consistently gives faster decays despite the smaller initial phase fluctuations, due to this being an excited state of the system.evolving a periodic region of volume V /ξ d = L/ξ = 500, and setting ϵ = 2.5 × 10 −3 and λ = √ 2 so that the false vacuum mass is m fv /m = 0.1.We additionally set the dimensionless modulation frequency to ωξ/c = 680, which is sufficiently large that the Floquet instability bands are above the Nyquist frequency for all of our simulations.This allows us to model the expected experimental situation where these instabilities are damped by the small-scale dynamics of the BEC, and do not affect the evolution of the IR modes; the actual experimental value of ω is unimportant so long as the Floquet instabilities are quenched.Our simulations use 2048 lattice sites and a timestep that is 1/16 times the modulation period 2π/ω, giving spatial and temporal resolution of ∆x/ξ ≈ 0.244 and c∆t/ξ ≈ 5.77 × 10 −4 , respectively.In Appendix B we show that our results are numerically converged at this resolution, and that the Noether charges of the cold-atom Hamiltonian (4) are conserved to within a few parts per billion.

B. Bubble nucleation rates
We extract decay rates for the analog system using ensembles of 1024 simulations, with each simulation corresponding approximately to a different possible classical history drawn from the path integral describing the full evolution of the many-body quantum state.We initialize each simulation as the homogeneous false vacuum φ = πφ 0 plus independent random draws of the vacuum fluctuations δ φ.We treat the latter as a zero-mean Gaussian random field with a power spectrum that (as shown in Fig. 3) interpolates between a relativistic spectrum in the IR and a white-noise spectrum in the UV.We have checked that this power spectrum remains statistically stationary over time by averaging over the ensemble of nondecayed trajectories, effectively testing that our initial state is indeed an eigenstate of the Hamiltonian near the false vacuum.
As well as the relative phase, we also initialize the relative density and the total phase and density using random draws from their corresponding vacuum spectra.It is crucial to initialize all four fields in this way to correctly capture the vacuum state.For example, neglecting the relative density fluctuations corresponds to initializing the effective Klein-Gordon field with zero momentum everywhere, when in fact this momentum field should also contain vacuum fluctuations.In practice, we initialize the total and relative atomic field modes in our code, which is equivalent at the linear level to working in terms of the density and phase fields.
We find that it is crucial that the positive-and negativemomentum Fourier modes ψ k and ψ −k are not treated as statistically independent random variables.Instead, one must draw the positive-and negative-momentum normal modes a k , a −k independently, and then obtain the Fourier modes of the atomic fields via a reverse Bogoliubov transformation.This induces a nontrivial correlation between ψ k and ψ −k that appropriately captures the quantum statistics of the false vacuum state.Failing to include these correlations in the initial conditions puts the system into an excited state that nucleates bubbles much more rapidly than the false vacuum state, and much more even than the white-noise state.
We truncate all of the fluctuation spectra at a maximum wavenumber of ξk UV ≈ 3.22, which is a factor of 4 smaller The initial fluctuations in each simulation are identical except for an overall ∼ n−1/2 scaling.The level of violation is roughly stationary throughout each simulation, and approaches zero for large n, despite the nonrelativistic behavior of the system on small scales.
than the Nyquist frequency of our simulations, ξk Nyq = πξ/∆x ≈ 12.9.Evidence from pure Klein-Gordon lattice simulations [87] suggests that changing this cutoff modifies the decay rate in a way that can be absorbed into a renormalization of the bare model parameters.We leave a detailed investigation of this effect in the analog system for future work, and here use a fixed UV cutoff for all of our simulations.The amplitude of the fluctuations relative to the homogeneous value of the field is set by the dimensionless number density n, which we scan over in the experimentally accessible range 10 ≤ n ≤ 50.We measure a decay rate from each ensemble of simulations by counting the number of nondecayed trajectories as a function of time, dividing by the total number of simulations to obtain an estimate of the time-dependent survival probability.In doing so, it is necessary to choose a definition for when an individual realization has decayed.We do this by setting a threshold on the volume average of the cosine of the relative phase, ⟨cos(φ/φ 0 )⟩ V .This quantity fluctuates near to −1 in the false vacuum, and grows rapidly after a bubble nucleates before saturating near +1 once the transition has percolated, as illustrated in Fig. 5.We compute the decay threshold separately for each ensemble as the lowest possible value of ⟨cos(φ/φ 0 )⟩ V for which no more than 1% of the simulations cross back below the threshold in any given timestep. 8  8 A more obvious choice would be to allow zero downward crossings through the threshold, as this would capture the notion that vacuum decay is an irreversible process.However, we find that enforcing zero downward crossings makes the algorithm easily confused by small fluctuations in ⟨cos(φ/φ 0 )⟩ V , and results in a choice for the threshold that is far too conservative.Manual Our resulting estimates of the survival probability are shown in the left panel of Fig. 6.As expected, the ensembles with smaller n, and therefore larger initial fluctuations, decay on much shorter timescales.After an initial transient, each ensemble reaches a regime of exponential decay, We fit a decay rate Γ to each curve, restricting the fit to survival probabilities between 50% and 1% in order to exclude the nonexponential regime at early times and noisy small-number statistics at late times, respectively.The resulting decay rates (in dimensionless units, and measured per unit volume) are shown in blue in the right panel of Fig. 6, and are well-described by an exponential scaling with respect to n, in qualitative agreement with the instanton prediction (22).
It is important to note however that the proportionality constant linking log(Γ/V ) and n does not agree with the instanton prediction; our simulations decay significantly faster than predicted in the instanton approach.This same behavior has been observed in pure Klein-Gordon lattice simulations [68], and is an expected consequence of performing instanton calculations using the bare lattice parameters, rather than the renormalized theory [87].It is also worth pointing out that our instanton calculations are based on the effective Klein-Gordon theory, rather than the full analog system, and therefore neglects effects such inspection of the results with a 1% allowance for downward crossings confirms that this accurately captures the common-sense notion of when the field has decayed (e.g.see Fig. 5).We have checked that varying this allowed fraction between 0.5% and 2% does not significantly impact our measured decay rates.
as the excess small-scale power identified in Sec.III.We plan to explore these issues in the context of the analog system in future work.
As well as our simulations using vacuum initial conditions, we carry out a suite of simulations using white-noise initial conditions.This corresponds to the nonrelativistic UV limit (21) of the full power spectrum derived from Bogoliubov theory, and matches the prescription used by several previous studies of vacuum decay in cold-atom analog systems [30,31,33,34].The resulting decay rates are shown in red in the right panel of Fig. 6.These are fit only to survival probabilities between 20% and 1%, as we find that it takes longer for these initial states to settle into a period of steady exponential decay.We see that, while the resulting decay rates also follow the expected exponential scaling with n, they are globally larger for white-noise initial conditions than for the vacuum case, despite the fact that the actual amplitudes of the fluctuations are smaller in the IR in the white-noise case (compare the blue and purple curves in Fig. 3).We interpret this as evidence that white-noise fluctuations correspond to an excited state of the analog system, and thus lead to faster decays, on average, than the vacuum initial conditions we have derived here.
Note that this does not imply that the white-noise spectrum is somehow unphysical.In fact, such a spectrum is the vacuum state for an alternative system with zero atomic scattering, g = 0.The enhanced decay rates shown in red in Fig. 6 can thus be interpreted as being due to a mismatch between the Hamiltonian describing the initial conditions and the Hamiltonian describing the time evolution.

C. Verifying Klein-Gordon behavior
While our results for the decay rates are in broad agreement with our expectations for relativistic vacuum decay, we can also directly test whether the relative phase field φ is indeed analogous to a relativistic Klein-Gordon field by computing the Noether charges for the corresponding Klein-Gordon theory, Since the Noether charges for the underlying nonrelativistic Hamiltonian are conserved with extremely high precision in our simulations (see Appendix B), any nonconservation of the Klein-Gordon charges (25) should be interpreted as being due to limitations of the relativistic analogy, rather than numerical errors.In Fig. 7 we show the violation of these charges for a series of simulations with a broad range of dimensionless number densities n.We find that violations in the Klein-Gordon energy and momentum are roughly stationary over time, and reach a regime where they scale like For our fiducial parameters this can be translated into a physical temperature using Eq. ( 27).The decay rates (extracted by fitting in the gray shaded region, shown here as dashed lines) are consistent with being temperature-independent up to roughly T ≈ 0.06; beyond this point, large fluctuations in the total modes couple to the relative modes and ruin the effective relativistic picture.
|∆H|/|H| ∼ n−1 and |∆P |/|P | ∼ n−1/2 respectively, so that in the limit of small fluctuations the analogy holds with high accuracy.However, in the experimentally accessible regime n ∈ [10,50] that we are interested in here, the violation is on the order of at least a few percent in the energy.In the momentum, the relative errors reach order unity, although this reflects the fact that the total momentum of the field averaged over the entire volume V is intrinsically close to zero.While we do not believe these errors invalidate the mapping onto the Klein-Gordon theory, further improvements in the accuracy of the analog may be possible.Specifically, so far we have ignored the backreaction of the fluctuations onto the mean-field dynamics, which would modify this mapping in a way that could plausibly be absorbed into a renormalization of the parameters of the effective Klein-Gordon theory.(Similar effects have recently been investigated in the case of pure Klein-Gordon theory [87].)This would be consistent with our finding that the level of charge violation scales with the fluctuation amplitudes.We conjecture that accounting for these corrections and identifying the appropriate Klein-Gordon parameters could substantially improve the level of charge violation over that shown in Fig. 7, and also bring our decay rates into closer quantitative agreement with the instanton prediction.We plan to explore this in detail in future work. .Enhancement in the fluctuation power spectra of the total and relative phonons as a function of temperature.The vertical axis shows the ratio between the finite-temperature and zero-temperature power spectra evaluated at the minimum wavenumber kIR = π/L ≈ 6.28 × 10 −3 ξ −1 , for which the enhancement is maximized.

VI. FINITE-TEMPERATURE EFFECTS
Thus far we have considered only zero-temperature states of the analog system.However, any realistic experiment will inevitably be at some finite temperature, and will therefore contain thermal as well as quantum fluctuations.These are potentially a nuisance factor in studying quantum vacuum decay, giving an excess contribution to the decay rate and altering the phenomenology of the nucleated bubbles [4].It is therefore valuable to estimate the temperature threshold at which these deviations from the zero-temperature case become significant, as this can then guide the development and interpretation of the analog experiments.
In the framework of the truncated Wigner approximation, we can model the thermal bath by including additional fluctuation power in our initial conditions. 9This amounts to replacing vacuum expectation values with traces over a thermal density matrix, resulting in a scaledependent enhancement to the relative phase power spectrum, as well as for the relative density, and the total phase and density.(Here coth x = (1 + e −2x )/(1 − e −2x ) is the 9 Other prescriptions and theoretical frameworks exist, including modeling the effects of the thermal bath by adding a stochastic driving term to the Gross-Pitaevskii equations [35,37,39,88].However, our treatment here allows us to model quantum and thermal fluctuations in a simple and conceptually unified way.A detailed comparison against alternative simulation methods would be interesting, but is beyond our present scope.
hyperbolic cotangent function.)It is convenient to work in terms of the dimensionless temperature, where the numerical value corresponds to our particular choice of experimental parameters (cf.Table I).Fig. 8 shows the survival probability in ensembles of simulations at various temperatures, with n = 40.For dimensionless temperatures T ≲ 0.06 we see that, notwithstanding some differences in the initial nonexponential transient phase, the exponential decay rates are all consistent with the zero-temperature result.At higher temperatures, rather than finding an enhanced rate of relativistic decays, we instead find that the exponential decay model becomes an increasingly poor fit to the empirical survival probabilities.We interpret this finding as indicating the breakdown of the relativistic analogy at high temperatures, and conjecture that this breakdown is due to the impact of thermal noise on the total phonon modes.In contrast to the relative modes, which have an effective mass m fv due to the potential barrier around the false vacuum, the total modes have a massless dispersion relationship ω k ≃ ck in the IR, allowing them to become excited to very large amplitudes by the thermal bath, as illustrated in Fig. 9.The coupling between the total and relative modes then becomes significant, and spoils the effective relativistic dynamics of the relative modes.As evidence for this interpretation, we note that the T ≲ 0.06 threshold determined empirically from our simulations is just below the theoretically-predicted threshold at which the total modes of this system should lose phase coherence, Tϕ = n/ L = 0.08 [35].
Our results show that dimensionless temperatures of T ≲ 0.06 should give us access to a setting closely resembling the zero-temperature dynamics of the analog vacuum decay process.This translates into physical temperatures of T ≲ 10.9 nK for our proposed parameters.Note that our interpretation in terms of the phase coherence temperature Tϕ = n/ L implies that this threshold should scale proportionally with the fluctuation-amplitude parameter n, so that the T ≲ 10.9 nK benchmark should be viewed as a minimal requirement, with lower temperatures giving us access to vacuum decay rates over a broader range of parameter space.This benchmark is readily accessible with current experimental setups, which routinely reach temperatures on the order of a few nK, and have even recorded temperatures as low as tens of pK [89].

VII. SUMMARY AND OUTLOOK
Quantum analog experiments present a powerful new tool for understanding relativistic vacuum decay.Here we have carried out a detailed study of one such proposed experimental setup, which uses a rapidly modulated coupling between two atomic Bose-Einstein condensates to engineer a metastable false vacuum state for the relative phase.We have derived the spectrum of quantum fluctuations around this state, and have shown that this spectrum asymptotically matches that of the effective Klein-Gordon field in the IR.
As well as providing further evidence for the suitability of the cold-atom analog for studying relativistic physics, this vacuum fluctuation spectrum is also a crucial input for semiclassical lattice simulations of this system.By carrying out a suite of such simulations, we have confirmed the key theoretical expectations for the analog false vacuum: that it undergoes exponential decay, at a rate that is exponentially sensitive to the amplitude of the vacuum fluctuations.We have also shown that using an alternative fluctuation spectrum -in this case, white noise, which has been used in several previous studies of this system -leads to an enhanced decay rate compared to the pseudorelativistic vacuum fluctuations, as this corresponds to putting the system in an excited initial state.
In carrying out these simulations, we have identified a realistic set of parameters that will allow us to study vacuum decay with current experimental capabilities.This includes a protocol for scanning over fluctuation amplitudes, and thus decay rates, while keeping all other natural scales of the system fixed, enabling detailed and controlled experimental studies of the decay rate.
As well as the zero-temperature fluctuation spectrum, we have derived the enhancement of the fluctuation power due to thermal noise at finite temperature.We find that, so long as the system is below a given temperature threshold (which we argue is set by the coupling between the total and relative phase degrees of freedom), the decay rate extracted from our simulations is consistent with that at zero temperature.For our proposed parameters, this threshold lies well within reach of current experiments, meaning that we should be able to empirically test the physics of quantum bubble nucleation in the near future.
Our results here rely on several simplifying assumptions, which we plan to relax in future work.In particular, we have treated the BEC system as periodic, neglecting boundary effects due to the external trapping potential.In a forthcoming companion paper, we will generalize our Bogoliubov analysis to derive the inhomogeneous vacuum fluctuations in a box trap, and investigate the impact of these boundary effects on the bubble nucleation rate.We have also neglected in our calculations the backreaction of the fluctuations onto the mean-field dynamics of the BEC, and corresponding renormalization of the bare parameters of the effective relativistic theory.Incorporating these effects should allow for a more precise understanding of the validity of the relativistic analogy, improve the initialization and interpretation of our lattice simulations, and enable more detailed comparisons with instanton predictions.These developments will enable the first experimental tests of relativistic vacuum decay. .Pointwise convergence of our numerical solutions for increasing spatial and temporal resolution (left and right panels, respectively) in simulations with n = 30.Each curve shows the maximum absolute pointwise difference in cos(φ/φ0) between one solution with the stated resolution and another with double the spatial or temporal resolution, starting from identical initial conditions.The vertical dotted lines show the time of bubble nucleation in the converged simulations.Note that the resolution used in our simulations discussed in Secs.V and VI corresponds to the red curves here.which we have split into a linear and a nonlinear piece.Each piece has its own chemical potential, which can be chosen for convenience -e.g., to minimize sinusoidal oscillations in the homogeneous mode of the total phaseas these have no effect on the relative phase φ.where the dimensionless coefficients a i , b i , (i = 1, . . ., k) are chosen such that the integrator is exact to order n in the small timestep δt.Integrators of this form are symplectic, in the sense that they exactly conserve phase space volume.We implement an efficient realization of this integrator from Yoshida [86], which uses k = 16 steps and is accurate to order n = 8.
In Fig. 10 we show convergence tests of our code for increasing spatial and temporal resolution, measuring numerical errors in terms of pointwise differences in the cosine of the relative phase field, cos(φ/φ 0 ).For the level of resolution used in our simulations in Secs.V and VI, we see that the maximum error is on the order of ∼ 10 −7 prior to bubble nucleation, and at most ∼ 10 −5 even long after bubble nucleation.This indicates that our simulations are numerically converged, even in the highly dynamical nonlinear regime.
We also test our code by checking for violations in conservation of the Noether charges associated with the cold-atom Hamiltonian (4), which correspond to the total number of atoms and the total momentum of the system, respectively [34].(Note that the total energy is not exactly conserved, due to the explicit time-dependence of the rf modulation term in the Hamiltonian.)As shown in Fig. 11, both charges are conserved to the level of a few parts per billion in simulations at our fiducial resolution.

Figure 1 .
Figure 1.Lattice simulation of vacuum decay in the 1D analog system.Nonlinear interactions between fluctuations around the false vacuum (blue) lead to the nucleation of a true vacuum bubble (red), which then expands relativistically.The simulation shown here corresponds to the blue curves in Fig. 7, and conserves the Hamiltonian of the effective relativistic theory to within ∼ 10% (see discussion in Sec.V C).

Figure 4 .
Figure 4.The three scattering lengths a11, a22, a12 of our proposed homonuclear 41 K mixture as a function of magnetic field strength.The quoted values and gray shaded region correspond to ±2 ppm ≈ ±1.4 mG either side of the zerocrossing, as given by Ref.[59].

Figure 5 .
Figure 5. Individual random realizations of vacuum decay in our n = 35 ensemble, showing how the volume-averaged cosine of the relative phase field evolves over time.Each curve corresponds to an independent simulation, which oscillates near ⟨cos(φ/φ0)⟩ V = −1 until a true vacuum bubble nucleates, at which point the trajectory grows until it saturates near ⟨cos(φ/φ0)⟩ V = +1.The colored curves are three randomly-selected trajectories, highlighted to illustrate the typical behavior.The black dotted line shows our empiricallydetermined decay threshold for this ensemble, as found using the procedure described in the main text.

Figure 6
Figure 6.Left panel : Survival probability for the false vacuum state as a function of time, as estimated using ensembles of 1024 simulations for each curve.We scan over the dimensionless number density n = ξn to probe a broad range of decay rates.The gray shaded region (0.01 ≤ Pr(survive) ≤ 0.5) is used to fit an exponential decay rate Γ for each curve (shown as dashed lines).Right panel : Dimensionless decay rate per unit volume as a function of n, computed for both vacuum and white-noise initial conditions.Both curves are well-described by a linear fit, as expected from instanton calculations.The white-noise case consistently gives faster decays despite the smaller initial phase fluctuations, due to this being an excited state of the system.

× 10 5 Figure 7 .
Figure 7. Fractional violation of the Klein-Gordon charges (25) as a function of the dimensionless BEC number density n.The initial fluctuations in each simulation are identical except for an overall ∼ n−1/2 scaling.The level of violation is roughly stationary throughout each simulation, and approaches zero for large n, despite the nonrelativistic behavior of the system on small scales.

Figure 8 .
Figure 8. Survival probability as a function of dimensionless temperature T = kBT /(gn), for n = 40.For our fiducial parameters this can be translated into a physical temperature using Eq.(27).The decay rates (extracted by fitting in the gray shaded region, shown here as dashed lines) are consistent with being temperature-independent up to roughly T ≈ 0.06; beyond this point, large fluctuations in the total modes couple to the relative modes and ruin the effective relativistic picture.

Figure 9
Figure 9. Enhancement in the fluctuation power spectra of the total and relative phonons as a function of temperature.The vertical axis shows the ratio between the finite-temperature and zero-temperature power spectra evaluated at the minimum wavenumber kIR = π/L ≈ 6.28 × 10 −3 ξ −1 , for which the enhancement is maximized.

Figure 10
Figure 10.Pointwise convergence of our numerical solutions for increasing spatial and temporal resolution (left and right panels, respectively) in simulations with n = 30.Each curve shows the maximum absolute pointwise difference in cos(φ/φ0) between one solution with the stated resolution and another with double the spatial or temporal resolution, starting from identical initial conditions.The vertical dotted lines show the time of bubble nucleation in the converged simulations.Note that the resolution used in our simulations discussed in Secs.V and VI corresponds to the red curves here.

Figure 11 .
Figure 11.Relative variations in the Noether charges (B6) for the same simulation shown in the red curves of Fig. 10.
We have performed an exhaustive search of other known Feshbach resonances in homonuclear mixtures of stable bosonic isotopes of the alkali metals ( 7 Li [31, 56], 23 Na [57, 58], 39 K [59, 60], 85 Rb [61, 62], 87 Rb [63, 64], and 133 Cs [65, 66]), and have not found any other interstate resonances where the condition a 11 ≃ a 22 [48]ω 2 ⊥ (y 2 + z 2 ).fv = 4ϵ(λ 2 − 1) m = 0.1 m TableI.List of fundamental and derived parameters for our proposed 1D cold-atom experiment.Here u = 1.661 × 10 −27 kg is the unified atomic mass unit and a0 = 5.292 × 10 −11 m is the Bohr radius.The scattering length a quoted here is the mean of the two intrastate scattering lengths; the difference is ∼ 1% (cf.Fig.4).The number density n and scattering strength g are scanned over by varying the number of atoms of each species N and the harmonic trap frequency ω ⊥ respectively, while holding the energy scale gn constant.fieldatthezero-crossing of a Feshbach resonance[48], but there is then no further freedom to tune a 11 and a 22 in order to set them equal to each other.Fortunately, as pointed out by Fialko et al. [31], 41 K (potassium-41) possesses a Feshbach resonance between the |F, m F ⟩ = |1, 0⟩ and |1, +1⟩ states with a zero crossing at B ≃ 675.256G, where the condition a 11 ≃ a 22 is realized naturally with a precision of ∼ 1%.
Evolution under each of these operators individually can be solved exactly; the nonlinear piece conserves the amplitude of each field and simply performs a local phase rotation, x→k represents a Fourier transform, and F −1 k→x its inverse.(Theseare implemented numerically as fast Fourier transforms, so that in practice Eq.(B4) is only exact under the assumption that the fields are band-limited with maximum wavenumber less than or equal to the Nyquist frequency on the lattice.)Whilethere is no exact solution for the evolution under O = O lin + O nlin from generic initial data, we can approximate this full evolution by chaining together a series of short steps with each of the individual operators,ψ(x,t 0 + δt) = e −ia1O lin δt ℏ e −ib1O nlin δt ℏ × • • • × e −ia k O lin δt ℏ e −ib k O nlin δt ℏ ψ(x, t 0 ) + O δt n+1 , × cos R(t, t 0 ) i sin R(t, t 0 ) i sin R(t, t 0 ) cos R(t, t 0 ) F x→k {ψ(x, t 0 )} , R(t, t 0 ) = ϵgn t − t 0 ℏ + λ ϵ/2[sin(ωt) − sin(ωt 0 )],(B4)where F