The fate of the false vacuum: Finite temperature, entropy and topological phase in quantum simulations of the early universe

Despite being at the heart of the theory of the"Big Bang"and cosmic inflation, the quantum field theory prediction of false vacuum tunneling has not been tested. To address the exponential complexity of the problem, a table-top quantum simulator in the form of an engineered Bose-Einstein condensate (BEC) has been proposed to give dynamical solutions of the quantum field equations. In this paper, we give a numerical feasibility study of the BEC quantum simulator under realistic conditions and temperatures, with an approximate truncated Wigner (tW) phase-space method. We report the observation of false vacuum tunneling in these simulations, and the formation of multiple bubble 'universes' with distinct topological properties. The tunneling gives a transition of the relative phase of coupled Bose fields from a metastable to a stable 'vacuum'. We include finite temperature effects that would be found in a laboratory experiment and also analyze the cut-off dependence of modulational instabilities in Floquet space. Our numerical phase-space model does not use thin-wall approximations, which are inapplicable to cosmologically interesting models. It is expected to give the correct quantum treatment including superpositions and entanglement during dynamics. By analyzing a nonlocal observable called the topological phase entropy (TPE), our simulations provide information about phase structure in the true vacuum. We observe a cooperative effect in which the true vacua bubbles representing distinct universes each have one or the other of two distinct topologies. The TPE initially increases with time, reaching a peak as the multiple universes are formed, and then decreases with time to the phase-ordered vacuum state. This gives a model for the formation of universes with one of two distinct phases, which is a possible solution to the problem of particle-antiparticle asymmetry.


I. INTRODUCTION
The evolution of the early universe described by inflation is now a standard model of cosmic evolution. This theory is widely accepted because of the observational evidence including the cosmic microwave background radiation (CMB) detected in all directions. This evidence points to a "Big Bang" origin as the beginning of the Universe. One of the building blocks is a quantum field theory (QFT) model that explains the origin of the early universe and the observed temperature fluctuations in the CMB. These effects originated in Coleman's theory of quantum tunneling of a scalar quantum field in an initially metastable vacuum [1,2]. The validity of the thin-wall approximations used in this theory and later variations are not verified. Quantum field models are exponentially complex and impossible to solve directly.
It is clearly not possible to repeat the start of the universe, so it has been proposed to verify the solutions experimentally. Such experiments would effectively be an analog quantum computer for the scalar field dynamics of the universe. Here we note that the geometry should be multi-mode and have no boundaries, to allow freespace nucleation. The proposal for a suitable quantum analog computer uses a two-species Bose-Einstein condensate (BEC) experiment in a one-dimensional uniform ring configuration, similar to those studied in several laboratories [3][4][5]. The Bose-Einstein condensate has a modulated coupling between two spin components, which creates a local minimum in the effective phase potential. This allows an experimental study of models of false vac-uum quantum tunneling in relativistic scalar field theories, with an engineered scalar field potential.
It is essential to also model the quantum simulator itself. Firstly, this provides insight into the quantum equations themselves, even if only in an approximation. Secondly, it is necessary to understand the performance of the BEC quantum simulator, and its realization in the laboratory as far as possible. In this paper, we employ numerical phase-space methods to simulate the quantum dynamics of the BEC quantum simulator under the conditions of laboratory experiments for up to 1024 modes. Earlier analyses have verified a false vacuum tunneling, but have not accounted for the effects of finite temperatures and Floquet instabilities in such experiments. Our results show that a momentum cutoff is essential, and that a low initial temperature is required.
A feature of this scalar field model is the existence of a spontaneously broken discrete phase symmetry. This leads to nonlocal topological effects, which are uncovered using phase-unwrapping image processing analysis on the data. There are two distinct types of quantum vacua created with opposite phases. These are locally identical but globally distinct. Similar models have been employed as possible explanations of particle-antiparticle asymmetry [6][7][8]. We show that such global topological effects can be quantified using the concept of an observational phase entropy, which can both increase and decrease with time.
Each true vacuum created from the decay of metastable vacua can expand into a separate "universe" with different topological phases, creating domains of vacuum with boundaries. The result of nucleation includes the tunneling rate, the fluctuations in density and arXiv:2010.08665v1 [quant-ph] 16 Oct 2020 temperature, and the collision of domain walls. All play an important role in the physical nature of the resulting universe formed. However, theoretical work so far relies on a number of assumptions which have not been tested in an experiment. False vacuum decay in zero space dimensions has been recently simulated using a quantum computer [9], but this uses many orders of magnitude fewer qubits than are needed in a spatial model. Here we simulate Coleman's original model, including multi-mode spatial effects, although without the gravitational effects required in a full cosmological theory.
The present paper describes a numerical simulation, whose purpose is to evaluate the feasibility of an experiment and to predict likely outcomes. We use the most successful dynamical phase-space representation for long times, which is the truncated Wigner (tW) approximation to QFT [10,11]. This uses a 1/N expansion for N bosons, and has had previous success in first-principles predictions of quantum dynamics in bosonic quantum fields. It gives correct predictions of tunneling in small bosonic quantum systems with shallow potentials, by comparison to exact number state and positive-P representation methods [12]. Cosmological models also employ relatively flat potentials [13], so this is not unrealistic.
Previous work analyzed quantum bubble nucleation using a 41 K BEC interferometer [14][15][16] at zero temperature. This proposed a condensate containing atoms of the same element with two spin components coherently coupled by a microwave field. The coupled BEC is initialized by a Rabi rotation into a metastable state, which decays into a stable state through quantum tunneling. The transition of the BEC from a metastable to a stable state is an ideal experiment for the investigation of Coleman's original idea. It simulates the relativistic scalar field as a relative phase between the spin components, with the speed of light represented by the speed of sound in the BEC. A related system has been studied using the presence of a seeded vortex within the condensate to initiate false vacuum decay [17].
Here we treat a finite-temperature theory, as will occur in a real experiment. We utilize a variation of the Bogoliubov method [18] for treating the quantum initial state, which employs a nonlinear chemical potential to eliminate divergences in the Bogoliubov theory [19]. We calculate the effects of thermal noise on vacuum tunneling using the truncated Wigner approximation, which has given successful quantum coherence predictions [10,[20][21][22][23][24]. We show that finite temperatures can enhance tunneling, and study how tunneling rates are modified at different initial laboratory temperatures. These results include a momentum cutoff to eliminate Floquet instabilities, and we show that a cutoff is essential by analyzing the effects of changing the cutoff.
We also treat the dynamics and time-evolution of the observable topological phase entropy, which has an intuitive, understandable interpretation. Tunneling events that form the vacuum are disordered, leading to a peak entropy as time evolves. The entropy nearly reaches the maximum possible for an appropriate choice of space and phase bins. As a result, there is a predicted dynamical reduction with time in the topological entropy, as true vacuum domains are formed. The final vacuum state is more ordered, since the false vacuum is unoccupied. Domainwall formation is minimized at low temperatures, which is important for cosmological interpretations [7], since it is known that high-temperature domain-wall formation can lead to anomalous CMB effects.

A. Field equations
Bosonic quantum fields with internal degrees of freedom [25][26][27] are used to describe the Higgs sector in the standard particle model. Global symmetries of the Hamiltonian are broken while creating the low-energy ground state or vacuum. The observation of the Higgs particle makes this an important fundamental concept.
The theory of a metastable or 'false' vacuum was developed by Coleman [1,2]. This used a simpler model, treating the fundamental quantum dynamics of a scalar quantum field with non-derivative self-interactions and a Lagrangian density (using a + − −− metric) of: The local field potential was defined to have two spatially homogeneous, locally stable equilibrium states, φ = φ + and φ = φ − . The first of these has a higher energy, with U (φ + ) > U (φ − ). This is unstable to quantum corrections, and is expected to decay to the true vacuum, φ − . A characteristic predicted to occur in such decays is that a true vacuum is formed by quantum tunneling at a particular space-time point, and subsequently grows at the speed of light.
The dynamics of the evolution of the system is described by Heisenberg field equations of form: The fact that these are operator equations makes them effectively insoluble, apart from approximations. Even if all the eigenstates were known, as in some onedimensional theories, there are exponentially many terms [28] in an expansion of generic initial states. Coleman analyzed a scalar quantum field theory of how such a true vacuum would arise, given an initial metastable state. This model also predicted the formation of individual early universe 'bubbles'. Such theories can be extended to include gravitational effects, and have been used as a theory of the early universe [29]. In these, the scalar field is renamed the inflaton field, and decay to a true vacuum causes an inflationary expansion [13], creating the 'Big Bang'. More recent cosmological studies often focus on post-inflationary events [30], which are less sensitive to quantum fluctuations.
The observation of the Higgs particle and evidence for CMB density fluctuations, confirms the importance of such quantum field theory models. Yet the original false vacuum energies are thought to be many orders of magnitude greater than any possible experiment, possibly ∼ 10 15 GeV [13]. The original event is also hidden from direct observation. In addition to such experimental problems, the quantum field theory itself is exponentially complex, and cannot be solved exactly. As a result, the theory has mainly been analyzed using classical or perturbative approximations [31]. The inclusion of general relativistic effects further complicates the analysis.
It is important to have a better understanding of at least the simplest quantum field theory models. A feature of the model we use is that it possesses a spontaneously broken discrete symmetry, which is known to provide a potential solution to the particle-antiparticle asymmetry problem [7]. Qualitative analysis of domain-wall formation at high initial temperatures has led to objections to this idea [8], related to CMB spectral observations. A full quantum dynamical treatment of the location of domain-walls is needed. Our simulations show that at high temperatures domain walls are prevalent. However, at low temperatures, domain walls are restricted to universe boundaries where they appear less likely to cause inhomogeneities in the CMB spectrum. Observing domain walls in an experimental setting would help to verify or refute this analysis.

B. Approximations and interpretation
Since digital quantum computers are orders of magnitude too small, we propose to solve the equations using an analog quantum computer: a laboratory BEC quantum simulator. In order to obtain insight into the expected performance of the simulator, the dynamics of the BEC system will be solved numerically, but with approximations. Here, we give an outline of those approximations, explaining where they will be expected to hold, and also discuss the interpretation of the simulations.
An interesting consequence of all quantum models for the universe is that the entire universe is described as a single quantum state: there is no external observer. Quantum measurement theory in the conventional Copenhagen model requires an observer to collapse the wave-function. As a result, there are foundational problems in interpreting the wave-function itself. This leads to the question of what can one identify in the simulation that will correspond to a universe?
The numerical simulations of this paper give predictions of the dynamics of the wave function for a model of the universe according to a Hamiltonian treatment. As such the averages taken over the simulations provide the ensemble predictions for the laboratory experiment if re-peated many times. An interesting question is whether or how a particular laboratory realization can relate to a particular dynamical trajectory.
We use a mapping of the wave-function to a Wigner field distribution but with a simpler, approximate timeevolution equation which ensures a positive Wigner distribution throughout the dynamics. In this truncated Wigner (tW) model, stochastic equations are written for complex amplitudes α k that represent modes k. These equations are solved numerically, the quantum noise being modeled stochastically. There is a direct correspondence between the observable experimental moments of the field quadratures and the moments of the real and imaginary parts of α k . The measured quantity of interest is the particle number which, once operator ordering is taken into account, corresponds to |α k | 2 up to an error of order ∼ 1. The values of |α k | 2 encountered in the simulations are macroscopic, and hence the difference between operator and simulation moments is negligible. In this sense, a probabilistic interpretation is possible, where the individual complex amplitudes trajectories correspond to an individual realization. Vacuum fluctuations may be considered as real events, and no additional collapse mechanism is required, as discussed in greater detail in Section (IV).
The dynamics of the tW distribution is approximate. Even though the local evolution errors when using the tW method are of order 1/N where N is the number of particles in each mode, these may grow during time-evolution to create macroscopic errors at later times [12,[32][33][34]. Such errors can increase during quantum tunneling. The tW method cannot describe the formation of macroscopic superposition states [35][36][37][38][39] , or predict certain macroscopic Bell violations [40], because such states cannot be described by a positive Wigner function. Quantum squeezing and entanglement can be described however.
As it is a feasibility study, we do not give exact results, because the quantum simulator experiment is intended to do this. Nevertheless, it is important to ask how reliable the tW phase-space method is. The tW distribution has been used to predict both squeezing and entanglement in systems of large particle number. Comparisons have been made of tW and exact positive-P methods for the dynamics of quantum squeezing in solitons [10,11,41], with excellent agreement between both methods and with experiments [20,21,42]. There is also good experimental agreement for large-scale quantum BEC interferometry, where fringe visibility decoherence times have been accurately predicted. In this regime, quantum entanglement in the form of Schrodinger's quantum steering has been inferred based on the simulations, for states of up to 40,000 atoms [24].
Other studies treated quantum tunneling in driven non-equilibrium systems, which showed agreement between tW and exact methods for shallow tunneling potentials [12,32]. This agreement disappeared for deeper potential wells, which is more likely to lead to macroscopic superposition states requiring negative Wigner distributions. Another quantum field system with metastable behavior is the quantum solitonic breather [28,41], for which tW methods have shown good agreement with conservation laws [23] and both exact positive-P representation and integrable methods [43,44] during the early stages of breather relaxation. While Coleman's original model proposed a deep potential, the models currently favored by many cosmologists do not. Instead, a shallow potential is thought to be more realistic in an inflationary universe [13]. As a result, it is reasonable to use a relatively flat quantum field internal potential, with a large particle number. In this regime, the numerical simulation methods used here appear reliable. Nevertheless, owing to the long time-scales involved -and possible error growth [45,46] -the main goal here is an experiment, regarded as an early universe quantum simulation.

C. Two-species Hamiltonian
For a BEC system having two occupied hyperfine levels with mass m, the Hamiltonian of the coupled field system includes an s-wave scattering potential [47]. It is important in our model that there is a strong mixing between the two spin species, without a phase separation. This requires that the inter-species interaction is minimized, and is assumed to be zero here. There will be losses due to spin-changing inelastic collisions, but these dissipative effects are neglected. The size of such effects is not known for the 41 K Feshbach resonance of interest. This provides a limitation on the accessible tunneling times, since the atoms must tunnel before they are absorbed.
A general two-species Hamiltonian includes both intra and inter-species scattering. The inter-species scattering length is often close to the intra-species one, since differences in nuclear spin orientation do not strongly perturb inter-atomic forces. However, this can change dramatically at a magnetic Feshbach resonance. The required tuning of the s-wave scattering interactions can therefore be achieved with the external magnetic field chosen so that cross-species scattering is suppressed. This is possible in 41 K, as well as in other isotopes like 7 Li.
In these cases, near the Feshbach resonance, one can write the Hamiltonian aŝ Here, writingΨ j ≡Ψ j (x, t) for brevity for the j-th Bose field, the individual spin-species Hamiltonians arê 4) and the microwave coupling Hamiltonian,Ĥ c , iŝ The componentsΨ j are the coupled field operators corresponding to different nuclear spin states, and the subscripts j, k = 1, 2 are the spin indices. These operators have commutation relations Ψ j (x, t) ,Ψ k (x , t) = is a restricted delta function [48] that includes a momentum cutoff, restricting the field to a lattice for numerical simulation.
The coefficient g is the s-wave scattering interacting strength between the atoms, which for a threedimensional system, g 3D , is given by: where a is the s-wave scattering length.
If the atoms are confined by a transverse harmonic trap with frequency ω ⊥ , where the transverse trapping energy ω ⊥ is much higher then the thermal energy, the Bose gas reaches a one-dimensional regime [49]. In this regime, the atoms are confined tightly within an effective s-wave cross-section A s = 2π(l ⊥ ) 2 , where l ⊥ = /mω ⊥ is the transverse harmonic oscillator length [50][51][52]. The coefficient g for a one-dimensional system is hence expressed as g = g 3D /A s , which gives: We neglect microwave spontaneous emission effects, as these are very weak for such microwave transitions.

D. Stephenson-Kapitza pendulum term
The coupling ν describes a microwave field that rotates the nuclear spin by resonantly coupling two hyperfine levels with a frequency separation of Ω HF in the external magnetic field. This is modulated in time in order to induce metastability using Stephenson's concept of a modulated pendulum [53][54][55], later popularized by Kapitza [55]. The depth of the field potential in the metastable state is determined by the dimensionless variable δ.
We work in a rotating frame such that this energy separation is removed from the Hamiltonian, using a different reference energy for each spin component. Here ν (t), the coupling strength between the spin components, is modulated [53][54][55][56][57] as a sinusoidal time-dependent variable with an additional modulation frequency ω, so that ν (t) = ν + δ ω cos ωt. (2.8) Provided the modulation is at a high frequency, this Hamiltonian is equivalent to the Coleman model of a relativistic scalar quantum field with an engineered quartic potential. Here the phonon velocity corresponds to the speed of light [14,16]. Modulation amplitudes are relatively large, so that the excitation is nearly bichromatic [58,59], or double-sideband.
The result of including this term is an engineered potential U (φ a ) for an effective scalar field φ a which physically is the phase difference between the two co-existing Bose-Einstein condensates with phases φ j . We define the atomic phase difference as φ a = φ 1 − φ 2 − π, so that the false vacuum is at φ a = 0 and the true vacua are at φ a = ±π. In the limit of strong particle-particle repulsion, this obeys the relativistic field equation (2.2), where the speed of light, c, is replaced by the phonon velocity, c = gρ 0 /m in the BEC.
The potential equation is given by [16]: where ω 0 = 2 √ νgρ 0 / and λ = δ 2gρ 0 /ν. This potential can be varied by the experimentalist, and is best understood in a dimensionless form described later. It is known that instabilities can form in this model at too low a modulation frequency. This is due to coupling between effects caused by the modulation frequency and high-frequency phonon modes [60,61]. Such effects therefore require the use of high enough modulation frequencies to move the instability region above any physical cutoff in momentum.
We will assume that there is a physical mechanism to remove high-momentum phonon modes. An example of this would be the use of a spatially modulated potential to introduce a band-gap. A second possibility is the use of a swept modulation frequency to reduce parametric gain by changing the unstable momenta. Ultimately, as pointed out by earlier workers, the inverse scattering length provides an intrinsic cutoff, ultimately at 1/a.
An experimental mechanism to achieve this is via an optically modulated trap potential. Instabilities were not observed in our previous numerical simulations, due to the use of a finite lattice that includes a momentum cutoff. An analysis of modulational instabilities in experiments with larger numbers of modes and higher phonon momenta is given in the Appendix, where typical parameters are presented.
We conclude that such instabilities are generally present, but can be suppressed in the proposed experiment by using high enough modulation frequencies combined with a momentum cutoff, as shown in Fig (1).

III. INITIAL STATE AT FINITE TEMPERATURE
Our initial state includes quantum and thermal fluctuations. We combine the Bogoliubov method [18] with the Wigner representation, so that both the initial thermal excitations and vacuum noise are taken into account. In this approximation, the system is assumed to have a macroscopic condensate mean population N c with density ρ c for one of the spin indices, say j = 1. Experimen- tally, this is produced using evaporative cooling methods [62][63][64].
The ground state field operatorΨ 1 is nearly equal to the square root of the condensate density, √ ρ c , which is assumed constant over the ring. This is valid for typical ultra-cold BEC experiments below threshold, provided trapping potential noise is small. There are zero momentum divergences in applying the Bogoliubov expansion to finite systems, which are removed here using a nonlinear chemical potential [19].
Initially, the microwave coupling ν is turned off, and the second spin species operatorΨ 2 is in the vacuum state. The BEC is prepared in a single species condensate with one spin component populated at a temperature T . While it is possible to achieve temperatures well below the condensate critical temperature, these are still not at zero temperature. The condensate has thermal phonon excitations, which can change the tunneling time. These also induce phase fluctuations and finite temperatures in the effective scalar field. This is expected to give modified tunneling compared to a metastable quantum field without extra noise. Since the exact quantum state prior to the Big Bang is not known precisely, our goal here is to determine the effect of thermal noise on a laboratory experiment. Even this may not capture the full effects of evaporative cooling, which is a complex dynamical process [65].

A. Nonlinear chemical potential
To model a finite temperature experiment, we assume that the initial density matrix is in a grand canonical ensembleρ GC at temperature T for j = 1, and a vacuum state for j = 2, so that: HereK is the 'Kamiltonian' , which includes a chemical potential to give a finite particle number in the thermal state, so thatK where µ N is the chemical potential for an initial pop-ulationN in one of the spin configurations. This can be any nonlinear function ofN [19], so that, including terms up to second order The chemical potential has no effect on the dynamics if ν = 0, since Ĥ , µ N = 0. Such a nonlinear chemical potential is useful for describing thermal fluctuations in a BEC. The utility of this approach is that on linearizing the total Kamiltonian,K, the zero frequency divergence of the Bogoliubov expansion [18] is eliminated with an appropriate choice of the quadratic coefficient µ 2 .
To allow an expansion around a well-defined coherent state phase, we suppose thatρ is an ensemble average of condensatesρ (φ c ) with a coherent phase φ c , so that: This ensemble corresponds to a number-averaged ensemble of states with Poissonian number fluctuations around N c , which is more realistic than using a fixed particle number. The actual number fluctuations in experiments are often super-Poissonian, and these additional number fluctuations can be included if required. For number-conserving interferometric measurements of relative phase, the initial coherent phase φ c is irrelevant. All phase choices give identical observables. It is therefore sufficient to include a single member of the phase ensemble with φ c = 0, following standard methods from laser physics [66][67][68]. There are other techniques for obtaining convergent Bogoliubov expansions [69,70]. These generally involve operator expansions having nonstandard commutators, and do not readily permit Wigner phase-space expansions.

B. Regularized Bogoliubov theory
The grand canonical Hamiltonian,K, is obtained from a Bogoliubov expansion [18,71,72] of the single-species field operator at time t = 0, assuming a condensate phase of φ 1 = 0. This is given bŷ Here ψ c = √ ρ c , which is the initial condensate density, and the field fluctuations are expanded as a sum over phonon momenta k: For a complete unitary transformation, all modes must be included, including the zero-momentum mode. Next, we introduce the condensate quadrature operatorP aŝ The number operatorN can be written to second order in the quantum field fluctuations, givinĝ Expanding the grand canonical Hamiltonian in δΨ 1 , the choice of µ 1 = 0 and µ 2 = g/L eliminates all first order terms as well as terms inP 2 , which would cause phase divergences and unphysical phase diffusion in equilibrium. The resulting convergent Bogoliubov expansion has u 0 = 1 and v 0 = 0 for the zero-momentum terms with k = 0, rather than the divergent expression found in the standard expansion.
The mode coefficients of the single-species field for k = 0 are expressed in terms of the excitation energy k and the free-particle energy E k , as where E k = 2 k 2 / (2m) is the free-particle energy and is the excitation energy. Due to periodic boundary conditions on a ring trap, the allowed values of momenta are for M momentum modes, assuming an odd mode number. The resulting effective Kamiltonian describing thermal excitations of the BEC is given bŷ Here, the phonon operatorsb † k ,b k describe the creation and annihilation of a quasi-particle in an excited state k. The resulting excitation in each mode with k = 0 is a propagating phonon. This expansion is particularly useful for the finite temperature system we are interested in here. Phonon excitations with k = 0 are populated according to the bosonic thermal distribution, For k = 0, we assume a vacuum state with n 0 = 0, since this operator cannot be coupled to an energy exchange process, owing to number conservation. This choice gives Poissonian number fluctuations. As discussed above it may be necessary in modeling experiments to include larger number fluctuations due to technical noise occurring in the evaporative cooling process.
After the initial preparation, a microwave pulse is used to rotate the Bose gas occupations so that the two spin species have an equal occupation, which is equivalent to a linear beam-splitter. We denoteΨ j (x, 0) as the initial quantum fields after the BEC is split into the two states. The two initial spin states after rotation have equal density ρ 0 = ρ c /2 and a relative phase of π.
This corresponds to a rotation matrix acting on the quantum fieldsΨ 1 andΨ 2 : where θ = π/2 and φ = −π/2. The system is then assumed to evolve according to the Hamiltonian (2.3) with a c.w. microwave coupling field present. This carrier has an appropriate phase relationship with the microwave pulse field so that the quantum system is initially in the metastable high energy state.

C. Dimensionless parameters
The equation of motion and quantum operators can be transformed into dimensionless form by introducing dimensionless time, distance, and frequency: The speed of sound in a weakly interacting BEC is given by c = gρ 0 /m, and the initial temperature of quantum degeneracy is T d [73]. We define a characteristic length, temperature and frequency as: The field amplitude in dimensionless coordinates is Ψ j =Ψ j √ x 0 , and the density in dimensionless form is given by ρ 0 = ρ 0 x 0 . This gives a characteristic energy scale that is the geometric mean of the mean field energy gρ 0 and coupling energy ν, as used previously. The six dimensionless parameters that define the physical system are therefore: where L is the dimensionless length of the simulation, T is the temperature of the initial BEC field with density ρ c = 2ρ 0 , and λ is the effective depth of the modulation. The corresponding dimensionless Hamiltonian, including the chemical potential, is: , the effective nonlinearity, depends on the other parameters. The chemical potential has no effect on dynamics, and is only required to remove singularities in the linearization of the initial state.

IV. PHASE-SPACE REPRESENTATIONS AND ENTROPY
To investigate vacuum nucleation at finite temperature, we take both quantum and thermal fluctuations into account. This is achieved by performing stochastic numerical simulations in the Wigner representation of the full quantum model of the two-component BEC system. These numerical simulations are not exact, and indeed there is no exact method known. Such large-scale quantum field calculations are exponentially complex. This lack of an exact solution is the motivation for a quantum simulation experiment. However, we can use stochastic methods to investigate realistic conditions and expected results.
The effect of quantum and thermal fluctuations results in a wide range of nucleation times. Numerical methods which are not limited to short simulation time are desirable to simulate the Coleman theory. Here we choose the truncated Wigner (tW) approach [10,11,74]. It has a sampling error that remains well-controlled over a long simulation times. This method gives the first quantum corrections in an M/N expansion [75], were M is the number of modes and N is the number of atoms. To give reliable results, it is therefore necessary that M/N << 1. This has been confirmed in comparisons with exact positive-P and complex-P simulations that do not use truncations [10,43].
There have been successful predictions of measured quantum squeezing in solitons [10,20,21,42,76], propagation effects in BEC lattices [77] and fringe visibility in finite temperature BEC interferometers [11,24]. These also indicate that the technique is reliable. Hence the tW method should be able to simulate these experiments given a large number of atoms. Truncation errors can build up at long times [33]. This means that the numerical simulations are expected to be correct if tunneling is not too slow. The ultimate goal is an experiment.
Tunneling phenomena are a stringent test of numerical quantum field simulations. Earlier comparisons of truncated Wigner quantum tunneling with exact methods have shown that it is correct for relatively shallow potentials [12,32]. This is also the regime of most interest in many cosmological models. As a result, we expect that this method will give good indications of the effects of thermal and quantum noise in the regimes of most interest.
Using this approach, the quantum state is represented by a stochastic phase space distribution of trajectories following the Gross-Pitaevskii equation. In a thermal state, the Wigner representation of the initial state has a complex Gaussian distribution. We perform our simulations for a one-dimensional system, whose equation of motion are obtained from (2.4). A typical example is shown in Fig (1).

A. Wigner representation
We require that the dynamics of the system is evolved quantum dynamically, hence quantum fluctuations must be taken into account. To do this, we transform the phonon operators using the Wigner representation correspondence [10,11,23,48]. Since the initial state is approximately Gaussian before phase-averaging, the corresponding initial Wigner representation is also Gaussian [77,78].
Each initial phonon mode is represented as a complex Gaussian variable, i.e.b k ∼ β k andb † k ∼ β * k . Thermal fluctuations are included in the modes with k = 0. A detailed explanation of how this is obtained using the nonlinear chemical potential method is explained elsewhere [19]. The Wigner representation of the initial quantum density matrixρ (0), after evaporative cooling, is a com-plex Gaussian distribution given by: (4.1) Here, W 1 [ψ 1 ] is a representation of a thermal state with finite temperature, and W 2 [ψ 2 ] is a vacuum state. These are positive distributions that can be sampled probabilistically, with samples given after evaporative cooling by: with α k and β k defined as independent complex Gaussian random variables in each vacuum mode and phonon mode respectively. These have mean values such that α 2 k = α k = 0, β 2 k = β k = 0. The only non-vanishing moments are

B. Reduction to dimensionless parameters
In dimensionless form, Eq (4.2) for the condensate after cooling and before rotation, is given by where | ψ c | 2 = ρ c = 2 ρ 0 is the mean field density of the single-species condensate before rotation. The first term is an inverse Fourier transform of the collective excitations in Wigner representation. The resulting values for u, v are where k = k / ω 0 is the Bogoliubov excitation energy in dimensionless form, so that, in dimensionless units:

C. Metastable state generation and detection
Assuming that the thermal phonons effectively behave as a canonical ensemble of free bosons, the thermal fluctuations for k = 0 are represented by the complex Wigner amplitude where η i, k is a complex Gaussian noise in dimensionless space, with η i, k η * j, k = δ ij δ k k . As a result, To create the metastable state described in Coleman theory in our proposed experiment, a radio frequency field with shifted phase of π/2 is applied to the singlespecies BEC. This prepares a superposition of initial states |1 and |2 , where the two-component condensate corresponds to the initial metastable state, together with finite temperature thermal fluctuations.
We denote ψ 1,0 and ψ 2,0 as the initial Wigner fields after the BEC is Rabi-rotated into the the two hyper-fine levels. These initial states are required to have equal density ρ 0 with a relative phase of π, which corresponds to a rotation matrix identical to that used for the Heisenberg fields, but now applied to the Wigner fields ψ 1 and ψ 2 : where θ = π/2 and φ = −π/2. A sample dynamical trajectory in the Wigner phasespace representation satisfies the equation: (4.10) Here we ignore the chemical potential term, which is identical for both components and has no effect on the relative phase dynamics. Transforming (4.10) into this dimensionless form, the time-evolution of the Wigner field trajectory is given by [16]: Increasing the modulation so that λ > 1 gives a local minimum in the corresponding effective potential. The corresponding dimensionless effective potential in the phase-difference φ a is (4.12) Using the fields ψ 1 and ψ 2 as initial conditions, one can propagate the Wigner fields in real time, using the equation of motion (4.11). The atomic relative phase then evolves approximately according to the relativistic field equation (2.2), so that: (4.13) The relative phase of the two spin components is dynamically evolved, which includes quantum tunneling effects.
Since we wish to evaluate the laboratory experiment, the full atomic equations are evolved, rather than just the reduced phase equations. We detect vacuum formation at a finite time by rotating to measure the relative phase from the resulting hyperfine populations, and the results are compared at different temperatures. A typical single trajectory is shown in Fig (1). This shows tunneling events occurring at isolated space-time points. As expected, the true vacuum regions grow at the speed of light (c = 1). Full details are given in Section (V).

D. Topological entropy
Entropy is an important phenomenon in all physical systems. It is one of the foundations of thermodynamics, and can be interpreted as a measure of disorder or randomness of a system. For a quantum system undergoing unitary evolution, such as the entire universe in this model, the von Neumann entropy is invariant. While von Neumann entropy can increase when the contributions for different entangled parts of the universe are summed, the overall von Neumann entropy for the universe is static. Alternative definitions of entropy have emerged that are based on measurement properties [79][80][81][82][83]. This leads to the question of which entropy measure can be used to quantify the disorder of an early universe simulator, and how it can be calculated and measured.
Here we investigate the disorder of the simulated universe using an observational macroscopic entropy that can be calculated and measured. This is based on a wellknown quantum entropy measure, the Wehrl entropy [84], where Q(α) = α|ρ |α /π M is the Husimi function [85]. The Q function is a positive, probabilistic representation, defined for all quantum states. It can be used to link the simulations and an interpretation of the quantum universe, based on the Q function [86]. In this interpretation, the universe simply corresponds to a particular sample of a Q-function probability. It is nontrivial to simulate the Q-function dynamics, since it does not satisfy a Fokker-Planck equation. Consequently, rather than solving for the Q-function directly -which would be equivalent to quantum field dynamics -we have utilize a closely related method, the truncated Wigner (tW) approximation [10,11,74]. This has much simpler dynamical equations.
Measurements that correspond to a Q-function trajectory are anti-normally ordered, so that averages over the symmetrically-ordered tW simulations do not directly correspond to those of a Q-function trajectory. These two distributions corresponding to two different representations, one approximate, the other accurate. In any experiment on a large BEC, the difference between the truncated Wigner and the Q distribution is microscopic, and has a negligible effect on macroscopic observables like the average phase difference. Since the ordering introduces differences of only a microscopic order, either can have the interpretation of macroscopic reality [87]. The retrocausal nature of the individual trajectories and their relationship to quantum measurement theory is treated in detail elsewhere [86]. In this interpretation, no additional collapse mechanism is required.
If we consider a mesoscopic observable, one can approximate the Q-function distribution by the Wigner function used here, since the two are related by a microscopic convolution of order [88]: However, even the Wehrl entropy measure, though simpler than the von Neumann entropy, is not readily measurable in a multi-mode system, as it requires an exponentially complex set of measurements. To resolve this problem, the idea of an observational entropy has been recently put forward, which uses a finite set of measurements to define entropy [89]. We use an observational version of the Wehrl entropy, which is a combination of the Wehrl and observational entropies. Each amplitude α, describing a possible universe in the Wigner or Q-representation, is reduced to a phase, measured and binned into a set S i which classifies phases into p distinct ranges in each of contiguous regions. This is a topological measurement, since the phase can only be established through a nonlocal phaseunwrapping algorithm [90], which allows one to distinguish topological phases of −π and π through continuity in space and time.
Given n i as the number of measured universes in the i th bin, from a total of n = p , we define a probability P i = n i /n, and a corresponding observational Wehrl entropy as (4.16) An early universe simulation must have a finite number of trajectory samples. Hence, we require an appropriate binning strategy to formulate sample probabilities, in which the total number of bins should be less than the number of samples, to give reliable estimates. The simplest strategy would be to use a binary binning with the relative phase at each point in space in either a false vacuum or in a true vacuum, thus ignoring topological effects. However, as shown in Fig (1), if we start with a false vacuum at φ = 0, we find two topologically distinct true vacuum states with relative phases of −π and π. These are distinguishable using nonlocal phase unwrapping methods in space or time. As a result, we require 3 phase bins in each spatial region, of −π ± π/2, 0 ± π/2 and π ± π/2 to capture the vacuum states in our early universe model. Since phase is measured in finite volume, we define it by averaging over a range of neighboring spatial lattice points. Therefore we identify in each space-time interval three phase bins, two for the true vacua, and one for the false vacuum. Compared to entropy as a microscopic quantity, this entropic measure is uniquely sensitive to macroscopic topological disorder. In fact, it is sensitive to disorder on the scale of the entire universe, or model universes in the case of the proposed laboratory quantum simulations using coupled Bose condensates. This has quite different properties to the microscopic von Neumann entropy.
Our simulations show a cooperative effect, where each vacuum bubble eventually becomes dominated by one or other of the topological phases. This provides a model for multiple universes with fundamentally different properties. An intriguing property of this type of scalar field symmetry breaking is that it is purely topological. There are no local measures that distinguish the different topological phases, although they are distinguishable using phase unwrapping.
It is speculated that discrete symmetry breaking could provide a mechanism for matter-antimatter asymmetry [7,8]. The basis for this is that scalar field behavior involves very high energies, with matter and antimatter being formed at much lower energies. As a result, small asymmetries could have a large influence at the lower energies of matter formation. The validity of this model hinges on the question of whether domain walls form, as these can alter the observed CMB spectrum. Our simulations indicate that domain wall formation is suppressed at low temperatures and confined to universe boundaries. Experimental evidence of domain wall formation and symmetry breaking can be found by measuring the topological entropy.

V. NUMERICAL RESULTS
This section will summarize the results of numerical studies of the effect of finite initial temperatures in the proposed BEC experiments. Our simulations use a dis-crete lattice, corresponding to a physical momentum cutoff at M = 256 modes. This is necessary to prevent modulational instabilities. The effect of removing the momentum cutoff by using a smaller lattice spacing with more modes is reported in the Appendix. It is experimentally challenging to measure relative phase unambiguously in a BEC, as this requires a simultaneous measurement of two complementary quadratures. Most of the numerical results presented here use the more accessible measure of relative number distribution, p z ∝ cos (φ a ), while we also present results for the relative phase φ a in the section on topological phase entropy.

A. Experimental parameters
To have a realistic numerical study, we have chosen possible parameters that correspond to a one-dimensional system of 41 K atoms in a ring trap near a Feshbach resonance. There are many choices of atomic species possible, including 7 Li, so this is only one scenario among many.
For the existence of quasi-particles in a finite temperature condensate, the circumference L of the trap should be shorter than the temperature-dependent phase coherence length l φ ≈ 2 ρo mk B T [73]. The restriction on the trap circumference L l φ limits the temperature T of the condensate in the ring trap, i.e. T T c = 2 ρc mk B L [16]. For condensates in a ring trap with a three-dimensional density ρ c,3D , assuming the condensate atoms are transversely confined in the effective s-wave scattering crosssection A s , the corresponding one-dimensional density is estimated to be ρ c = N c /L ≈ A s ρ c,3D . Following the suggested parameters in [15,16], the parameters of the proposed experiments are listed in Table I.
The partial differential equations (4.11) are solved using an interaction picture fourth-order Runge-Kutta (RK4) method with the extensible open source MATLAB package xSPDE [91]. From the experimental parameters listed in I, the corresponding typical dimensionless parameters used in the numerical simulations are listed in Table II.

B. Observational criteria
In order to convert the relative phase of the two species into number density distribution, a π/2 radio frequency pulse can experimentally be applied to the coupled fields. Vacuum nucleation can be observed from the relative number density distribution, where ρ 1 ( x) and ρ 2 ( x) are the number density of the two species respectively, after applying the second Rabi rotation.    Figure 2 shows a single trajectory example of onedimensional false vacuum dynamics. The simulation starts with thermal states of a two-component condensate at a low reduced temperature of τ = 1 × 10 −5 . The coupled field system is in the metastable state initially, with p z ∼ 1 at time t = 0, indicated by the yellow contour. The system starts to decay into a stable true vacuum state with p z ∼ −1, indicated by the blue contour at times t 2. In this example, five bubbles are formed of true vacua. These bubbles expand until they either meet at continuous domain walls of false vacuum (at x ≈ −43 and x ≈ 25), which correspond to topologically distinct phases, or else form localized oscillons (at x ≈ −30, 6 and 46).
From the simulations of our model at finite temperature, the thermal energy introduces extra thermal fluc-  As can be seen in Figure 3, no stable domain walls or periodic oscillons are formed as the temperature increases, and no true vacuum bubbles survive on long time scales. As quantum tunneling is enhanced at higher temperatures, more bubbles are formed in the true vacua, but most of them are short-lived. Tunneling is accelerated at high temperatures, which leads to strong fluctuations between the true vacua and the false vacua, and many relatively unstable domain walls.
To quantify tunneling events, we can examine the average relative phase of the coupled fields along the axial coordinate, where As illustrated in Figure 4, the relative average phase cosφ a = 1 corresponds to an initial false vacuum. As the tunneling starts and the false vacuum decays to the true vacua, the average relative phase is expected to gradually decrease from 1 to −1 in a complete transition. At very low temperatures, the presence of the true vacuum bubbles is noticeable with cosφ a < −0.5. However, at higher temperatures, the presence of the true vacuum bubbles is less noticeable due to the influence of the thermal fluctuations, and cosφ a only goes to just below 0. We define a threshold value cosφ a = 0.9 as the initiation of the false vacuum tunneling event.  Recalling that the tW method provides quantum estimation from a set of stochastic trajectories, hence to investigate the thermal effect on the tunneling rate at finite temperature conditions, we repeat the single trajectory simulation and determine the probability of bubble creation ( cosφ a ≤ 0.9) over time P ( t).
We then obtain the survival probability of the false vacuum F given that F = 1 −´ t 0 P ( t )d t . The tunneling rate Γ for long time scales is calculated from a linear fit in log scale given that F = exp(−Γ t) [92]. We performed simulations over a range of coupling strengths ν and reduced temperatures τ . Results are presented in Figure 5, 6 and 7. Previous work showed that larger couplings ν extend the tunneling time [15,16], which confirmed an expected slowing-down of quantum tunneling with dissipation [93]. To illustrate thermal effects, we compare the Bogoliubov thermal state results with the coherent state results which thermal noise is neglected. Although the coherent initial state is not a true Bogoliubov ground state, it has very similar behavior to the low temperature ground state. In the Wigner representation, the coherent initial state is represented by adding the quantum noise in each vacuum mode (4.4) to the classical false vacuum state. In dimensionless from, the coherent Wigner fields after the BEC is Rabi-rotated are Gaussian random variable as already described in (4.7).

C. Thermally induced changes in decay rates
How does the finite initial temperature affect the decay of the false vacuum? From all three figures ( Figure  5 to 7), the tunneling rates Γ determined from both the coherent state and the thermal states at all tested temperatures τ show a power law dependence on ν. The gradients of log Γ for each effective modulation depth λ are similar. If we compare the change of the tunneling rate Γ at a fixed modulation depth λ, from each of Figure 5, 6 and 7, one can see that the rate of tunneling is increased as the temperature τ increases. For the case of very low temperatures, τ = 1 × 10 −5 , as the modulation depth λ increases the tunneling rate is reduced more significantly than at higher temperatures. On the other hand, in the case of the highest reduced temperature studied (τ = 3 × 10 −4 ), this reduction of the tunneling rate due to an increase in λ is less significant  than at lower temperatures. The tunneling of the false vacuum is less restricted by λ at higher temperatures. This result suggests that at high temperatures, the effect of thermal fluctuations dominates the quantum decay of the false vacuum. For a fixed coupling strength ν, an increase of temperature increases the probability of penetrating the modulation depth barrier λ, and hence increases the formation rate of the true vacuum. The correlated ground state with τ = 10 −5 has a slightly lower tunneling rate than a coherent state, as it is stabilized by the Bogoliubov correlations. The qualitative effects of bubble nucleation and growth are similar in both cases.

D. Topological entropy results
The topological entropy allows us to quantify the phase disorder caused by the formation of true vacua and domain walls. Results for the topological entropy with 10 4 distinct realizations are shown in Fig (8). Initially the coarse grained observational entropy is nearly zero, since all phase-space coordinates are in the metastable false vacuum. The observational entropy initially increases with time evolution, in contrast to the von Neumann quantum entropy which is constant with time. The entropy reaches a maximum as tunneling occurs, giving a nearly maximally disordered state with entropy S < 8 ln 3 = 8.79, using = 8, and p = 3. The entropy then reduces as the true vacua grow, reducing phase disorder. The entropy function is shown to reach a maximum as a result of tunneling and after t 40, it stabilizes to near S ∼ 8 ln 2, as the false vacuum is eliminated. Between τ = 10 −6 and τ = 10 −4 shown by black and red curves, there is a decrease in the maximal entropy S max and lifetime of the peak. From τ = 10 −5 to 10 −4 , as S max decrease the lifetime of the peak increase. Beyond τ = 10 −4 , the system thermalizes, as shown by the green curve.
At very low temperatures, the final state has almost no false vacuum. As a result, S ≤ 5.55 as seen in the simulations of Figure 8. At higher temperatures the false vacuum never disappears completely, due to unstable domain wall formation, thus increasing the final state disorder. This is shown by the red and black scattered lines in Figure 8 and leads to the prediction of stable steadystate entropy values S ∼ 5.3. The remaining disorder arises from the randomness due to the two topological phases of the true vacua and the remaining false vacuum domain wall or oscillating oscillons. The amount of disorder possible is limited by the spatial extent chosen for phase bins, and the number of samples used. showing two different colors that depict the two different topological phases for the true vacuum. The true vacuum bubbles expand in space until they meet after a time t 15 at x ≈ −19 and after a time t 24 at x ≈ 46. The neighboring true vacuum regions with distinct topological phases are separated by a domain wall of false vacuum. On long time scales, Figure 1 clearly depicts the formation of an asymmetric structure of the topological phase in the true vacuum.
In Figure 9 we display the evolution of the entropy function S( t) for different values of ν using 10 4 tW trajectories, and a comparison plot for the evolution of average relative phase cos(φ a ) of the corresponding single trajectory examples. For larger couplings ν, the maximum entropy S max is reduced and the lifetime of the peak is extended. As the coupling ν is increased, the tunneling initiation is generally delayed and the tunneling time is extended.
The single trajectory results presented in Figure 9 (b) are consistent with the expectation of a slowing-down of tunneling. We find similar behavior for the modulation depth λ in the regime of metastability (for λ > 1). In the case of large λ, bubble nucleation is delayed, with more results given in the Appendix.
The time evolution of the entropy function S( t) for different numbers of spatial modes M using 10 4 trajectories is shown in Figure 10. The results are compared with the corresponding average relative phase cos(φ a ) shown in Figure 10 (b). The evolution of S( t) and cos(φ a ) for higher spatial mode numbers are specified by black and red scattered lines in Figure 10. Systems with higher mode numbers have a more chaotic phase fluctuation, whereas the system with lower mode numbers specified by the blue dashed line show a smooth tunneling to the true vacuum.
The behavior is reflected in the entropy plots in Figure 10, where the narrow peaks of the entropy function for black and red scattered lines indicate a more chaotic phase ordering. On the other hand, the flattened peak of the entropy function for blue scattered line shows a relatively stable phase. In our simulations, the rapid transition to chaotic fluctuations at higher spatial modes can be minimized by increasing the characteristic frequency to ω = 200 to achieve a stable tunneling to true vacuum. This behavior is caused by Floquet mode instabilities, and is an artifact of the microwave modulation frequency, as explained in greater detail in the Appendix.

VI. CONCLUSION
Prior to the present paper, the proposed experiment on the decay of a false vacuum had only been studied previously in the zero-temperature limit [14,16,60]. As a result the effects of finite temperature thermal noise on vacuum tunneling and nucleation using BECs with two spin components were not understood. This mixture of two BEC spin components can be used as the relativistic analogous quantum field, where the relative phase of the two species corresponds to the metastable state (phase π) and stable state (phase zero) respectively. The components can be coupled via a microwave field, which creates an unstable vacuum.
To create a metastable vacuum from this unstable vacuum, one can make use of the classical concept of a modulated pendulum [53][54][55][56], where a modulated amplitude of microwave coupling allows one to engineer the metastable vacuum potential. We consider a symmetric intra-component case in this work for its simplicity. In previous work, we suggested the use of a pair of Zeeman states of 41 K where |1 = |F = 1, m F = 1 and |2 = |F = 1, m F = 0 as the two spin states in the proposed BEC experiment. This Zeeman pair has a symmetrical intra-component Feshbach resonance where the inter-state scattering lengths between the two components are zero [94].
We have shown in this paper that quantum vacuum nucleation is accelerated at finite temperature, although the bubbles formed of true vacua may be short-lived at high temperatures due to thermalization of the BEC. The formation of the true vacuum bubbles at finite temperature generally follows the expected behavior, apart from an accelerated tunneling rate. Clearly, lower temperatures move one further into the true quantum regime.
Our results also show that higher oscillator frequencies ω can remove short-wavelength instabilities, provided that there is a high-momentum cutoff present. This increases the feasibility of a false vacuum BEC experiment. The proposed table-top experiment using Zeeman states of 41 K may help to test vacuum tunneling theory under real experimental conditions where thermal effects are unavoidable. This provides insight into early universe models with high temperatures, which is present in some early universe models [7].
Such an experiment can detect a unique topological feature of the scalar field vacuum. Unlike the von Neumann quantum entropy, which is time-invariant, the topological phase entropy is predicted to reach a maximum at the time when most tunneling occurs. Intriguingly, phase entropy can decrease with time, because the final vacuum state is more ordered than the state occurring while tunneling. This is a true topological feature of the present model, since the relative phase is only uniquely defined after a nonlocal unwrapping.
Such measurements may be useful in investigating proposals in which a discretely broken symmetry provides a model for particle-antiparticle symmetry. We find that domain wall formation is prevalent at high temperatures, in agreement with qualitative predictions [7,8]. However, at low temperatures domain walls are restricted to universe boundaries where they would be far removed from having a direct influence on CMB inhomogeneities, which was thought to be a possible problem with such theories.
Alternative implementations include a homogeneous, two dimensional simulation. This allows even more complex topological vacuum structures to form. Such experiments are realizable in micro-gravity with a shell geometry [95,96], for example in the NASA CAL space-station environment. There are many other possible approaches using different quantum technologies including superconducting circuits or discrete Bose-Hubbard lattices, which are not analyzed here.

ACKNOWLEDGMENTS
The authors acknowledge helpful discussions with Andrei Sidorov. This research has been supported by the Australian Research Council Discovery Project Grants schemes under Grants DP180102470 and DP190101480 .

APPENDIX: MODULATIONAL INSTABILITIES
Modulational instabilities can occur if the microwave modulation frequency ω is too low relative to the momentum cutoff. In such a case, the microwave modulation cannot be adiabatically eliminated, and parametric instabilities occur. The effect of such Floquet modes was studied by Braden et al. [60,61]. In the work above, we have included a momentum cutoff to prevent this, as explained in the main text.
In this Appendix, we analyze the effects of modulational instabilities at finite temperatures with a higher momentum cutoff. The value of cos (φ a ) of the relative phase of such a system evolves to around ∼ 0, indicating that the transition from the false vacuum to a true vacuum is inhibited, as the relative phase becomes completely randomized. This is studied numerically by choosing a smaller lattice spacing and lartger M , to resolve the high wave-number unstable modes.
The unstable modes occur in a narrow band centered at wave-numbers as derived in [60,61]. The dimensionless critical unstable wavenumber is given by: where σ = cosφ a = ±1 is the relative phase of the fields.  Figure 11 shows the cutoff wave-numbers k nyq in the simulations of previous sections (V B) and (V C) in comparison with the corresponding critical wave-numbers k c at various values of ν. The Nyquist wavenumber determined by the lattice spacing ∆ x, is k nyq = π/∆ x. The unstable wavenumber k c is the wavenumber of the highest gain Floquet mode.
As shown in Figure 11, such effects are excluded in the M = 256 simulations presented in sections (V B) and (V C) at all values of ν. Therefore, the short-wavelength fluctuations between the true vacua and the false vacua found at higher temperature ( Figure 3) are purely due to the thermalization of the condensate. Thermal effects were excluded in earlier work [61]. The dynamics at finite temperatures presented in Figure 3 shows similarities to the dynamics in the presence of the Floquet modes in [61]. In the following we will investigate the combined effect of thermal fluctuations and Floquet instabilities on the dynamics, by setting k nyq above the critical k c . Figure 12 shows our results when both thermal effects and unstable Floquet modes are included. Here we reduce the lattice spacing ∆ x in the simulations by increasing the number of simulation modes to M = 1024. This lattice spacing corresponds to a Nyquist wavenumber of k nyq ≈ 32.14, well above the critical wavenumber k c ≈ 15.39 for ν = 7 × 10 −3 .
The most significant effect of the unstable Floquet modes is that true vacua formed at finite temperature are gradually destroyed, and the system is eventually dominated by chaotic fluctuations. At a reduced temperature of τ = 1 × 10 −4 , the presence of Floquet modes causes fluctuations with wavelengths shorter than the dominant thermal fluctuations.    In Figure 2, 3 and 4 we show the effect of thermal fluctuations independently by setting k nyq < k c to exclude Floquet modes. The true vacua in the absence of the Floquet modes have a metastable structure with average relative phase cosφ a acquiring a non-zero constant value ( Figure 4). However, this behavior is different in the presence of the Floquet modes, where the structure of the true vacua is short-lived.
In the example presented in Figure 12, most of the vacuum bubbles only survive over a time duration 5 t 30 (which corresponds to a real experimental time duration ∼ 20ms). At t 30, one can expect the average relative phase cosφ a to be around ∼ 0 as the chaotic fluctuations dominate. This is confirmed by the averaged result using 8000 trajectories shown in Figure 13.  Figure  12. The errors in cosφ a are less than 1%.
The modulation depth λ is known to correspond to the strength of the unstable Floquet modes [60]. In the absence of thermal effects, earlier studies [61] showed that increasing the modulation depth λ can reduce the timescale of the Floquet modes, and results in the stabilization of long-wavelength structure in the decay of false vacuum. This effectively delays the nucleation of the true vacuum bubbles. For systems at finite temperature where short-wavelength thermal fluctuations coexist with the true vacua, we found that this delay of bubble nucleation is valid. Figure 14 shows the examples on the effect of increasing λ at a fixed reduced temperature τ = 1×10 −4 .
At time t 10, our results show that the systems with larger λ (λ = 1.3, 1.5) break through the threshold value of the tunneling initiation cosφ a = 0.9 later than the λ = 1.1 system, indicating a delay of the bubble nucleation. In the presence of Floquet instabilities, increasing λ results in a significant reduction of the average relative phase. The peak value of | cosφ a | drops from ∼ 0.35 to ∼ 0.1 as λ is increased from 1.1 to 1.5. For decay at finite temperature, we have shown in the main text Figure 4 that this phase signal is also weakened by the influence of thermal fluctuations in the absence of the Floquet instabilities.
Increasing λ in the presence of Floquet modes further reduces the phase signal, which increases the difficulty of measuring true vacua in an experiment. Therefore, it is expected that true vacuum bubbles are destroyed by Floquet instabilities regardless of λ. At time t 10 in Figure 14, the average relative phase of all three systems reaches ∼ 0 eventually. The true vacuum bubbles in the system with the lowest λ survive longer due to their larger sizes. In the case where Floquet instabilities are not negligible, the tuning of the modulation depth λ plays a role in the balance between time duration and strength of true vacua signals.
In order to remove the unstable Floquet modes from the system, one can increase the oscillation frequency ω to shift k c to higher wave-numbers above the increased momentum cutoff used in this Appendix. Figure 15 illustrates the simulated dynamics of the false vacuum at a fixed temperature τ = 1 × 10 −4 , and increasing ω.
In Figure 15, the results show that the stabilization of true vacua is sensitive to an increase of ω. Both the lifetime and the phase signal of the true vacua are enhanced. In comparison with Figure 12 ( ω = 50 at τ = 1 × 10 −4 ), increasing ω to 150 ( Figure 15) extends the sizes of the true vacuum bubbles and the chaotic fluctuations are suppressed in the true vacua. The survival time of the true vacuum bubbles before domination by fluctuations is also improved slightly when ω is increased, from roughly t 30 for ω = 50 to t 40 for ω = 150.
As mentioned, the simulated system with ω = 50 in Figure 12 includes the Floquet modes by setting k nyq ≈ 32.14 well above the critical wavenumber k c ≈ 15.39; while for the system with ω = 150, the Floquet modes are partially removed as the critical wavenumber is increased to k c ≈ 28.79. This partial removal of the Floquet modes reduces the chaotic fluctuations and results in the partial stabilization of true vacua.
If we further increase the modulation frequency to ω = 200, the critical wavenumber is shifted to k c ≈ 33.57 which is slightly higher than k nyq . In this case, the Flo-  quet modes are almost completely excluded. Figure 16 shows that in the dynamics of the false vacuum with ω = 200, one can see that the size of the true vacuum bubbles are extended, and a break-down of the vacuum bubbles is not observed in the simulation.  Other parameters are as in Figure 15.
Increasing ω to suppress the chaotic fluctuations induced by the presence of the Floquet instabilities has a significant effect on the stabilization of true vacua. Figure 17 shows that the suppression of the chaotic fluctuations enhances the average relative phase cosφ a . The relative phase can reach cosφ a ≈ −0.75 even at the high reduced temperature τ = 1 × 10 −4 , which is comparable to a measurement at a lower temperature τ = 1 × 10 −5 ( cosφ a ≈ −0.8 ∼ −0.9, in Figure 4).   Figure 15. The errors in cosφ a are less than 1%.
In summary, bubble nucleation and break-down are delayed by increasing ω to remove the Floquet modes. From the averaged results using 8000 trajectories, true vacua reach their peak size at a time t ≈ 10 for ω = 50 and t ≈ 20 for ω = 150, and are later destroyed by chaotic fluctuations due to Floquet modes (i.e. cosφ a ≈ 0). For ω = 200, the time of the peak relative phase is delayed to t ≈ 35 and cosφ a remains ∼ −0.8 over the remaining simulation time duration. The high dimensionless frequency ω = 200 in the simulations corresponds to an oscillation frequency ω = 2π × 38.24kHz, which is achievable in current experiments. We emphasize that this still requires an overall momentum cutoff, although at a higher wavenumber.