Nonlinear Dynamics of the Cold Atom Analog False Vacuum
Abstract
We investigate the nonlinear dynamics of cold atom systems that can in principle serve as quantum simulators of false vacuum decay. The analog false vacuum manifests as a metastable vacuum state for the relative phase in a twospecies BoseEinstein condensate (BEC), induced by a driven periodic coupling between the two species. In the appropriate low energy limit, the evolution of the relative phase is approximately governed by a relativistic wave equation exhibiting true and false vacuum configurations. In previous work, a linear stability analysis identified exponentially growing shortwavelength modes driven by the timedependent coupling. These modes threaten to destabilize the analog false vacuum. Here, we employ numerical simulations of the coupled GrossPitaevski equations (GPEs) to determine the nonlinear evolution of these linearly unstable modes. We find that unless a physical mechanism modifies the GPE on short length scales, the analog false vacuum is indeed destabilized. We briefly discuss various physically expected corrections to the GPEs that may act to remove the exponentially unstable modes. To investigate the resulting dynamics in cases where such a removal mechanism exists, we implement a hard UV cutoff that excludes the unstable modes as a simple model for these corrections. We use this to study the range of phenomena arising from such a system. In particular, we show that by modulating the strength of the timedependent coupling, it is possible to observe the crossover between a second and first order phase transition out of the false vacuum.
Keywords:
false vacuum, phase transitions, analog cold atom systems, early Universe1 Introduction
False vacuum decay, the quantum firstorder phase transition out of a metastable “false vacuum” state, plays an important role in current models of the early Universe. This includes theories of false vacuum eternal inflation (see e.g., ref. Guth:2007ng ()), symmetry breaking phase transitions in GUT theories Guth:1981uk (), and perhaps even the Standard Model of particle physics Ellis:2009tp (). The dynamics of these transitions has implications for anthropic solutions to the cosmological constant problem, the observability of false vacuum eternal inflation in the multiverse Aguirre:2007an (); Feeney:2010dd (); Feeney:2010jj (), the future fate of the Higgs vacuum Ellis:2009tp (), and primordial gravitational wave signals Kosowsky:1991ua () within the band of future experiments such as the Laser Interferometer Space Antenna (LISA) (see e.g., ref. Axen:2018zvb ()).
Further, as an intrinsically quantum mechanical phenomenon (i.e., requiring ), false vacuum decay touches on fundamental issues in quantum field theory, including the notions of measurement and the emergence of classicality. Therefore, a complete computational and conceptual framework to understand false vacuum decay has implications both for our understanding of the cosmos and for the foundations of quantum field theory.
In the standard view, false vacuum decay is essentially a field theory version of tunnelling through a barrier familiar from quantum mechanics Coleman:1977py (); Callan:1977pt (); Coleman:1980aw (). However, an important distinction between quantum mechanics and quantum field theory is that the former involves only a single degree of freedom, while the latter involves interactions between infinitely many degrees of freedom Banks:1973ps (). This has important consequences, both conceptually and computationally. In the case of quantum mechanics, it is straightforward to evolve the wavefunction using the Schrödinger equation, thus obtaining a complete numerical solution to the tunneling problem. However, an analogous approach is infeasible in quantum field theory due to the exponential complexity of the state space. In the standard formalism Coleman:1977py (); Callan:1977pt (), false vacuum decay is reformulated as a quantum mechanics problem, resulting in a drastic dimensional reduction of the underlying phase space. Bubble nucleation is then interpreted as a quantum tunneling event, with no realtime description available. We recently presented an alternative description of vacuum decay based on a careful study of dynamics in the full field space, considering the cooperative dynamics of many field degrees of freedom Braden:2018tky (). Using this approach, we were able to reproduce the expected decay rates of the standard formalism, while also giving a realtime description of the decay process that does not rely on quanutm tunneling. Given the differing approximations and subsequent interpretations of these two approaches, it would be greatly informative to study a physical system that undergoes false vacuum decay, where Nature automatically includes all quantum effects.
One proposal to build such a system uses a twocomponent dilute gas cold atom BoseEinstein condensate (BEC) to simulate a relativistic scalar field with a false vacuum potential Fialko:2014xba (); Fialko:2016ggg (). In this proposal, the false vacua are generated from unstable local maxima by periodically modulating the direct conversion between the two components. From the field theory perspective, this amounts to modulating the overall amplitude of the potential, which can lead to a stabilization of the longwavelength modes as in the Kapitza pendulum. However, as recently shown by Braden et al. Braden:2017add (), this simultaneously destabilizes shorter wavelength modes, whose exponential growth is expected to lead to a breakdown of the analogy with the relativistic scalar field system, and a classical destabilization of the false vacuum. Therefore, the feasibility of such an experiment rests on a physical mechanism to remove these dynamically unstable shortwavelength modes, leaving only the longwavelength degrees of freedom for which the false vacuum description is valid. Another recent work explored the role of impurities (vortices) in the condensate in seeding vacuum decay Billam:2018pvp ().
In this paper, we explore the viability of the analog false vacuum in detail through the use of nonlinear simulations of cold atom BECs. Our primary focus is on the viability of generating a metastable state within the BEC, rather than providing a detailed dictionary between the BEC system and a relativistic scalar field. First, we verify the destabilization of the condensate by the parametrically excited shortwavelength instabilities, confirming the results obtained in the linear analysis of ref. Braden:2018tky (). We also briefly explore the subsequent nonlinear dynamics. In this case, not only is the false vacuum interpretation of the experiment disrupted, but the entire relativistic scalar field analogy breaks down. In order to restore the false vacuum interpretation, corrections to the GPE description must appear at wavenumbers around the Floquet band to either damp out or otherwise remove these exponentially unstable modes. Simultaneously, these corrections should be negligible for the longer wavelength modes needed to form the bubbles. In order to be viable, such corrections must be motivated by real physical effects, rather than numerical tricks such as a convenient choice of grid spacing. We briefly outline some physically plausible modifications to the shortwavelength dynamics of the condensate. Since modifications to the GPE at shortwavelengths will generically leak into the dynamics of modes required to form bubbles, the impact of these modifications must be accounted for when making quantitative predictions. A detailed study of these subtleties is beyond the scope of this work and will be explored in future publications. Regardless of the scalar field reinterpretation, the BEC system still represents a metastable state that decays as a result of quantum effects (if the system is safe from the unstable Floquet modes). Therefore, even in the absence of a direct link to early Universe models, such systems can be used to investigate various fundamental issues in quantum field theory.
The remainder of this paper is organized as follows. First, section 2 presents the theoretical description of dilute gas cold atom BoseEinstein condensates, modeled by the GrossPitaevskii equation (GPE). Next, section 3 briefly reviews the connection between the GPE and the dynamics of a relativistic scalar field, in particular how to generate a false vacuum potential minimum. Our main results are contained in section 4 and section 5. In section 4, we explore the manifestation of Floquet instabilities in the full nonlinear condensate dynamics, verifying that they lead to a breakdown of the effective relativistic field theory description of the condensate fluctuations. We briefly comment on some physical effects that could tame these instabilities and restore the analog false vacuum interpretation. Meanwhile, section 5 explores the range of nonlinear field theory dynamics that can arise once the Floquet modes have been exorcised, showing how we can slowly deform a secondorder phase transition into a firstorder phase transition through the tuning of an experimental parameter. Finally, we conclude in section 6. Some more technical aspects of the analysis are presented in a set of appendices. In appendix A we briefly outline our semiclassical lattice approach—the truncated Wigner approximation—used for numerical simulations. Appendix B introduces our dimensionless units and equations of motion, and provides a dictionary between our choices and the previous literature. Some further evidence that we are seeing the Floquet instability and details of the deformed linear dispersion relationship induced by finitedifferencing stencils are given in appendix C. Finally, tests of our numerical code are presented in appendix D, including direct convergence tests and a demonstration that relevant Noether charges are conserved.
2 2Component Cold Atom BoseEinstein Condensates in the Dilute Gas Limit
The analog false vacuum system under investigation is a 2component BEC. There are currently two different configurations that can in principle be utilized as false vacuum decay simulators. The first one is a singlespecies BEC in a doublewell potential in the socalled tightbinding regime, exhibiting symmetric and antisymmetric singleparticle states Albiez:2005 (); Bucker:2011 (). False vacuum decay requires a modulation of the tunnelcoupling via the trapping parameters. At present, this configuration is capable of mimicking false vacuum decay in +dimensions. To avoid dimensional restrictions altogether, we focus on a different configuration: a two species BEC of atoms with mass , e.g., K or Rb atoms in two different hyperfine states with a modulated linear coupling between the two species. In this paper we will consider the case of a condensate confined to a ringtrap, which we model as a +dimensional condensate. The analogy between the two species system and false vacuum decay has been derived in detail in ref. Braden:2017add (), although the derivation applies to the doublewell setup as well. In the following, we will summarize the key results of ref. Braden:2017add (), which motivate the nonlinear numerical simulations presented in the subsequent sections.
We consider an ultracold highly diluted 2component system consisting of two singleparticle states ( and ) of mass . Within this setup the atomatom interactions can be approximated by wave contact potentials: , and . Atoms in different hyperfine states have angular momentum and hence a slightly different energy. It is possible to drive transitions between two hyperfine states with a given transition rate through the utilization of external fields, e.g., a radiofrequency field. When the BEC forms, the collective longwavelength excitations of the condensed atoms (i.e., the condensate) acquire a vacuum expectation value (vev), which we denote by the complex numbers . The Hamiltonian density for the two interacting condensates is given by
(1) 
Here we are dropping the effects of the external confining potentials on the dynamics, and we have defined
(2) 
The combined evolution of the 2component system is given by
(3a)  
(3b) 
where we have assumed the interatomic scattering strengths are of the form and defined . The effects of quantum fluctuations are then approximated by generating realizations of complex Gaussian random fields, and superimposing these on the appropriate background homogeneous values of . These serve as the initial conditions for subsequent evolution of the coupled GPEs. The entire procedure is known as the truncated Wigner approximation. Further details about the truncated Wigner approximation, as well as our numerical approach to solve the coupled GPEs and generate initial quantum fluctuations, are given in appendix A. For the remainder of this work, we consider only the limit and , which we refer to as a symmetric theory. To generate an effective false vacuum potential, we further consider a periodically modulated coupling,
(4) 
In order to connect with the dynamics of relativistic scalar fields, it is convenient to consider an alternative set of Hermitian canonical coordinates, the condensate densities and phases , such that
(5) 
In terms of these, the coupled GPEs (3) become
(6a)  
(6b) 
To map the coupled 2component BEC onto the relativistic false vacuum decay Braden:2017add (), it is natural to instead use the total density and phase, and the relative density and phase. For notational convenience we will denote these by
(7a)  
(7b) 
respectively. We will refer to the pairs of canonically conjugate variables and as the total and relative phonons, respectively.
3 Effective Relativisitic Scalar Field Dynamics
We now briefly review the effective relativistic scalar field theory lurking within the dynamics of the cold atom BECs. For a detailed path integral derivation and additional details, the interested reader can refer to refs. Braden:2017add (); Fialko:2014xba (); Fialko:2016ggg (). Analogous results can also be obtained by working directly with the equations of motion. In order to streamline the presentation, here we provide the key results needed to interpret the remainder of the paper. Throughout, we assume that the GPE description of the previous section holds for all relevant modes of the condensates.
To relate the BEC dynamics to that of a relativistic scalar field, we first consider the evolution of a purely homogeneous background field configuration. We denote solutions in the homogeneous approximation by an overbar . Potential vacua are identified by looking for eigenstates of the Hamiltonian, with the eigenvalue () playing the role of the chemical potential and sets .^{1}^{1}1Including fluctuations changes the time evolution of the spatial average of the total phase compared to the evolution in the homogeneous background, so that . This can be interpreted as a correction to the chemical potential from the fluctuations. We verified that this effect appears in our nonlinear simulations. For , we find stationary solutions when , with the state having the higher energy and thus representing a potential false vacuum. We next consider fluctuations around these stationary solutions to the homogeneous background equations:
(8a)  
(8b)  
(8c)  
(8d) 
where an overbar indicates a solution to the coupled GPEs in the homogeneous limit.^{2}^{2}2Because of the nonlinear relationship between the Cartesian condensate variables () and the density () and phase variables (), the volume averaged density and the background density differ by a UV divergent quantity . The equality holds only when the fields are exactly spatially homogeneous. The density fluctuations are then integrated out, and only the modes well below the healing length are considered. The result is a (nonlinear) effective Lagrangian for the phases that has the form of a pair of relativistic scalar fields.
As explicitly laid out in Braden et al. Braden:2017add (), a number of conditions must hold for an effective relativistic scalar field theory for the fluctuations to emerge from the procedure outlined above. These include:

the evolution of the coupled condensates are given by the GPE;

the background dynamics of the condensate are well approximated as homogeneous;

the local fluctuations in the particle densities are small, and can be treated quadratically in the action;

the shorttime dynamics of the fluctuation modes are simple and can be integrated out to obtain an effective timeaveraged description; and

the relative and total phonons can be treated as effectively decoupled, leading to a truncation rather than integration of the total phase fluctuations.
The first four of these are generic requirements, independent of the particular couplings in the GPE, while the fifth assumption (the decoupling between total and relative phonons) is specific to fluctuations in theories with expanded around a background state of equal densities in the two condensates. This latter condition holds for the stationary solutions with . In more general cases, the background condensates can be made spatially inhomogeneous. By modifying the kinetic term in the resulting effective Lagrangian, under certain conditions this leads to an interpretation of a scalar field evolving in a classical background (such as a gravitational field).
Under these assumptions, the relative and total phases of the condensates behave approximately as decoupled relativistic scalar fields. The relative phase obeys
(9) 
which is the equation of motion for a field with potential
(10) 
In the above we have defined a mean particle density
(11) 
and the parameter
(12) 
which controls the shape of the effective potential, and can be experimentally tuned. For , this reduces to the sineGordon model, while for , the potential for the effective relativistic scalar develops a series of local minima, as illustrated in fig. 1.
Meanwhile, the Lagrangian has a global symmetry associated with an overall phase rotation of both condensates, so the total phase is a Goldstone mode and does not develop a potential. We can therefore identify the relative phase phonons as the appropriate variable in which to study metastability, even in the nonlinear regime.
4 Floquet Instabilities and Nonlinear Dynamics of Coupled Condensates
Above we briefly outlined the connection between dilute gas cold atom BECs, described by coupled GrossPitaevskii equations, and relativistic scalar fields. For simplicity, we only considered the symmetric case, resulting in equal densities for each condensate in the homogeneous false vacuum state.^{3}^{3}3Only is required to obtain equal condensate densities, not the additional requirement that we have included in our definition of a symmetric theory. The equal background densities, in turn, ensured that the total and relative phonons decouple from each other in the linear regime. Ignoring backreaction and rescattering of the inhomogeneities, the background dynamics is independent of the fluctuations and acts as an external input to the fluctuation equations. Therefore, the choice ultimately allows us to treat the linear dynamics of the total and relative phonons independently. As listed in the previous section, a number of assumptions about the condensate dynamics must hold in order for the false vacuum derivation to be valid. In this section, we will use nonlinear simulations to test these assumptions. Specifically, we will assume the condensates evolve according to the coupled GPEs, and that initially they are wellapproximated as spatially homogeneous. Given this, the emergence of an effective relativistic scalar field theory posessing an analog false vacuum rests on two key assumptions about the resulting fluctuation dynamics:

perturbations in the local number densities are small and can be integrated out at quadratic order in the action, and

all relevant dynamical condensate modes only feel the timeaveraged effects of the modulations.
Further, reducing to a single effective scalar requires a decoupling between the total and relative phase phonons. The validity of these assumptions for the linearized fluctuation dynamics was studied in detail by Braden et al. Braden:2017add (). We showed that for , the introduction of periodic modulation in of frequency induces a band of exponentially unstable modes with wavenumber centered at
(13) 
where . This exponential growth occurs in both the relative phase and relative density fluctuations, which if left unchecked will lead to a breakdown of the first two assumptions above. The emergence of these exponentially growing modes required only that the condensates:

followed the coupled GPEs (3),

were subjected to a periodic modulation in ,

could be treated as nearly spatially homogeneous, and

small initial inhomogeneous fluctuations were present.
Subject to these assumptions, the presence of the Floquet band is unavoidable. The third condition is not strictly necessary to excite Floquet instabilities, but consideration of an inhomogeneous background will significantly modify the scalar field interpretation, so we will not pursue it here. The fourth condition simply provides the seed for exponential linear growth, and can be fulfilled either by requiring that the initial fluctuations satisfy the uncertainty principle, or by populating the Floquet band through nonlinear interactions described by the coupled GPEs.
The first obstacle to realizing an analog false vacuum is to ensure that either the Floquet instabilities are under control or that their presence does not qualitatively alter the dynamics of the longer wavelength modes essential for bubble nucleation. One possibility is that the full nonlinear dynamics sequesters the effects of linear instability from the IR modes required to form bubbles. Below we demonstrate that, by itself, nonlinear evolution with the coupled GPEs does not provide a mechanism to isolate the effects of the Floquet instability from the IR modes.
Therefore, the existence of the analog false vacuum depends on the existence of a physical mechanism in the UV to prevent the Floquet modes from appearing. Some examples of these corrections include:

corrections to the wave scattering approximation from atoms retaining memory of the internal structure of the interaction potential between collisions;

corrections to the shortwavelength dynamics of the condensate from the finite mode occupancy;

the discrete nature of the atoms leading to a breakdown of the continuum field description;

interactions of the collective phonon modes with the uncondensed thermal modes of the gas; and

corrections to the onedimensional approximation from excitations of transverse degrees of freedom associated with the trapping potential in the transverse directions.
If the fourth effect is large, then we expect the decay dynamics to be driven by thermal fluctuations, rather than vacuum fluctuations. Meanwhile, the final effect will lead to strong modifications to the effective relativistic field theory description. Therefore, the first three effects are the most promising from the viewpoint of analog false vacuum decay. Although all three of these effects will most strongly modify the dynamics at large wavenumbers, they may also change the IR dynamics of the bubbles. Therefore, detailed analysis is required to understand any modifications to the false vacuum decay interpretation that results from these corrections to the truncated Wigner description.
There are two potentially relevant UV scales: the healing length of the condensate (the length scale over which the condensate will respond to a defect or boundary) and the interatomic separation. For symmetric condensates, the wave numbers associated with the healing length and interatomic separation are
(14) 
The validity of the dilute gas BEC approximation requires , making the healing length the first relevant UV scale.
Determining if any of these mechanisms is active and sufficient to remove the Floquet modes requires detailed physical modelling of any proposed experimental setup, which we leave to future work. An example of a relevant effect in another experimental setup can be found in ref. 2018arXiv180707564N (). Referring to (13), we also see that a degree of experimental tunability exists, since either increasing the modulation frequency or decreasing the condensate density will move the instability to higher wavenumbers relative to the healing length. However, one must also beware that increasing the driving frequency could cause additional modes (such as transverse excitations in the trap or additional internal excited states of the atoms) to become excited. If this happens, then the dynamics will differ dramatically from the coupled GPEs discussed here.
4.1 Nonlinear dynamics
We now study the full nonlinear dynamics of the coupled condensates, properly accounting for the presence of the Floquet band. For concreteness, throughout this section we choose parameters
(15) 
which were suggested as reasonable values by Fialko et al. Fialko:2014xba (). The initial fluctuation amplitudes for the nonlinear simulations are set by the dimensionless particle density
(16) 
again as suggested in ref. Fialko:2014xba (). The maximal Lyapunov exponent for fluctuations around the fiducial false vacuum as a function of dimensionless wavenumber are shown in fig. 2. We see a narrow band of unstable wavenumbers centered around , as predicted by (13). Note that for these choices of parameters, this is slightly above the UV scale associated with the healing length, so it is not unrealistic to expect corrections to appear. However, since we wish to understand the implications of assuming the coupled GPEs hold, we proceed without considering any such corrections. Fig. 2 provides a sharp prediction for the wavenumbers at which the GPE description must break down in order to maintain the analogy with false vacuum decay.
Since the Floquet modes represent a piece of physics beyond that of the timeaveraging analysis used to obtain the false vacuum, we expect that the full dynamics will undergo a drastic change as numerical simulation parameters are tuned to either include or exclude their effects. In fig. 3 we illustrate the full nonlinear dynamics of the relative condensate phase using the numerical procedure described in appendix A. To demonstrate the crucial role the Floquet band has on the dynamics, we have utilized a series of numerical simulations with varying lattice spacings (i.e., Nyquist wavenumbers) to isolate the effects of the exponentially unstable modes. The two panels of fig. 3 show the evolution of the relative phase for one of these simulations. The left panel has so that , and the lattice cutoff is just below the start of the instability band illustrated in fig. 2. As a result, the unstable Floquet modes cannot be resolved by our simulation, and we expect the timeaveraging analysis to be valid. Indeed, we see the dynamical nucleation and subsequent expansion of domain wallantiwall pairs (i.e., bubbles in one dimension) in the relative phase, indicating a firstorder phase transition. By contrast, the right panel takes and , which is above the upper edge of the instability band.^{4}^{4}4To ensure that the dynamics of the shortwavelength modes is properly resolved, we have chosen the lattice cutoff small enough to saturate the convergence in both and (see appendix D), and ensure the Noether charges are preserved to the level. To demonstrate the change in behavior occurs from effects localized in the Floquet band, we also ran a simulation with and . This grid just resolves the full Floquet band, and is only a % change in the grid spacing compared to the left panel of fig. 3. The qualitative behavior matches that of the right panel of fig. 3, with the false vacuum state being lost due the Floquet modes. As expected, we see a drastic change in the resulting dynamics of the relative phase. During the initial stages (while the Floquet modes still have a small amplitude), the evolution matches that of the left panel. However, at the exponentially growing linear modes become and undergo strong modemode coupling. The evolution is subsequently dominated by shortwavelength modes of large amplitude, competely erasing the firstorder phase transition dynamics. The Floquet modes therefore qualitatively change the condensate dynamics from that of an effective relativistic scalar field trapped in a false vacuum. In each simulation, we used the same realization of the initial fluctuations. In particular, since the lattice cutoff is below the Floquet band in the left panel, this means we did not initially populate the exponentially unstable modes. Although physically unrealistic since it violates the uncertainty principle, this was done to demonstrate that the nonlinear dynamics of the coupled GPEs will populate these modes. Since fluctuation power is being dynamically generated rather than input as an initial condition, it also demonstrates the need to have a physical model for the behavior of these modes. Here we have taken this model to be the coupled GPEs, and shown that this assumption is inconsistent with the existence of a false vacuum decay.
This is further demonstrated in fig. 4, where we show the lattice average of for the simulations in fig. 3. At the false vacuum , while in the true vacuum , so that if the firstorder phase transition interpretation holds, the expectation value should go from approximately to . Corrections to the pure behavior seen in fig. 3 arise from domain walls and fluctuations present in the condensates. However, we clearly see that in the absence of Floquet modes this intuition holds, while when the Floquet modes are present the transition to is lost and instead we end up with . This latter behaviour is what we expect if the relative phase is randomly distributed.
Having examined the evolution of the relative condensate phase, we now turn to the local density fluctuations. This will allow us to study two interesting points regarding the validity of the relativistic scalar field interpretation of the condensate dynamics. First, as outlined in section 3, the scalar field description rests on integrating out the density perturbations to quadratic order to generate a kinetic term for the relative phase phonons. This step will be invalid if the density fluctuations become large as a result of the Floquet instability. In the linear regime, relative density perturbations grow commensurate with the relative phase perturbations for the exponentially unstable modes, and we therefore expect that this assumption will be badly violated in the case where Floquet instabilities remain. Second, in the linear regime the total and relative phase perturbations decouple (when expanded around either the false or true vacuum). This decoupling property was essential in truncating the total phase dynamics to obtain an effective action for the relative phase alone.^{5}^{5}5As mentioned above, the emergence of the relative phase as the key tunneling variable comes directly from the Hamiltonian term proportional to , and thus holds to nonlinear order as well. Therefore, the decoupling of the relative and total phase phonons allows for a simpler treatment of single field tunneling, as opposed to a more complex two field problem. As we show below, in the absence of the Floquet modes, both the small amplitude and decoupling approximations hold throughout the simulation. However, they are both broken as soon as the dynamics of the Floquet band is considered.
Fig. 5 shows the evolution of the relative density between the two condensates for a simulation without (left) and a simulation with (right) the unstable Floquet band present. The simulations used correspond to the left and right panels of fig. 3. As expected, in the absence of the Floquet instability, the fluctuations remain small throughout the evolution. We also note that the locations of the bubble walls and collision regions evident in the relative phase also manifest as coherent fluctuations in the relative density. Meanwhile, when the Floquet modes are present the evolution is dramatically different. At early times, the evolution matches that of the simulation without the Floquet band. However, at , fluctuations of short wavelengths dominate the evolution, indicating a breakdown of the interpretation in terms of a relativistic scalar field.
In fig. 6 we instead plot the evolution of the total condensate densities. As with the relative density fluctuations, we see stable behavior in the absence of the Floquet modes. Unlike the case of the relative density fluctuations, we see no evidence for the bubbles in the evolution of the total density. This indicates that, in the absence of the exponentially growing modes, the relative and total phonons remain decoupled. However, as soon as the Floquet dynamics are correctly captured by the simulations, large fluctuations in the total particle density develop as well. Since the total density phonons are stable in linear perturbation theory, this indicates that the growth of fluctuations is driven by nonlinear scattering with the growing Floquet modes of . Thus, the decoupling between the total and relative phonons is destroyed by the presence of the Floquet resonance.
Finally, fig. 7 shows the RMS fluctuations in the relative and total density fluctuations, normalized to the (constant) mean total density over the lattice. This provides an alternative visualization of the breakdown of the effective scalar field description, and subsequent loss of decoupling between the relative and total phonons. In particular, in the absence of the Floquet instability, the fluctuation amplitude remains stable throughout the course of the simulations. However, as soon as the Floquet band is present, we see a rapid initial growth in the fluctuations of both the relative and total energy densities. The growth slows somewhat at , consistent with partial quenching of the linear instability by nonlinear modemode coupling as seen in the real space evolution above. Further, the growth of fluctuations in the total density slightly lags the growth in the relative density. Such a lag is consistent with total density fluctuations being induced by nonlinear scattering with the exponentially growing relative density fluctuations. This further shows that the loss of decoupling between the total and relative phonons occurs because of nonlinear interactions arising from the rapid growth of Floquet modes in the linear regime.
In this section, we have investigated the effects of Floquet instabilities arising in the analog false vacuum decay proposal of Fialko et al. Fialko:2014xba (); Fialko:2016ggg (), which were originally pointed out in Braden et al. Braden:2017add (). The existence of these instabilities leads to a complete breakdown of the analogy between the BEC evolution and a relativistic scalar field, and thus spoils the false vacuum analogy. Therefore, the assumption that the condensates obey the coupled GPEs at all scales is inconsistent with the creation of a false vacuum by external periodic modulation of a coupling constant. Meanwhile, in cases where the unstable modes are artificially removed from the dynamics, the false vacuum interpretation is restored, at least at a qualitative level. Therefore, realizing false vacuum decay experiments within the current proposal requires the presence new shortwavelength physics not captured by the GPE to enter at length scales longer than . This new physics must remove the exponentially growing Floquet modes, while simultaneously exerting a minimal influence on the dynamics of the bubbleforming modes. Two examples of how this may occur are: damping of the growing modes, and an effective projection of the modes to remove them from longwavelength condensate dynamics. However, we expect the effects of realistic physical corrections will also leak into the dynamics of the longer wavelength modes needed to nucleate the bubbles. Therefore, detailed analysis of the nonlinear dynamics is required to understand how the corrections needed to tame the Floquet instabilities modify the interpretation of the bubble nucleation dynamics. We leave such an exploration to future work.
5 First and Second Order Phase Transitions of the Analog False Vacuum
In the previous section we confirmed the presence of linear Floquet instabilities within the GPE description of the analog false vacuum cold atom BEC proposal. Using nonlinear simulations, we further demonstrated the need to consider physical corrections beyond the use of coupled GPEs and the truncated Wigner approximation in order to preserve the false vacuum interpretation. We left the complex issues of how these modifications may arise, and how they impact the interpretation in terms of a relativistic scalar to future work. To help motivate these future studies, we now investigate the nonlinear dynamics under the assumption that we can find an experimental realization where the Floquet modes have been removed and the truncated Wigner approximation is valid, in which case the effective scalar field description should hold. A range of dynamical symmetry breaking phenomena can be arranged, including both second and firstorder phase transitions. In each case, the transition is driven by the dynamics of quantum fluctuations, opening the possibility of experimentally testing many fundamental aspects of quantum field theory.
For simplicity, throughout this section we remove the Floquet modes by setting a hard lattice cutoff so that the Nyquist frequency is below the Floquet band. However, one should be aware that the introduction of the finite lattice cutoff can lead to both a numerical pileup of fluctuation power at the grid scale, and aliasing of short wavelength power into longer wavelengths. Each of these can significantly modify the nucleation rate of bubbles, and in the former case can lead to a spurious nucleation of bubbles. This is especially a problem when the cutoff wavenumber is below that required to properly capture the full profile of a nucleated bubble. To overcome this limitation, in this section we also take the frequency of the external driver to be sufficiently large that the Floquet band occurs at higher wavenumbers than those needed for a convergence of the IR modes. The grid spacing for all simulations is , which is sufficiently converged that the simulation results are indistinguishable by eye as the grid spacing is varied. This is simply a numerical trick. In a real experiment one must ensure that choosing a large driving frequency does not excite additional degrees of freedom, such as internal states of the trapped atoms or higher harmonics of the trap. However, our goal is to highlight some of the interesting behavior that can be obtained in a pristine theoretical environment, rather than deal with the many subtleties and additional effects that must be accounted for in an actual experiment.
As shown in refs. Fialko:2014xba (); Fialko:2016ggg (); Braden:2017add () and outlined in section 3, the temporal modulation of creates an effective potential for the longwavelength modes of
(17) 
illustrated in fig. 1, where the parameter
(18) 
is experimentally tunable. For , is a local maximum, while it is a local minimum for . In the (timeaveraged) scalar field interpretation of , increasing adjusts the field from sitting on top of a hill to sitting in a local minimum. For initial states localized around the false vacuum , we expect two distinct mechanisms by which the false vacuum state can decay. The value of determines which mechanism is relevant.
For , modes with will experience a tachyonic instability, We confirmed this at the level of linear perturbation theory in ref. Braden:2017add (), where we also numerically obtained and dependent corrections. Based on this picture, we expect that for , falls off the top of the hill through a spinodal instability. This leads to the rapid creation of a set of domains (with the field localized near either the or vacuum) connected by a network of domain walls. Meanwhile, for the local maxima become local minima of the potential, corresponding to false vacuum states. In this regime, we obtain a firstorder phase transition and the false vacuum decays via the nucleation of bubbles. Of course, we do not necessarily expect a sharp distinction between the secondorder and firstorder phase transitions to occur at . In particular, when , the barrier between the false and true vacuum is narrow and shallow. Therefore, we expect that in some spatial regions the initial condensate fluctuations will probe over the barrier of the potential and fall towards the true vacuum. For very shallow barriers, these transitions happen in a significant fraction of the volume, and rather than having a few welldefined bubble nucleations, many regions of the true vacuum appear simultaneously. Such dynamics is intermediate between the tachyonic growth of longwavelength modes and the regime of rare welldefined bubble nucleations.
Fig. 8 illustrates the nonlinear stabilization of the low tachyonic modes by increasing the modulation amplitude of the interspecies conversion rate. The top two panels have . We observe the rapid emergence of a dynamically evolving domain wall network, as expected from the tachyonic growth of the longwavelength modes. The transitionary behaviour between first and secondorder phase transitions is seen in the middle panels of fig. 8. Finally, welldefined bubble nucleations appear for sufficiently larger than one, as seen in the bottom panels of the figure.
Fig. 8 also illustrates the dependence of the false vacuum decay rate on the experimentally tunable parameter . This effect appears in the figure as a delay in the time at which bubbles nucleate as the value of is increased,^{6}^{6}6This is clear in fig. 8 since we start from the same realization of the initial fields in each case. In the more general case, where the initial conditions are also sampled, then this effect would only be evident upon considering an ensemble of timeevolutions, rather than direct comparison of individual simulations. and allows for a tuning of typical decay time relative to the expected experimental lifetime of the condensate. Additionally, the decay rate is sensitive to the number density of condensed atoms, which changes the amplitude of fluctuations in relative to the absolute mean. Since the RMS amplitude of fluctuations is set by through the uncertainty principle, this is analogous to adjusting the value of . This may provide a window to adjust the expected decay rate while holding the effective scalar potential (i.e., the scalar field theory) fixed.
5.1 Validity of Scalar Field Interpretation and Phonon Decoupling
As we saw in section 4, when Floquet instabilities induced by the external driving of the system are present both the phase and local number density fluctuations grow rapidly. This invalidates the assumption that density fluctuations are small, which is needed to derive the effective scalar field description. We now investigate whether or not a similar breakdown occurs when the shortwavelength Floquet modes have been excised from the evolution, but an exponential instability remains in the longwavelength modes. In particular, in the regime of spinodal instability (), the exponentially growing longwavelength modes that destabilize the local maximum are simply another band of exponentially unstable linear modes. It is therefore important to explicitly check the validity of various approximations in this regime.
Fig. 7 shows the evolution of the RMS of both the total () and relative () energy densities for the simulation. Although initially the relative RMS of the perturbations grows, unlike in the previous section it saturates before becoming order one. Further, the fluctuations in the total density remain fairly stable throughout the simulation. It thus appears that the effective relativistic scalar field interpretation of these simulations will continue to hold (at least at leading order), and that the total and relative phonons remain decoupled for the entire evolution.
This is further illustrated by the evolution of the relative and total densities within the same simulation, seen in fig. 10. Comparing to fig. 3, we see that the density perturbations remain significantly smaller than in the case of excitations in the higher order instability band. Additionally, large coherent structures in the relative density appear, presumably due to the numerous coherent domain walls in the relative phase. However, no such structures appear in the total density fluctuations. Combined with the lack of growth in these fluctuations overall, this suggests that the total phase phonons remain decoupled from the relative phase phonons in this case as well.
6 Conclusions
In this paper, we investigated the nonlinear stability of analog false vacuum decay experiments in cold atom BoseEinstein condensates. In particular, we first showed that using the coupled GrossPitaevskii equations as a physical model for the evolution of the condensates is inconsistent with the generation of an effective false vacuum through modulation of the interspecies conversion rate. To do this, we confirmed the presence of the linear Floquet instabilities studied in ref. Braden:2017add () in full nonlinear simulations. We then demonstrated that excitation of these modes leads to a complete breakdown of the effective scalar field description, and therefore the false vacuum decay analogy. Since the exponentially unstable modes appear at nonzero wavenumber, it is possible that their effects can be removed from the system if the GrossPitaevskii equations used to model the condensates break down at wavelengths above the instability scale. However, this requires modeling of the condensates beyond the coupled GPEs, and such corrections must be analyzed on a casebycase basis to determine whether a given experiment is viable. Further, the original analogy between the relativistic scalar field and cold atom dynamics was obtained assuming the coupled GPEs were the correct description. Therefore, the mapping of any corrections to the coupled GPEs into the scalar field dynamics must also be understood in order to make quantitative comparisons. Given these complexities, we did not pursue a full study here. Rather, our results should be used as an indication that care must be exercised when developing analog cosmological systems using cold atoms.
If feasible, an analog false vacuum system would provide a unique opportunity to experimentally study fundamental issues in both cosmology and relativistic quantum field theory. Under the assumption that Floquet modes can be removed, we investigated the variety of behavior accessible to experiments by tuning the amplitude of the oscillating interspecies conversion rate. This parameter controls the shape of the potential for the relative phase, taking a local potential maximum to a local potential minimum. In the absence of temporal modulation, if the field begins at the local potential maximum, we demonstrated that the longwavelength modes of the condensates experience tachyonic growth, and the initial homogeneous state decays via a spinodal instability. Meanwhile, well above this critical amplitude we observed the stabilization of the longwavelength tachyonic modes and the decay of the metastable state through the formation of domain wallantiwall pairs (i.e., onedimensional bubbles). The transition between these two regimes is not sharp, but instead near the critical amplitude the decay occurs in an intermediate regime posessing partial characteristics of both the spinodal and bubble nucleation instabilities. In each of these regimes, we found that removing the shortwavelength Floquet instability was sufficient to ensure the density perturbations remain small and the approximate relativistic scalar field description remains valid. The study of this crossover behavior is one example of the type of physical phenomena accessible to the cold atom systems studied in this paper.
Our results show the need for detailed physical modelling of the condensate dynamics in order to find an experimental system to realize analog false vacuum decay. In future work, we will investigate corrections to the GPE in specific experimental scenarios in order to determine if the necessary conditions to remove the Floquet instabilities can in principle be met. We will also determine the detailed mapping (should one exist) to the analog scalar field theory. Theoretical work is also necessary to narrow down the specific observables within the cold atom system (e.g., correlation functions) that will be most useful for shedding light on theoretical aspects of vacuum decay. Vacuum decay is both strongly nonlinear and strongly quantum mechanical, making it highly likely that we will learn something truly new from analog false vacuum decay. We believe this is strong motivation for future investigation.
JB and HVP were supported by the European Research Council (ERC) under the European Community’s Seventh Framework Programme (FP7/20072013)/ERC grant agreement number 3006478CosmicDawn. JB is supported by the Simons Foundation Origins of the Universe program (Modern Inflationary Cosmology collaboration). The work of JB, HVP, and AP was partially performed at the Aspen Center for Physics, which is supported by the National Science Foundation grant PHY1607611. The participation of JB, HVP, and AP at the Aspen Center for Physics was supported by the Simons Foundation. MCJ is supported by the National Science and Engineering Research Council through a Discovery grant. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science. AP was supported by the Royal Society. SW acknowledges financial support provided by the Royal Society University Research Fellow (UF120112), Royal Society Enhancement Grants (RGF/EA/180286,RGF/EA/181015), EPSRC Project Grant (EP/P00637X/1), and partial support from STFC consolidated grant No. ST/P000703/. This work was partially enabled by funding from the University College London (UCL) Cosmoparticle Initiative.
References
 (1) A. H. Guth, Eternal inflation and its implications, J. Phys. A40 (2007) 6811 [hepth/0702178].
 (2) A. H. Guth and E. J. Weinberg, Cosmological Consequences of a First Order Phase Transition in the SU(5) Grand Unified Model, Phys. Rev. D23 (1981) 876.
 (3) J. Ellis, J. R. Espinosa, G. F. Giudice, A. Hoecker and A. Riotto, The Probable Fate of the Standard Model, Phys. Lett. B679 (2009) 369 [0906.0954].
 (4) A. Aguirre, M. C. Johnson and A. Shomer, Towards observable signatures of other bubble universes, Phys. Rev. D76 (2007) 063509 [0704.3473].
 (5) S. M. Feeney, M. C. Johnson, D. J. Mortlock and H. V. Peiris, First Observational Tests of Eternal Inflation: Analysis Methods and WMAP 7Year Results, Phys. Rev. D84 (2011) 043507 [1012.3667].
 (6) S. M. Feeney, M. C. Johnson, D. J. Mortlock and H. V. Peiris, First Observational Tests of Eternal Inflation, Phys. Rev. Lett. 107 (2011) 071301 [1012.1995].
 (7) A. Kosowsky, M. S. Turner and R. Watkins, Gravitational radiation from colliding vacuum bubbles, Phys. Rev. D45 (1992) 4514.
 (8) M. Fitz Axen, S. Banagiri, A. Matas, C. Caprini and V. Mandic, Multiwavelength observations of cosmological phase transitions using LISA and Cosmic Explorer, Phys. Rev. D98 (2018) 103508 [1806.02500].
 (9) S. R. Coleman, The Fate of the False Vacuum. 1. Semiclassical Theory, Phys. Rev. D15 (1977) 2929.
 (10) C. G. Callan, Jr. and S. R. Coleman, The Fate of the False Vacuum. 2. First Quantum Corrections, Phys. Rev. D16 (1977) 1762.
 (11) S. R. Coleman and F. De Luccia, Gravitational Effects on and of Vacuum Decay, Phys. Rev. D21 (1980) 3305.
 (12) T. Banks, C. M. Bender and T. T. Wu, Coupled anharmonic oscillators. 1. Equal mass case, Phys. Rev. D8 (1973) 3346.
 (13) J. Braden, M. C. Johnson, H. V. Peiris, A. Pontzen and S. Weinfurtner, A New Semiclassical Picture of Vacuum Decay, 1806.06069.
 (14) O. Fialko, B. Opanchuk, A. I. Sidorov, P. D. Drummond and J. Brand, Fate of the false vacuum: towards realization with ultracold atoms, EPL 110 (2015) 56001 [1408.1163].
 (15) O. Fialko, B. Opanchuk, A. I. Sidorov, P. D. Drummond and J. Brand, The universe on a table top: engineering quantum decay of a relativistic scalar field from a metastable vacuum, J. Phys. B50 (2017) 024003 [1607.01460].
 (16) J. Braden, M. C. Johnson, H. V. Peiris and S. Weinfurtner, Towards the cold atom analog false vacuum, JHEP 07 (2018) 014 [1712.02356].
 (17) T. P. Billam, R. Gregory, F. Michel and I. G. Moss, Simulating seeded vacuum decay in a cold atom system, 1811.09169.
 (18) M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani and M. K. Oberthaler, Direct observation of tunneling and nonlinear selftrapping in a single bosonic josephson junction, Phys. Rev. Lett. 95 (2005) 010402.
 (19) R. Bücker, J. Grond, S. Manz, T. Berrada, T. Betz, C. Koller et al., Twinatom beams, Nature Physics 7 (2011) 608 EP .
 (20) N. Navon, C. Eigen, J. Zhang, R. Lopes, A. L. Gaunt, K. Fujimoto et al., Synthetic dissipation and cascade fluxes in a turbulent quantum gas, arXiv eprints (2018) [1807.07564].
 (21) C. W. Gardiner, J. R. Anglin and T. I. A. Fudge, The stochastic grosspitaevskii equation, Journal of Physics B: Atomic, Molecular and Optical Physics 35 (2002) 1555.
 (22) G. N. Felder and I. Tkachev, LATTICEEASY: A Program for lattice simulations of scalar fields in an expanding universe, Comput. Phys. Commun. 178 (2008) 929 [hepph/0011159].
 (23) A. V. Frolov, DEFROST: A New Code for Simulating Preheating after Inflation, JCAP 0811 (2008) 009 [0809.4904].
 (24) R. Easther, H. Finkel and N. Roth, PSpectRe: A PseudoSpectral Code for (P)reheating, JCAP 1010 (2010) 025 [1005.1921].
 (25) Z. Huang, The Art of Lattice and Gravity Waves from Preheating, Phys. Rev. D83 (2011) 123509 [1102.0227].
 (26) J. Sainio, PyCOOL  a Cosmological ObjectOriented Lattice code written in Python, JCAP 1204 (2012) 038 [1201.5029].
 (27) J. Braden, J. R. Bond and L. MersiniHoughton, Cosmic bubble and domain wall instabilities I: parametric amplification of linear fluctuations, JCAP 1503 (2015) 007 [1412.5591].
 (28) J. Butcher, Implicit RungeKutta Processes, Math. Comp. 18 (1964) 50.
 (29) J. P. Boyd, Chebyshev and Fourier Spectral Methods, Dover Books on Mathematics. Dover Publications, Mineola, NY, second ed., 2001.
Appendix A Stochastic Semiclassical Simulations: The Truncated Wigner Approximation
In this appendix we briefly outline our methodology to solve for the dynamics of a twocomponent cold atom BEC described by the coupled GPEs. The method, when applied to the GPE, is widely used in the cold atom community and known as the truncated Wigner approximation Gardiner:2002 (). More cosmologically oriented readers will recognize the spiritual similarity with semiclassical scalar field lattice simulations used to study preheating Felder:2000hq (); Frolov:2008hy (); Easther:2010qz (); Huang:2011gf (); Sainio:2012mw (). The basic idea is to model quantum fluctuations by sampling realizations of the initial quantum state as encoded in the Wigner functional. The subsequent complex dynamical evolution is then treated by solving the classical equations of motion. This amounts to propagating a classical probability distribution functional on (very highdimensional) field space via the method of characteristics. Given this, it is natural to interpret the results of a single simulation as the outcome of a single experimental realization. This approach is known to be exact when the field dynamics is linear and the initial state is Gaussian. In more complicated nonlinear situations, it can sometimes capture the crucial elements of the dynamics, although the validity must be checked on a casebycase basis.
Initially, the condensates are assumed to be of the form
(19) 
where the are background solutions, here taken to be stationary points of the coupled GPEs in the spatially homogeneous limit.^{7}^{7}7The total phase can experience a linear timeevolution at the stationary points, but we include the background chemical potential to remove this. Quantum corrections to the “classical” background are modeled by a realization of a Gaussian random field
(20) 
where the are draws of uncorrelated complex Gaussian random deviates with unit variance , and is the simulation volume. We generate our samples using the BoxMueller transform
(21) 
where and are uniform random deviates on the interval . Since the fluctuations are assumed to be Gaussian, the spectra should match the initial quantum 2point statistics^{8}^{8}8An obvious extension allows this to be extended to incorporate correlations between the condensates.
(22) 
In an experiment, this will depend on the preparation of the initial state; while in a purely theoretical study it should be adjusted to accurately reflect the physics under investigation. In particular, these initial fluctuations must be mapped into those of the the relativistic scalar fields in order to correctly interpret the results in the effective scalar field theory. Some care is need to choose the correct initial fluctuations to reproduce the initial Minkowski state for the effective scalar field . To make contact with the previous literature, here we follow the approach of Fialko et al. Fialko:2014xba (); Fialko:2016ggg () and choose for constant , resulting in a filtered white noise spectrum for . Considering the nonlinear transformation into the density and phase variables, we see that this does not correspond precisely to the standard Minkowski vacuum of the effective relativistic scalar field. We leave consideration of how this impacts the decay rate, and the appropriate initial state to properly simulate the false vacuum, to future work.
From this sampled initial configuration, the fields are evolved using the coupled GPEs. For our numerical scheme, we use a tenthorder accurate GaussLegendre integrator (see e.g. refs. Braden:2014cra (); Butcher:1964 ()) for the temporal evolution, and a collocation based Fourier pseudospectral approach boyd01 () for spatial discretization. When useful, we also compare with a secondorder accurate finite differencing approach for the spatial discretization, with
(23) 
and
(24) 
for an arbitrary scalar function . In all cases, we enforce periodic boundary conditions, explicitly in the case of a finitedifference stencil, and implicitly when using a Fourier pseudospectral method.
Appendix B Dimensionless Variables
In this appendix we briefly present the dimensionless form of our evolution equations and summarize the most important dynamical timescales. It is convenient to define dimensionless time, space, and condensate variables (here denoted by a bar),
(25) 
through the introduction of three numerical scales , and with dimensions of , , and respectively. It is also convenient to introduce the numerical sound speed () and energy (),
(26) 
Expressed in these variables, the dimensionless Hamiltonian is
(27) 
The (dimensionless) combination sets the overall scale of the numerical Hamiltonian, and thus does not enter directly into the evolution equations. It does, however, set the overall scale of the initial dimensionless fluctuations . Meanwhile, the numerical energy can be used to transform each of the remaining dimensionful couplings into a dimensionless coupling constant.
Alternatively, working directly with the equations of motion, we have
(28) 
We have introduced the dimensionless numerical coefficients
(29) 
For convenience, we now specialize to the two condensate case, and define
(30) 
For our numerical simulations, we scale the condensates by their average background densities , set the numerical sound speed to be the propagation speed of the relative phase phonons, and the dimensionless nonlinear potential coupling to unity
(31) 
We immediately see
(32) 
so that the condensates are measured in units of their background density (11), and positions in units of the healing length (see (14)). The resulting dimensionless equations of motion suitable for numerical simulations are
(33a)  
(33b) 
The commutator for the scaled field variables is
(34) 
which we can write on a discrete lattice with lattice sites
(35) 
For our particular normalization of the fields (32), the relative amplitude of dimensionless fluctuations in the condensate scales with the total number of condensed particles as ,
(36) 
Since the commutator sets the overall amplitude of vacuum fluctuations through the uncertainty relation, the relative amplitude of density fluctuations to the density in the zero mode scales as .
b.1 Comparison with Dimensionless Units of Fialko et al.
To ease comparison with previous work, here we also present the dimensionless units used in Fialko et al. Fialko:2014xba (); Fialko:2016ggg () and translate into the notation of our paper. To avoid notational confusion, we denote the various normalization constants of Fialko et al. with superscript and the corresponding dimensionless quantites by . As above, the dimenionless units in our paper are denoted by an overbar . In those papers, they continue to set the numerical sound speed to the propagation speed of the relative phase fluctuations. However, they measure time in units of the oscillation frequency of the effective scalar field associated with the relative phase, rather than choosing to measure position in units of the healing length. The condensate wavefunctions are further normalized to the inverse numerical volume scale . The corresponding numerical scales are related as
(37) 
Finally, they measure in units of instead of , so that
(38) 
Appendix C Demonstration of Linearity of Floquet Modes and the Effective Discrete Lattice Wavenumber
In section 4 we demonstrated the presence of new shortwavelength dynamics that were neglected in previous studies of analog false vacuum decay. This new dynamics destroyed the effective scalar field description, and also resulted in a loss of the phase transition as measured by the mean value of . Since this dramatic change occurred at precisely the wavenumbers predicted by a Floquet analysis of the linear perturbations, we concluded that this was driven by the Floquet instability. We now explicitly demonstrate that the new shortwavelength dynamics are indeed a result of a linear instability. To do this, we rerun the simulations in fig. 3 using a secondorder accurate centered finitedifferencing approximation for the Laplacian instead of a Fourier pseudeospectral approximation. The resulting evolution of the relative phase is shown in fig. 11.
As outlined in the subsection below, the effective wavenumber felt by linear fluctuations is determined by the choice of discrete Laplacian. In general it differs from the continuum value as shown in fig. 12. In particular, for the Nyquist mode, the effective linear wavenumber for the secondorder accurate Laplacian stencil is times the continuum wavenumber. Since Floquet instabilities arise from linear dynamics, on the lattice the Floquet modes will feel this effective wavenumber. In contrast, nonlinear interactions not involving derivatives are sensitive to the continuum wavenumber of the fluctuations, and therefore they are not directly influenced by the choice of spatial discretization. Therefore, if the new dynamics seen with decreasing lattice spacing were associated with nonlinear interactions amongst the fluctuations, it would appear at the same Nyquist frequency, regardless of the choice of Laplacian stencil. For linear fluctuations, the dynamics will instead appear at different wavenumbers determined by the effective linear wavenumber associated with the choice of discrete Laplacian. As seen in fig. 11, the emergence of the smallscale instability is sensitive to the effective lattice wavenumber, indicating they arise from linear perturbation dynamics. Moreover, they appear at precisely the wavenumber predicted by (linear) Floquet theory. This provides very strong evidence that the effects we are seeing are a well understood physical effect associated with evolution by the coupled GPEs.
c.1 Effective Lattice Wavenumber for Discrete Laplacian Stencils
In this subsection we briefly derive the linear dispersion relationship associated with a discrete approximation to the Laplacian operator. The distortion of the effective wavenumber for linear fluctuations from the continuum limit by the use of finite differencing stencil is wellknown, but is presented here for completeness.
Consider a dimensional discrete lattice with sites labelled by the vector index . We will denote function values at position by . For simplicity, assume that the lattice sites are arranged in a rectangular grid with uniform lattice spacing , so that . Consider finitedifference stencils for the Laplacian operator of the form
(39) 
with the vector indices . The choice of coefficients specify the stencil. Now consider a linear operator of the form
(40) 
with a constant. For simplicity, we consider only a single function defined on our lattice. The generalization to many linearly coupled fields is straightforward. Taking a discrete Fourier transform of (40), we obtain in onedimension
(41) 
with
(42) 
In the continuum, the RHS of (41) would be , so we identify an effective wavenumber associated with the Laplacian stencil as a function of the continuum wavenumber
(43) 
where is the Nyquist wavenumber. To avoid spurious numerical damping associated with asymmetric stencils, we use a symmetric stencil so that the imaginary contribution vanishes, and we have
(44) 
In fig. 12, we show this effective wavenumber for several loworder symmetric onedimensional Laplacian stencils.
We also note one further complication to the picture above. Our code evolves the condensate fields , and the numerical Laplacian is thus defined to act directly on these fields rather than the densities and phases . However, our analytic derivation of the Floquet instability was done in the density and phase variables. Because of the nonlinear transformation between these sets of variables, additional distortions to the required numerical stencils appear in the equations of motion for and . However, these distortions vanish in the limit of strictly linear inhomogeneities. The excellent match between our analytic predictions for the location of the Floquet band and the finitedifferencing results is thus even further evidence that the new physics seen in our nonlinear simulations is simply the manifestation of the linear Floquet instability.
Appendix D Numerical Convergence and Conservation Tests
In this appendix we present a variety of convergence and consistency tests to illustrate the precision of our numerical approach. We provide two tests, including: a direct test of numerical convergence with variation of either the time step or grid spacing ; and a demonstration that relevant conserved charges are timeinvariant. These tests will demonstrate that we are correctly solving for the condensate dynamics under the assumption they are described by the coupled GPEs. Of course, more precise physical modeling will result in corrections to the GPE description. However, these corrections are not modeled by the errors introduced by approximate numerical methods.
Throughout this appendix, we consider convergence properties for two fiducial simulations presented in the main text: the simulation in the top left panel of fig. 8 (with and no Floquet modes), and the simulation shown in fig. 3 (with and Floquet modes). For reference, the parameters are
(45) 
and
(46) 
respectively. The initial conditions for each simulation (including the particular realization of the fluctuations) were the same. The former simulation doesn’t possess a false vacuum, and thus experiences a spinodal instability as shown in section 5 rather than nucleating bubbles. However, by removing the explicit timedependence, the energy of the system will be conserved and thus provide us an additional diagnostic tool for testing our numerics. As well, the Floquet instabilities will not be present in this case, allowing spatial convergence to be tested in a much simpler way. In order to allow a direct comparison between simulations with different time steps or grid spacings, the initial conditions are identical for each simulation. There are several dynamical timescales that must be resolved by our choice of timestep . In the linear regime, these are roughly the driving frequency , the linear frequency of the Nyquist mode for the total () and relative () phonons, and the growth rate of the Floquet modes.^{9}^{9}9Technically, the oscillation frequencies should be obtained from the imaginary parts of the Floquet exponents, but and can be used as reasonable proxies. The free oscillation frequencies of the Nyquist modes can be easily obtained (see Braden et al. Braden:2017add ())
(47a)  
(47b) 
Since , we have . To obtain spatial resolution, we will be interested in the limit where the Floquet band is below the Nyquist mode, which occurs roughly when . Therefore, for spatially resolved simulations, we expect that and will dominate the temporal convergence properties. Of course, once the fields become strongly nonlinear, new timescales may emerge, and the interactions between different modes may modify the effective oscillation frequencies much like plasma effects. However, we will take the linear scales given above as a guideline.
d.1 Direct convergence tests
First we show direct pointwise convergence tests with respect to changes in both the time step and grid spacing . These demonstrate the rapid convergence displayed by our GaussLegendre time stepping and pseudospectral discretization, respectively. To explore pointwise convergence, we consider the following norm between two simulations and
(48) 
which measures the maximal pointwise difference between the two simulations over the entire simulation volume. Here, the superscript labels the simulations, and is a collection of spatial points at which to evaluate our functions. In our convergence testing, we will consider differences between simulations in which either the time step or grid spacing differ by a factor of , while holding the grid spacing or time step fixed respectively. When only the time step is varied, the spatial grids match exactly, and we simply have to ensure that we compare simulations at the same time. For our equally spaced Fourier collocation grid, each subsequent spatial grid refinement is a superset of the previous one, so the comparison can be done directly on the coarser grid. Although we expect to see rapid convergence because of our use of highly accurate numerical methods, we also expect the presence of exponentially growing Floquet modes to slowly degrade the quality of the convergence for two reasons. First, our temporal integrator is symplectic, and thus preserves phase space volumes. For linear fluctuations, this corresponds to accurately tracking the overall amplitude of oscillating Fourier modes. However, the numerical tradeoff for this is a small error in the oscillation phase, and a corresponding error in pointwise comparisons at a fixed time. For exponentially growing modes, this error will grow with time as the modes grow in amplitude. For highly nonlinear fluctuations, there will generally be localized stuctures that will propagate through the simulation volume. Small temporal phase errors lead to small changes in the shapes and velocities of these structures. Over time, the same structure present in two simulations with different choices of may follow a slightly different spacetime path, leading to a growing pointwise difference between the simulations.
Below we show the maximal pointwise norm (48) between simulation pairs with either (fig. 13) or (fig. 14) varied. As the time step is decreased, initially we see the pointwise error decrease by a factor of , as expected for a tenthorder accurate integrator. This eventually saturates at a level of for the simulations without Floquet modes present, while the saturation level increases with time in the presence of Floquet modes. As explained above, this is not unexpected due to the exponential growth of the Fourier modes, and the tradeoff made in temporal phase and amplitude accuracy described above. To further identify the exponentially growing Floquet modes as the root cause of the timedependent saturation of the temporal accuracy, we also performed the same temporal convergence tests for and for the simulation, corresponding to a gridspacing that just excludes and just includes the unstable Floquet band. We found a growing saturation level for the latter case, while the convergence plateaued in a manner similar to the simulation for the former case. As a test that these errors did not arise from aliasing into longwavelength fluctuations, we performed the corresponding convergence tests with a finitedifferencing Laplacian stencil (which does not lead to aliasing) and found the same growth in the saturation level.
Meanwhile, as the spatial grid spacing is decre