Nonlinear Dynamics of the Cold Atom Analog False Vacuum

TL;DR

This paper explores the nonlinear dynamics of cold atom systems as quantum simulators for false vacuum decay, revealing how unstable modes influence the system's evolution and phase transition behavior.

Contribution

It provides a numerical analysis of the nonlinear evolution of unstable modes in a cold atom analog false vacuum, and demonstrates how coupling modulation can induce different phase transitions.

Findings

Unstable short-wavelength modes threaten false vacuum stability.

Short-scale modifications to the GPE can stabilize the false vacuum.

Modulating coupling strength induces crossover between phase transition types.

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 two-species Bose-Einstein 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 short-wavelength modes driven by the time-dependent coupling. These modes threaten to destabilize the analog false vacuum. Here, we employ numerical simulations of the coupled Gross-Pitaevski equations (GPEs) to determine the non-linear 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 time-dependent coupling, it is possible to observe the crossover between a second and first order phase transition out of the false vacuum.

Equations

$$\mathcal{H}=\frac{{\hslash }^2}{2m}{\left|\nabla {\psi }_i\right|}^2+\frac{g_{ij}}{2}{\left|{\psi }_i\right|}^2{\left|{\psi }_j\right|}^2-\frac{\nu }{2}{\sigma }_{ij}^{\mathrm{x}}\left({\psi }_i{\psi }_j^{\ast }+{\psi }_i^{\ast }{\psi }_j\right).$$

$$\sigma^{\rm\bf x}=\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right)\,.$$

$$i\hslash {\dot{\psi }}_1$$

$$i\hslash {\dot{\psi }}_2$$

$$\nu (t)={\nu }_0+\delta \hslash \omega \cos (\omega t).$$

$${\psi }_i=\sqrt{{\rho }_i}{e}^{i{\varphi }_i}\text{mbox}and{\psi }_i^{\ast }=\sqrt{{\rho }_i}{e}^{-i{\varphi }_i}.$$

$$\hslash {\dot{\rho }}_i$$

$$\hslash {\dot{\varphi }}_i$$

$$\vartheta$$

$$\phi$$

$$\varrho (x,t)$$

$$\vartheta (x,t)$$

$$\epsilon (x,t)$$

$$\phi (x,t)$$

$$\ddot{\phi }-{\nabla }^2\phi +\frac{4{\nu }_0\bar{n}\left(g_{\mathrm{s}}-g_{\mathrm{c}}\right)\left(1-2{\delta }^2\right)}{{\hslash }^2}\left(\sin \phi +\frac{{\lambda }^2}{2}{\sin }^2(\phi )\right)=0,$$

$$V(\phi )=\frac{4{\nu }_0\bar{n}\left(g_{\mathrm{s}}-g_{\mathrm{c}}\right)\left(1-2{\delta }^2\right)}{{\hslash }^2}\left(-\cos \phi +\frac{{\lambda }^2}{2}{\sin }^2\phi \right).$$

$$\bar{n}=\frac{\bar{\varrho }}{2},$$

$${\lambda }^2=\frac{2{\delta }^2(g_{\mathrm{s}}-g_{\mathrm{c}})\bar{n}}{{\nu }_0},$$

$$\frac{{\hslash }^2{k}^2}{\bar{n}m(g_{\mathrm{s}}-g_{\mathrm{c}})}\approx 2\left(\sqrt{1+\frac{{\hslash }^2{\omega }^2}{4{\bar{n}}^2(g_{\mathrm{s}}-g_{\mathrm{c}}{)}^2}}-1\right)-4\sigma \frac{{\nu }_0}{\bar{n}(g_{\mathrm{s}}-g_{\mathrm{c}})},$$

$$k_{\mathrm{heal}}=\frac{\sqrt{g_{\mathrm{s}}\bar{n}m}}{\hslash },k_{\mathrm{atom}}=\frac{8{\pi }^2{\hslash }^2}{g_{\mathrm{s}}m}.$$

$$\lambda =\delta \sqrt{\frac{2g_{\mathrm{s}}\bar{n}}{{\nu }_0}}=1.3,\frac{{\nu }_0}{g_{\mathrm{s}}\bar{n}}=2\times 10^{-3},\frac{\hslash \omega }{g_{\mathrm{s}}\bar{n}}=100\sqrt{\frac{{\nu }_0}{g_{\mathrm{s}}\bar{n}}},$$

$$\frac{\hslash \bar{n}}{2\sqrt{m{\nu }_0}}=10^3,$$

$$V_{\mathrm{eff}}(\phi )=V_0\left(-\cos \phi +\frac{{\lambda }^2}{2}{\sin }^2\phi \right),$$

$$\lambda =\delta \sqrt{\frac{2g_{\mathrm{s}}\bar{n}}{{\nu }_0}}$$

$${\psi }_i(\mathbf{x})={\bar{\psi }}_i+\delta {\hat{\psi }}_i(\mathbf{x}),$$

$$\delta {\hat{\psi }}_i=\frac{1}{\sqrt{V}}\sum _k\sqrt{{\mathcal{P}}_i(k)}{\hat{a}}_k{e}^{i\mathbf{k}\cdot \mathbf{x}},$$

$${\hat{a}}_k=\sqrt{-\ln ({\hat{A}}_k)}{e}^{2\pi i{\hat{P}}_k},$$

$$⟨\delta {\hat{\psi }}_i(\mathbf{k})\delta {\hat{\psi }}_j({\mathbf{k}}^{\prime }){⟩}_{\mathrm{Q}}={\mathcal{P}}_i(k){\delta }_{ij}\delta (\mathbf{k}-{\mathbf{k}}^{\prime }).$$

$${\nabla }^{2}f(x_i)\to L[f](x_i)\equiv \frac{1}{2d{x}^2}\left(f(x_{i+1})-2f(x_i)+f(x_{i-1})\right)$$

$$(\nabla f(x_i){)}^2\to G[f](x_i)\equiv \frac{1}{d{x}^2}\left[{\left(f(x_{i+1})-f(x_i)\right)}^2+{\left(f(x_{i-1})-f(x_i)\right)}^2\right]$$

$$\bar{t}={\omega }_{\mathrm{p}}t,\bar{x}={\kappa }_{\mathrm{p}}x,\mathrm{and}\bar{\psi }=\sqrt{\Lambda }\psi ,$$

$$c_{\mathrm{p}}\equiv \frac{{\omega }_{\mathrm{p}}}{{\kappa }_{\mathrm{p}}}\mathrm{and}{E}_{\mathrm{p}}=\hslash {\omega }_{\mathrm{p}}.$$

$$\bar{H}=\frac{1}{\Lambda {\kappa }_{\mathrm{p}}^d}\int {\mathrm{d}}^d\bar{x}\frac{1}{\hslash {\omega }_{\mathrm{p}}}\left(\frac{(\hslash {\kappa }_{\mathrm{p}}{)}^2}{2m_i}{\left|\bar{\nabla }{\bar{\psi }}_i\right|}^2+\frac{g_{ij}}{2\Lambda }{\left|{\bar{\psi }}_i\right|}^2{\left|{\bar{\psi }}_j\right|}^2-\frac{{\nu }_{ij}}{2}\left[{\bar{\psi }}_i{\bar{\psi }}_j^{\ast }+{\bar{\psi }}_i^{\ast }{\bar{\psi }}_j\right]\right).$$

Full text

1 Introduction

False vacuum decay, the quantum first-order 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 $\hbar\neq 0$), 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 tunneling 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 real-time 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 real-time 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 two-component dilute gas cold atom Bose-Einstein 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 long-wavelength 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 short-wavelength modes, leaving only the long-wavelength 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 short-wavelength 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 short-wavelength dynamics of the condensate. Since modifications to the GPE at short-wavelengths 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 Bose-Einstein condensates, modeled by the Gross-Pitaevskii 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 second-order phase transition into a first-order 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 semi-classical 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 finite-differencing 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 2-Component Cold Atom Bose-Einstein Condensates in the Dilute Gas Limit

The analog false vacuum system under investigation is a 2-component BEC. There are currently two different configurations that can in principle be utilized as false vacuum decay simulators. The first one is a single-species BEC in a double-well potential in the so-called tight-binding regime, exhibiting symmetric and anti-symmetric single-particle states Albiez:2005 ; Bucker:2011 . False vacuum decay requires a modulation of the tunnel-coupling via the trapping parameters. At present, this configuration is capable of mimicking false vacuum decay in $1$+$1$-dimensions. To avoid dimensional restrictions altogether, we focus on a different configuration: a two species BEC of atoms with mass $m$, e.g.,  41K or 87Rb 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 ring-trap, which we model as a $1$+$1$-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 double-well 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 ultra-cold highly diluted 2-component system consisting of two single-particle states ($|1\rangle$ and $|2\rangle$) of mass $m_{1}=m_{2}=m$. Within this setup the atom-atom interactions can be approximated by $S$-wave contact potentials: $g_{11}$, $g_{22}$ and $g_{12}=g_{21}$. 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 $\nu$ through the utilization of external fields, e.g.,  a radio-frequency field. When the BEC forms, the collective long-wavelength excitations of the condensed atoms (i.e.,  the condensate) acquire a vacuum expectation value (vev), which we denote by the complex numbers $\psi_{i}$. The Hamiltonian density for the two interacting condensates is given by

[TABLE]

Here we are dropping the effects of the external confining potentials $V_{\rm ext,i}$ on the dynamics, and we have defined

[TABLE]

The combined evolution of the 2-component system is given by

[TABLE]

where we have assumed the inter-atomic scattering strengths are of the form $g_{11}=g_{22}=g_{\rm s}$ and defined $g_{12}=g_{\rm c}$. 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 $\psi_{i}$. 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 $g_{11}=g_{22}=g_{\rm s}$ and $g_{\rm c}=0$, which we refer to as a symmetric theory. To generate an effective false vacuum potential, we further consider a periodically modulated coupling,

[TABLE]

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 $\rho_{i}$ and phases $\phi_{i}$, such that

[TABLE]

In terms of these, the coupled GPEs (3) become

[TABLE]

To map the coupled 2-component 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

[TABLE]

respectively. We will refer to the pairs of canonically conjugate variables $(\vartheta,\varrho)$ and $(\varphi,\varepsilon)$ 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 $\bar{\cdot}$. Potential vacua are identified by looking for eigenstates of the Hamiltonian, with the eigenvalue ($\mu_{\rm bg}$) playing the role of the chemical potential and sets $\bar{\vartheta}=-\frac{\mu_{\rm bg}}{\hbar}t$.111Including fluctuations changes the time evolution of the spatial average of the total phase $\langle\vartheta\rangle_{\rm V}$ compared to the evolution in the homogeneous background, so that $\frac{{\rm d}}{{\rm d}t}\langle\vartheta\rangle_{\rm V}\neq\frac{{\rm d}}{{\rm d}t}\bar{\vartheta}=-\frac{\mu_{\rm bg}}{\hbar}$. 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 $g_{11}=g_{22}$, we find stationary solutions when $\cos\bar{\varphi}=\pm 1$, with the $\cos\bar{\varphi}=-1$ 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:

[TABLE]

where an overbar $\bar{\cdot}$ indicates a solution to the coupled GPEs in the homogeneous limit.222Because of the nonlinear relationship between the Cartesian condensate variables ($\psi_{i}$) and the density ($\rho_{i}$) and phase variables ($\phi_{i}$), the volume averaged density and the background density differ by a UV divergent quantity $\langle\varrho\rangle_{\rm V}=\sum_{i}\left|\bar{\psi}_{i}\right|^{2}+\langle\left|\delta\psi_{i}\right|^{2}\rangle_{\rm V}\neq\sum_{i}\left|\bar{\psi}_{i}\right|^{2}$. 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; 2. 2.

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

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

the short-time dynamics of the fluctuation modes are simple and can be integrated out to obtain an effective time-averaged description; and 5. 5.

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 $g_{11}=g_{22}$ expanded around a background state of equal densities in the two condensates. This latter condition holds for the stationary solutions with $\cos\bar{\varphi}=\pm 1$. 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

[TABLE]

which is the equation of motion for a field with potential

[TABLE]

In the above we have defined a mean particle density

[TABLE]

and the parameter

[TABLE]

which controls the shape of the effective potential, and can be experimentally tuned. For $\lambda=0$, this reduces to the sine-Gordon model, while for $\lambda>1$, the potential for the effective relativistic scalar develops a series of local minima, as illustrated in fig. 1.

Meanwhile, the Lagrangian has a global $U(1)$ 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 Gross-Pitaevskii 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.333Only $g_{11}=g_{22}$ is required to obtain equal condensate densities, not the additional requirement $g_{\rm c}=0$ 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 $g_{11}=g_{22}$ 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 well-approximated 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:

1. A.

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

all relevant dynamical condensate modes only feel the time-averaged effects of the $\nu$ 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 $\delta\ll 1$, the introduction of periodic modulation in $\nu$ of frequency $\omega$ induces a band of exponentially unstable modes with wavenumber centered at

[TABLE]

where $\sigma=\cos{\bar{\varphi}}=\pm 1$. 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), 2. 2.

were subjected to a periodic modulation in $\nu$, 3. 3.

could be treated as nearly spatially homogeneous, and 4. 4.

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 $S$-wave scattering approximation from atoms retaining memory of the internal structure of the interaction potential between collisions; 2. 2.

corrections to the short-wavelength dynamics of the condensate from the finite mode occupancy; 3. 3.

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

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

corrections to the one-dimensional 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

[TABLE]

The validity of the dilute gas BEC approximation requires $k_{\rm heal}\ll k_{\rm atom}$, 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. However, the discrete nature of the atoms provides a promising mechanism to remove the Floquet modes. The continuum description must break completely at $k_{\rm atom}$, and hence modes above this scale cannot exist within the condensate. As well, we expect that the finite mean free path of the atoms will induce additional viscous corrections at wavenumbers somewhat below $k_{\rm atom}$, which will preferentially damp large wavenumbers. Since $k_{\rm heal}$ encodes the scale above which the dispersion relationship changes to that of free particles, these additional viscous corrections may enter at scales not too far removed from the healing length. Another example of a relevant effect in another experimental setup can be found in ref.  2018arXiv180707564N , where modifications of the trap geometry are used to induce “synthetic dissipation” of high-k modes and the ability to tune the dissipation scale. Referring to (13), we also see that a degree of experimental tunability exists, since either increasing the modulation frequency $\omega$ or decreasing the condensate density $\bar{n}$ 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

[TABLE]

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

[TABLE]

again as suggested in ref.  Fialko:2014xba . The maximal Lyapunov exponent for fluctuations around the fiducial false vacuum as a function of dimensionless wavenumber $\bar{k}=\frac{\hbar k}{\sqrt{g_{\rm s}\bar{n}m}}$ are shown in fig. 2. We see a narrow band of unstable wavenumbers centered around $\frac{\hbar k}{\sqrt{g_{\rm s}\bar{n}m}}\approx 1.705$, 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 time-averaging 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 $dx$ (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 $\frac{\sqrt{g_{\rm s}nm}}{\hbar}dx\approx 1.89$ so that $\frac{\hbar k_{\rm nyq}}{\sqrt{g_{\rm s}nm}}\approx 1.66$, 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 time-averaging analysis to be valid. Indeed, we see the dynamical nucleation and subsequent expansion of domain wall-antiwall pairs (i.e.,  bubbles in one dimension) in the relative phase, indicating a first-order phase transition. By contrast, the right panel takes ${\frac{\sqrt{g_{\rm s}\bar{n}m}}{\hbar}dx\approx 0.12}$ and ${\frac{\hbar k_{\rm nyq}}{\sqrt{g_{\rm s}\bar{n}m}}\approx 26.53}$, which is above the upper edge of the instability band.444To ensure that the dynamics of the short-wavelength modes is properly resolved, we have chosen the lattice cutoff small enough to saturate the convergence in both $dt$ and $dx$ (see appendix D), and ensure the Noether charges are preserved to the $\mathcal{O}(10^{-15})$ level. To demonstrate the change in behavior occurs from effects localized in the Floquet band, we also ran a simulation with $\frac{\sqrt{g_{\rm s}nm}dx}{\hbar}\approx 1.79$ and $\frac{\hbar k_{\rm nyq}}{\sqrt{g_{\rm s}nm}}\approx 1.75$. This grid just resolves the full Floquet band, and is only a $\mathcal{O}(10)$% 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 $\frac{g_{\rm s}\bar{n}}{\hbar}t\sim 100$ the exponentially growing linear modes become $\mathcal{O}(1)$ and undergo strong mode-mode coupling. The evolution is subsequently dominated by short-wavelength modes of large amplitude, competely erasing the first-order 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 $\cos\varphi$ for the simulations in fig. 3. At the false vacuum $\cos\varphi_{\rm fv}=-1$, while in the true vacuum $\cos\varphi_{\rm tv}=1$, so that if the first-order phase transition interpretation holds, the expectation value should go from approximately $-1$ to $1$. Corrections to the pure $\pm 1$ 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 $\langle\cos\varphi\rangle_{\rm V}\approx 1$ is lost and instead we end up with $\langle\cos\varphi\rangle_{\rm V}\approx 0$. 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.555As mentioned above, the emergence of the relative phase as the key tunneling variable comes directly from the Hamiltonian term proportional to $\nu$, 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 $\frac{g_{\rm s}\bar{n}}{\hbar}t\sim 100$, $\mathcal{O}(1)$ 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 $\varrho$ fluctuations is driven by nonlinear scattering with the growing Floquet modes of $\varepsilon$. 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 $\frac{g_{\rm s}\bar{n}}{\hbar}t\sim 100$, consistent with partial quenching of the linear instability by nonlinear mode-mode 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. Although this presents an obstacle to successfully realizing an analog false vacuum decay experiment, it is important to remember that the coupled GPEs are an approximation to the full dynamics of the BECs. In an actual experiment, additional dynamics will be present to modify the GPE description. From the viewpoint of the successfully realizing false vacuum decay in the current proposal, we require new short-wavelength physics not captured by the GPE to enter at length scales longer than $\lambda_{\rm unstable}=2\pi k_{\rm unstable}^{-1}$. This new physics must remove the exponentially growing Floquet modes, while simultaneously exerting a minimal influence on the dynamics of the bubble-forming 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 long-wavelength condensate dynamics. A number of physical mechanisms can lead to these types of effective modifications to the coupled GPEs. Due to the finite interatomic spacing, the continuum description assumed by the GPEs fails completely for wavelengths of order this spacing, guaranteeing all shorter wavelength modes will effectively be projected from the system. Further, we expect the effects of this spacing and the finite mean free path of the atoms to induce viscous like corrections at somewhat longer wavelengths, similar to the case of a normal fluid, leading to the damping of high-wavenumber modes. Similarly, interactions between the condensate and the cloud of uncondensed particles will induce damping of condensate modes, which may help to counteract the Floquet instability. Finally, a particularly intriguing possibility is the induce synthetic dissipation into the system, allowing for a tunable cutoff in the GPE description 2018arXiv180707564N . 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 briefly indicated some physical mechanisms that may act to effectively remove the instabilities from the system. We left the complex issues of detailed modeling of these mechanisms, 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 first-order 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 $\frac{\hbar\omega}{g_{\rm s}\bar{n}}=3200\sqrt{\frac{\nu_{0}}{g_{\rm s}\bar{n}}}\approx 143.11$ 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 $\frac{\sqrt{g_{\rm s}\bar{n}m}}{\hbar}=\frac{25}{2048}\sqrt{\frac{g_{\rm s}\bar{n}}{\nu_{0}}}$, 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 $\nu$ creates an effective potential for the long-wavelength modes of $\varphi$

[TABLE]

illustrated in fig. 1, where the parameter

[TABLE]

is experimentally tunable. For $\lambda\leq 1$, $\varphi=\pi$ is a local maximum, while it is a local minimum for $\lambda>1$. In the (time-averaged) scalar field interpretation of $\varphi$, increasing $\lambda$ adjusts the field from sitting on top of a hill to sitting in a local minimum. For initial states localized around the false vacuum $\varphi=\pi$, we expect two distinct mechanisms by which the false vacuum state can decay. The value of $\lambda$ determines which mechanism is relevant.

For $\lambda<1$, modes with ${k\lesssim 2\frac{\sqrt{\nu_{0}g_{\rm s}\bar{n}}}{\hbar}\sqrt{1-\lambda^{2}}}$ will experience a tachyonic instability, We confirmed this at the level of linear perturbation theory in ref.  Braden:2017add , where we also numerically obtained $\frac{\nu_{0}}{g_{\rm s}\bar{n}}$ and $\delta$ dependent corrections. Based on this picture, we expect that for $\lambda<1$, $\varphi$ 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 $\varphi=0$ or $\varphi=2\pi$ vacuum) connected by a network of domain walls. Meanwhile, for $\lambda>1$ the local maxima become local minima of the potential, corresponding to false vacuum states. In this regime, we obtain a first-order phase transition and the false vacuum decays via the nucleation of bubbles. Of course, we do not necessarily expect a sharp distinction between the second-order and first-order phase transitions to occur at $\lambda\approx 1$. In particular, when $0<\lambda-1\ll 1$, 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 well-defined bubble nucleations, many regions of the true vacuum appear simultaneously. Such dynamics is intermediate between the tachyonic growth of long-wavelength modes and the regime of rare well-defined bubble nucleations.

Fig. 8 illustrates the nonlinear stabilization of the low-$k$ tachyonic modes by increasing the modulation amplitude $\delta$ of the interspecies conversion rate. The top two panels have $\lambda<1$. We observe the rapid emergence of a dynamically evolving domain wall network, as expected from the tachyonic growth of the long-wavelength modes. The transitionary behaviour between first- and second-order phase transitions is seen in the middle panels of fig. 8. Finally, well-defined bubble nucleations appear for $\lambda$ 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 $\lambda$. This effect appears in the figure as a delay in the time at which bubbles nucleate as the value of $\lambda$ is increased,666This 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 time-evolutions, 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 $\psi_{i}$ relative to the absolute mean. Since the RMS amplitude of fluctuations is set by $\hbar$ through the uncertainty principle, this is analogous to adjusting the value of $\hbar$. 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 short-wavelength Floquet modes have been excised from the evolution, but an exponential instability remains in the long-wavelength modes. In particular, in the regime of spinodal instability ($\lambda<1$), the exponentially growing long-wavelength 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. 9 shows the evolution of the RMS of both the total ($\varrho$) and relative ($\varepsilon$) energy densities for the $\lambda=0$ 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 $\lambda=0$ 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 Bose-Einstein condensates. In particular, we first showed that using the coupled Gross-Pitaevskii 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 Gross-Pitaevskii equations used to model the condensates break down at wavelengths above the instability scale. Some of these possible mechanisms were briefly outlined in section 4. However, this requires modeling of the condensates beyond the coupled GPEs, and such corrections must be analyzed on a case-by-case 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 long-wavelength 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 long-wavelength tachyonic modes and the decay of the metastable state through the formation of domain wall-antiwall pairs (i.e.,  one-dimensional 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 short-wavelength Floquet instability was sufficient to ensure the density perturbations remain small and the approximate relativistic scalar field description remains valid. The study of this cross-over 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 non-linear 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.

Acknowledgements.

JB and HVP were supported by the European Research Council (ERC) under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 3006478-CosmicDawn. 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 PHY-1607611. 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.

Appendix A Stochastic Semiclassical Simulations: The Truncated Wigner Approximation

In this appendix we briefly outline our methodology to solve for the dynamics of a two-component 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 semi-classical 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 high-dimensional) 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 case-by-case basis.

Initially, the condensates $\psi_{i}$ are assumed to be of the form

[TABLE]

where the $\bar{\psi}_{i}$ are background solutions, here taken to be stationary points of the coupled GPEs in the spatially homogeneous limit.777The total phase can experience a linear time-evolution at the stationary points, but we include the background chemical potential $\mu_{\rm bg}$ to remove this. Quantum corrections to the “classical” background $\bar{\psi}_{i}$ are modeled by a realization of a Gaussian random field

[TABLE]

where the $\hat{a}_{k}$ are draws of uncorrelated complex Gaussian random deviates with unit variance $\left\langle\left|\hat{a}_{k}\right|^{2}\right\rangle=1$, and $V$ is the simulation volume. We generate our samples $\hat{a}_{k}$ using the Box-Mueller transform

[TABLE]

where $\hat{P}_{k}$ and $\hat{A}_{k}$ are uniform random deviates on the interval $[0,1]$. Since the fluctuations are assumed to be Gaussian, the spectra $\mathcal{P}_{i}$ should match the initial quantum 2-point statistics888An obvious extension allows this to be extended to incorporate correlations between the condensates.

[TABLE]

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 $\varphi$. To make contact with the previous literature, here we follow the approach of Fialko et al. Fialko:2014xba ; Fialko:2016ggg and choose $\mathcal{P}(k)=\mathcal{A}\Theta(k_{\rm cut}-\left|\bm{k}\right|)$ for constant $\mathcal{A}$, resulting in a filtered white noise spectrum for $\delta\hat{\psi}_{k}$. 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 $\psi_{i}$ are evolved using the coupled GPEs. For our numerical scheme, we use a tenth-order accurate Gauss-Legendre 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 second-order accurate finite differencing approach for the spatial discretization, with

[TABLE]

and

[TABLE]

for an arbitrary scalar function $f$. In all cases, we enforce periodic boundary conditions, explicitly in the case of a finite-difference 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),

[TABLE]

through the introduction of three numerical scales $\omega_{\rm p}$, $\kappa_{\rm p}$ and $\Lambda$ with dimensions of $T^{-1}$, $L^{-1}$, and $L^{d}$ respectively. It is also convenient to introduce the numerical sound speed ($c_{\rm p}$) and energy ($E_{\rm p}$),

[TABLE]

Expressed in these variables, the dimensionless Hamiltonian is

[TABLE]

The (dimensionless) combination $\Lambda\kappa_{\rm p}^{d}$ 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 $\delta\bar{\psi}_{i}$. Meanwhile, the numerical energy $\hbar\omega_{\rm p}$ 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

[TABLE]

We have introduced the dimensionless numerical coefficients

[TABLE]

For convenience, we now specialize to the two condensate case, and define

[TABLE]

For our numerical simulations, we scale the condensates by their average background densities $\bar{n}$, set the numerical sound speed to be the propagation speed of the relative phase phonons, and the dimensionless nonlinear potential coupling to unity

[TABLE]

We immediately see

[TABLE]

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

[TABLE]

The commutator for the scaled field variables is

[TABLE]

which we can write on a discrete lattice with lattice sites $\bar{x}_{m}$

[TABLE]

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 $N^{-1/2}$,

[TABLE]

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 $N^{-1}$.

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 F and the corresponding dimensionless quantites by $\bar{\cdot}^{\rm F}$. As above, the dimenionless units in our paper are denoted by an overbar $\bar{\cdot}$. 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 $\omega_{0}=2\frac{\sqrt{\nu_{0}g_{\rm s}\bar{n}}}{\hbar}$ 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 $\left(\kappa_{\rm F}^{\rm F}\right)^{-d}$. The corresponding numerical scales are related as

[TABLE]

Finally, they measure $\nu_{0}$ in units of $g_{\rm s}\bar{n}$ instead of $\hbar\omega_{\rm P}^{\rm F}$, so that

[TABLE]

Appendix C Demonstration of Linearity of Floquet Modes and the Effective Discrete Lattice Wavenumber

In section 4 we demonstrated the presence of new short-wavelength 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 $\cos\varphi$. 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 short-wavelength dynamics are indeed a result of a linear instability. To do this, we rerun the simulations in fig. 3 using a second-order accurate centered finite-differencing 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 second-order accurate Laplacian stencil is $\frac{2}{\pi}$ 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 small-scale 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 well-known, but is presented here for completeness.

Consider a $d$-dimensional discrete lattice with sites $\bm{x}_{\vec{{\rm m}}}$ labelled by the vector index $\vec{{\rm m}}=(m_{1},m_{2},\dots,m_{d})$. We will denote function values at position $\bm{x}_{\vec{{\rm m}}}$ by $f_{\vec{{\rm m}}}\equiv f(\bm{x}_{\vec{{\rm m}}})$. For simplicity, assume that the lattice sites are arranged in a rectangular grid with uniform lattice spacing $dx$, so that $\bm{x}_{\vec{{\rm m}}}=\vec{{\rm m}}dx$. Consider finite-difference stencils for the Laplacian operator of the form

[TABLE]

with the vector indices $\bm{\alpha}\in\{(\alpha_{1},\alpha_{2},\dots,\alpha_{d}),\alpha_{i}=-\frac{N}{2},-\frac{N}{2}+1,\dots,\frac{N}{2}-1\}$. The choice of coefficients $c_{\bm{\alpha}}$ specify the stencil. Now consider a linear operator $\hat{\mathcal{O}}$ of the form

[TABLE]

with $\mathcal{C}$ a constant. For simplicity, we consider only a single function $f$ defined on our lattice. The generalization to many linearly coupled fields is straightforward. Taking a discrete Fourier transform of (40), we obtain in one-dimension

[TABLE]

with

[TABLE]

In the continuum, the RHS of (41) would be $(k^{2}+C)\tilde{f}_{k}$, so we identify an effective wavenumber $k_{\rm eff}(\bm{k})$ associated with the Laplacian stencil as a function of the continuum wavenumber $\bm{k}$

[TABLE]

where $k_{\rm nyq}=\frac{\pi}{dx}$ is the Nyquist wavenumber. To avoid spurious numerical damping associated with asymmetric stencils, we use a symmetric stencil $c_{\bm{\alpha}}=c_{-\bm{\alpha}}$ so that the imaginary contribution vanishes, and we have

[TABLE]

In fig. 12, we show this effective wavenumber for several low-order symmetric one-dimensional Laplacian stencils.

We also note one further complication to the picture above. Our code evolves the condensate fields $\psi_{i}$, and the numerical Laplacian is thus defined to act directly on these fields rather than the densities $\rho_{i}$ and phases $\phi_{i}$. 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 $\rho_{i}$ and $\phi_{i}$. 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 finite-differencing 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 $dt$ or grid spacing $dx$; and a demonstration that relevant conserved charges are time-invariant. 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 $\lambda=0$ and no Floquet modes), and the simulation shown in fig. 3 (with $\lambda=1.3$ and Floquet modes). For reference, the parameters are

[TABLE]

and

[TABLE]

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 time-dependence, 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 $dt$. In the linear regime, these are roughly the driving frequency $\omega$, the linear frequency of the Nyquist mode for the total ($\omega_{\rm tot}$) and relative ($\omega_{\rm rel}$) phonons, and the growth rate of the Floquet modes.999Technically, the oscillation frequencies should be obtained from the imaginary parts of the Floquet exponents, but $\omega_{\rm tot}$ and $\omega_{\rm rel}$ can be used as reasonable proxies. The free oscillation frequencies of the Nyquist modes can be easily obtained (see Braden et al. Braden:2017add )

[TABLE]

Since $\frac{\nu_{0}}{g_{\rm s}\bar{n}}\ll 1$, we have $\omega_{\rm rel}\approx\omega_{\rm tot}$. To obtain spatial resolution, we will be interested in the limit where the Floquet band is below the Nyquist mode, which occurs roughly when $\omega=2\omega_{\rm rel}$. Therefore, for spatially resolved simulations, we expect that $\omega_{\rm tot}$ and $\omega_{\rm rel}$ 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 $dt$ and grid spacing $dx$. These demonstrate the rapid convergence displayed by our Gauss-Legendre time stepping and pseudospectral discretization, respectively. To explore pointwise convergence, we consider the following norm between two simulations $\psi^{(n)}$ and $\psi^{(n+1)}$

[TABLE]

which measures the maximal pointwise difference between the two simulations over the entire simulation volume. Here, the superscript (n) labels the simulations, and $\mathbb{L}$ 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 $dt$ or grid spacing $dx$ differ by a factor of $2$, 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 $dt$ 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 $dt$ (fig. 13) or $dx$ (fig. 14) varied. As the time step is decreased, initially we see the pointwise error decrease by a factor of $2^{10}\sim 1000$, as expected for a tenth-order accurate integrator. This eventually saturates at a level of $\mathcal{O}(10^{-13})$ 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 time-dependent saturation of the temporal accuracy, we also performed the same temporal convergence tests for ${\frac{\sqrt{g_{\rm s}\bar{n}m}}{\hbar}dx=1.89}$ and ${\frac{\sqrt{g_{\rm s}\bar{n}m}}{\hbar}dx=1.79}$ for the ${\lambda=1.3}$ simulation, corresponding to a grid-spacing 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 ${\lambda=0}$ simulation for the former case. As a test that these errors did not arise from aliasing into long-wavelength fluctuations, we performed the corresponding convergence tests with a finite-differencing Laplacian stencil (which does not lead to aliasing) and found the same growth in the saturation level.

Meanwhile, as the spatial grid spacing is decreased we see an exponential improvement in the pointwise convergence owing to our use of a pseudospectral scheme. In particular, for the simulations shown here, an increase in the number of grid points by a factor of $8$ leads to over a ten order of magnitude improvement in convergence without the chaotic behaviour induced by the Floquet band. As with the temporal convergence, we see the saturation level increase with time for the simulation including Floquet modes.

D.2 Conserved charge preservation

The previous subsection explicitly demonstrated the rapid convergence properties of our temporal and spatial discretization schemes. We also observed a saturation of the pointwise accuracy of our simulations, which we believe arises from the exponential linear growth of fluctuation amplitudes and the resulting chaotic dynamics in the full field space. Since we lack an exact analytic solution to compare to, technically these results demonstrate that we are converging to a solution, but not necessarily the correct solution of the equations of motion. To demonstrate that we are indeed solving the coupled GPEs, we now show that our numerical scheme preserves the various conserved quantities of the continuum equations. As emphasized above, the goal of this work is to explicitly demonstrate the need for modeling of the short wavelength behavior of the condensate, not to explore how these corrections modify the false vacuum decay picture. Therefore, we work under the assumption that the coupled GPEs are indeed the correct evolution equations. Of course, the various phenomenological corrections briefly mentioned in section 4 will lead to a violation of these conserved charges, but this will be the subject of future work.

When $\delta=0$, we have the following three locally conserved charges

the total condensate particle density $\varrho\equiv\sum_{i}\left|\psi_{i}\right|^{2}$, 2. 2.

the total momentum density $\mathcal{P}\equiv\frac{i}{2}\sum_{i}\left(\nabla\psi_{i}^{\dagger}\psi_{i}-\psi_{i}^{\dagger}\nabla\psi_{i}\right)=\sum_{i}\rho_{i}\nabla\phi_{i}$, and 3. 3.

the condensate energy density $\mathcal{H}\equiv\sum_{i}\left(\frac{\hbar^{2}}{2m}\left|\nabla\psi_{i}\right|^{2}+\frac{g}{2}\left|\psi_{i}\right|^{4}\right)-\nu\left(\psi_{1}^{\dagger}\psi_{2}+\psi_{1}\psi_{2}^{\dagger}\right)$.

These are readily identified, respectively, with invariance of the action under a global phase rotation of the $\psi_{i}$’s, spatial translations, and time translation. When $\delta\neq 0$ we lose time translation invariance, and the energy density is no longer conserved. Although not considered here, if we further consider inhomogeneous condensates evolving in an external potential, then the spatial translation invariance will be broken and momentum no longer conserved. However, even in this limit the total particle number will be conserved, which is the origin of the ungapped linear fluctuations associated with the total phase phonons.101010In the limit $\nu=0$, our coupled GPEs decouple into two independent condensates, each satisfying the nonlinear Schrodinger equation (NLS). Therefore, in this limit the individual condensate densities, momenta, and energies are conserved. Further, in $1$+$1$-dimensions the NLS is integrable, with an infinite hierarchy of additional conserved quantities. We will not consider this limit here, although how the introduction of $\nu$ leads to the breaking of these conserved quantities is interesting. In fig. 15, we demonstrate the preservation of these conserved quantities for our fiducial convergence testing simulations. To make a fair comparison, we choose ${\frac{\sqrt{g_{\rm s}\bar{n}m}}{\hbar}dx=\frac{50\times 2}{4720}\frac{\nu_{0}}{g_{\rm s}\bar{n}}\approx 0.19}$ and ${\frac{g_{\rm s}\bar{n}}{\hbar}dt=\frac{2\pi}{256}\frac{g_{\rm s}\bar{n}}{\hbar\omega}\approx 5.49\times 10^{-3}}$ for both simulations. In all cases, we obtain excellent preservation of all of the Noether charges for the simulations both with and without the Floquet instability. The results are particularly striking for the mean total density, where discrete jumps are visible associated with the finite precision of machine arithmetic. Similar discrete jumps also occur for the energy density for the case where the Hamiltonian is time-independent. When $\nu$ oscillates in time and the Floquet modes are present in the system, the external driver $\nu$ drives exponential growth of the unstable modes. As a result, energy is injected into the system, leading to a growth in the mean energy density with time. We verified that this is the case, but do not include the results in this appendix as here we want to focus on only the conserved charges.

For simplicity, here we have shown the conservation properties for a single choice of spatial and temporal resolution. However, the accuracy of the conservation laws is sensitive to the choice of $dx$ and $dt$. Although not presented here, we also investigated how well the conserved charges were preserved for the same suite of simulations used in the direct convergence testing shown in fig. 13 and fig. 14. From these investigations, we found that the preservation of the total density $\varrho$ was primarily sensitive to the choice of time-step $dt$, with little sensitivity to the choice of grid-spacing $dx$. Meanwhile, the preservation of the total momentum $\mathcal{P}$ had a saturation level with decreasing $dt$ that was sensitive to the choice of lattice spacing. Finally, the energy density $\mathcal{H}$ displays convergence properties intermediate between the total density and momentum, with more sensitivity to the grid-spacing $dx$ than the number density, but spatial saturation occuring at a larger grid spacing than for the momentum.