Inflationary preheating dynamics with ultracold atoms

TL;DR

This paper explores how dynamical instabilities in ultracold atom systems can amplify quantum loop corrections, drawing parallels to early-universe particle production, and proposes experimental protocols for observing these effects.

Contribution

It provides a detailed analytical and numerical study of tachyonic instabilities in ultracold atoms, linking macroscopic growth to microscopic quantum processes, and suggests optimized experimental methods.

Findings

Amplification of loop corrections via dynamical instabilities in ultracold gases.

Establishment of a correspondence between growth rates and Feynman diagrams.

Proposed protocols for experimental realization of instability-driven amplification.

Abstract

We discuss the amplification of loop corrections in quantum many-body systems through dynamical instabilities. As an example, we investigate both analytically and numerically a two-component ultracold atom system in one spatial dimension. The model features a tachyonic instability, which incorporates characteristic aspects of the mechanisms for particle production in early-universe inflaton models. We establish a direct correspondence between measureable macroscopic growth rates for occupation numbers of the ultracold Bose gas and the underlying microscopic processes in terms of Feynman loop diagrams. We analyze several existing ultracold atom setups featuring dynamical instabilities and propose optimized protocols for their experimental realization. We demonstrate that relevant dynamical processes can be enhanced using a seeding procedure for unstable modes and clarify the role of initial quantum fluctuations and the generation of a non-linear secondary stage for the amplification of modes.

Equations

$$H$$

$$+\frac{g_{12}}{2}\int dz\left\{({\phi }_1^{†}{)}^2{\phi }_2^2+4|{\phi }_1{|}^2|{\phi }_2{|}^2+({\phi }_2^{†}{)}^2{\phi }_1^2\right\},$$

$$H_{3D}=\int d^{3}r\left[\frac{{\left|\nabla \psi \left(\mathbf{r}\right)\right|}^2}{2M}+V(\mathbf{r}){\left|\psi \left(\mathbf{r}\right)\right|}^2+\frac{g}{2}{\left|\psi \left(\mathbf{r}\right)\right|}^4\right],$$

$$\psi (\mathbf{r})={\phi }_1(z){\psi }_1(x,y)+{\phi }_2(z){\psi }_2(x,y)+\dots ,$$

$$h_{\text{SC}}\sim \left[({\phi }_1^{†}{)}^2{\phi }_2^2+\text{h.c.}\right].$$

$${\tilde{h}}_{\text{SC}}$$

$$=\frac{1}{2}\left[{\left({\tilde{\phi }}_S^{†}\right)}^2{\tilde{\phi }}_0^2+{\left({\tilde{\phi }}_A^{†}\right)}^2{\tilde{\phi }}_0^2+\text{h.c.}\right],$$

$${\mu }_2{\phi }_2=\left[-\frac{{\nabla }_z^2}{2M}+V(z)+g_2|{\phi }_2{|}^2\right]{\phi }_2.$$

$$H$$

$$+\frac{g_{12}}{2}\int dz\left[({\phi }_1^{†}{)}^2{\phi }_{2,s}^2+({\phi }_{2,s}^{†}{)}^2{\phi }_1^2\right].$$

$$i{\partial }_t{\phi }_1$$

$$h_0[{\phi }_{2,s}]$$

$$h_1[{\phi }_{2,s}]$$

$${\phi }_1(t,z)$$

$$\left(\begin{array}{cc}{h}_0& h_1\\ -h_1& -h_0\end{array}\right)\left(\begin{array}{c}{u}_k(z)\\ v_k(z)\end{array}\right)={\omega }_k\left(\begin{array}{c}{u}_k(z)\\ v_k(z)\end{array}\right).$$

$$\frac{1}{2R_{TF}}\int dz{e}^{iqz}{\phi }_{2,s}^{(b)}(z)\to \delta (q)$$

$${\omega }_k^2={\left(\frac{k^2}{2M}+\frac{2g_{12}{\mu }_2}{g_2}-{\mu }_1\right)}^2-{\left(\frac{g_{12}{\mu }_2}{g_2}\right)}^2.$$

$$\pm k_{\ast }=\pm \sqrt{2M\left({\mu }_1-{\mu }_2\frac{2g_{12}}{g_2}\right)}$$

$$\delta k=\sqrt{2M}\left(\sqrt{{\mu }_1-{\mu }_2\frac{g_{12}}{g_2}}-\sqrt{{\mu }_1-{\mu }_2\frac{3g_{12}}{g_2}}\right)$$

$$\frac{1}{2}⟨\left\{{\phi }_1^{†}(t,k_{\ast }),{\phi }_1(t,k_{\ast })\right\}⟩\simeq e^{2{\gamma }_{\ast }t},$$

$${\bar{\gamma }}_1\lesssim {\gamma }_{\ast }=\frac{g_{12}}{g_2}{\mu }_2.$$

$$i{\partial }_t{\phi }_1$$

$$+g_{12}\left({\phi }_1^{\ast }{\phi }_2^2+2{\phi }_1|{\phi }_2{|}^2\right),$$

$$i{\partial }_t{\phi }_2$$

$$+g_{12}\left({\phi }_2^{\ast }{\phi }_1^2+2{\phi }_2|{\phi }_1{|}^2\right).$$

$$W({\phi }_1,{\phi }_2)=\frac{4}{{\pi }^2}\exp \left(-2|{\phi }_1{|}^2-2|{\phi }_2-{\phi }_{2,s}{|}^2\right).$$

$$N=\sum _k\left[n_1(k,t)+n_2(k,t)\right].$$

$$F_{\alpha \beta }(t,z;t^{\prime },z^{\prime })$$

$$-⟨{\phi }_{\alpha }(t,z)⟩⟨{\phi }_{\beta }(t^{\prime },z^{\prime })⟩,$$

$$n_1(k,t)+2\pi \delta (0)|{\bar{\varphi }}_1^2|=\int dz{e}^{-ikz}{F}_{12}(t,z;t,0),$$

$$n_2(k,t)+2\pi \delta (0)|{\bar{\varphi }}_2^2|=\int dz{e}^{-ikz}{F}_{34}(t,z;t,0),$$

$${\rho }_{\alpha \beta }(t,z;t^{\prime },z^{\prime })$$

$${\mathcal{D}}_{\alpha \gamma }(t,z)F_{\gamma \beta }(t,z;t^{\prime },z^{\prime })$$

$${\mathcal{D}}_{\alpha \gamma }(t,z){\rho }_{\gamma \beta }(t,z;t^{\prime },z^{\prime })$$

$$F_{\alpha \beta }(t,t^{\prime },k)\sim e^{\gamma (k)(t+t^{\prime })}\equiv f_{\ast }(t,t^{\prime },k),(\alpha ,\beta =1,2)$$

Full text

I Introduction and overview

Ultracold quantum gases provide a unique opportunity for studying nonequilibrium phenomena that are otherwise very difficult to access experimentally. Ensembles of trapped atoms can be largely isolated from the environment, such that they offer the possibility to address fundamental aspects of quantum many-body systems governed by unitary time evolution. In addition, cold-atom setups provide a very flexible testbed, with tunable interactions or with reduced dimensionality by shaping the confining optical potential, realizing a multitude of different relevant physical situations.

It has been pointed out that characteristic aspects of the evolution of the early universe may be addressed with table-top experiments Schmiedmayer and Berges (2013), such as acoustic oscillations in the time evolution of an ultracold quantum degenerate gas Hung et al. (2013), or inflationary quasiparticle creation and thermalization dynamics in coupled Bose-Einstein condensates Opanchuk et al. (2013); Posazhennikova et al. (2016). The latter may involve paradigmatic dynamical phenomena such as nonequilibrium instabilities seeded by quantum fluctuations Felder et al. (2001a), which trigger important nonlinear phenomena such as enhanced particle creation from secondary instabilities Berges and Serreau (2003) and the subsequent approach to nonthermal fixed points with Bose condensation Berges et al. (2008) long before thermalization sets in.

In this work we analyze how two-component bosonic cold-atom systems in one spatial dimension may be used to address relevant aspects of the nonequilibrium dynamics of tachyonic instabilities arising in the context of inflationary models Felder et al. (2001a, b). Experimentally, we have in mind the twin-atom beam experiments of the Schmiedmayer laboratories Bücker et al. (2011), and the spin-changing collision dynamics arising in bosonic spin-one systems as employed by the Oberthaler group Linnemann et al. (2016). We identify the relevant dynamical processes of the cold-atom system and describe them using quantum-statistical field theory. In this way we establish a one-to-one correspondence between macroscopic observables of the Bose gas and microscopic processes in terms of Feynman loop diagrams for correlation functions.

More precisely, we first analyze the generation of instabilities that are seeded from initial quantum fluctuations in the Bose gas. These primary instabilities trigger an exponential increase of the occupation numbers of characteristic unstable modes with momenta $\pm k_{*}$. This early stage is followed in time by a secondary growth period, which describes the nonlinear enhancement of occupancies in a wider momentum range. In FIG. 1 we sketch the initial state and the subsequent processes for the Bose gas with components $\varphi_{1}$ and $\varphi_{2}$. The initially unstable state $\varphi_{2}$ decays via the excitation of characteristic modes of $\varphi_{1}$ at momenta $\pm k_{*}$. As these modes become highly occupied, their interactions lead to the excitation of modes with higher momenta of multiples of $k_{*}$. Remarkably, macroscopically measurable growth rates and their characteristic momenta give, in principle, direct information about the topology of the Feynman diagrams that describe the underlying microscopic processes.

While for typical inflaton models the secondaries have a significant impact on the dynamics, their role can be restricted by total particle number limitations in cold-atom experiments. In particular, we demonstrate that for the existing twin-atom beam setup Bücker et al. (2011) the experimental detection of the relevant higher-loop processes may be difficult. To overcome this limitation, we propose to seed the primary instabilities: Instead of preparing the whole condensate in the ‘excited’ state $\varphi_{2}$, one keeps a fraction of all atoms in the ‘lower’ state $\varphi_{1}$ and imprints a modulation to transfer them into the most unstable modes $\pm k_{*}$. We perform classical-statistical (Truncated Wigner) simulations to numerically support our findings.

This publication is organized as follows: In section II we discuss our model system and relate it to possible experimental realizations. Subsequently, we analytically determine the early-time dynamics of this many-body system including the primary instability in section III. In the following section IV we demonstrate the existence of secondaries being present in this system numerically and give a simple heuristic argument for their generation. In section V we complement the observations with a more refined explanation in terms of Feynman diagrams, which allows us to estimate the secondaries’ growth rates. Furthermore, we propose a seeding strategy to make these non-linear phenomena experimentally better accessible in section VI. To get a better understanding of the situation with and without seed, we consider a classical toy model that encompasses the necessary ingredients for the generation of secondaries in section VII. Finally, we conclude in section VIII.

II Two-component Bose gas

The inflationary early universe undergoes very rapid expansion, after which all energy ends up in a vacuum-like state without entropy or particles. Typical models describe this state in terms of a coherent Bose field, the inflaton, whose decay finally leads to the observed particle content of today’s universe. Very efficient mechanisms for particle production from inflaton decay involve nonequilibrium instabilities, such as parametric resonance Kofman et al. (1994, 1997) or tachyonic instabilities Felder et al. (2001a, b). The latter can be particularly fast and nonlinear many-body phenomena beyond simple classical-field or Bogoliubov–type approximations play a crucial role.

A prominent class of inflationary scenarios – hybrid inflation – involve multiple scalar fields Linde (2008). In this work, we consider a model with a bosonic complex field operator $\varphi_{a}$ with two components $a=1,2$. This model allows us to address essential features of instability dynamics and the build-up of non-linearities, where we will concentrate on tachyonic instabilities associated to imaginary parts in the dispersion relation. While the couplings to other degrees of freedom of the Standard Model of particle physics play an important role at later times, we focus on the bosonic self-interactions that dominate the dynamics at earlier times of the nonequilibrium evolution. It should be emphasized that, even though typical couplings in cosmological inflaton models are usually very small, the dynamics become strongly correlated because of the instabilities that amplify small nonlinear corrections exponentially fast.

Although the time scales depend on the specifics of any microscopic model chosen, essential aspects of the dynamics can be be studied in a wide range of different many-body systems. This is based on the observation that there is an effective loss of details about the initial conditions and model parameters as time proceeds. In recent years this culminated in the discovery of new universality classes far from equilibrium Berges et al. (2008), which characterize a significant part of the systems’ dynamical evolution in terms of scaling exponents and scaling functions that are the same for a wide range of different relativistic models Berges et al. (2008); Berges and Hoffmeister (2009); Berges and Sexty (2011); Gasenzer et al. (2012); Berges and Sexty (2012); Gasenzer et al. (2014); Orioli et al. (2015) as well as non-relativistic systems of ultracold quantum gases Scheppach et al. (2010); Nowak et al. (2011, 2012); Berges and Sexty (2012); Orioli et al. (2015); Schachner et al. (2016); Karl and Gasenzer (2016).

In an ultracold atom setup with different field components $\varphi_{a}$, the index $a$ typically labels internal states of the bosonic atoms such as hyperfine states or eigenstates within an external trapping potential. The quantum many-body Hamiltonian of our model system is given by

[TABLE]

where we abbreviated $|\varphi_{a}|^{2}=\varphi^{\dagger}_{a}\varphi_{a}$ and $|\varphi_{a}|^{4}=\varphi^{\dagger}_{a}\varphi^{\dagger}_{a}\varphi_{a}\varphi_{a}$. The interaction constants $g_{a}$ and $g_{12}$ are related to intra- and interspecies scatterings, respectively, and we only consider the repulsive case $g_{a},g_{12}>0$. $M,V(z)$ and $\mu_{a}$ are the atomic mass, the one-dimensional trapping potential and two chemical potentials, respectively. The latter controls the total particle number $N$ and the energy gap $\Delta E=\mu_{1}-\mu_{2}$ between the two components.

The system with Hamiltonian (1) is invariant under a global $\mathcal{O}(2)$ transformation, $(\varphi_{1},\varphi_{2})\rightarrow(\varphi_{1},\varphi_{2})e^{i\theta}$, which implies total particle number conservation. For $g_{12}=0$, the symmetry group is enlarged to $\mathcal{O}(2)\times\mathcal{O}(2)$, i.e. $(\varphi_{1},\varphi_{2})\rightarrow(\varphi_{1}e^{i\theta_{1}},\varphi_{2}e^{i\theta_{2}})$ corresponding to individual number conservation. For the special case of $g_{1}=g_{2}=g_{12}$, the interaction part is invariant under $\mathcal{O}(2)\times\mathbb{Z}_{2}$ due to the symmetry $(\varphi_{1},\varphi_{2})\rightarrow(\varphi_{2},\varphi_{1})$.

The above Hamiltonian can be realized in various ways. We specifically have in mind two different experimental setups that rely on external or internal degrees of freedom and should allow the observation of the dynamical features to be discussed in the following sections. The first experimental setup concerns the twin-atom beam experiment Bücker et al. (2011). The authors realized an effectively one-dimensional Bose gas by a magnetic atom-chip trap and applied sophisticated optimal-control techniques to access specific excited states of the radial trapping potential, which in turn serve as the required two components. More precisely, the full three-dimensional system is described by the Hamiltonian

[TABLE]

where $V(\mathbf{r})$ is the trapping potential in three spatial dimensions (3D) and $g$ denotes the interaction constant. The fields $\varphi_{a}(z)$ depending only on the (1D) spatial $z$-direction arise as coefficients in an expansion of the 3D field operator, taken as

[TABLE]

where $\psi_{a}(x,y)$ denote the two lowest lying eigenfunctions of $V(\mathbf{r})$. Performing the $x,y$ integrals in the grand-canonical Hamiltonian corresponding to $H_{3D}$ and truncating the above expansion yields a Hamiltonian of the form $\eqref{eq:Hamiltonian1}$. This truncation is expected to be a valid approximation if the radial potential is chosen such that the energetically higher lying states are not excited. The precise form of the radial potential then determines the parameters entering the two-component model in consideration (cf. appendix B).

As an alternative to exciting states of the external potential, one may utilize internal degrees of freedom of spinor Bose gases Kawaguchi and Ueda (2012). In such setups, the gas is usually rendered effectively one-dimensional by optical dipole traps. We emphasize that our model system is closely related to spin-1 systems but can not be mapped onto so-called pseudo-spin-1/2 systems. The latter have been extensively studied, e.g. as binary mixtures in the context of the miscible-immiscible phase transition Sabbatini et al. (2011). It is necessary to consider larger spins because the Hamiltonian (1) differs from a pseudo-spin-1/2 system by the “spin-changing” terms

[TABLE]

These terms turn out to play a crucial role for the dynamical instabilities considered in this work. Experimentally, spin-changing collisions arise naturally in bosonic spin-1 systems and have been utilized, e.g., for the creation of entangled twin-atom states Gross et al. (2011). In such a setup, the bosonic gas is confined to a hyperfine state $F=1$ with three magnetic quantum numbers $m_{F}=0,\pm 1$. The spin changing terms are then given by

[TABLE]

where $\tilde{\varphi}_{m_{F}}$ label the spinor components and we have introduced $\tilde{\varphi}_{S}=\frac{1}{\sqrt{2}}\left(\tilde{\varphi}_{1}+\tilde{\varphi}_{-1}\right)$ and $\tilde{\varphi}_{A}=\frac{i}{\sqrt{2}}\left(\tilde{\varphi}_{1}-\tilde{\varphi}_{-1}\right)$. Thus $m_{F}=0$ corresponds to $\varphi_{2}$, and $\varphi_{1}$ would be doubly degenerate as the (anti-) symmetric combinations of $m_{F}=\pm 1$. In the presence of the quadratic and negligibly small linear Zeeman effect, the resulting Hamiltonian mainly differs from (1) by this degeneracy and we restrict ourselves to the model containing only two components for simplicity. Nevertheless, we will argue that our results only depend on a few generic ingredients, such that many aspects also apply to the system with three components. For completeness, we note that the full spin-1 Hamiltonian can be found, e.g., in Ref. Uchino et al. (2010).

III Primary Instabilities

Before we turn to simulations in section IV, we first obtain an analytic understanding of primary instabilities. This is done by linearizing the equations of motion in fluctuations around the initial condensate. We identify instabilities as exponentially growing solutions characterized by a dispersion relation with non-vanishing imaginary part.

In recent twin beam experiments Bücker et al. (2011); Wasak et al. (2014) it was possible to transform all atoms from the lower state $\varphi_{1}$ to the energetically higher state $\varphi_{2}$. To lowest order in fluctuations, this configuration is considered to be a stationary state which we denote by $\varphi_{2,s}$. The latter is the solution of the equation

[TABLE]

In order to detect the dynamical instability, we describe the evolution of $\varphi_{1}$ by using the approximate Hamiltonian

[TABLE]

The Heisenberg equations of motion for the atoms then read

[TABLE]

where we have abbreviated

[TABLE]

The approximate dynamical equations for $\varphi_{1}$ can be solved by a Bogoliubov transformation

[TABLE]

which leads to the Bogoliubov-de-Gennes equations

[TABLE]

We refer to the appendix A for details of the Bogoliubov transformation in the presence of unstable modes. The eigenvalues $\omega_{k}$ of (12) and the Bogoliubov mode functions $u_{k}$ and $v_{k}$ determine the low-energetic excitations and thus the quasi-particle spectrum. We proceed by specifying the trapping potential as either $a)$ harmonic $V(z)=\frac{m}{2}\omega^{2}z^{2}$, or $b)$ a box potential $V(z)=V_{\infty}\Theta\left(|z|-\frac{L}{2}\right)$ with $V_{\infty}\rightarrow\infty$. In the former case, we use the stationary solution $\varphi_{2,s}^{(a)}=\sqrt{\frac{\mu_{2}}{g_{2}}\left(1-\frac{z^{2}}{R_{\text{TF}^{2}}}\right)}$ within the Thomas-Fermi (TF) approximation with the TF radius given by $R_{\text{TF}}=\sqrt{\frac{2\mu_{2}}{Mw^{2}}}$. The normalisation of $\varphi_{2}$ relates the particle number $N_{2}$ and the chemical potential. For case $b)$, we obtain $\varphi_{2,s}^{(b)}=\sqrt{\frac{\mu_{2}}{g_{2}}}$ with $\mu_{2}=n_{2}g_{2}$ where $n_{2}=\frac{N_{2}}{L}$ is the homogeneous particle density. Going from case $a)$ to case $b)$ corresponds to the replacement

[TABLE]

in the equations of motion in momentum space. Solving these equations, we arrive at the squared dispersion relation

[TABLE]

The system becomes unstable for $\omega_{k}^{2}<0$ with the two most unstable modes being

[TABLE]

as determined from the maximal $-\omega_{k_{*}}^{2}\equiv\gamma_{*}^{2}=\left(\frac{g_{12}}{g_{2}}\mu_{2}\right)^{2}$. Because of the symmetry $k\rightarrow-k$ in the Bogoliubov spectrum, one produces two atom beams in opposite directions with momentum $k_{*}$. Henceforth, these beams are termed twin beams. The bandwidth of the instability is given by

[TABLE]

Because of the property $\omega_{k}^{2}<0$ in the instability regime, they are called tachyonic instabilities in a cosmological context Felder et al. (2001a, b). The initial growth of the occupation of atoms associated with the most unstable mode is

[TABLE]

where $\varphi(t,k)$ is the momentum mode appearing in the Fourier expansion $\varphi(t,x)=\frac{1}{\sqrt{L}}\sum_{k}e^{-ikx}\varphi(t,k)$. Note that the above growth rate has been calculated for the homogeneous case $b)$ and serves as a reference value for the trapped case $a)$, where nearby modes will also contribute to the instability (cf. eq. (13)). Since these generally grow slower, we expect that the true most unstable growth rate $2\bar{\gamma}_{1}$ will effectively be damped, i.e.

[TABLE]

IV Secondaries

After having determined the dynamics at early times, which are governed by the dispersion relation (14), we proceed with a numerical simulation. This allows us to access subsequent nonlinear dynamics including secondary growth rates and to validate our analytical calculations for the primary instabilities.

We employ the classical statistical or truncated Wigner approach (TWA) to get more insight into the fluctuations of the atoms in the excited as well as the ground state. To this end, we consider the classical equations of motions

[TABLE]

We choose a coherent field in the excited state corresponding to the solution $\varphi_{2,s}$ and initialize the lower state with quantum noise. Our main observable will be the mode occupation number defined as the ensemble-averaged absolute value of the Fourier transformed fields, i.e. n_{a}(k,t)=\braket{\left missing}{\varphi_{a}(k,t)}{{}^{2}}_{W}. The subscript denotes averaging of the classically evolved fields with respect to different initial conditions, which are sampled according to the initial Wigner distribution

[TABLE]

The mode occupation number is connected to the conserved particle number

[TABLE]

The classical-statistical approximation is essentially valid for high typical occupation numbers. For a detailed investigation in the case of one-dimensional bosonic gases, see Ref. Berges and Gasenzer (2007).

Figure 2 presents results of a simulation with a harmonic trapping potential (case $a)$). The detailed parameter sets for each simulation can be found in appendix B. Shown is the occupation number of $\varphi_{1}$ as a function of momentum in units of the estimated most unstable mode $k_{*}$. The corresponding dispersion is plotted in the inset.

The figure gives results at different times. At initialisation, there are only quantum fluctuations present. However, the system features a dynamical instability which manifests itself in the growth of modes close to $k_{*}$. At early times the momentum range of the instability is in excellent agreement with our analytical estimate as can been seen from the squared dispersion relation (14) displayed in the inset. At some intermediate time, additional modes around $3k_{*}$ become occupied. These momentum modes are initially stable, and they start growing during a secondary stage of amplification after the primary instability led to sizeable nonlinearities.

The process of the generation of these secondaries is depicted in figure 1 and can be understood as follows: The system is prepared in an unstable state, where the whole condensate is transferred to the excited state $\varphi_{2}$. Consequently, the condensate decays into the energetically favourable state $\varphi_{1}$. By symmetry this decay occurs as the emission of atom pairs in opposite directions, i.e. twin beams with momentum $\pm k_{*}$. Since the occupation number of the unstable modes grows exponentially, they become highly occupied and the system enters a regime where their self-interaction is no longer negligible. At this point the Bogoliubov approximation breaks down and the stage of secondary amplification sets in. Processes involving quartic self-interactions between modes then give rise to a particularly fast growth of the mode with momentum $3k_{*}$, since the process $\left(k_{*},k_{*}\right)\rightarrow\left(-k_{*},3k_{*}\right)$ involves the maximum number of already highly occupied modes allowed by momentum conservation.

Furthermore, the Hamiltonian (1) contains another possibility: Two twin beam atoms with momentum $k_{*}$ can interact with the condensate mode in $\varphi_{2}$ which is expected to increase the occupation number of the $2k_{*}$ mode in $\varphi_{2}$. This is indeed the case as can be seen in figure 3, where we have plotted the relevant mode occupation numbers as a function of time.

We emphasize that the occupation numbers are obtained from an equal-time correlation function of two fields. In frequency space such a quantity involves an integral over all possible frequencies and the non-equilibrium processes are not restricted by “on-shell” energy conservation. Though we are not dealing with a homogeneous system, and momentum conservation can hold at most approximately, corresponding observations are also made for the spatially translation invariant case relevant for cosmology Berges and Serreau (2003).

The simulations reveal that the condensate stays approximately constant for early times, which justifies the Bogoliubov approximation of section III to understand the initial instability. As expected, the primary unstable mode grows exponentially and we give a more refined explanation for the secondary growth rates in terms of Feynman diagrams in the next section V. For now, we note that the growth rates are higher than those from the primary instability – the secondaries start later but speed up a lot. As soon as the growth of all modes stops, one enters a nonperturbative regime where the dynamics is no longer dominated by single scattering processes. In figure 2, we can see that the clear peak structure of the spectrum is washed out in this regime. One also observes from the growth of the zero-momentum mode $n_{1}(0)$ an indication for the onset of the out-of equilibrium formation of a quasi-condensate in $\varphi_{1}$ (see figure 3).

V Diagrammatic Interpretation

In order to make analytical progress in the non-linear regime, we consider the time evolution of the so-called statistical two-point correlation function

[TABLE]

where $\alpha,\beta=1,2,3,4$ label $\varphi_{1},\varphi_{1}^{\dagger},\varphi_{2},\varphi_{2}^{\dagger}$, respectively. Here $\{A,B\}\equiv AB+BA$ denotes the anti-commutator, which is the appropriate object to consider in the quantum theory because expectation values calculated from the classical-statistical simulations always concern symmetrized field operator products Berges and Gasenzer (2007). Furthermore, the anti-commutator is experimentally accessible since it is related to our main observables $n_{a}(k,t)$ via

[TABLE]

where $\bar{\phi}_{a}$ denotes the condensate part and $n_{a}$ the respective non-condensate fractions. Introducing also the commutator or so-called spectral fucntion as

[TABLE]

one can derive a set of quantum evolution equations for the nonequilibrium time-evolution of two-point correlation functions (see e.g. Ref. Berges (2015)). For given symmetrized part of the self-energy, $\Sigma^{F}(F,\rho)$, and anti-symmetrized part, $\Sigma^{\rho}(F,\rho)$, they read

[TABLE]

where the differential operator $\mathcal{D}_{\alpha\gamma}$ is given in appendix C along with a more detailed discussion of the equations. The self-energies $\Sigma^{F}$ and $\Sigma^{\rho}$ contain all interaction contributions such as direct scattering processes or decays. For known self-energies the equations (26) and (27) are exact and equivalent to Kadanoff-Baym or Schwinger-Dyson equations Berges (2015).

The primary instability results in an exponential growth of the statistical two-point function in Fourier space as $F_{\alpha\beta}(t,t^{\prime},k_{*})\sim e^{\gamma_{*}(t+t^{\prime})}\quad(\alpha,\beta=1,2)$. Thus we can estimate nonlinear effects from a resummed loop expansion of the self-energies appearing in the quantum evolution equations (26) and (27). For the non-relativistic Bose gases in one dimension, this is an expansion in the dimensionless diluteness parameter $\zeta=\sqrt{na_{s}}$ with the scattering length $a_{s}$ and the average density $n$. In the present case, we have $\zeta\sim{\mathcal{O}}(${10}^{-5}\text{,}$)$ depending on the precise parameters of the scenarios discussed in this work and for sufficiently short times, before the typical occupancies become of order $\sim 1/\zeta$, it represents an accurate approach for the description of the dynamics.

To explain the phenomenon of secondary instabilities, it is sufficient to consider the series up to two loops. In the following, we show how the secondaries can be given an interpretation in terms of Feynman diagrams. The analysis follows along the lines of reference Berges and Serreau (2003) for the inflaton decay in a cosmological context, which is also reviewed in Berges (2015). We emphasise that the analytical results we present will only provide a “leading-log” estimate for the exponential growth rates. A quantitative analysis along these lines including also field amplitudes would require a numerical treatment.

To determine the secondary growth rates it is sufficient to consider the primary growth of the $1,2$ components of the statistical two-point function in Fourier space as

[TABLE]

with $\gamma(k)=\sqrt{\gamma_{*}^{2}-\left(\frac{k^{2}-k_{*}^{2}}{2M}\right)^{2}}$, see (86). The leading one-loop correction to the evolution equation of $F_{\alpha\beta}$ with $\alpha,\beta=3,4$, i.e. the excited field, comes from the self-energy $\Sigma^{F}$. Thus we neglect all terms of quadratic order in $\rho$, which is justified for $F^{2}\gg\rho^{2}$ for typical modes. In momentum space, the corresponding contribution to the RHS of (26) is given by the memory integral

[TABLE]

Due to the exponential growth, the latest times dominate the time integral. Consequently, we approximate $f_{*}(t,s,p)\approx f_{*}(t,t^{\prime},p)$ and

[TABLE]

where we have used the equal time commutation relations. Furthermore, we simplify the memory integral to include only recent times with $\int_{t_{0}}^{t^{\prime}}\mapsto\int_{t^{\prime}-c/\mu_{2}}^{t^{\prime}}$. For an appropriate choice of $c$, the resulting error for the growth rates will be $\sim\log c$. With this approximation we obtain

[TABLE]

Turning to the momentum integral, we note that it is dominated by $\gamma(k_{*})=\gamma_{*}$ and consider a saddle-point approximation using

[TABLE]

where we assume $k_{*}>0$ and $\frac{k_{*}^{2}}{2M}>\gamma_{*}$ for simplicity. At $k=2k_{*}$ we obtain

[TABLE]

which is valid for $\gamma_{*}(t+t^{\prime})\gg\left(\frac{M\gamma_{*}}{k_{*}^{2}}\right)^{2}$. Reinstating all constants, we find the approximation

[TABLE]

which yields the growth rate estimate

[TABLE]

In a similar fashion, we can estimate the leading contribution to the $\alpha,\beta=1,2$ components to be

[TABLE]

The general interpretation is depicted in figure 4, where we show the two- and one-loop diagrams that can be associated with the generation of secondaries. In this simple picture, the number of internal propagator lines encodes which modes become unstable as well as how fast they grow.

Since for systems out of equilibrium different loop diagrams can become important at different times, we consider a dynamical power counting and give parametric estimates for the time-scales of the different diagrams next. A general loop diagram containg $n$ vertices, $m$ propagator lines and $k$ condensate fields will scale parametrically as $\sim\zeta^{n-k}F^{m}$. In accordance with the above estimates, we consider diagrams that contain the most unstable mode as internal propagator lines, since these contributions are expected to dominate the dynamics.

At first loop-order, we have the diagram

for the lower field $\varphi_{1}$. This will play an important role in the dynamics as soon as it becomes parametrically of order one. We can estimate the time-scale by $1\simeq\zeta e^{2\gamma_{*}t}$, resulting in

[TABLE]

This is a conservative lower bound since we have assumed the most unstable exponential growth for all modes which is somewhat overestimating the effect. Around $t_{3}$, we find sizeable deviations from the behaviour of the most unstable and of the condensate mode. As a consequence, there are corrections of order one which come from diagrams with an arbitrary number of loops. Thus the dynamics is no longer characterized in terms of the small expansion parameter $\zeta$. This signals a nonperturbative regime as indicated above (see figure 2).

At second order in $\zeta$, we have the setting-sun diagram

for the lower field $\varphi_{1}$ shown in figure 4$a)$. Proceeding as for the previous diagram, we use $1\simeq\zeta^{2}e^{6\gamma_{*}t}$ to find the parametric estimate for the time-scale

[TABLE]

Remarkably $t_{2}<t_{3}$, i.e. the two-loop correction becomes important earlier than the previous one-loop diagram. This two-loop process can be associated with the generation of secondaries. Employing momentum conservation, three propagator lines with $k_{*}$ combine to $3k_{*}$ and give rise to a secondary growth rate of $6\gamma_{*}$.

Similarly, we can estimate the one-loop correction

for the excited field $\varphi_{2}$ to become important at the time

[TABLE]

which follows from $1\simeq\zeta^{2}\left(\zeta^{-1/2}\right)^{2}e^{4\gamma_{*}t}$. It is important to note that the diluteness parameter enters twice here, once for each vertex and as $\sim\zeta^{-1/2}$ for the non-vanishing field, such that the diagram is effectively of order $\zeta$. This process can be associated with the secondary growth at $2k_{*}$ in $\varphi_{2}$ with rate $4\gamma_{*}$ calculated above (see figure 4b).

Our findings for the time-scales $t_{1}<t_{2}<t_{3}$ agree well with the time evolution of the statistical two-point function plotted in figure 5. The secondary at $t_{1}$ in $\varphi_{2}$ starts before the one at $t_{2}$ in $\varphi_{1}$ and finally the system becomes non-perturbative at $t_{3}$. The observed growth rates are lower than predicted because nearby modes with different rates $\gamma(k)<\gamma_{*}$ contribute and effectively damp already the primary instability. In order to test our estimates for the secondary rates, we correct for this damping by choosing the largest pointwise derivative for the primary rate and find $\bar{\gamma}_{1}\approx 0.94\times\gamma_{*}$, which is slightly smaller than $\gamma_{*}$ as expected due to the harmonic trapping potential. Similarly, the secondary rates are estimated from the numerics as $\bar{\gamma}_{2}\approx 1.86\times\gamma_{*}\approx 1.98\times\bar{\gamma}_{1}$ and $\bar{\gamma}_{3}\approx 2.80\times\gamma_{*}\approx 2.98\times\bar{\gamma}_{1}$. Thus we find a quantitative agreement of the estimated secondary rates at the onset of the exponential growth after correcting for the systematic error of the primary instability. At later times, the interaction with other modes becomes non-neglible and consequently the approximated loop calculation does not describe the dynamics anymore. Instead the growth eventually stops. As a result, one can only observe the secondaries for one or two orders of magnitude in growth for the statistical two-point correlator.

The above reasoning is very general and also explains why secondaries are a generic feature of seemingly unrelated systems, such as ultracold atoms and cosmological inflation. The important ingredients are a primary instability, such as the tachyonic one considered in this work, and a four-vertex of the model Hamiltonian. In this sense the present system provides a minimal model to study non-linear phenomena related to secondaries. In particular, the mechanism is insensitive to the system being non-relativistic or relativistic as in the case of inflation.

VI Secondary Instabilities with Seed

In principle, the process of secondaries repeats itself and leads to the occupation of even higher momentum modes. It is this amplification mechanism that is proposed to enhance inflationary particle production over a wide range of momenta. It would be fascinating to test this non-linear aspect of inflaton dynamics in a real cold-atom experiment. Unfortunately, our simulations suggest that it is difficult to detect this effect in a realistic setup without further optimization.

To be specific, we have simulated the twin beam experiment of Bücker et al. (2011) for characteristic parameters and plotted the result in figure 6. The spectrum of $\varphi_{1}$ exhibits twin beams which have also been observed in experiment. The green dashed-dotted line shows the range of initial quantum fluctuations which we have cut off next to the primary instability. We have explicitely checked that the twin beams are insensitive to the choice of the cutoff as long as the fluctuations are present within the range of unstable modes. Our choice reveals the possible excitation of secondaries at the expected $\pm 3k_{*}$. They have not been observed in the experiment, probably because of the small occupation number compared to fluctuations which is due to particle number limitations ($N\approx 800$). The naive ansatz of increasing the number of available atoms does not solve the problem because it increases the width $\delta k$ of the primary instability. At the same time $k_{*}$ shifts to lower momenta and the instability eventually vanishes (see equations (16) and (15)). This results in a non-trivial optimization problem limited by accessible setup parameters.

The above results suggest that the most important ingredient for the generation of secondaries is a high occupation of the primary instability. Thus, we propose to seed the initial instability significantly. Instead of preparing the whole condensate in the excited state $\varphi_{2}$, we keep a fraction of all atoms in the lower state $\varphi_{1}$ and imprint a modulation to transfer them into the most unstable modes $\pm k_{*}$, i.e. we choose a coherent initial field value of

[TABLE]

where the fraction $0<\beta<1$ controls the seeding strength. Tuning the seed, we expect to effectively skip the regime of primary growth and immediately observe secondaries. Furthermore we have often employed momentum conservation in the previous section. Strictly speaking this is not given in the presence of a harmonic trapping potential (case $a)$). Recently, almost perfect box potentials have been achieved experimentally Gaunt et al. (2013). Thus we may eliminate the systematical error of not having momentum conservation by choosing box-type potentials (case $b)$) in the following. Figure 7 shows the direct comparison of the two scenarios with and without seed.

The situation without seed is similar to the twin beam case. Comparing the early occupation (in red) to the squared dispersion relation (inset) again demonstrates the validity of the Bogoliubov approximation. However, this time the primary instability is very wide which is due to a different parameter set and mainly results from a larger particle number. At later times, the twin beams wash out and we find again the spectrum that we have previously associated with non-perturbative dynamics.

This is to be compared to the situation with seed: The primary instability is highly occupied early on and consequently the Bogoliubov dispersion does not dominate the dynamics. Instead we observe clear peaks at odd multiples of the seeded mode. They are created from a cascade of scattering events as explained in section IV. Moreover the underlying dynamics is much faster than the primary growth such that the latter is suppressed in the experiment. One could argue that we are merely observing a tree level scattering process here. That we are dealing with an amplified nonlinear process can be verified via the growth rate of the secondary instabilities, which are shown in figure 8. Due to the presence of a large initial one-point function in the lower field $\varphi_{1}$, the seed shifts the dominant dynamics from the two-point $F$ to the macroscopic field. In fact, we have plotted the absolute value of the latter, where we clearly observe a higher growth rate similar to the unseeded case of the previous section. Moreover, the rate is on the order of $6\gamma_{*}$ which we have estimated in the previous section. We postpone the discussion of the observed exponential growth of the $5k_{*}$ mode until the end of the next section. Furthermore, we observe that the $2k_{*}$ and $3k_{*}$ modes show a very peculiar early-time behaviour that was not present in the unseeded case. A doubly logarithmic plot (not shown) reveals that time-dependence follows a power law $\propto t^{2}$. Since the nonlinear corrections become important only at sufficiently late times $\gamma_{*}t\gg 1$, a naive explanation for this effect is the following: If the fields admit a Taylor expansion with vanishing constant term, i.e.

[TABLE]

with $C=\varphi(0)=0$ and $C^{\prime}=\frac{\varphi^{\prime}(0)}{\gamma*}$, then $\left|\varphi(t)\right|^{2}\propto t^{2}$ for $\gamma_{*}t\ll 1$. We will verify this reasoning in a classical toy model in the next section.

VII Simplified Four-Mode Model

In the previous section, we have seen that the seeding procedure simplifies the dynamics in the sense that the evolution gets dominated by one-point functions. This motivates us to study the coupled equations (19a) and (19b) in more detail. We have also seen that the dynamics involve mainly a few modes, namely odd multiples of $k_{*}$ of the lower field $\varphi_{1}$ and even multiples of $k_{*}$ of the excited field $\varphi_{2}$. Consequently, we transform the coupled equations for box-type potentials to Fourier space and restrict them to $\varphi_{mk_{*}}^{(1)}$ with $m$ odd and $\varphi_{nk_{*}}^{(2)}$ with $n$ even. We further restrict ourselves to the first possible secondaries at $2k_{*}$ and $3k_{*}$, i.e. $m=\pm 1,\pm 3$ and $n=0,\pm 2$. Employing the symmetry $k\rightarrow-k$, we are left with a system of four coupled non-linear differential equations for the complex functions $\varphi_{1}\equiv\varphi^{(1)}_{k_{*}}(t),\varphi_{3}\equiv\varphi^{(1)}_{3k_{*}}(t),\varphi_{0}\equiv\varphi^{(2)}_{0}(t),\varphi_{2}\equiv\varphi^{(2)}_{2k_{*}}(t)$ of the form

[TABLE]

The precise form of $\mathcal{H}_{j}$ is given in appendix B. The model has four independent parameters: Three interaction constants $c_{1},c_{2},c_{12}$ and the dimensionless combination $\frac{k_{*}^{2}}{M\mu_{2}}$. Most importantly, it contains nonlinear terms which introduce momentum-conserving scattering processes between the different modes in accordance with the original Hamiltonian (1). We do not expect to find any numerical agreement between this highly simplified model and the full quantum dynamics. However, we suspect that in the presence of seeds the feature of secondaries can be qualitatively understood with this rather simple differential equation.

The initial value problem that we have considered before corresponds to a high initial occupation of $\varphi_{0}(t_{0})=\sqrt{\frac{\mu_{2}}{c_{2}}}$, which is determined by the stationary solution of the equations of motion with all other modes set to zero. All other initial values we set to very small values $\varphi_{1}(t_{0}),\varphi_{2}(t_{0}),\varphi_{3}(t_{0})\ll\varphi_{0}(t_{0})$. Then a linearised analysis for the most unstable mode yields

[TABLE]

which for the corresponding mode occupation number is solved by $\left|\varphi_{1}\right|^{2}=Ae^{2\gamma_{*}t}+Be^{-2\gamma_{*}t}$. Here we recover again the Bogoliubov result and notice that a non-vanishing initial occupation is necessary to trigger the primary instability. In the quantum system, these initial fluctuations are always present as the “quantum-half”. To keep the discussion general, we consider an arbitrary initial occupation $\left|\varphi_{1}(t_{0})\right|^{2}=A+B=\beta^{2}\left|\varphi_{0}(t_{0})\right|^{2}$. In accordance with our previous simulations, we choose the lower and upper fields to be oscillating in phase, i.e. $A-B=0$. Consequently, the primary instability is given by

[TABLE]

with a small (seeding) parameter $0<\beta<1$, which justifies the linearization at early times.

If the system shows a clear separation of scales, then the zero mode will stay relatively highly occupied, while the most unstable mode grows exponentially. Thus, we consider only terms with the highest powers of the most unstable mode in the remaining two equations of motions, i.e. we approximate

[TABLE]

Of course, these equations can only be valid for a certain intermediate regime, where other contributions can be neglected. As we will see below, they lead to exponential growth, such that terms involving the fields themselves will gain importance over time and invalidate the approximation at late times. On the other hand, the initial values of the fields have to be sufficiently small compared to $\varphi_{k_{*}}^{(1)}(t_{0})$ in order for the approximation to be valid early on. The precise regime of validity is rather difficult to estimate. Nevertheless, we will use this approximation to determine the exponential growth rates and justify the approximation afterwards numerically.

In fact, we can directly integrate equations (46) as

[TABLE]

where $F(x|m)$ denotes the elliptic integral of the first kind with parameter $m$. Expanding the solutions at late times $\gamma_{*}t\gg 1$, we recover the results of the loop calculation,

[TABLE]

To substantiate our analytical predictions, we show the results of a numerical simulation of the coupled evolution equations in figure 9. We observe very similar dynamics compared to the statistical propagator of the first TWA simulation (cf. figure 5). The primary $\varphi^{(1)}_{k_{*}}$ exhibits to very good accuracy the expected exponential growth rate of $2\gamma_{*}$ and even the growth rates of the secondaries $\varphi^{(2)}_{2k_{*}}$ and $\varphi^{(1)}_{3k_{*}}$ agree well with the estimates $4\gamma_{*}$ and $6\gamma_{*}$, respectively. We conclude that the simplified four-mode model encompasses the basic ingredients for the generation of secondaries, i.e. a primary instability, the necessary non-linear scattering terms and one highly occupied condensate mode $\left(\varphi^{(2)}_{0}\right)$ serving as a particle bath.

We now turn to the unexpected short-time dynamics that we have observed in the seeded case. Expanding the solutions (47) and (49) at early times $\gamma_{*}t\ll 1$ yields

[TABLE]

which explains the observation of a power law $\propto t^{2}$ seen in figure 8. To verify the analytics, we have plotted the results of another simulation in figure 10, where we have increased the initial occupation of the primary unstable mode in order to compare to the situation with seed. We cleary observe that the secondaries show initially a growth $\propto t^{2}$. Additionally, one can still fit an exponential growth to the $3k_{*}$ mode (not shown in the plot), but not to the $2k_{*}$ mode, which agrees with the observation of the seeded TWA simulation. The strong seed assures an initially small occupation of the secondary modes compared to the primary, which allows us to make the approximation of equation (46) at early times. In conclusion, the naive interpretation of the previous section, given in terms of a Taylor expansion with vanishing time-independent constant, indeed explains the behaviour of the secondaries at early times.

For values of $\beta$ close to one, i.e. the strongly seeded situation, the linearised analysis, which we employed to determine the primary instability, may be questioned. Since we cannot solve the system analytically, we can only give some heuristic comments on the effect of seeding: A high initial occupation of the most unstable mode breaks the translational invariance of the linearised equations of motion and leaves a discrete symmetry. One can then apply Bloch’s theorem: We consider e.g. the secondary growth $\varphi^{(1)}_{3k_{*}}(t)\propto e^{3\gamma_{*}t}$ observed in figure 8. In the spirit of the simplified few-mode model, this field will couple e.g. to the Fourier mode $k=9=3+3+3$ giving rise to the rate $\varphi^{(1)}_{9k_{*}}(t)\propto e^{9\gamma_{*}t}$, which follows from the nonlinear term $\left(\varphi^{(1)}_{(9-3-3)k_{*}}\right)^{*}\varphi^{(1)}_{3k_{*}}\varphi^{(1)}_{3k_{*}}$ by momentum conservation. However, in the seeded case, this growth rate will already manifest itself in a lower mode because we may identify modes in different Brillouin zones as $k\equiv k\mod 2k_{*}$. Of course, this influence is suppressed, but this suppression might be compensated by the high occupation numbers. In particular, the nonlinear term $\left(\varphi^{(1)}_{(5-3-3)k_{*}}\right)^{*}\varphi^{(1)}_{3k_{*}}\varphi^{(1)}_{3k_{*}}$ can be turned into $\left(\varphi^{(1)}_{-3k_{*}}\right)^{*}\varphi^{(1)}_{3k_{*}}\varphi^{(1)}_{3k_{*}}$ by replacing $-k_{*}\rightarrow-3k_{*}$ only once. Thus, also the mode $5k_{*}$ can show the high growth rate of $9\gamma_{*}$, which may explain the observation in figure 8 and concludes our analysis of secondaries.

VIII Conclusion and Outlook

In the present study we have demonstrated that the generic feature of the non-linear secondary growth of unstable modes, which was first predicted in the context of inflationary particle production, prevails in various cold-atom setups. We have analyzed possible choices of realistic parameter sets and indicated why recent twin beam experiments have not detected this effect. The main obstacle concerns particle number limitations. In order to make the secondary growth of unstable modes experimentally accessible, we have proposed a seeding procedure which effectively amplifies the underlying nonlinear corrections.

We have furthermore demonstrated that the growth rates of secondary amplifications allow an interpretation in terms of loop diagrams, which reflect the structure of interactions of the underlying Hamiltonian. This is a direct macroscopic visualization of (quantum) loop processes. In turn, one may deduce essential aspects of the interaction structure of an unknown Hamiltonian from measured instabilites. This provides a striking example where macroscopic observables determine fundamental microscopic properties of a many-body system.

In this work, we have focused on the early-time quantum dynamics including primary and secondary growth. Furthermore, we have pointed out a highly simplified model that essentially captures the same dynamical features. Our simulations and analytical estimates suggest that the full quantum model subsequently enters a non-perturbative regime characterized by the importance of diagrammatic contributions with an arbitrary number of loops. This regime cannot be simulated by a simplified few-mode model due to the strongly correlated nature of many interacting degrees of freedom. The $\mathcal{O}(N)$-symmetric inflaton model - showing the corresponding behaviour at early times - exhibits universal scaling dynamics associated to an attractor far from equilibrium during a similar non-perturbative regime Berges et al. (2008); Orioli et al. (2015). It will be interesting to study the possibility of observing a so-called nonthermal fixed point in a realistic cold-atom experiment. However, one spatial dimension is special due to phase space restrictions and near-integrability. The model considered in this work encorporates two-body scattering exchanging particle types between the two fields. The presence of this interaction can be thought of as arising from breaking the symmetry to a smaller subgroup $\mathcal{O}(4)\rightarrow\mathcal{O}^{\prime}(2)$, where $\mathcal{O}^{\prime}(2)$ denotes a combined rotation of the two fields $\phi_{1/2}$. This might facilitate the presence of a universal attractor far from equilibrium. In that context, we would also like to study the influence of the dimensionality of space-time for the universality class in more detail. We reserve a detailed analysis of the long-time evolution of the discussed two-field system for future studies.

Acknowledgements.

We thank T. Gasenzer, M. Oberthaler, A. Piñeiro Orioli, M. Prüfer, and J. Schmiedmayer for fruitful discussions. Significant parts of this work are taken from the master’s thesis of T.V.Z at Heidelberg University. This publication is part of and supported by the DFG Collaborative Research Centre ”SFB 1225 (ISOQUANT)”.

Appendix A General Bogoliubov Transformation

In this appendix we summarize the inhomogeneous Bogoliubov transformation for unstable systems (see also the appendix of Garay et al. (2001)). We consider a general $d-$dimensional Hamiltonian of the quadratic form

[TABLE]

with complex bosonic Heisenberg field operators satisfying canonical commutation relations

[TABLE]

and we require $H_{0}^{\dagger}=H_{0}$ for hermiticity. The general Bogoliubov transformation reads

[TABLE]

where $a_{n},a_{n}^{\dagger}$ are Schrödinger creation and annihilation operators spanning the Fock space in consideration. In general, they do not fulfill canonical commutation relations. Note also, that the eigenergies $E_{n}$ are in general complex numbers. To determine these, we consider the Heisenberg equations of motion in compact matrix form as

[TABLE]

This is a non-hermitian eigenvalue problem for the doublet of mode functions $y_{n}=(u_{n},v_{n})$ known as the Bogoliubov-de Gennes equations. Using the hermitian conjugate of this equation we derive the general orthogonality relation between two eigenfunctions $y_{n},y_{m}$ as

[TABLE]

where we have defined the inner product

[TABLE]

By direct computation, we observe that for an eigenvector $y_{n}$ with eigenvalue $E_{n}$, there exist further eigenvectors with related eigenvalues

[TABLE]

Labeling the corresponding quantum numbers as $n^{*},\tilde{n}$ and $\tilde{n}^{*}$, respectively, we note that the Bogoliubov expansion contains terms like

[TABLE]

showing a clear redundancy in the description. Redefining $a_{n}+a^{\dagger}_{\tilde{n}^{*}}\rightarrow a_{n}$ and dropping the corresponding $\tilde{\;}$ labels in the sum lifts this redundancy and we are left with two orthonormality relations

[TABLE]

which hold for every eigenfunction $y_{n}(x)=\begin{pmatrix}u_{n}(x)\\ v_{n}(x)\end{pmatrix}$ that solves the BdG equations. Now the eigenvalues appear in pairs $n,n^{*}$ with $E_{n^{*}}=E^{*}_{n}$. In the case of real $E_{n}$, we have $n=n^{*}$. We can use these relations to invert the Bogoliubov transformation as

[TABLE]

and determine the commutation relations as

[TABLE]

We emphasize that these are not canonical commutation relations. In general, it is impossible to diagonalise an arbitrary quadratic Hamiltonian while retaining canonical commutation relations. Finally applying the Bogoliubov transformation, we can rewrite the Hamiltonian in diagonal form

[TABLE]

where we have dropped an irrelevant constant. Note that the sum runs over all quantum numbers $n$ including their dual $n^{*}$ (for complex $E_{n}$).

Appendix B Simulation Details and Numerical Parameters

For all simulations, we have considered Rubidium atoms with atomic mass $M=$1.44\text{\times}{10}^{-25}\text{,}\mathrm{kg} and scattering length $a_{s}=$5.3\text{\,}\mathrm{nm}. The interaction strength is given by $g=\frac{4\pi\hbar^{2}a_{s}}{M}$. We have used the classical-statistical or TWA approach to simulate the quantum dynamics numerically. For completeness, we also show the occupation of the excited field in figures 11, 12 and 13.

B.1 Proof of Principle: Generation of Secondaries

The first simulation is based on the following parameters used to demonstrate the principle of secondaries instabilities. We used the particle number $N=$6\text{\times}{10}^{6}\text{,}, the harmonic frequency $\omega=$20\text{\,}\mathrm{Hz} of the longitudinal trapping potential, $600$ lattice points with equal spacing $a=$0.33\text{,}\mathrm{\SIUnitSymbolMicro m} and assumed a cylindrical trap in orthogonal direction with radial size $R=$2.7\text{\,}\mathrm{\SIUnitSymbolMicro m}. We have reduced the atomic interaction as $\tilde{g}=g/500$. The two fields $\varphi_{1,2}$ arise as coefficients in an expansion of the 3D field into radial eigenfunctions as

[TABLE]

where $\psi_{1}(\rho,\phi)=\psi_{0,1}(\rho,\phi)$ and $\psi_{2}(\rho,\phi)=\frac{1}{\sqrt{2}}\left[\psi_{1,1}(\rho,\phi)+\psi_{-1,1}(\rho,\phi)\right]$ are the energetically lowest lying states with $\psi_{l,n}$ denoting the single particle eigenfunctions of the cylindrical box. For simplicity, we implicitely assume a slight breaking of the cylindrical symmetry in order to lift the degeneracy of the first excited state. Then the interaction strengths are given by the following overlap integrals

[TABLE]

The chemical potentials are fixed by the total particle number and the energy difference between the one-particle states in the cylindrical potential, i.e.

[TABLE]

We have seeded the lower state $\varphi_{1}$ with quantum fluctuations and used the first 150 eigenfunctions of the harmonic oscillator as a controlled UV-cutoff. One simulation consists of 50 samples of initial values and their time evolution. We employed fixed boundary conditions and used a total simulation time of $T=25\frac{\hbar}{\mu_{2}}$ and readout steps of $dt=0.1\frac{\hbar}{\mu_{2}}$.

B.2 Analysis of the real Twin-Beam Experiment

The comparision to the twin-beam experiment was done with the following parameters: $N=800$ particles, harmonic frequency $\omega=2\pi\times$16.3\text{,}\mathrm{Hz} and harmonic radial traps with $\omega_{x}=2\pi\times$1.83\text{\,}\mathrm{kHz}$,\omega_{y}=2\pi\times$2.5\text{\,}\mathrm{kHz}, giving the oscillator lengths $a_{i}=\sqrt{\frac{\hbar}{m\omega_{i}}}$. We used $600$ lattice points with equal spacing $a=$0.15\text{,}\mathrm{\SIUnitSymbolMicro m}$$. This time the fields $\varphi_{1,2}$ are defined via

[TABLE]

with $\psi_{n}$ denoting harmonic oscillator eigenfunctions which results in the following interaction constants

[TABLE]

The chemical potentials are fixed as $\mu_{2}=$492\text{,}\mathrm{Hz}$\times h$ and $\mu_{1}=\mu_{2}+\hbar\omega_{x}$. The total simulation time is $T=10\frac{\hbar}{\mu_{2}}$ in steps of $dt=0.1\frac{\hbar}{\mu_{2}}$. We used the first $151$ eigenfunctions to seed the initial fluctuations in $\varphi_{1}$ for 100 samples and employed fixed boundary conditions.

B.3 Proposal of Seeding Procedure

For the direct comparison of the situation with and without seed, we used the following parameter set: $N_{2}=30000$ particles in $\varphi_{2}$, a system length of $L=$80\text{,}\mathrm{\SIUnitSymbolMicro m} (no longitudinal trap!) and a cylindrical trap size of $R=$3\text{\,}\mathrm{\SIUnitSymbolMicro m}. The same reasoning as for the first simulation gives the effective interaction constants $g_{1},g_{2},g_{12}$ from overlap integrals (without tilde, i.e. we used the true $g$ instead of $\tilde{g}$). We used $300$ lattice sites with equal spacing $a=$0.13\text{,}\mathrm{\SIUnitSymbolMicro m}$$. The chemical potentials read $\mu_{2}=g_{2}N/L$ and $\mu_{1}$ is fixed as in the first simulation. The total simulation time is $T=10\frac{\hbar}{\mu_{2}}$ in steps of $dt=0.05\frac{\hbar}{mu_{2}}$. We employed periodic boundary conditions, which is justified for sufficiently large homogeneous systems. One simulation consists of 100 samples of initial fluctuations.

Without seed: Only initial fluctuations in $\varphi_{1}$.

With seed: Fluctuations on top of seed modulation at $k_{*}\approx 15\times\frac{2\pi}{a}$ in $\varphi_{1}$ with fraction $\beta=\frac{1}{2}$ of $\varphi_{2}$.

B.4 Classical Toy Model

The simple four-mode model featuring secondaries derived from the coupled GP equations is given by

[TABLE]

with $\mathcal{H}_{j}\left[\varphi_{0},\varphi_{1},\varphi_{2},\varphi_{3}\right]=\mathcal{H}_{\text{lin}}^{(j)}[\varphi_{j}]+\mathcal{H}^{(j)}_{\text{NL}}\left[\varphi_{0},\varphi_{1},\varphi_{2},\varphi_{3}\right]$. The linear and nonlinear terms are given by

[TABLE]

and

[TABLE]

For the simulation we used units of $\mu_{2}$, which sets the physical time scale $\bar{t}=\frac{\hbar}{\mu_{2}}$. We consider a weak dimensionless coupling $\frac{c_{2}}{\mu_{2}}=$1\text{\times}{10}^{-10}\text{,} and $\frac{c_{1}}{c_{2}}=11$, $\frac{c_{12}}{c_{2}}=7$ which are related to physical couplings via $g_{X}=c_{X}L$ with a ficticious system length $L$ and $X=1,2,12$. Note that the growth rate $\gamma_{*}$ is independent of $L$. Furthermore, we chose $\frac{k_{*}^{2}}{2M\mu_{2}}=3$. The initial values of the non-condensate modes are all equal and suppressed with respect to the condensate mode by $\beta=${10}^{-5}\text{\,}.

Appendix C Quantum Evolution Equations

Rephrasing the two-field system of (1) in terms of effective action techniques, we consider the corresponding classical action

[TABLE]

where $\phi$ denotes classical fields in contrast to the Heisenberg operators $\varphi$. The dynamics of the quantum theory are encoded in the generating functional $Z[J,R]=\int\mathcal{D}\phi\;e^{i\left(S[\phi]+J\cdot\phi+\phi\cdot R\cdot\phi\right)}$ and $\mathcal{C}$ is the Schwinger-Keldysh contour. We employ the 2PI effective action $\Gamma\left[\bar{\phi},G\right]$, which is defined as the Legendre transform of $W[J,R]=-i\log Z[J,R]$ w.r.t. the sources $J,R$. In this formulation the quantum equations of motion for the macroscopic field $\bar{\phi}$ and the full propagator $G$ are given by

[TABLE]

In addition to the macroscopic field

[TABLE]

it is convenient to introduce the statistical propagator $F$ and the spectral function $\rho$ as

[TABLE]

which are the symmetric $(F)$ resp. anti-symmetric $(\rho)$ parts of $G$ w.r.t. $\mathcal{C}$. Note that here and in the following $\alpha,\beta=1,2,3,4$ label $\varphi_{1},\varphi_{1}^{\dagger},\varphi_{2},\varphi_{2}^{\dagger}$. In this basis, one can rewrite the quantum evolution equations as

[TABLE]

supplemented by an evolution equation for the macroscopic field, which is not relevant for the present dicussion. The differential operator $\mathcal{D}$ has a contribution from the inverse classical propagator

[TABLE]

where in $4\times 4$-matrix notation the operator

[TABLE]

and $\sigma^{j}$ denote Pauli matrices and we have abbreviated

[TABLE]

and the operator

[TABLE]

As such, the equations (26) & (27) are exact. In order to make progress one has to specfiy an approximation for the self-energies $\Sigma$, which are split into local $\Sigma^{(0)}$ and non-local parts $\Sigma^{\rho/F}$. The superscripts $F$ resp. $\rho$ denote symmetric resp. anti-symmetric contributions to $\Sigma$ w.r.t. the closed time-contour.

Setting $\Sigma=0$ (tree-level) and fixing the macroscopic fields $\bar{\phi}_{2}=\sqrt{\frac{\mu_{2}}{g_{2}}}$, $\bar{\phi}_{1}=0$, we find

[TABLE]

where $H^{\prime}=H_{0}+2g_{12}\frac{\mu_{2}}{g_{2}}-\mu_{1}$ and similiar for $F_{21},F_{22}$ from which we recover the Bogoliubov result (14). As long as $\bar{\phi}_{1}=0$, the $1,2$ components decouple from the $3,4$ components of $F$ and e.g. $F_{31}\equiv 0$. We can derive a similar equation,

[TABLE]

from which one can read off the dispersion relation for the excited field as

[TABLE]

which is manifestly stable.