Towards the cold atom analog false vacuum

Analog condensed matter systems present an exciting opportunity to simulate early Universe models in table-top experiments. We consider a recent proposal for an analog condensed matter experiment to simulate the relativistic quantum decay of the false vacuum. In the proposed experiment, two ultra-cold condensates are coupled via a time-varying radio-frequency field. The relative phase of the two condensates in this system is approximately described by a relativistic scalar field with a potential possessing a series of false and true vacuum local minima. If the system is set up in a false vacuum, it would then decay to a true vacuum via quantum mechanical tunnelling. Should such an experiment be realized, it would be possible to answer a number of open questions regarding non-perturbative phenomena in quantum field theory and early Universe cosmology. In this paper, we illustrate a possible obstruction: the time-varying coupling that is invoked to create a false vacuum for the long-wavelength modes of the condensate leads to a destabilization of shorter wavelength modes within the system via parametric resonance. We focus on an idealized setup in which the two condensates have identical properties and identical background densities. Describing the system by the coupled Gross-Pitaevskii equations (GPE), we use the machinery of Floquet theory to perform a linear stability analysis, calculating the wavenumber associated with the first instability band for a variety of experimental parameters. However, we demonstrate that, by tuning the frequency of the time-varying coupling, it may be possible to push the first instability band outside the validity of the GPE, where dissipative effects are expected to damp any instabilities. This provides a viable range of experimental parameters to perform analog experiments of false vacuum decay.


Introduction
False vacuum minima are ubiquitous in relativistic field theories with complex interaction potentials. Although classically stable, large spatial regions where the fields are trapped in such false vacua can undergo decay via the nucleation of bubbles (either through quantum [1,2] or classical statistical fluctuations [3]). In cosmology, such transitions are a crucial ingredient in many incarnations of the multiverse (see e.g. [4] for a review), and may have also played a role in the early history of our own Hubble volume if symmetry breaking phase transitions occurred at some sufficiently high temperature [5]. Since cosmological evolution is intrinsically time-dependent, it is of great interest to have a time-dependent understanding of the phenomenology of bubble nucleation. Unfortunately, standard descriptions of this process [1,2] cannot properly tackle the time-dependent nature of the problem. Instead, the temporal dependence is introduced in an ad hoc manner with bubble nucleation events assumed to be independent and distributed in the spacetime according to a Poisson process. While some analytic progress has been made in understanding a Lorentzian picture of quantum tunnelling [6,7], here we focus on a possible laboratory experiment that can be used as a quantum computer for simulating vacuum decay in ensembles of "Universes" in order to understand the true quantum dynamics of first order phase transitions. Should such experiments be feasible, they would provide insight into fundamental open questions such as: • In the standard instanton description of bubble nucleation [1,2], one imagines that a spherically symmetric bubble with a critical radius "appears" via quantum tunneling, and eventually is well-described by the classical field configuration of a bubble of true vacuum expanding with a constant acceleration. What precisely happens before and during the "appearance" of the bubble? Are there any precursors to bubble nucleation, corresponding to intermediate field configurations, before a final transition to the expanding bubble? How should we understand the transition to the classical limit at late times?
• Instanton techniques do not give a complete description of multi-bubble configurations, and nucleation events are considered to be completely independent. In the full dynamics, are there non-trivial correlations between nucleation events? If so, there are potentially important implications for if/how percolation occurs in quantum first order phase transitions in quantum field theory.
• A deeper understanding of percolation could also inform our understanding of eternal inflation, where accelerated expansion of space can prevent percolation. Are there important correlations that give rise to clusters of bubbles, with different properties than the clusters that appear for uncorrelated nucleation events [8,9]? If so, there are potentially observable consequences from the collisions between bubbles in such clusters [10][11][12][13].
• Since the bubble interior has the symmetries of an open Friedman-Robertson-Walker (FRW) spacetime at late times, bubble nucleation represents the quantum creation of a Universe. Are there any novel observable signatures of the quantum origins of the Universe stored in the resulting field configurations?
Recently, Fialko et al [14,15] proposed such a condensed matter system to simulate the false vacuum decay of a relativistic scalar field in the laboratory. This is one amongst a variety of analog condensed matter systems that have been proposed to study various elements of relativistic quantum field theory. An interesting aspect of the false vacuum decay proposal is that it purports to simulate fully nonlinear quantum dynamics of an interacting field theory. This is in contrast to many proposed analog gravity systems (such as those involving analog Hawking and Unruh radiation [16][17][18], superradiance [19], and cosmological particle production [20,21]), which consider the kinematics of small amplitude non-interacting fluctuations propagating on an externally controlled classical background. Similarly, proposals to study the dynamics of solitons [22] and nonlinear structures such as oscillons (see e.g. [23,24]) typically focus on classical evolution in a nonlinear system [25]. In the linearized cases, the connection to relativistic fields is made through a recreation of the linear equations of motion for the fluctuations, rather than through nonlinear phenomena induced by mode-mode coupling.
In its simplest form, the proposed system for false vacuum decay consists of a pair of ultra cold gases with 2→2 intra-and inter-particle scattering interactions. The two particle species are then coupled through a radio-frequency field, which leads to a bilinear coupling describing the direct conversion of one species into the other. Below the critical temperature, the atoms condense into a Bose-Einstein condensate (BEC) phase of long-wavelength collective motion, whose dynamics can be modeled by the coupled Gross-Pitaevskii equations (GPEs); see e.g. [26] for a pedagogical review. This allows for a quantitative treatment of the dynamics of the long wavelength mean field condensate.
The (complex) condensates are conveniently rewritten in terms of the squares of their moduli (interpreted as a local condensate particle density) and their phases (whose gradients are proportional to the local condensate velocities). By integrating out fluctuations in the number density, the effective dynamics of the two phase variables is approximated by a pair of relativistic scalar fields. In the limit that the bilinear coupling goes to zero, the Hamiltonian is independent of the phase of each condensate, resulting in two global U (1) symmetries. The introduction of a non-zero coupling breaks the symmetry associated with a global change to the relative phase, while preserving the symmetry under a simultaneous rotation of both condensates. This generates a sinusoidal potential for the relative phase, while the total phase remains massless.
The final step to obtain an analog false vacuum is to modulate the bilinear coupling periodically in time. For sufficiently long wavelengths, the effect of this temporal modulation is to convert the local maxima of the original sinusoidal potential into a series of local minima, thus leading to an effective scalar field with a series of degenerate false vacua. However, the same oscillating coupling that is introduced to stabilize the zero mode of the system simultaneously induces additional instabilities at higher wavenumbers that are not accounted for in the time-averaged effective potential. This is an example of parametric resonance, which has been well-studied in the context of early Universe cosmology [27][28][29] (see [30] for a recent review), and has been experimentally verified in a single-component BEC [20]. If these wavenumbers are part of the long-wavelength collective modes described by the coupled GPEs, then the exponential growth of these modes will eventually destabilize the condensate as a whole, simultaneously invalidating the conclusion that the field is sitting in a false vacuum. Although we have described a particular implementation, the key feature for false vacuum decay is that the GPEs describe the dynamics of the longwavelength condensates, but not the specific details of the atom-atom interactions at very short wavelengths.
In this paper, we perform a detailed linear stability analysis of the proposed analog false vacuum setup, assuming the initial state can be described as a homogeneous background with small inhomogeneous perturbations on top. We assume the GPE provides an adequate description of the condensate dynamics, and perform the calculation in the full "UV"theory on the condensed matter side, rather than in the effective scalar description. We demonstrate the stabilization of the zero-mode within linear perturbation theory as we increase the amplitude of the modulation of the bilinear coupling. For modes below the healing length, we also show that the growth rates and transition to linear stability agree well with the predictions of the effective action for the effective relativistic scalar field. As anticipated above, while the zero-mode is stabilized through the increasing amplitude of the modulations, we simultaneously observe the emergence of exponentially unstable modes with nonzero wavenumbers.
We identify the oscillation frequency of the modulations as the key parameter, with the characteristic wavenumber of the resonant instabilities scaling approximately as the square root of the driving frequency on scales shorter than the healing length. If these exponentially unstable modes are correctly modeled by the GPE, they will lead to a catastrophic fragmentation of the proposed false vacuum. For sufficiently short wavelengths we expect viscous and stochastic corrections to the GPE will appear. By pushing the unstable modes into the viscous regime, we thus expect they will damp away, thus preserving the metastability of the system. However, the required frequency must not be too large, or additional degrees of freedom (e.g. radial modes for condensates in a ring trap, or mediating transitions for hyperfine splitting) can be excited. Therefore, to fully assess the feasibility of studying false vacuum decay in the lab, the corrections to the GPE at short wavelengths must be modeled, and the corresponding impact of these corrections on the dynamically induced instabilities must be understood. We do not undertake this comprehensive analysis here. However, the theoretical link between the GPE and a relativistic scalar undergoing false vacuum decay is interesting in its own right, and allows for a new theoretical toolbox to be applied to the problem. Therefore, even absent a complete understanding of the experimental feasibility, theoretical investigation of the proposed analog system is of interest.
In order to consider the simplest physical situation, and avoid additional technical complications, in this paper we assume that the two condensates have equal effective masses, equal self-interaction scattering rates, and equal background densities (which we refer to as the symmetric limit). This symmetric setup illustrates the basic mechanisms at play for more general choices of condensate parameters. Of course, in a real experiment these symmetry assumptions may not be realized, so it is also of interest to perturb off the symmetric limit (while still maintaining the assumption of a homogeneous background). The linearized analysis in these asymmetric cases will be presented in a forthcoming publication, where we describe the rich additional dynamics that emerge.
The remainder of this paper is organised as follows. In section 2 we briefly review the physics of dilute gas Bose-Einstein condensates (BECs), present the equations of motion describing the condensates, and briefly comment on the regime of validity of the equations. We then connect the BEC dynamics to relativistic scalar field theory in section 3, and explicitly demonstrate how a false vacuum effective potential can be engineered for the relativistic scalar degree of freedom. Next, we obtain suitable homogeneous background solutions around which to perform our fluctuation analysis in section 4. The main results of the paper are contained in section 5, where we analyze the linear fluctuation dynamics around the proposed false vacuum using the machinery of Floquet theory. Finally, we conclude with a summary and some future directions in section 6. To streamline the main presentation, technical details regarding the emergence of a relativistic scalar field (appendix A) and the derivation of the linearized perturbation equations (appendix B) are contained in a pair of appendices.

Coupled Bose Condensates
The condensed matter system under consideration is a 2-component BEC, an ultra-cold gas of N atoms containing two single-particle states |1 and |2 such that N 1 +N 2 = N . An example for such a system is 87 Rb in two different hyperfine states. If the gas is ultra-cold and highly diluted, the atom-atom interaction potential can be approximated by an S-wave contact potential. With this approximation, there are three kinds of atom-atom contact interactions present in the system. They are parameterized by three coupling constants: g 11 and g 22 within each component, and g 12 = g 21 across components. Each hyperfine state has a different total angular momentum, and hence slightly different energies. One can employ external fields to set up three-level scattering processes to drive transitions between the hyperfine states with transition rate ν.

Model Hamiltonian and Gross-Pitaevskii Equations
We first derive the relevant equations for the 2-component BEC. The boson field annihila-tionΨ i and creation operatorsΨ † i for the single-particle states satisfy In terms of these, the Hamiltonian for the BECs in external trapping potentials V ext,i (x) is given byĤ = d d xĤ, with the Hamiltonian densitŷ The coupled Gross-Pitaevskii equations for the 2-component BEC follow from the Heisenberg equations of motion, When working with the path integral in the next section, it will be convenient to have an explicitly real Hamiltonian density. Integrating by parts, we can re-express (2.3) aŝ up to a boundary term. In our analysis, we will ignore the effects of the external potentials V ext,i . We could imagine, for example, that the condensate is trapped in a ring with strong confining potentials in the two orthogonal directions; see [21] for such a setup. The transverse excitations can then be integrated out leaving a theory for the longitudinal modes with no effective trapping potential. However, one needs to be careful not to excite radial eigenmodes as exhibited in [21]. Setting V eff,i = 0, it is then consistent to consider dynamics around a homogeneous background. Including the external potential results in a spatially inhomogeneous ground state that must be perturbed. This significantly complicates the linear stability analysis carried out here, so we leave investigation of this case to future work.
If cooled down to sufficiently low temperatures the system undergoes spontaneous symmetry breaking through the formation of a background condensate. It is convenient to redefine our field operatorsΨ (2.7) The Gross-Pitaevskii equations are now obtained by replacing the quantum operatorsΨ i in (2.5) with This obviously neglects nonlinear corrections from interactions amongst the quantum fluctuations around ψ i . Through direct transformation of the resulting equations we obtain so that in the absence of inhomogeneities, iρ i = 0 at every point in space. Accounting for inhomogeneities which follows from the Gross-Pitaevskii equation regardless of the number of condensates. We recognize as the total particle current for the condensates. In the absence of loss of particles through the boundary, the total particle number is conserved. For simplicity, we redefine the coupling constants As stated in the introduction, in this paper we restrict our considerations to the case of symmetric condensates, with m 1 = m 2 ≡ m and g 11 = g 22 . We will assume these conditions hold in the remainder of the paper, while the effects of relaxing these assumptions will be explored in an upcoming publication. Further, it is convenient to work with the following alternative variables: which are (up to a scaling) canonical coordinates. We can directly express the original variables ψ i as Substituting into the Hamiltonian density and accounting for scaling, we obtain the Hamiltonian density K for the variables , ε, ϑ, ϕ: We introduce K to distinguish it from the function obtained by substituting (2.15) directly into (2.6). The overall factor allows K be varied directly with respect to (2.14) in Hamilton's equations.
The remainder of our analysis assumes that the coupled GPEs provide an adequate theoretical description of the two-component BEC. In reality, we expect that at sufficiently short wavelengths additional corrections must be included, with a complete breakdown of the GPE description occurring eventually. To estimate where this occurs, recall that the four particle contact interactions ∼ g ij |Ψ i | 2 |Ψ j | 2 arise as the S-wave scattering limit of the full interaction rate, with the potential interactions between the particles (not to be confused with the external trapping potential) approximated as a delta function. Assuming g 11 = g 22 = g and g 12 = 0, so that there is only a single scattering length a s , the S-wave scattering length a s is defined by g = 4π 2 a s m , (2.18) and the dilute gas approximation requiresna 3 s 1 wheren = n 1 +n 2 2 is the mean particle density.
There are two relevant length scales to consider. The first is the length scale on which fluctuations can probe the internal structure of the interparticle interactions, which occurs when the wavelength of the perturbation is of order the S-wave scattering length. This defines a wavenumber, The second important length scale is the "healing length", which characterizes the length scale over which the condensate can adjust in the presence of a boundary or defect. The wavenumber k heal associated with the healing length L heal can be estimated by comparing the gradient energy ∼ 2 k 2 heal /m to the interaction term ∼ gn in the Hamiltonian, yielding Note that in the dilute gas approximation, k atom is parametrically larger than k heal . When fluctuations probe the internal structure of the interparticle interactions, we cannot replace these interactions by an effective delta function potential, and the original Hamiltonian operator (2.3) from which the GPE descends is no longer valid. Therefore, the GPE is not a good description of the system when fluctuations can probe the internal structure of the interparticle interactions. Meanwhile, the GPE is expected to be a good description on scales below the healing length. We therefore conclude that the breakdown of the GPE will occur somewhere in the range k heal < k < k atom .

Connection to false vacuum decay
We now demonstrate how the cold atom system described above can be connected to a seemingly unrelated physical phenomenon -relativistic false vacuum decay. To make this connection, we first demonstrate that under suitable conditions, the phases ϑ and ϕ of the condensates behave as effective relativistic scalar field degrees of freedom. In the limit ν = 0, the Lagrangian is invariant under independent global rotations of each phase. The formation of the condensates ψ i leads to a spontaneous breakdown of these symmetries, resulting in the presence of two Goldstone modes. The introduction of ν explicitly breaks the U (1) (i.e. shift) symmetry associated with the relative phase, generating a sinusoidal potential in the presence of the condensate. Under a suitable set of experimental conditions (see below), the dynamics of the relative phase variable is described by the sine-Gordon Lagrangian. At this point, the potential for the relative phase posesses a series of stable degenerate global minima, and unstable degenerate local maxima. Each of these local extrema is a stationary point in the strictly homogeneous limit, but small perturbations around the local maxima experience tachyonic instabilities, indicating the maxima are not vacua of the system. The final step to connect to false vacuum decay is to create a series of local potential minima to act as metastable states. The key insight (c.f. Fialko et al [14,15]) is that in the long-wavelength limit the local maxima of the original sinusoidal potential can be converted into local minima through periodic temporal modulation of the coupling ν. Explicitly, as we demonstrate below, the relative phase ϕ has a leading order effective Lagrangian of the form where c s , V 0 , and λ are constants whose explicit form will be given below. Due to the complexity of intermediate expressions, in the main text we sketch the derivation of the effective action (3.1). We explicitly indicate the key features that lead to a relativistic scalar field description, and the required approximations to obtain the sine-Gordon model.
A detailed derivation is given in appendix A. Finally, we show how modulating ν converts the sine-Gordon model into a model suitable for studying false vacuum decay.

Effective Relativistic Field Dynamics from Bose Condensates
It is convenient to apply the path integral formalism for the fluctuation dynamics in the 2component system. The astute reader will undoubtedly notice the close parallel between the derivation in this section and the usual transformation from the Hamiltonian to Lagrangian formalism in the path integral language, which illustrates the generality of the construction. The propagator has the form where the Lagrangian density is Here, H is the real Hamiltonian density (2.6) defined in section 2 and µ is the chemical potential, which is a Lagrange multiplier associated with the conservation of total particle number. Although we work directly in Lorentzian signature, we note that after appropriate Wick rotation we obtain the same results in Euclidean time via the grand partition function.
To obtain an effective action for the phase variables ϑ and ϕ, we must integrate out the density variables and ε. With this in mind, we expand where we note that the time-independence ofn follows from the coupled GPEs. If we assume the inhomogeneous fluctuations in the local number density are small, we can expand (3.3) to second order in δ and δε. We will denote the ith order in this expansion by L In order to simplify the result while maintaining the essential aspects of the derivation, we now introduce the additional symmetry assumptionε = 0. In the next section we will see this constraint is closely related to the assumption g 11 = g 22 . Expanding the Lagrangian, we obtain Since the Lagrangian density is quadratic in δ and δε, we can perform the Gaussian integrals to obtain an effective Lagrangian for the phase variables ϑ and ϕ. From the viewpoint of obtaining a relativistic scalar field theory, the key observation is that the θ δ and φδε terms lead to contributions ∼θ 2 and ∼φ 2 in this effective Lagrangian density. Combining this with the gradient and potential energy contributions already present in L (0) δ , we see that modulo various environmental corrections and distortions to the dispersion relationship, we will obtain a Lagrangian for a relativistic scalar field. To obtain the sine-Gordon model for ϕ, we now make the following set of approximations (details can be found in appendix A): 1. we keep only the leading order gradient energies in (∇ϕ) 2 , (∇ϑ) 2 , and ∇ϕ · ∇ϑ, 2. we assume the experimental parameters satisfy ν/(g ± g c )n 1, 3. we work in the long-wavelength limit and take 2 k 2 /(g c ± g)mn 1 (e.g. on scales larger than the healing length), 4. we ignore backreaction of the fluctuations onto the homogeneous background, and 5. we ignore corrections induced by interactions between ϕ and ϑ by truncating rather than integrating out the dynamics of ϑ.
Combined with the assumptions of a homogeneous background,ε = 0, and smallness of δ and δε, we obtain the following effective Lagrangian for ϕ: which is of course the Lagrangian for the sine-Gordon model. In order to determine the precise relationship between the analog false vacuum dynamics and the BEC experiment, one should carefully account for the various correction terms as well as integrate out the dynamics of ϑ. We leave a detailed investigation of the effects of the environmental corrections and comparison to numerical simulations to future work. Given the many approximations made to establish this analogy, we will be able to use the 2-component BEC not only to test the universality of the false vacuum decay, but also its robustness. At the same time one needs to be careful when extracting quantitative information, e.g. decay rates, from this system. In particular, even if the approximations hold for short time-scales, errors may accumulate over time if the decay rates are too small. Detailed analysis is required to compare the time-scale for this accumulation to the typical lifetime of the condensate, which is beyond the scope of this paper.

Generation of False Vacuum Minima Through Modulated Coupling
From the above analysis, we see that the introduction of the particle conversion parameterized by ν induces an effective sinusoidal potential into the dynamics of the relative phase degree of freedom. Up to higher derivative corrections (not included above), non-canonical contributions to the kinetic term (which are small in the long-wavelength regime and for ν (g − g c )n), and some (possibly nonlocal) corrections from interaction with the environment, we showed that the relative phase variable ϕ has the effective Lagrangian of the sine-Gordon model. Following [14,15], we now demonstrate that the local maxima of cosine potential can be transformed into local minima for the long-wavelength degrees of freedom through an appropriate modulation of the coupling ν: (3.7) Rather than integrate out this time-dependence directly in the propagator, we will first derive an effective time-averaged Hamiltonian for the full system, and use this to derive modifications to the quadratic Lagrangian presented above [31]. Given a Hamiltonian density of the form H = H 0 + ∆H cos(ωt) , (3.8) in the absence of instabilities in the short-wavelength modes that can destabilize the system, modes with characteristic time-scales much less than ω −1 will evolve in the effective Hamiltonian density At each subsequent order of the expansion in ω −1 , an additional commutator with ∆H is introduced. Heuristically, the new factor of ∆H introduces a factor ofn ω, while the action of the commutator differentiates each sinusoidal function and reduces the power of n by one. This means that for the parameterization (3.7) the expansion in ω −1 becomes an expansion in δ. For the purposes of illustration, here we ignore issues of operator ordering and use the correspondence between quantum commutators and classical Poisson brackets whereĤ indicates the quantum operator and H the corresponding classical function. For the case considered here, we have where we have included the to account for the fact that our preferred variables are scaled canonical coordinates. In the Poisson bracket approximation, the contribution from the gradient terms vanish. Therefore, we consider only the potential contributions. A straightforward calculation gives which can be easily (but tediously) verified by working directly with the variables ψ i and ψ * i instead. We now consider how this correction to H modifies the effective action presented above. As before, we assumeε = 0 and m 1 = m 2 . Further, in order to most directly relate to the previous analysis, we make the same sequence of approximations as above. The time-averaged analog of (3.6) is where we have usedL to indicate the effects of the time-averaging. We thus see that the introduction of high frequency oscillations in ν has led to the creation of a stabilizing effective potential for the slow dynamics of the system. Comparing with (3.1) we identify the sound speed (3.14) Note that our expansion is only valid for δ 1, so that (3.14) does not imply the oscillations can induce an imaginary sound speed. At tree level, the potential has the form where we have defined The effective potential (3.15) is illustrated in figure 1. In the remainder of the paper we will explore the validity of this approximation, specifically the assumption that the high temporal frequency dynamics is stable and can be suitably averaged to obtain the effective action (3.13). Unfortunately, some of the elegance of the above derivation is spoiled when we include the higher derivative corrections. Further complications arise if we also break the symmetry between the two condensates, by allowing eitherε = 0 or m 1 = m 2 . We leave exploration of the precise connection between relativistic field theory and symmetry properties of the BECs to future work.

Effective Action for Linear Phase Perturbations
The effective Lagrangians (3.6) and (3.13) are valid for the nonlinear dynamics of ϕ subject to the assumptions outlined in appendix A and summarized above. Due to its nonlinear nature, we can use (3.13) to infer the existence of a metastable false vacuum in the effective relativistic dynamics of ϕ.
However, to connect with the linear fluctuation analysis presented in the remainder of the paper, it is useful to instead obtain the effective Lagrangians for linear fluctuations in ϕ and ϑ. We therefore decompose and work to quadratic order in the full set of fluctuation variables. As will be seen in section 4, we are interested in states with cosφ = ±1 ≡ σ, so we impose this additional constraint. In contrast to (3.6), we make no assumptions about the magnitude of ν/(g − g c )n or 2 k 2 /4mn(g ± g c ). The details of this expansion are provided in appendix A, where we see that the total and relative phase fluctuations decouple providedε = 0, g 11 = g 22 and m 1 = m 2 . In particular, we obtain the following effective Lagrangian for δϕ As we will see below, the equation of motion derived from this effective Lagrangian matches that obtained from direct manipulation of the equations of motion (5.10), indicating that our effective action approach reproduces the correct behavior in the regime of linear fluctuations.
We can similarly incorporate (3.12) to obtain an effective Lagrangian for linear perturbations of the relative phase evolving with the effective time-averaged Hamiltonian which we will use to make analytic predictions for the growth rates of the homogeneous mode below.

Homogeneous Background Solutions
Before proceeding with the analysis of fluctuations, we must first find a background to expand around. In order to find suitable states, we search for eigenstates of the (nonlinear) We can absorb the constant µ into a global phaseψ i = ψ i e iµt witḣ ψ i = 0. Alternatively, we can evolve the condensates with a modified Hamiltonian and search for zero eigenmodes of F .
For simplicity, we assume a homogeneous condensate, which is consistent with the equations of motion in the absence of a trapping potential in the GPE as we have assumed here. It is convenient to define an effective potential where V is the spatial volume of the system. The variables (2.14) introduced in section 2 provide a convenient basis to study the homogeneous dynamics. The Hamiltonian K for these variables is Here,θ andn form a canonically conjugate pair, as doε andφ, and the stability properties are easily analyzed from the Hessian of (4.4). For notational simplicity, it is convenient to definē as well as the dimensionless dynamical parameter x ≡ε n . (4.6) To simplify some notation, we now arrange our variables into the vector It is easily verified that the equations of motion can be obtained by extremizing U bg , so that stationary solutions occur when ∇ q U bg = 0, where the derivative is with respect to our variables (4.7). Stationarity ofε follows from extremization with respect toφ, ε = ∂U bg ∂φ = 0 = 2ν sinφ n 2 +ε 2 . (4.8) Concentrating on the states where both condensates are occupied, we thus see that at the stationary points cosφ = ±1 ≡ σ . (4.9) Next, enforcing stationarity inφ gives us an equation for the relative density differenceε: Finally, the chemical potential µ follows from the stationarity ofθ: To determine the stability of these solutions, consider small deviations δq around a point q 0 which satisfy with the matrix Since ∂θU bg = 0, we can consider (ε,φ) in isolation, as explicitly seen from the characteristic equation for D ij : (4.14) Unfortunately, due to the degeneracy of the Hamiltonian, the simplistic approach given here cannot address the full stability of the system (since the η = 0 eigenvalues are degenerate). However, these eigenvalues are associated with the (n,θ) dynamics, which are seen to possess no dynamical instabilities in the homogeneous limit. The required components of the Hessian are: Thus around our states of interest with sinφ = 0, we have ∂φεU bg = 0. As outlined in section 3, we are interested in cases when the false vacuum minima forφ arise from the local potential maxima after the introduction of temporal modulations in ν. We are therefore interested in unstable states (i.e. η 2 > 0) with ∂φφU bg < 0. The latter condition immediately givesν bg σ < 0 , (4. 16) while the former tells us ∂ 2 U bg ∂ε 2 > 0 . In theε = 0 state, the density and chemical potential are then related by with σ = cosφ = ±1. σν > 0 corresponds to the true vacuum, and σν < 0 to the false vacuum.
directly with the full-time dependent equations of motion, rather than those derived from the time-averaged Hamiltonian (3.19). In appendix B we derive the equations of motion for linear perturbations about an arbitrary homogeneous background in a symmetric two condensate experiment.
Here we want to study fluctuations around the proposed false vacuum. From section 4 this impliesε = 0, sinφ = 0, and ν cosφ < 0. Imposing the former two background conditions in (B.3), we see that the perturbation variables δ and δϑ decouple from δε, and δϕ. It is also convenient to work with the dimensionless density perturbations The equations for the total density (δξ) and phase (δϑ) perturbations are From this we immediately read off the dispersion relationship This of course matches the equation of motion derived from our second order effective Lagrangian (A.12). For g + g c < 0 (i.e. an imaginary sound speed), there is a tachyonic instability for the modes with 2 k 2 ≤ 4mn(g + g c ).
The situation is more interesting for the relative density (δζ) and phase (δϕ) phonons, as they couple to each other through ν: In this case, the dispersion relationship is For ν = 0, the low-k modes thus propagate with sound speed In the absence of ν, we see that the effect of a non-zero cross-coupling g c is to split the total and relative perturbations into two distinct groups characterized by their sound speeds. On the other hand, for g − g c > 0 the σν > 0 states become gapped, with an effective zero-mode frequency while the σν < 0 states instead develop a tachyonic instability for 2 k 2 < 4m |ν|. Finally, if g − g c < 0, an additional tachyonic instability appears for modes with 2 k 2 < 4mn |g − g c | − 4mσν in the σν > 0 phases. However, in the σν < 0 phases, the modes with 2 k 2 < −4mσν are stabilised, while the modes with −4mσν < 2 k 2 < −4mσν + 4mn |g − g c | develop a tachyonic instability.
In order to most directly connect with the analog scalar field dynamics, it is convenient to convert these coupled first order equations into a second order equation for δϕ. Allowing ν = ν(t) while maintaining cosφ = σ = ±1, we have This is equivalent to the Euler-Lagrange equations (A.13) for the effective Lagrangian (A.11) obtained from working to quadratic order in all of our perturbation variables. At linear order, the dynamics of δϕ thus differ from that of a free massive relativistic scalar field in Minkowski space through: the presence of a friction term ∂ t δϕ, a modified dispersion induced by the ∇ 4 term, and the time-dependence of the Laplacian term through ν(t). Furthermore, owing to the presence of the inverse Laplacian, the equations written in second order form are in general nonlocal and depend on the boundary conditions. All three of these corrections can ultimately be traced back to the non-canonical multiplier of the field kinetic energy term in (3.18). In taking the standard sine-Gordon limit, we assumed that the gradient energy of the density fluctuations could be neglected, and that ν 0 /n(g − g c ) 1. The effect of these two approximations is to discard the non-canonical multiplier, so that these corrections do not appear in the equations derived from (3.6). We expect the linear dynamics to be different from a relativistic field in Minkowski space because of the explicit time-dependence of the system at hand. In fact, in the weak resonance regime and in the absence of a potential, we can engineer an effective field in an FRW spacetime, see for example [32,33].

Floquet Formalism
In the previous subsection we looked at the linear stability of perturbations around homogeneous states withφ = iπ (i ∈ Z) that apply when ν is time independent, in which case it is sufficient to consider the dispersion relationship. For ν = 0, we saw that no instabilities were present provided g ± g c > 0. On the other hand, for σν < 0 a tachyonic instability is induced in the coupled fluctuations between the relative density and relative phase perturbations. Now we analyze the effects of introducing a periodic modulation ν = ν(t) = ν 0 + δ ω cos(ωt) . (5.11) As seen above, around the stationary background solutions, δ and δϑ are decoupled from the dynamics of δε and δϕ. Furthermore, only δε and δϕ feel the presence of ν. It is therefore sufficient to consider the perturbations δε and δϕ in isolation. For convenience, we introduce the dimensionless parameters . (5.12) Assuming the experiment can be treated as spatially periodic, we obtain the following equations for the Fourier modes (δζ k and δφ k ) in terms of the time coordinate τ = ωt: (5.13) We perform our analysis using (5.13), but to connect with the usual discussions of parametric resonance in field theory, we note that this can be rewritten as the following damped wave equation for the phase alone: The equations of motion (5.13) take the forṁ a periodic operator with period T = 2πω −1 . From Floquet theory, solutions take the form where p(t) is a real function satisfying p(t + 2T ) = p(t). The complex constants Λ are known as Floquet exponents, with the solution y satisfying (5.21) known as the corresponding Floquet mode. For (5.13), y ∈ R 2 and Tr L = 0. It follows that for each choice of parameters k s , δ,ν s ,ω s there are two Floquet exponents satisfying Λ 1 = −Λ 2 , with the two corresponding Floquet modes forming a basis for solutions to the equation. Furthermore, sinceL is real, the Floquet exponents must be either purely real or purely imaginary. If (Λ) > 0, the mode grows exponentially in time, and if the initial condition has a nonzero projection onto this mode, the solution will grow exponentially in time. The exception to the simple form (5.21) occurs for degenerate Floquet exponents, where an additional power law prefactor is obtained. However, these degenerate modes generally form a measure zero subset in the space of parameters.
Since we want to study the stability of the proposed false vacuum, we focus on the Floquet exponent with the largest positive real component, which we denote With our convention (5.21), Λ is dimensionless and normalized so that after one period of L (i.e. at t = 2πω −1 ), the corresponding Floquet mode has grown by a factor e (µ) and acquired a phase e i (Λ) . In order to connect to the either the time-scales of the internal BEC dynamics or the effective relativistic scalar field, we will often find it convenient to work instead withω s Λ max orν −1/2 sωs Λ max respectively.

Analytic predictions of time-averaged analysis
The detailed instability structure of the Floquet modes as a function of experimental parameters arises from a complex interplay between the timescale of the external driver ω −1 and the timescale for the undriven dynamics of the modes. For sufficiently long-wavelengths, the internal time-scale for evolution of the modes will be much longer than ω −1 . In this regime, we therefore expect the time-averaging analysis will provide a suitable approximation to the dynamics, and we can use it to derive analytic predictions for the Floquet exponents. For shorter wavelengths, this separation of scales will no longer hold, and we expect the time-averaged analysis to fail, an intuition which we will confirm below. Similarly, for sufficiently large driving amplitudes, we expect that the dynamics of the external driver will become dominant, again invalidating the time-averaging analysis.
We begin with the time-averaged effective Lagrangian for linear perturbations (3.19) in the limit k → 0, and define a new field perturbation where λ andν s are defined in (3.16) and (5.12), respectively. The corresponding potential for this canonically normalized field is The equation of motion for δϕ c (k = 0) is We immediately read off the prediction for the Floquet exponents: If we take the limitν s 1 at fixed λ (which here is equivalent to dropping the non-canonical corrections to the kinetic term), we obtaiñ Below, we compare these predictions of the time-averaging formalism with the full solutions to the linear perturbation equations. We will see that the time-averaged Hamiltonian analysis accurately reproduces the results of a full Floquet analysis of the long-wavelength modes as long asν s 1 andω s 5 (in a way that will be made more precise below). Further, it correctly captures the suppression of the zero-mode instability for increasingν s . However, it does not perfectly capture the shift in the stabilizing value of λ that occurs asν s is increased. This is not surprising, since the calculation of the effective Hamiltonian proceeds through an expansion in δ as explained in section 3. Whenν s ∼ 1, we have λ 2 ∼ δ 2 , and thus expect corrections in the intercept which occurs at λ 2 ∼ 1.

Floquet Analysis: δ = 0
With the basic formalism in place, we now compute the Floquet exponents for the linear system (5.13). Examining the equations of motion, we see that they depend onν s andk s only through the combination Therefore, we first explore the linear instabilities of the system in terms of the three parameters κ 2 ,ω s , and σδ. The results for the maximal Floquet exponent as function of κ 2 and σδ at various fixed choices ofω s are given in figure 2 and figure 3. In a given experimental realization of the BEC dynamics, the values of σδ will be fixed, while κ 2 ∈ σν s ,k 2 max 4 + σν s , wherek max is the dimensionless wavenumber above which the GPE description begins to break down. The inhomogeneous modes properly captured by our analysis would thus live on a vertical line in the Floquet charts ( figure 3). By increasing δ, the vertical line shifts to the right. Increasingν s either shifts the line down (if σν s < 0) or shifts the line up (if σν s > 0). If this vertical line intersects a region of instability, we expect the false vacuum decay interpretation of the experiment will be spoiled.
Examining (5.13), one can show the instability charts described in the previous paragraph possess the following symmetry properties for a fixed value ofω s : the instability charts for κ 2 > − 1 2 , σδ ≥ 0, andω s ≥ 0. We note that the last condition is the physically meaningful one, as long as g − g c > 0. Some representative examples for various choices ofω s are shown in figure 3, where we plot Λ max for a range of values of κ 2 and σδ. As is typical for Floquet theory, there exist a series of instability bands emerging from the |δ| = 0 axis interspersed with regions of stability. The behavior of some representative mode functions are shown in figure 4.
Examining the various charts in figure 3, we see they share a number of common morphological features at each value ofω s > 0. The first important feature is the large instability band centered at κ 2 = − 1 2 and δ = 0 with width ∆κ 2 = 1 at δ = 0. As seen from the top left panel of figure 4, the modes within this band experience rapid tachyonic growth with a small superimposed oscillation. Physically, the modes within this large tachyonic region are well approximated by the time-averaged Hamiltonian analysis presented section 3. In addition to this large tachyonic band, there are additional instability bands that appear at larger values of either κ 2 − 1 2 or |σδ|. This rich instability structure is not properly captured (even qualitatively) by the time-averaged analysis, since the bare frequencies ω k,rel of these modes are comparable to the driving frequency ω. The width of these higher order bands vanishes as |δ| → 0, which is the expected behaviour since δ parameterizes the strength of the periodic perturbations. These higher-order bands are separated by regions of stability, where Λ max = 0. For the range of κ 2 and σδ values in the plot, the modes within the non-tachyonic instability bands experience weak parametric resonance, undergoing oscillations with an exponentially growing envelope as seen in the upper right panel of figure 4. With our definition of κ 2 these higher order bands move towards the origin and increase in width asω s is decreased.
Alternatively, if we instead parameterize our equations withκ 2 =ω −1 s κ 2 , the locations of the higher order bands are relatively stable as we varyω s . However, the vertical height of the tachyonic band now shrinks asω s increases, with a width 2ω −1 s centered atκ 2 = −ω −1 s .
For |δ| 1, the center of the first non-tachyonic instability band follows from energymomentum conservation for fluctuations produced by a homogeneous source with energy ω [34]. To leading order in the fluctuation amplitude, the lowest order process in coupling constants produces a pair of modes with equal and opposite wavenumbers. At this order the modes are produced on shell and the external coupling acts as a source, so energy conservation yields ω = 2 ω k,rel = 2 k 2 2m + 2σν 2 k 2 2m + 2σν + 2(g − g c )n . from which we see it is convenient to work with the variableκ 2 = 2ω −1 s κ 2 when considering the higher-order instability bands. In the limitsω s 1 andω s 1 we have, respectively We illustrate this in figure 5, where we see that (5.31) coincides with the center of the first instability band. Furthermore, for the rangeω s 1 relevant to analog false vacuum decay, the asymptotic limit provides an accurate approximation. Similar considerations give the locations of the higher order resonance bands.
We now estimate the driving frequency ω required to push this instability band out of the regime of validity of the GPE, assuming for simplicity g c = 0. Let us denote the Here m pr is the proton mass and a 0 the Bohr radius. We have assumed our atoms have nucleon number Z so m ≈ Zm p and the S-wave scattering length is given by a s = ra 0 . To make a very conservative estimate, we should require that the mode at the center of the instability band k floq satisfies k floq = k atom . Denoting this frequency by ω atom we immediately find We can expect values of Zr 2 ∼ 10 3 − 10 4 , and therefore the necessary driving frequency is in the GHz. This is somewhat problematic, as hyperfine transitions are associated with  energies of order MHz. A less conservative estimate is to require that the instability is at least pushed above the wave number associated with the healing length, which should require somewhat lower modulating frequencies. By definitionk heal = 1, so we obtain This is parametrically smaller by the factorna 3 s . For example, for densities of order 10 20 m −3 and S-wave scattering lengths a s ∼ 10a 0 ∼ 1nm, this yields modulating frequencies a factor ∼ 10 −3 times smaller than ω atom , on the order of 100 kHz to MHz. This is a more consistent modulation for radio frequencies associated with energies of hyperfine transitions. In conclusion, it appears challenging to push the first instability band all the way to k atom .
Having obtained the stability charts for the parameters κ 2 ,ω s , and σδ; and briefly commented on their interpretation for the stability of the BECs, we now quantitatively relate the stability of various modes to the parameters describing the symmetric BECs, namelyk s ,ν s ,ω s , and δ.
In the context of an analog false vacuum decay experiment, the motivation for introducing the time-dependent oscillations into ν is to provide a stabilization mechanism for the unstable fixed point of the background dynamics withφ = π. Therefore, we first demonstrate that this long-wavelength stabilization mechanism operates by studying the dynamics of the k = 0 mode. This also allows us to make a quantitative comparison with the predictions of the time-averaged Hamiltonian analysis, from which the emergent false vacua were derived. In figure 6 we plot Λ max (k = 0) as a function of λ for a range ofν s values withω s held fixed and vice versa. For comparison, we also plot the prediction (5.27) as dashed lines (left panel) or dots (right panel). From the left panel, we see forν s 1 the leading order time-averaging analysis provides an accurate description of the linear instability of the k = 0 mode. However, for sufficiently largeν s , we see strong corrections to this relationship emerge. From the right panel, we see that for |ν s | 1 the shape of the lowest instability band is essentially independent ofω s (when scaled appropriately). However, as we takeω s → 1 additional instabilities present at larger values of λ move towards smaller values of the stabilization parameter λ. In figure 7 we further illustrate the convergence properties of the solution ofω s is varied. We see that even for the large valueν s = 0.5, the zero-mode dynamics has converged forω s 5, indicating that a time-averaged analysis of the long-wavelength modes is sufficient. Figure 6 demonstrates the stabilisation of the homogeneous degrees of freedom, as   predicted by a time-averaged analysis for the zero mode in section 3. However, there are additional fluctuations with nontrivial spatial structure whose stability must also be considered. The internal time-scale for the dynamics of these modes decreases with wavenumber, and for sufficiently short wavelengths the time-averaged analysis is insufficient. In figure 8 we show the stability of the k = 0 modes as we vary the parameters of the system. Since our interest in this paper is the case where the homogeneous modes have been stabilized, we choose parameters such that λ = 1.3 as a representative example. Even if the zero mode is stabilized, if the range ofk s that are correctly described by the GPE happen to intersect one of the higher order Floquet bands, then the system as a whole will be unstable to a statistically homogeneous amplification of inhomogeneities. Of course, if the time-scale for these instabilities is much longer than the lifetime of the experiment, then they may be neglected as subleading corrections on the experimental dynamics. From figure 8 we see that bands of higher k modes do indeed experience exponential instabilities. The location of the bands with respect tok s has very little sensitivity toν s or δ (i.e. λ). The primary effect of increasingν s is to increase the width of the instability bands and make the higher bands more prominent. Meanwhile, the primary effect of increasing δ (or λ) is to increase the peak strength of the instability bands. Therefore, the net effect of increasing either of these parameters will be to reduce the time-scale for any unstable Floquet modes present in the system to destabilize the condensate. However, they cannot be adjusted in such a way as to remove a pre-existing band of unstable modes from the system. By contrast, in figure 9, we see that the primary effect of increasingω s is to push the instabilities to higher wavenumbers. In addition, the width of the bands with respect tõ k s narrows, while roughly maintaining its peak amplitude. This second set of properties is robust as long asν s 1. We thus identifyω s as the key control parameter to push unstable inhomogeneous modes to shorter wavelengths, where the GPE description eventually fails and we need additional physics to stabilize them to avoid spoiling the false vacuum analog. This is, of course, entirely expected sinceω −1 s measures the ratio of timescales associated with the internal scattering of the individual atoms to that of the external modulation. By increasingω s , these scattering interactions must occur on much shorter timescales to respond toω s , corresponding to much larger wavenumbers. 1, the location of the instability band (ink s ) is essentially independent of both λ andν s . In these units, increasing λ both the stabilizes the zero mode and increases the strength of the higher order Floquet bands. Meanwhile, increasingν s primarily increases the width of the instability band, while slightly decreasing its amplitude.

Conclusions and Future Work
We undertook a thorough analysis of the analog relativistic false vacuum decay proposal of Fialko et al [14,15]. If experimentally feasible, this provides the exciting possibility of testing strongly nonlinear and nonperturbative quantum field dynamics in a fully timedependent setting. The specific phenomena under consideration, relativistic false vacuum decay, has direct application to multiverse models of cosmology. Due to the exponential complexity of quantum field theory, no tools exist to treat the decay problem exactly. Furthermore, extant approximations do not rigorously incorporate the time-dependent nature of the decay, nor can they be used to explore correlations between tunnelling events. When applied to cosmology, these correlations can significantly influence the distribution of fundamental constants seen by observers in the multiverse. Furthermore, they can have important implications for the probability of observing signatures of the multiverse through bubble collisions. Reaching a deeper understanding of these issues can thus shed light on some of the most fundamental issues in modern cosmology. In order to assess the feasibility of the proposal, we first carefully derived the relationship between the theoretical description of the BEC and the emergent relativistic scalar field with a false vacuum potential. We cataloged the various corrections to the simple interpretation as an effective canonical scalar field, and indicated the required conditions for these corrections to be small. Since the only aspect of the BEC physics used in the derivation was the applicability of the GPE, this illustrates the generality of the mechanism beyond any specific experimental proposal. We then explored the linear stability properties of the proposed false vacua using Floquet theory. The introduction of external modulations in the interspecies conversion rate (which was the key ingredient to generate effective false vacuum minima) simultaneously induces new bands of instabilities in the shorter wavelengths. Furthermore, for a sufficiently strong modulation amplitude, the long-wavelength modes can again be destabilized.
The key control parameter determining the wavelengths of the newly destabilized modes is the frequency of the modulations in the interspecies conversion rateω s . We demonstrated that the corresponding instability wavenumber scales as √ω s for the relevant regimeω s 1. Since the GPE description does not hold to arbitrarily short wave-lengths, it is possible to push the higher order Floquet instabilities out of the regime of validity of the coupled GPEs studied here. Physically, we expect the discreteness of the particles will induce viscosity and stochasticity into the equations at short wavelengths, which will act to damp the exponential growth of the modes. In the pessimistic scenario where modes up to k atom obey the GPE, driving frequencies of order of tens of GHz are required. In the extremely optimistic scenario where modes beyond the healing length are damped, frequencies in the 100 kHz range are instead needed. However, one must also ensure that the increased modulation frequency does not allow additional modes of the system to be excited, such as the radial eigenmodes of a ring trap. For the particular setup described in [15] based on radio-frequency fields, ensuring the required frequency is below the MHz scale requires that the linear dynamics of modes a factor of a few beyond the healing length must be modified from the GPE description. Whether this condition holds in reality requires further investigation. In order to illustrate the key dynamics with minimal technical complications, we made two key assumptions in our analysis: that the background condensate is spatially homogeneous and that the two particle species are symmetric (with equal masses and interspecies scattering rates). On technical grounds, the first assumption allows for a diagonalization of the fluctuation modes through the Fourier transform, explicitly decoupling the equations of motion and transforming a PDE into a collection of uncoupled ODEs. Physically, one could imagine the condensates are trapped in a ring with strong trapping potentials in only two directions. In this case the effective one-dimensional dynamics should be well approximated by the spatially homogeneous and periodic setup. However, one needs to be careful not to excite the radial eigenmodes of the system through modulations in the transition rate. Alternatively, accounting for the inhomogeneity of the background induced by the external trap, one can imagine scenarios where there is a narrow boundary layer near the edges of the trap, and then a nearly homogeneous plateau in the middle. Our analysis should then account for the dynamics of perturbations with wavelengths much smaller than the width of the plateau. An interesting extension is to include the dynamical effects of the trapping potentials for the condensate, which will introduce inhomogeneity into the background. Assuming there is still a false vacuum interpretation in this case, the fluctuations on this inhomogeneous background can then be studied using the methods developed in [35].
The symmetric assumption (i.e. g 11 = g 22 and m 1 = m 2 ) along withε = 0 was primarily adopted for technical simplicity and to create the most direct connection between the cold atom system and the relativistic false vacuum decay system. Relaxing these conditions generally leads to a mixing between the total and relative phonons, so that the effective field theory description describes two coupled relativistic fields. Further, the sound speeds between the modes need not be the same, so the relativistic interpretation can be spoiled. From the viewpoint of linear stability, the coupling between the modes can lead to the emergence of a second unstable mode, modify the stabilization condition for the zero mode, and shift the value ofk s at which the first non-tachyonic instability band emerges. The rich resulting dynamics will be explored in a future publication.
As well as these extensions, an interesting possibility is to consider the interactions between many condensates, rather than just two. This could allow for an experimental realization of a complex potential landscape, analogous to those believed to descend from low-energy string theory constructions and that form a fundamental ingredient in many current multiverse paradigms for cosmology. We do not expect the linear stability analysis to differ dramatically with the addition of additional condensates. On the contrary, the nonlinear dynamics resulting from the decay of the false vacuum may display considerably more variability with the introduction of additional decay and entropy production channels.
Since this paper focuses solely on the linear stability analysis of the "UV complete" theory of the Bose-Condensates, we leave these interesting possibilities to future work.

A Effective Action for Condensate Phases
In this appendix we derive the effective action for the total and relative phases under the assumption of symmetric condensates (i.e. m 1 = m 2 , g 11 = g 22 ) and a homogeneous background around which to expand our fluctuations. Furthermore, we assume that we can fix the dynamical background variableε = 0. This is a common approximation in the literature, but we note that a full consideration of the field dynamics may break the assumption ε = 0. We explicitly indicate the required approximations to arrive at an action describing a relativistic field with canonical kinetic term, as well as provide some arguments for why these approximations may hold. A fully rigorous justification requires tracking the dynamics of the homogeneous background, including the asymmetry between the background condensate densitiesε. However, even with the restriction to symmetric condensates, the technical details are already rather cumbersome, so we leave an investigation of the effects of introducing asymmetries between the condensates to future work.
We first split the density fields into a homogeneous background plus an inhomogeneous fluctuation: As written, this does not uniquely define the perturbations due to ambiguity in the k = 0 mode amplitude. We fix this ambiguity by forcing δ = δε = 0, where · represents a spatial average. However, we note that in some cases it is convenient to demand that n andε (along with the corresponding averages for the angles) follow the homogeneous background equations of motion as presented in section 4. The first step is to expand the Lagrangian density to quadratic order in the density perturbations δ and δε. Since we are interested in the nonlinear dynamics of the phase variables, we do not expand them in fluctuations. Denoting the ith order of this expansion by L (i) δ , a straightforward calculation gives At the next order in the expansion, we obtain terms schematically of the form where i = 1, 2 labels the individual condensate species and δρ represents fluctuations in either condensate. Smallness of δε and δ thus ensures that we are keeping the leading terms involving derivatives of the density perturbations. In cases wherenδφ δρ, (with δφ representing fluctuations in any of the phase variables) consistency also requires us to drop the terms involving two derivatives of the phase variables in L (1) δ . An example of such a state is one where the condensate fluctuations δψ i are complex Gaussian random fields. For general states this condition need not hold. However, comparing to the O(∇ 2 ) terms in L (0) δ , we see that the contributions will be suppressed by the smallness of the density perturbations. In the interests of simplicity, we will therefore neglect the O(∇ 2 ) contributions to L (1) δ . Before proceeding, we also divide the total phase into a homogeneous background and perturbation ϑ(x, t) =θ(t) + δϑ(x, t) . (A.4) As with the density perturbations, it may be convenient to define the backgroundθ via either removal of the tadpole δϑ = 0, or to defineθ(t) to be a solution of the homogeneous background equations. We now introduce the quantity E[θ,φ,n] ≡ − θ − 2(g + g c )n + 2ν cosφ + 2µ (A.5) which is the homogeneous background equation of motion forθ whenε = 0. By integrating the Gaussian density perturbations δ and δε, we obtain an effective Lagrangian for the phase variables L ϑ,ϕ eff = − nθ − (g + g c )n 2 + 2µn + 2νn cosφ + 1 4(g + g c ) E 2 (∇δϑ) 2 − nδθ (∇ϕ) 2 + 2νn (cos ϕ − cosφ) − δθ + 2ν(cos ϕ − cosφ) Here we have ignored spatial dependence of the functional determinant to put it back inside the integral. Since we will not consider the effects of this term in the remainder of the paper, this approximation has no impact on our results. We also assume we are working in an experimental regime where ν/(g − g c )n 1, and restrict ourselves to the long-wavelength modes with 2 k 2 /4(g ± g c )mn 1. The effective action (A.6) then simplifies L ϑ,ϕ eff = − nθ − (g + g c )n 2 + 2µn + 2νn cosφ + 1 4(g + g c ) E 2 + 2 4(g + g c ) δθ 2 − 2n 4m (∇δϑ) 2 − nδθ + 2 4(g − g c )φ 2 − 2n 4m (∇ϕ) 2 + 2νn (cos ϕ − cosφ) (A.7) + 2E 4(g + g c ) − δθ + 2ν(cos ϕ − cosφ) + ν 2 (g + g c ) (cos ϕ − cosφ) 2 − ν (g + g c ) δθ (cos ϕ − cosφ) .
Finally, if we ignore the backreaction of the fluctuations on the homogeneous dynamics, then the barred quantities will follow the equations of motion for a homogeneous background as presented in section 4. From this, we see that E = 0 on background solutions withε = 0, so we are left with L ϑ,ϕ eff = − nθ − (g + g c )n 2 + 2µn + 2νn cosφ δθ (cos ϕ − cosφ) (∇ϕ) 2 + 2νn (cos ϕ − cosφ) (A.8) Since we have assumed the fluctuations do not backreact or rescatter off the background, the first line of (A.8) acts as an irrelevant constant for the dynamics of the perturbations. Under the sequence of approximations given above, we end up with an effective Lagrangian for the phase variables that is quadratic in the total phase perturbations δϑ. Therefore, we can also integrate out the perturbations δϑ to obtain an effective action for the relative phase variable ϕ alone. This is most easily done by going to Euclidean time τ E = it, and expanding δϑ is Matsubara modes. We will neglect these additional environmental contributions here, but we mention in passing a subtlety that appears if the Euclidean time average of ∂ τ E δϑ is non-vanishing. In this case, the zero mode of δϑ is not representable as part of the τ E -periodic Matsubara sum and must be treated separately, generating additional environmental contributions. Truncating the ϑ dynamics and working to leading order in ν/(g + g c )n, we obtain the effective sine-Gordon Lagrangian (3.6) presented in section 3.

A.1 Effective Action for Linear Perturbations
To directly connect to the linear fluctuation analysis presented in this paper, it is also of interest to obtain the effective Lagrangian for linear perturbations to the phase variables around our fiducial background false vacuum states. We demonstrate that under suitable approximations, the effective action obtained this way produces equivalent results to direct manipulation of the equations of motion. As above, due to the complexity of the resulting equations, we will restrict ourselves to the case of symmetric condensates here. Furthermore, we assume the background is fixed, with cosφ = ±1 ≡ σ.
Since we now wish to concentrate on the equations for linear perturbations around a homogeneous background, we write ϑ(x, t) =θ(t) + δϑ(x, t) (A.9a) ϕ(x, t) =φ(t) + δϕ(x, t) (A.9b) with d d xδϕ = 0 = d d xδϑ. Expanding the Lagrangian density to quadratic order in δϑ and δϕ (in addition to δ and δε), we see that the O(∇ 2 ) contributions to L (1) δ disappear, and that L (2) δ becomes independent of the phase fluctuations. Therefore, we can begin directly with (A.6), except with the fluctuation determinant on the fourth line set to a constant. Further, we perturb around a homogeneous background solution so E = 0. Expanding (A.6) to quadratic order in δϑ and δϕ we see the dynamics of the total and relative phonon modes decouple. We are left with the following pair of effective Lagrangian densities 2(g + g c )n 2n L ϑ lin = 1 2 δθ 1 + 2 ∇ T · ∇ 4mn(g + g c ) −1 δθ −n (g + g c ) m (∇δϑ) 2 2 (A.10) and 2(g − g c )n 2n L ϕ lin = In (A.10) we have dropped a total time derivative. Using the Euler-Lagrange equations on each of these effective Lagrangians, we obtain δθ −n (g + g c ) m 1 − 2 ∇ 2 4mn(g + g c ) ∇ 2 δϑ = 0 (A.12) and d dt in agreement with (5.5) and (5.10) respectively.

(B.1d)
We now convert to our preferred variables, the relative and total phonon modes. To simplify some expressions, we recall the dimensionless fractional difference in the background densities defined in section 4,