Macroscopic Quantum Tunneling Escape of Bose-Einstein Condensates

TL;DR

This paper investigates the non-exponential decay in macroscopic quantum tunneling of Bose-Einstein condensates, combining experiments with theoretical models to reveal the role of interactions and many-body effects.

Contribution

It introduces a modified WKB model accounting for mean-field effects and compares mean-field and many-body dynamics using different simulation methods.

Findings

Non-exponential decay observed in experiments

Modified WKB model aligns with experimental data

Many-body effects suggest limitations of mean-field theory

Abstract

Recent experiments on macroscopic quantum tunneling reveal a non-exponential decay of the number of atoms trapped in a quasibound state behind a potential barrier. Through both experiment and theory, we demonstrate this non-exponential decay results from interactions between atoms. Quantum tunneling of tens of thousands of 87 Rb atoms in a Bose-Einstein condensate is modeled by a modified Jeffreys-Wentzel-Kramers-Brillouin model, taking into account the effective time-dependent barrier induced by the mean-field. Three-dimensional Gross-Pitaevskii simulations corroborate a mean-field result when compared with experiments. However, with one-dimensional modeling using time-evolving block decimation, we present an effective renormalized mean-field theory that suggests many-body dynamics for which a bare mean-field theory may not apply.

Figures

Tables

Table 1. Table 1: Experimental and Estimated Parameters. From the weak configuration, with peak height, saddle height, horizontal distance from local potential minimum to saddle point ( x 0 subscript 𝑥 0 x_{0} ) and initial total particle numbers trapped behind the barrier ( N 0 subscript 𝑁 0 N_{0} ).

Peak [nK]	Saddle [nK]	$x_0$ [$\mu m$]	$N_0/{10}^3$
$146\pm 10$	$46\pm 3$	$16\pm 1$	$423\pm 21$
$170\pm 11$	$52\pm 10$	$17\pm 2$	$548\pm 37$
$190\pm 13$	$58\pm 4$	$18\pm 1$	$559\pm 33$
$221\pm 16$	$69\pm 5$	$20\pm 1$	$630\pm 26$
$260\pm 18$	$81\pm 6$	$21\pm 1$	$583\pm 32$
$299\pm 21$	$92\pm 6$	$23\pm 1$	$663\pm 54$

Equations

$$V(x,y,z)$$

$$V_{\mathrm{b}}(x,y,z)$$

$$\omega (z)$$

$$\Gamma ={\Gamma }_{\mathrm{bg}}+\exp (\alpha +\beta \mu ).$$

$$V_{\mathrm{tri}}(x)=\{\begin{array}{cc}-\frac{V_{\mathrm{s}}}{w}|x\pm x_0|+V_{\mathrm{s}},& |x\pm x_0|<w\\ 0,& \text{otherwise}\end{array}$$

$$P_{\mathrm{tri}}$$

$${\tau }_{\mathrm{tri}}$$

$${\Gamma }_{\mathrm{tri}}=\frac{P_{\mathrm{tri}}}{{\tau }_{\mathrm{tri}}}.$$

$$\dot{N}(t)=-({\Gamma }_{\mathrm{tri}}+{\Gamma }_{\mathrm{bg}})N(t)=-\Gamma N(t),$$

$$\mu =b{N}^{1/3},$$

$$\displaystyle\begin{split}&V_{\mathrm{eff}}(x,t)=\left\{\begin{matrix}-\frac{V_{\mathrm{mf}}}{w}|x\pm x_{0}|+V_{\mathrm{mf}},&|x\pm x_{0}|<w,\\ 0,&\text{otherwise},\\ \end{matrix}\right.\\ &V_{\mathrm{mf}}=V_{\rm s}\bigl{[}1+a\frac{N(t)}{N_{0}}\bigr{]}.\\ \end{split}$$

$$P_{\mathrm{eff}}$$

$${\tau }_{\mathrm{eff}}$$

$$\hat{H}=-J\sum _{i=1}^{L-1}({\hat{b}}_{i+1}^{†}{\hat{b}}_i+h.c.)+\sum _{i=1}^L[\frac{U}{2}{\hat{n}}_i({\hat{n}}_i-\hat{1})+V_i^{\mathrm{ext}}{\hat{n}}_i].$$

$$i\hslash \dot{{\psi }_i}=-J({\psi }_{i+1}+{\psi }_{i-1})+g|{\psi }_i{|}^2{\psi }_i+V_i^{\mathrm{ext}}{\psi }_i.$$

Full text

I Introduction

Quantum tunneling, first studied in $\alpha$-decay, is one of the most significant and earliest effects observed in quantum mechanics Gurney and Condon (1929). Tunneling provided a basis for the theory of molecular spectra in the 1920s, when it was demonstrated that transitions between chiral isomers occurred at measurable rates Hund (1927). This applies to biological systems, where the transition rate between isomers is slow enough to allow stable life, and even to cosmology, where the transition for polar molecules such as ammonia shows a measurable rate Anderson (1972). Tunneling has applications in many different systems, from two-proton decay and double beta decay in nuclear physics to the interdisciplinary study of tunneling in enzymes Rommel et al. (2012), with roots in biology, chemistry, and physics. As electronic devices reach the nanoscale, quantum tunneling will play a larger and larger role in understanding and developing nanoelectronics, such as tunneling diodes Choi et al. (2010). Moreover, with the advent of Josephson junctions, the epitome of macroscopic quantum devices based on tunneling, we can now measure voltage with unprecedented accuracy Dresselhaus et al. (2011).

The definitions and details of macroscopic quantum tunneling (MQT), quantum tunnelling at macroscopic scales, is discussed in Section II.1. Studying MQT in the context of Bose-Einstein condensates (BECs) presents many advantages in terms of both fundamental explorations and future MQT device design. First, BECs offer a high degree of controllability: interactions can be tuned over seven orders of magnitude using Feshbach resonances Pollack et al. (2009). Second, experimental advances in radio-frequency (RF) magnetic traps Fernholz et al. (2007), as well as optical trapping Ashkin et al. (1986), allow for greater access to controllable experiments needed to study MQT. Third, BECs enable manipulation of many-body states Bakr et al. (2010) that are inaccessible in other experimental settings. Many-body simulations, which we discuss next, elucidate the importance of this point. Fourth, BECs have controllable statistics and spatial dimensions. Fifth, BECs allow site-resolved microscopy in the context of optical lattices in 1, 2, and 3 dimensions. Sixth, atomic interferometry now permits observation of up to tenth order correlators for quasi-1D systems Langen et al. (2015).

To complement experimental advances, numerically exact simulations of bosonic Josephson junctions using the Bose-Hubbard model demonstrate a substantial deviation from mean-field theory Sakmann et al. (2009). However, until recently MQT has largely been treated with semiclassical approximations, such as in a double well Albiez et al. (2005a); Raghavan et al. (1999a), in Landau-Zener tunneling Cristiani et al. (2002), and in escape tunneling Moiseyev et al. (2004); Carr et al. (2004); Dekel et al. (2007). Theoretical methods for these semiclassical estimations include Jeffreys-Wentzel-Kramers-Brillouin (JWKB) and instanton approaches, along with the nonlinear Schrödinger equation (NLS) and hydrodynamic formulations thereof. Beyond these mean-field and instanton techniques, matrix product state methods such as time-evolving block decimation (TEBD) and multi-configurational time-dependent Hartree for bosons (MCTDHB) theory both provide numerical solutions of many-body dynamics. For example, the use of TEBD to simulate the Bose-Hubbard model for superfluid decay, has confirmed numerical limits on instanton computations Danshita and Polkovnikov (2010, 2012). Non-Hermitian quantum mechanics, which is frequently used in scattering problems and now applied to tunneling problems as an effective model, is broadening the view of tunneling phenomena, such as its use in asymmetric tunneling and interchain pair tunneling Santos and da Silva (2014); Zhang et al. (2017). Another possible many-body method, the time-adaptive MCTDHB, was used in our Josephson example above Sakmann et al. (2009). MCTDHB has also been applied to the quantum escape problem Beinke et al. (2015); Lode et al. (2012a); though this work examined depletion, the method has not yet produced predictions for entanglement. Measures such as entropy, entanglement, and correlations help illustrate when semiclassical or lower-order mean-field approximations fail Diego A. Alcala (2016).

In this Article, we present a combined theoretical and experimental study of MQT. Single-particle quantum tunneling can be modeled with the Schrödinger equation, with well-known exponential decay in the number of atoms trapped over time Gurney and Condon (1929). However, we have performed a macroscopic experiment exhibiting non-exponential decay of BEC tunneling from a single trapping well to unbound space Potnis et al. (2017). Figures LABEL:sub@fig:Torontoexp:a and LABEL:sub@fig:Torontoexp:b sketch the experimental process and the single well trap. Figure LABEL:sub@fig:Torontoexp:c displays raw data from the experiment, presenting an example of the observed non-exponential decay. The experimental data demonstrates that atomic interactions have participated in the tunneling process. Inspired by this result, we develop an alternative theoretical model in which the interactions cause the barrier to change dynamically, leading the decay to deviate from the single-particle exponential case. Subsequently we suggest an alternate interpretation in terms of an effective or renormalized mean-field theory drawn from TEBD simulations and accounting for the effects of condensate fragmentation and depletion, which can be tested in future MQT experiments.

This Article is organized in the following manner. We offer a brief discussion covering the nuances and regimes of quantum tunneling in Section II, including four subtopics: statistical properties, the role of interactions, type of trapping potential, and dimension of the system. The aim of this discussion is to briefly touch on the vast tunneling landscape and spark interest and ignite ideas for other physicists intrigued by the many open questions in MQT research. Next, the details of the MQT experiment are covered in Section III, with its experimental settings, results, and a 3D mean-field simulation. Furthermore, we present a case study analysis that exhibits excellent agreement with the experimental non-exponential decay result and further verifies the assumption of mean-field dynamics. In Section IV we then address the question of simpler effective 1D models. The experiment proceeds through three distinct regimes: initial transient classical spilling, quantum tunneling, and decay dominated by background loss; we model the tunneling and decay regimes, the main subject of this Article. We use a modified JWKB method in which the inter-particle interaction is taken into account via an effective mean-field interaction parameter, and we include the background loss in the tunneling rate. The model contains two fitting parameters which help to illuminate the experimental findings. In Section V, we propose an explanation for the effective mean-field used in the JWKB model. Motivated by advances in nonlinear optics, we demonstrate how a mean-field model can effectively reproduce many-mode many-body dynamics of the quantum tunneling process of a meta-stable state. Finally, we summarize our conclusions and future research avenues in Section VI.

II Regimes of macroscopic quantum tunneling

There are several factors that can affect the behavior of quantum tunneling, which may push the tunneling process into distinct tunneling regimes. These factors include, but are not limited to, the following. First, the statistical properties of the trapped particles are a primary concern, e.g. whether exploring bosonic, fermionic, or anyonic tunneling Zürn et al. (2013); del Campo et al. (2006), and whether the system can be treated as few-body or macroscopic. Second, interaction between particles – varying from zero to non-zero Yu et al. (2012); Abbarchi et al. (2013); Lee et al. (2003); Raghavan et al. (1999b); Shelykh et al. (2008), weak to strong Salasnich et al. (2001); Carr et al. (2005); Streltsov et al. (2008) and attractive to repulsive Ueda and Leggett (1998); Diego A. Alcala (2016) – may suppress or enhance tunneling, or even reform the barrier Aghamalyan et al. (2013). Third, the shape of the potential well and assorted lattice and other geometries induce and alter quantum phases, such as in spin chains and disordered systems Lewenstein et al. (2007); Bloch et al. (2008). Fourth, the dimension of the system can give rise to extra degrees of freedom which can for instance introduce chaos in the tunneling dynamics Mitchell and Steck (2007); Steck et al. (2002).

These factors, together with the notion of macroscopicity, delineate the extensive territory of MQT into distinct tunneling behaviors. The term macroscopic may connect to an intuitive concept of “large” objects. In fact the “macroscopic” in MQT has abundant aspects beyond one’s first impression. In the following, we lay out some of these possibilities, as a motivation for our own particular investigation, and to inspire future experiments and theory on the plethora of MQT phenomena, mostly still unexplored.

II.1 Macroscopicity, Statistical Properties, and Quasiparticles

Macroscopicity, which is a prominent characteristic of MQT, features at least four facets. First, macroscopicity separates phenomena in a semiclassical phase space, but these phenomena are not necessarily describable by classical mechanics. However, sometimes an extension of the semiclassical limit is apt. For example, the phase-space analogy for a driven double-well BEC in a quantum two-mode approximation suggests that it can be mapped to a macroscopic superposition state of two pendulum rotor states Mahmud et al. (2005). Second, macroscopicity can refer to a large quantity of particles, can activate many degrees of freedom, and can lead to emergent behavior. One example is ultracold bosons on a ring, where from the weakly interacting semiclassical limit in terms of dark solitons to the strongly interacting limit in terms of Yrast states, phase coherence can break down and phase slip enables continuous winding and unwinding of the system Kanamoto et al. (2008, 2009, 2010). Third, macroscopicity can be associated with the notion of massive objects, where quantum mechanics meets gravity and may lead toward quantum gravity. One theoretical avenue is the Schrödinger-Newton equation, the Schrödinger equation with an additional gravitational potential term Diósi (1984); Penrose (1996); SNe (1998). Fourth, macroscopicity may exemplify complexity, which is garnering wider attention. For instance, recent studies of complex networks relying on quantum mutual information examine critical points for transverse Ising and Bose-Hubbard models Vargas and Carr (2015).

Tunneling bears a close relation to transport. Loosely speaking, transport is above the barrier (often due to a “push”) whereas tunneling is through the barrier. Thus, one can consider tunneling of mass, charge, spin, etc. Quantum tunneling is applicable far beyond the original concept of particle tunneling. For instance, we may speak of tunneling of magnetization in ferromagnetic films, or in a spin-1 BEC. Studies of the transverse spin wave, with and without superflow, show different conditions for a certain anomalous tunneling behavior where tunneling occurs without reflection Watabe and Kato (2009); Watabe et al. (2011). We can also consider the tunneling between molecular states, such as ammonia and other pyramidal molecules Anderson (1972); Tserkis et al. (2014); Jona-Lasinio et al. (2002), and the Josephson junction serves as the simplest possible model of such processes Zapata et al. (1998).

The statistical properties of particles can be divided into fermions and bosons, and in lower dimensions include anyons; whether they contain only a single component or are a mixture of components; and whether the basic object exhibiting MQT is best described as emergent quasi-particles like vortices, solitons, or skyrmions. Most such possibilities are accessible in cold quantum gases. The studies of few-particle system have already paved the way to some of these distinctions, especially between bosons and fermions. Bose and Fermi Tonks-Girardeau gases show different deviations from exponential decay at short times as a direct consequence of different ground state energies del Campo et al. (2006), shown in Fig. LABEL:sub@fig:schematic:b. In few-fermion tunneling, pairing needs to be considered, a distinction from the boson case Zürn et al. (2013). The influence of the inter-atomic interaction in these systems is also investigated, details in Section II.2. Non-exponential decay of ultracold single-atom tunneling is expected to occur in certain parameter regimes García-Calderón and Romo (2016).

The tunneling of a BEC mixture or multi-component BEC also broadens the possibilities for MQT. A mixture can be divided into several cases: several different atomic species, each Bose-condensed; different internal states of the same atomic species, such as different hyperfine states of 87Rb; or use of different external states of a trap, in for instance MQT of quantum vortices García-March et al. (2011). For example, a study of MQT describes the mixing of two weakly-linked superfluids of interacting fermions, making it possible to obtain atomic Josephson junction equations describing the system as a whole Salasnich et al. (2008). In general the MQT of domain walls, skyrmions etc. in such mixtures may carry different information and different decoherence properties from a bulk BEC.

Understanding the nature of the particles involved is a vital step in understanding any tunneling regime. Calculations can be remarkably simplified, for instance, if a quasi-particle representation is valid, or if a many-body system can be analyzed on a single-particle basis or a single-particle-like basis. Sometimes, a quasiparticle also induces new phenomena. For instance, the tunneling of nonequilibrium quasiparticles through a Josephson junction brings decoherence and produces energy decay in superconducting qubits and resonators Martinis et al. (2009); Catelani et al. (2012). However, strongly-interacting systems may not exhibit quasiparticle representations, as we will discuss in the next section.

II.2 The Role of Interactions

Inter-atomic interaction poses crucial considerations for quantum tunneling. A well-known example is Josephson dynamics Abbarchi et al. (2013); Lee et al. (2003); Raghavan et al. (1999b); Shelykh et al. (2008), which occurs in systems of weakly coupled macroscopic quantum states. For weak interaction strengths, tunneling dominated effects such as the AC and DC Josephson effects Albiez et al. (2005b); S. Levy and Steinhauer (2007) and coherent temporal oscillations Yu et al. (2002) occur, where condensates tunnel between two wells. When the interaction exceeds a critical value, the populations become self-trapped Raghavan et al. (1999b), where tunneling is suppressed and condensates are mainly located in one potential well. Similar effects also happen in exciton and polariton condensates Shelykh et al. (2008). However, by biasing a more strongly-interacting system it is still possible to reduce tunneling times by tens to hundreds of orders of magnitude Dounas-Frazer et al. (2007); Carr et al. (2010).

Repulsive and attractive inter-atomic interactions can lead to different tunneling dynamics. Again, few-body tunneling presents several precursors to these effects in many-body systems. Fermionization of two distinguishable fermions occurs during the tunneling dynamics of a repulsively interacting system Zürn et al. (2012); Rontani (2012); Gharashi and Blume (2015). Sequential single-particle tunneling is found in both repulsive and attractive systems Lundmark et al. (2015), while pairing phenomena can be investigated in the strongly attractive interaction region Zürn et al. (2013); Lundmark et al. (2015); Gharashi and Blume (2015). As inter-particle interaction varies from strongly attractive to strongly repulsive, tunneling rates diverge within a wide range of orders of magnitude Lundmark et al. (2015); Gharashi and Blume (2015).

The strength of the interaction in a system determines the appropriate tunneling theory. For instance, bosons with weak interactions can often be described by mean-field theory, where the many-body wavefunction is dominated by a semi-classical field, i.e., a complex scalar field Ueda and Leggett (1998); Salasnich et al. (2001); Carr et al. (2005). Even weak interactions produce an effective barrier, which deforms as the wavefunction escapes, as depicted by the blue curve in Fig. LABEL:sub@fig:schematic:a, where repulsive atomic interactions change the barrier the wavefunction encounters, resulting in a non-exponential decay.

Increased interaction strength can lead to correlations, fluctuations, and entanglement Osterloh et al. (2002), which render mean-field theory ineffective. We include here a non-exhaustive list of resulting effects on MQT. First, the BEC may be fragmented or depleted. When the energy of the system exceeds a threshold as a consequence of strong interactions, fragmentation – macroscopic occupation of more than one mode of the single-particle density matrix – is induced Streltsov et al. (2008). The tunneling of an initially parabolically-trapped ultracold Bose gas in a coherent state into open space develops fragmentation after some propagation time Lode et al. (2012b); Diego A. Alcala (2016); the fragmented components are not said to be phase-coherent, as the notion of phase is tied to a single macroscopically occupied mode. Thus semiclassical wave theory becomes less and less applicable, and tunneling does not proceed in the same manner as the original concept. Nevertheless, particles with an energy below that of a potential barrier tunnel, so the term, even in this non-semiclassical context, still applies. Apart from fragmentation, condensates can also deplete, where there is non-macroscopic occupation of many modes Diego A. Alcala (2016). Second, interactions can enhance or decrease tunneling rates. The tunneling rate of a quasibound many-body state is sped up (slowed down) by repulsive (attractive) interactions Diego A. Alcala (2016). Third, fluctuations affect the tunneling rate. Quantum fluctuations of the Josephson-Leggett mode in a Josephson junction drastically enhances the MQT escape rate Ota et al. (2011); Asai et al. (2014). Finally, dissipation may also alter the tunneling behavior. In the Josephson-Leggett mode, quantum dissipation suppresses MQT Asai et al. (2014).

For strong interactions, 1D bosons can show similar transport properties as noninteracting fermions, due to the boson-fermion dual representation in 1D Qian et al. (2011). A famous example is a Tonks-Girardeau gas, a strongly interacting system, where bosonic atoms exhibit fermion-like behavior Girardeau (1960). A non-exponential decay regime arises in a bosonic Tonks-Girardeau gas for short times, which exhibits few-body decay features. A fermionic Tonks-Girardeau gas, on the other hand, shows bosonization and deviations from exponential decay at long times del Campo et al. (2006). As an extreme example of strong interactions, consider the unitary Fermi gas, and its holographic dual in a weakly curved gravitational representation Adams et al. (2012). The MQT of a unitary Fermi gas remains an exciting and untested prediction of dynamical holographic approaches Adams et al. (2014).

II.3 Trapping Potential

The potential dictates the tunneling environment, and plays an essential role in characteristic tunneling times. An important case across many fields of physics, chemistry, etc. is tunneling from one or a few single-well discrete modes to continuous modes in free space, i.e., the quasi-bound or escape problem. This is exactly the system studied in this Article, where in our case tunneling escape from a single well is assisted by inter-atomic interactions Potnis et al. (2017). For a harmonic trap, power law behavior in MQT of BECs is observed near the critical point of collapse, while for an anharmonic trap, there is no power law behavior Haldar et al. (2013).

Double well tunneling stands in strong contrast to the quantum escape problem, since for instance a small bias in the system can lead to an exponential suppression of tunneling. Both wells contain discrete modes for atoms, thereby the tunneling process describes motion from discrete modes in one well to discrete modes in the second well. Josephson effects Abbarchi et al. (2013); Lee et al. (2003); Raghavan et al. (1999b); Shelykh et al. (2008); Albiez et al. (2005b), including the well-known AC and DC Josephson effects Albiez et al. (2005b); S. Levy and Steinhauer (2007), which arise in double-well-type systems, epitomize this form of MQT. Related contexts in this direction range from adiabatic transport of BECs  Nesterenko et al. (2009) to interaction supported transport of BECs Nesterenko et al. (2014), and polarized fermion tunneling in 3D Macrì and Trombettoni (2013). Under common assumptions, only a few modes in each well are required for consistent analysis with experiments.

In addition to these examples, other potentials to consider include periodic potentials, like lattices. In a static periodic lattice potential, the expansion of matter waves is quadratic for short times Creffield et al. (2010). Other than potential shapes, barrier materials also modify the tunneling process. Comparing two different barrier materials, pillared silicon and aerogel, quantum reflection is suppressed differently by mean-field interactions at low velocity Pasquini et al. (2006). In addition to the bare potentials above, there are also dressed potentials, such as radio frequency dressed potentials Zobay and Garraway (2001). Studies show that on-site interactions in neutral atoms are dramatically enhanced in radio frequency dressed traps, which can affect the tunneling behavior Shotter et al. (2008).

Beyond static potentials, dynamic potentials, such as those directly influenced by a time-dependent force, provide a unique milieu for tunneling dynamics. A ring trap with a driven space-time reflection symmetry breaking (PT symmetry) potential can induce chaos in the system Heimsoth et al. (2013); Valdez et al. (2016). Shaking a lattice can reduce the tunneling rate, and even completely suppress it for certain values of the shaking or driving parameter Lignier et al. (2007); Creffield (2007); Sias et al. (2008); Creffield et al. (2010); the same phenomenon also occurs in a driven double well system Grossmann et al. (1991). In a single potential well with a periodically modulated amplitude, dynamical tunneling is suppressed by nonlinear interaction, in comparison to the non-interacting case Wüster et al. (2012). A less obvious dynamic barrier can form indirectly by other means. For instance, interactions between particles can deform the barrier, rendering a previously static barrier dynamic, as shown in Fig. LABEL:sub@fig:schematic:a.

II.4 Dimension of the System

Tunneling usually occurs at the weakest points of a trap, where the escape rate is the largest and most favorable. For instance, in this Article, tunneling proceeds through the weakest points of the potential well, the two saddle points as depicted schematically in Fig. LABEL:sub@fig:Torontoexp:b. This is similar to a dam breaking; water begins to spill out from the structurally weakest points. In this view, tunneling can be recast as a 1D or quasi-1D problem. But, higher dimensions can elicit new phenomena. One such phenomenon is chaos, which can significantly complicate the dynamics Mitchell and Steck (2007). The cluttered paths in Fig. LABEL:sub@fig:schematic:c sketch possible chaos in the experimental trap, and dashed lines at the saddle points indicate MQT. Chaos-assisted tunneling sometimes gives rise to tunneling oscillation Steck et al. (2002), and produces irregular fluctuations in the tunneling rate Mouchet et al. (2001); Tomsovic and Ullmo (1994). Another consequence, due to the presence of extra degrees of freedom, is incoherent oscillations of a large number of polarized fermions in a 3D double well Macrì and Trombettoni (2013).

III A Macroscopic Quantum Tunneling Experiment

We have made the first step in reporting a non-exponential decay in a single-well tunneling experiment due to inter-atomic interactions Potnis et al. (2017). In this section, we offer a description of the experiment, present previously unpublished findings from both the experiment and supporting 3D mean-field simulations showing reasonable agreement with the data, and finally exhibit a case study of the experiment using an effective 1D modified JWKB model derived in detail in Section IV.

III.1 Experimental Design

The experiment studied tunneling of a 87Rb BEC prepared in a quasi-bound state with repulsive inter-atomic interaction. The trapping potential had harmonic confinements in the $x$ and $z$ directions with trapping frequencies $\omega_{x}$ and $\omega_{z}$. Due to gravity and a magnetic field gradient, there was a tilt in the $y$ direction (vertical direction) with a constant acceleration $a$. A slice through the $y$–$z$ plane of the potential well is shown in Fig. LABEL:sub@fig:Torontoexp:b, with complete trapping potential as follows:

[TABLE]

Here ${V_{0}}$ is the peak height of the barrier, $z_{\rm R}=8\,\mu$m is the Rayleigh range, $\omega_{0}$ is the barrier waist, $\omega_{0}=1.3(0.1)\,\mathrm{\mu m}$, and $\omega(z)$ is the Gaussian beam waist.

The experiment uses two trapping configurations: weak and tight. The parameters of the trap in the weak configuration are $\omega_{x}=2\pi\times 32.7(0.24)\rm{\,Hz}$, $\omega_{z}=\omega_{x}/2$, and $a=2.08(0.04)\mathrm{\,m/s^{2}}$, while for the tight configuration, $\omega_{x}=2\pi\times 86.6(0.6)\rm{\,Hz}$, $\omega_{z}=\omega_{x}/2$, and $a=8.40(0.06)\mathrm{\,m/s^{2}}$. As a result of a tighter confinement in the latter, the initial total number of atoms trapped in the potential is around one fifth of that in the weak configuration. This decrease of particle number affects the validity of the Thomas-Fermi approximation in our previous mean-field calculation Potnis et al. (2017), and requires an offset of the kinetic energy term. The details of this correction is described in Section IV.2. Although the first report of our experiment focused on the tight trapping configurations Potnis et al. (2017), our theory starts with the weak configuration and then covers the tight configuration afterward, expanding the reach of our study to different regimes.

Figure LABEL:sub@fig:Torontoexp:a depicts the time sequence of the experimental procedures. First, a cloud of 87Rb atoms in the $\mid F=2,m_{F}=2\rangle$ ground state is loaded into the trapping potential from a hybrid trap. Next, evaporative cooling techniques lower the temperature of the atoms until they form a BEC. Then, the barrier height is ramped down non-adiabatically in 20 ms (5 ms) for the weak (tight) configuration. After ramping procedure, the condensate is held in the trap for a variable time from 0.1 ms to 1.2 s. Finally, the trapping potentials are abruptly turned off and the cloud is imaged after time-of-flight (TOF) expansion.

Table 1 shows major parameters in the weak configuration. The relative saddle heights are about one-third of their corresponding peak heights, where we note that MQT occurs through the two saddles or weak points of the potentials, not the peak. All the plots in this Article reference peak heights. In the trap, $x_{0}$ is the horizontal distance from the bottom (local minimum point) to one saddle point. We will see in Section IV that our theoretical fit parameters agree well with this parameter. $N_{0}$ is the total particle number in the trap at the beginning of experimental observation of the dynamics, after the non-adiabatic ramp down of the barrier. Figure LABEL:sub@fig:Torontoexp:c shows a typical raw data set from the experiment. We find clear non-exponential decay in the number of atoms. Further presentation and analysis of the data is deferred to the following two sections.

III.2 Experimental Data and 3D Mean-Field Model

Data for the chemical potential as a function of atom number and time dependence of the number of atoms trapped in the quasi-bound state are shown in Fig. 3 and 4 for both weak and strong trapping configurations. Mean-field simulations are then performed for comparison, using the split-step operator method for the 3D Gross-Pitaevskii Equation (GPE); imaginary time evolution is used to find the ground state. Measured trap parameters and initial atom numbers are used as input parameters to the simulations without any free parameters. To mimic the experimental procedure in Fig. LABEL:sub@fig:Torontoexp:a, the barrier height is linearly ramped down to a final value. Absorbing boundary conditions are introduced to avoid reflections of escaped atoms from the edge. Finally, the chemical potential extracted from the experiment uses a Thomas-Fermi approximation Potnis et al. (2017). This leads to good agreement in the weak trapping case in Fig. LABEL:sub@fig:mu3D:a but a systematic error in Fig. LABEL:sub@fig:mu3D:b where the Thomas-Fermi approximation is less applicable due to a larger kinetic or zero-point energy. In contrast, in Fig. 4 the mean-field time sequence is qualitatively close but is up to a factor of 2 off in both trapping configurations. These studies suggest a mean-field model can reproduce the main features of this experiment, but may bear some correcting either due to characterization of the trapping potential or due to other beyond mean-field effects, as we will explore in Sec. V.

Both experimental data and mean-field predictions clearly show non-exponential decay in both trapping configurations. Figure LABEL:sub@fig:Nt3D:a covers the weak configuration with barrier heights $V_{0}$ = 260(18) nK (blue circles), $V_{0}$ = 190(13) nK (purple diamonds) and $V_{0}$ = 170(11) nK (teal squares), and Fig. LABEL:sub@fig:Nt3D:b the tight configuration with barrier heights $V_{0}$ = 330(35) nK (red circles), $V_{0}$ = 290(30) nK (blue diamonds) and $V_{0}$ = 240(25) nK (green squares). From bottom to top, 3D GPE simulations are shown by the solid gray lines with barrier heights for Fig. LABEL:sub@fig:Nt3D:a: 120, 130, 140, 210, and 220 nK; and in Fig. LABEL:sub@fig:Nt3D:b: 230, 240, 290, 300, 340 and 350 nK. The barrier heights used in the simulations were chosen to match the results from the experiment. In the case of the tight configuration (Fig. LABEL:sub@fig:Nt3D:b), we see that the barrier heights for which the simulation closely matches the data are consistent with the experimentally measured barrier height. However, in the case of the weak configuration (Fig. LABEL:sub@fig:Nt3D:a), there is a discrepancy. This could be possibly be due to aberrations developed in the barrier beam off-axis, which results in an incorrect estimation of the barrier height. Quantum tunneling starts at around 40 ms (20 ms) in weak (tight) configuration. The non-exponential decay feature of this quantum tunneling region (shown in Fig. LABEL:sub@fig:Nt3D:a and Fig. LABEL:sub@fig:Nt3D:b) will be emphasized in Fig. 5.

We further consider the relation between decay rate and chemical potential in Fig. LABEL:sub@fig:Nt3D:c for the tight trapping configuration, with barrier heights of $V_{0}$ = 240(25) nK (green squares), $V_{0}$ = 290(30) nK (blue diamonds), and $V_{0}$ = 330(35) nK (red circles). Their corresponding 3D mean-field simulations are shown for barrier heights of 230 and 240 nK, 290 and 300 nK, 340 and 350 nK. The 3D mean-field simulations fit the tight configuration well, which did not exhibit strong classical spilling in the initial stages of the dynamics, but not the weak trapping configuration, which had well-identified distinct classical tunneling and quantum tunneling regimes, and is discussed further in Sec. IV. We observe that the decay rate can be fitted with a simple exponential function of the chemical potential of form

[TABLE]

Here $\Gamma_{\rm bg}$ is the background loss rate, and $\Gamma_{\rm bg}=0.31(0.02)\,\mathrm{Hz}$ in the experiment. These fits are included in Fig. LABEL:sub@fig:Nt3D:c as dashed black lines.

III.3 Effective 1D JWKB Description: An Experimental Case Study

In the weak configuration, there is spilling dominated dynamics for approximately the first 40 ms, quantum tunneling from about 40 ms to about 1 s, and background loss dominates thereafter. While in the tight configuration, there is spilling dominance until around $0\sim 20\,\mathrm{ms}$, quantum tunneling from 20 ms to 1.1 s, and experimental background loss thereafter. During the escape, trapped atoms tunnel from a quasi-bound state into the continuum via the weakest points, the saddle points, as the rate is largest and favorable. Figure LABEL:sub@fig:Torontoexp:b schematically shows the corresponding representative trajectories. An intriguing question is whether or not, in this rather complicated potential possibly supporting chaotic semiclassical trajectories, one can capture the basic dynamics of MQT with the well-known JWKB model in an effective 1D picture. To explore such an idea, as illustrated in Figure 1, we choose the most probable direct path for MQT, from the local minimum behind the barrier to the saddle points.

Figure 5 presents one case study with a barrier height of 190(13) nK in the weak configuration and serves to convey the details of our approach. We divide the dynamical process into three main parts: transient spilling, mean-field assisted MQT, and background loss dominated dynamics, division indicated by the vertical dashed lines. In our analysis we discard the classical spilling transients. MQT is defined as beginning when the chemical potential is equal to $V_{\rm s}$, the difference in potential energy between the saddle points and the local minimum of the potential. For the weak configuration, this is at about $t=20$ to 40 ms, beginning region A. In this region, the decay process slows down and shows a non-exponential decay feature; an exponential fit (solid black line) has $\chi^{2}=3.32$ as opposed to our model fit $\chi^{2}=1.21$. This non-exponential decay is caused by the mean-field effect, or atomic interactions, as confirmed by the effective JWKB model fit and the 3D GPE simulations. As we show in Section IV, this fit requires the interactions to match the data. As the decay process progresses, the number of atoms and the chemical potential decreases, and so does the mean-field effect. Thus, we obtain an effectively dynamical barrier height, which decreases with time. Combined, these factors produce a faster decay at the beginning of region A and a much slower decay by the end of the region. Finally, in region B, the decay process is dominated by the background loss. The blue curve is our theoretical fit and red points with error bars are the experimental data. The error region for our fit, from the uncertainty in our resulting fit parameters ($a$ and $w$ as discussed in Section. IV.1), is shown in green. The error region shown in yellow is the combined error including uncertainty in (i) fit parameters, (ii) atom number (error bars for red points), and (iii) experimental parameters; errors added in quadrature. The major contributions to the experimental error are the standard deviations of initial particle numbers ($\delta N_{0}$) at the start of the tunneling regime, and uncertainty in the peak barrier height ($\delta V_{\mathrm{s}}$), about 6% (11%) in the weak (tight) configuration. For early times, $t<100\mathrm{\mu s}$, both $\delta N_{0}$ and $\delta V_{\mathrm{s}}$ contributed about $O(10^{3}\sim 10^{4})$ to the total error envelope. After this time, the contribution from $\delta N_{0}$ decreases 1-2 orders of magnitude, while $\delta V_{\mathrm{s}}$ contribution remains about the same. All other error contributions are generally at least one order of magnitude smaller than $\delta V_{\mathrm{s}}$. We found similar results for the tight configuration, but the smaller particle numbers resulted in larger error envelopes relative to particle numbers; see Section. IV.2 for details. We also calculate the reduced chi squared of our model for all the data sets in both configurations, which confirms that our theoretical fits are well within experimental error.

In fact, this pattern in Fig. 5 occurs in all the experimental data sets, in both weak and tight configurations, as we discuss in the following section.

IV Effective 1D JWKB Model of Macroscopic Quantum Tunneling

In this section, we derive and explore in detail an effective 1D modified JWKB model, comparing quantitatively to experimental results. The term “modified” refers to inclusion of mean-field effects not normally considered in based JWKB analysis, which turn out to be key to the non-exponential decay observed in the experiment.

IV.1 Modified JWKB Model

Our initial hypothesis was that we could capture non-exponential tunneling simply through barrier shape. A triangle barrier is the simplest case with a barrier width which increases as the particles escape and the chemical potential decreases, leading to a slow-down in tunneling. The triangle also has the advantage of being analytically tractable. We found that a square barrier does not reproduce the experimental data, while a Gaussian increases the analytical difficulty without improving the accuracy of the model. We take the triangle height to be the saddle height $V_{\mathrm{s}}$ and triangle centers at $x=\pm x_{0}$, both as determined in the experiment, and single free parameter, the full width at half maximum (FWHM), $w$ (i.e., its slopes are $\pm V_{\mathrm{s}}/w$):

[TABLE]

as shown in Fig. 6. Using Eq. (3), the JWKB tunneling probability is given by

[TABLE]

the semi-classical oscillation time in the potential by

[TABLE]

and the tunneling rate through the saddle by

[TABLE]

The experimentally measured background loss rate is $\Gamma_{\mathrm{bg}}\approx 0.31(0.02)\,\mathrm{Hz}$. The number of atoms remaining behind the barrier as a function of time is given by a decay equation

[TABLE]

We further assume our system can be described by the Thomas-Fermi picture, which justifies using the analytical approximation for the chemical potential Potnis et al. (2017), from which we obtain

[TABLE]

with constant of proportionality $b=1.15(8)$ nK ($3.03(5)$ nK) as measured for the weak (tight) configuration. Note that any dependence on barrier height is negligible within experimental and model fitting parameter error.

The next logical step to adapt JWKB to this context is to replace $E$ in Eq. (4) and Eq. (5) with the appropriate single-particle energy in the presence of the effective potential and interactions from Eq. (8). However, this approach alone fails to accurately reproduce the experimental results over the entire experimental time-window: tunneling stops too soon and the dynamics of the final half of the tunneling process are not accounted for. Instead, along with replacing $E$ with $\mu$ in Eq. (8), the explicit inclusion of mean-field effects is also required to fit the data, which makes the tunneling rate $\Gamma$ in Eq. (7) time dependent, similar to results from the complex scaling method Schlagheck and Wimberger (2007). In order to take into account the time-dependent mean-field effects, we add time dependence to the height of our triangle potential, capturing a mean-field effective potential from the GPE of form $V_{\mathrm{eff}}(x)=V(x)+g|\psi(x,t)|^{2}$, as is also necessary to be self-consistent with the Thomas-Fermi approximation underlying Eq. (8); see also Moiseyev et al. (2004); Carr et al. (2004). The simplest model capturing such effects is a potential of form

[TABLE]

Here, $N_{0}$ is the initial particle number in the trap when quantum tunneling regime begins ($t_{\mathrm{QT}}$), $a>0$ is a unitless mean-field fit parameter, and $V_{\mathrm{mf}}$ is the time-dependent potential pre-factor that captures the mean-field. The effect of $a$ in $V_{\mathrm{mf}}$ is to introduce a dynamic effective barrier which decreases to the bare potential as the atoms escape. A schematic of the bare triangle barrier Eq. (3) and effective mean-field triangle Eq. (9), at the beginning of the tunneling regime $V_{\mathrm{eff}}(x,t_{\mathrm{QT}})$, are shown in Fig. 6. Note, mean-field interaction would necessarily alter the ground-state wavefunction of the system; the JWKB model used here is a minimal model which is able to reproduce the experimental findings with only two free parameters, $w$ and $a$, and may hint at many-body dynamics as will be discussed in Section V. Finally, we make the appropriate replacements of Eqs. (8) and  (9) into Eqs. (4) and (5) to give

[TABLE]

the modified JWKB tunneling probability and semi-classical oscillation time, which are subsequently substituted into Eqs. (6) and (7).

Figure 7 contains the fitted potential FWHM parameter $w$, the saddle height $V_{\rm s}$, and the initial ($t=t_{\mathrm{QT}}$) mean-field saddle height $V_{\mathrm{mf}}=V_{\rm s}(1+a)$. The increasing trend in $w$ reflects the fact that the experimental trap widens with increasing height, as shown in Table 1. The trend in $V_{\mathrm{mf}}$ overall increases with increasing peak height, leading to a deeper potential well in Fig. LABEL:sub@fig:WKBpara:a. Deeper wells are able to hold larger number of atoms, and the more atoms the stronger the mean-field effect (larger barrier). These larger barriers cause slower maximal tunneling rates as shown in Fig. LABEL:sub@fig:loggammamu:b, except for 260 nK. The case of peak height 260 nK has a smaller increase from bare saddle height $V_{\rm s}$ to $V_{\mathrm{mf}}$ than expected; $V_{\mathrm{mf}}$ for 260 nK is close to that of 221 nK even though the bare height $V_{\rm s}$ is more than 10% larger, indicating additional unaccounted for systematic error in the experimental setup for this barrier height. The initial number of atoms for $V_{0}=260\,\mathrm{nK}$ ($N_{0}\approx 583,000$) is much closer to $V_{0}=190\,\mathrm{nK}$ ($N_{0}\approx 560,000$) than $V_{0}=221\,\mathrm{nK}$ ($N_{0}\approx 630,000$). This is further supported by the results in the tight configuration runs, see Section IV.2.

In the experiment, we were able to fit a nearly exponential relation between the experimental tunneling rate $\Gamma_{\rm exp}$ and the chemical potential, Eq. (2), for weak (Fig. LABEL:sub@fig:loggammamu:a) and tight (Fig. LABEL:sub@fig:Nt3D:c) configurations. Figure LABEL:sub@fig:loggammamu:b shows both experiment and 3D mean-field simulations for the weak configuration with barrier heights of $V_{0}$ = 170(11) nK (teal squares), $V_{0}$ = 190(13) nK (purple diamonds) and $V_{0}$ = 260(18) nK (blue circles). From the modified JWKB fits, we calculate the instantaneous decay rate $\Gamma$ as a function of chemical potential $\mu$, and use Eq. (2) to fit to this relation, Fig. LABEL:sub@fig:loggammamu:b; in light of the discussion in Section IV.2, we omit tight configuration data. By examining Fig. 8, one finds that the fit captures the gross features of both the experimental data and modified JWKB curves, especially when comparing curves which correspond with the same barrier height (teal, purple, and blue curves between both subfigures. In Fig. LABEL:sub@fig:loggammamu:b, for barrier heights $V_{0}\leq 221\,\mathrm{nK}$, it can be seen that quantum tunneling is still contributing, as the tails for the fits strongly deviate from the modified JWKB curves.

Finally, we demonstrated the power of the resulting JWKB prediction for $N(t)$ in our case study in Fig. 5. Figure 9 shows $N(t)$ predictions from the theoretical fitting results from Fig. 7 for the other five sets of data in the weak configuration, with the same trends of non-exponential MQT behavior. For brevity, we do not illustrate complete error bar analyses here, but they show similar trends to the case study.

In Section. III.3 we put forth the question of whether or not an effective 1D JWKB model could reproduce the experimental findings, and in this section we have demonstrated that it is indeed possible. Two underlying assumptions were used, that we could simplify to a 1D model and that the JWKB was applicable. JWKB is applicable when the spatial derivative of the de Broglie wavelength $\lambda_{\mathrm{dB}}=\hbar/p(x)$, is small, $d\lambda_{\mathrm{dB}}/dx\ll 1$, where $p(x)$ is the semiclassical momentum. Clearly this is the case here since our model fits the data well. How can this be true for a 3D system? It turns out that the 3D potential used here can be efficiently treated by averaging over all semi-classical paths, and using the JWKB approximation for tunneling through the barriers, as established in Makri and Miller (1989). Our effective JWKB potential, Eq. (9), uses a flat potential inside the well, while a true 1D path would necessarily include the harmonic and Gaussian contributions from Eq. (1). However, this allowed the model to capture, on average, the classical oscillation frequency with the barriers caused by the complicated semi-classical trajectories in the true 3D well. The key point is that our JWKB approach is an effective model, not a first principles calculation, and the results fit the data well.

IV.2 The Tight Configuration

According to our calculations from the 3D GPE simulations for the kinetic energy term, the chemical potential needs to be modified for the tight configuration potential. This is because we estimate the chemical potential with the Thomas-Fermi approximation Potnis et al. (2017). But, for the tight configuration, due to the tighter confinement, we have a smaller particle number in the potential, which is almost one-fifth of that in the weak configuration. This makes the kinetic energy term larger in the tight configuration where the acceleration is large, which leads to a modification of chemical potential of about 6 to 21 nK for different barrier regimes in the tight configuration. Figure 10 shows the kinetic energy term in the weak and tight configuration. These are calculated by numerically solving the full 3D GPE equation using imaginary time propagation. For the weak configuration, the kinetic term is much smaller, Thomas-Fermi is an excellent approximation and indeed it fits the data well. While on the contrary, the kinetic term becomes larger and can’t be neglected in the tight configuration. All the plots which include the chemical potential in this Article are already fixed with its kinetic energy term according to Fig. 10.

The tight configuration has larger fluctuations and the traps deviate from their ideal shape in Eq. (1), an effect to which our theoretical model is sensitive; $V_{ts}$ in the tight configuration had about twice as much uncertainty as the weak configuration. Similar to the weak configuration, the largest contribution to the error envelopes was the uncertainty in barrier height and the fluctuation in the particle number, with all other errors typically 1-2 orders of magnitude smaller. Particle number error dominated early times ($t<100$ ms) with errors up to $O(10^{3}\sim 10^{4})$ before dropping 1-2 orders of magnitude, and the barrier error dominated thereafter with errors up to $O(10^{3}\sim 10^{4})$. So, although the error contributions were not significantly larger in the tight configuration, due to the fewer total particles involved in tunneling these experimental uncertainties resulted in larger error envelopes relative to the number of particles.

All fits for the tight configuration fall into two trends. Figure 11 shows two representative fits for peak heights of 290 nK and 330 nK. Both plots distinguish the error in the fit due to uncertainty in our fit parameters (green envelope), and the total uncertainty including experimental uncertainty (yellow envelope); contrast with weak configuration in Fig. 5, in which the total uncertainty is not much larger than the data error bars. Quantum tunneling is only observable for the first 600 ms, with background losses quickly dominating the dynamics thereafter, in contrast to the weak configuration where MQT is observable for about twice as long. Data sets for the tight configuration either had many points with large experimental error bars and large fluctuations like Fig. LABEL:sub@fig:WKB2ndfit:a, or had few data points with smaller error bars but still large fluctuations like Fig. LABEL:sub@fig:WKB2ndfit:b; all but $V_{0}=240\,\mathrm{nK}$ had reduced chi-squared values for $N(t)$ of $O(0.05\sim 0.10)$. Although the total uncertainty for the tight configuration did not allow for physical insight into the modified JWKB free parameters, the model was still able to capture the overall trend in Fig. 11.

V Effective Mean-Field

This section is outlined as follows. We first present the literature which motivates using a renormalized mean-field parameter for depleted and fragmented BECs. Next, we present the many-body Hamiltonian and a corresponding mean-field equation. Following that, we describe the method by which we create a meta-stable state to study quantum tunneling. We show the failure in a straightforward mean-field application, and how a modified interaction parameter captures the low-order many-body dynamics from TEBD. Finally, we draw conclusions from the MQT experiment, and suggest future studies to distinguish between true mean-field and many-body effects in MQT experiments.

Experimentally, BECs have overall been very well understood with mean-field theory in the form of the GPE. However, the presence of many-body processes can induce fluctuations, fragmentation, and depletion, thereby rendering mean-field models inaccurate or even ineffective. Many techniques and approaches are used to move beyond the GPE Proukakis (2008); of particular interest here are those using an effective interaction parameter Stoof et al. (1996); Stoof (1997); Proukakis et al. (1998); Shi and Griffin (1998); Morgan (2004, 2005). We will demonstrate how, similar to multi-component optical systems being well modeled by an effective scalar GPE-like equation derived from a multi-mode or large vector NLS, fragmented and depleted BECs described by a many-body Hamiltonian can, for the purposes of MQT, be well described by a mean-field model with an effective or renormalized interaction parameter.

Advances in nonlinear optics and strong connections to the NLS motivate the effective interaction parameter used in our JWKB model. A deep theoretical understanding of the propagation and self-focusing of partially incoherent beams in nonlinear media, which can lead to spatial incoherent solitons, has been developed through several equivalent methods Christodoulides et al. (2001): an infinite set of coupled nonlinear Schrödinger equations (coherent density approach) Christodoulides et al. (1997a, b); Coskun et al. (1998); Eugenieva et al. (2000), propagation equation for mutual coherence function Nayyar (1997); Hasegawa (1975); Snyder and Mitchell (1998); Shkunov and Anderson (1998), and self-consistent multimode theory Mitchell et al. (1997); Christodoulides et al. (1998); Akhmediev et al. (1998); Królikowski et al. (1999); Sukhorukov and Akhmediev (1999). Similarities in the propagation equation for the mutual coherence method to the NLS allow for analytical techniques to be extended for partially incoherent regimes; e.g., derivation of an analytical expression for the collapse threshold of spatially partially coherent beans in inertial bulk Kerr media Bang et al. (1999). It has also been shown that BECs at finite temperature have analogous behavior to incoherent light in nonlinear media Buljan et al. (2005), further suggesting that analogies to nonlinear optics can offer useful insight.

For the many-body dynamics, consider $N$ bosons at zero temperature in the canonical ensemble in a quasi-bound state, one suitable for quantum tunneling as laid out in this Article. One appropriate model to study many-body dynamics for this system is the Bose-Hubbard Hamiltonian (BHH), which can be invoked using an optical lattice of $L$ sites with deep enough sites for a tight binding and single band approximation, or alternately taking it as a discretization scheme in an appropriate limit:

[TABLE]

In Eq. (12), $U$ determines the on-site two-particle interactions and $J$ is the energy of hopping. An external box trapping potential, such as in schematic plot Fig. LABEL:sub@fig:schematic:a, is given by $V_{i}^{\mathrm{ext}}$ with height $h$. The field operator $\hat{b}_{i}^{\dagger}$ ($\hat{b}_{i}$) creates (annihilates) a boson at the $i\mathrm{th}$ site, satisfying the usual commutation relation $[\hat{b}_{i},\hat{b}_{j}^{\dagger}]=\delta_{ij}$, and $\hat{n}_{i}\equiv\hat{b}_{i}^{\dagger}\hat{b}_{i}$. We will work in hopping units: energies are scaled to $J$ and time $t$ to $\hbar/J$. To simulate the many-body dynamics of the BHH, we use TEBD, a matrix-product-state method which is able to efficiently simulate one dimensional many-body quantum systems, and allows access to a wide variety of many-body quantities like fluctuation and entanglement MPS ; Wall and Carr (2012); Schollwock (2005).

To describe the system from a mean-field perspective, the discrete nonlinear Schrödinger equation (DNLS) may either be obtained via discretization of the NLS or from a mean-field approximation of the BHH Mishmash et al. (2009):

[TABLE]

In Eq. (13), the condensate order parameter, $\psi_{i}$, is normalized to the number of atoms, $N=\sum^{L}_{i=1}|\psi_{i}|^{2}$, and $g\equiv U$ from the BHH. Note however, $g$ in the DNLS is related, but not equivalent, to the interaction parameter in the GPE. To be exact, $g=g^{(1)}\int dx|W^{(0)}(x)|^{4}$, where $g^{(1)}$ is the quasi-1D interaction strength, proportional to scattering length, and $W^{(0)}(x)$ is the lowest order Wannier state; see Mishmash (2012) for details. Mean-field simulations are performed using a fourth-order Runge-Kutta adaptation of Eq. (13). The BHH approaches the DNLS equation in the mean-field limit $N\to\infty$, $U\to 0$, $NU=\mathrm{const.}$, which we take in units of $J$, $NU/J=\mathrm{const.}$ Both the BHH and the DNLS are single band models, valid when the many-body wavefunction covers many sites and has variations larger than the lattice constant.

We use box boundary conditions and initially set the barrier, $V^{\mathrm{ext}}_{i}$, to zero over the first 15 lattice sites, and $h=0.15$ over all the rest; the lattice size $L$ is large enough that the particles don’t reach the end of the grid within simulation time, mimicking escape into open space. After propagation in imaginary time, the wavefunction is mostly trapped in the first 15 sites. We then reduce the external barrier to be $h=0.10$ from sites 16 to 19, thus rendering the state meta-stable, and evolve in real time. We reduce the height to induce larger fragmentation and depletion of the many-body wavefunction, but show how the renormalized mean-field still captures the key features of tunneling dynamics.

In Fig. LABEL:sub@fig:MF_Plot:a we plot the normalized number of trapped atoms as a function of time for the many-body simulations with $NU/J=1.00$, mean-field with $g=NU/J=1.00$, and effective mean-field with two values of $g_{\mathrm{eff}}$. The values of $g_{\mathrm{eff}}$ were found by sweeping over $g$ values until appropriate upper and lower bounds on the many-body dynamics were found. In Fig. LABEL:sub@fig:MF_Plot:b, we plot the two largest $\{\lambda_{\ell}\}$, the eigenvalues of the single particle density matrix $\langle\hat{b}_{i}^{\dagger}\hat{b}_{j}\rangle$, and $\lambda_{\geq 3}$, the sum of all but the two largest eigenvalues. We calculate eigenvalues for the single particle density matrix over the entire system ($\lambda_{l}^{\mathrm{sys}}$, solid curves in Fig. LABEL:sub@fig:MF_Plot:b), as well as focusing in on the single-particle density matrix in the trapping potential ($\lambda_{l}^{\mathrm{trap}}$, dashed lines in Fig. LABEL:sub@fig:MF_Plot:b). Information about the whole system $\lambda_{i}^{\mathrm{sys}}$ indicates a large degree of fragmentation, with up to 30% occupation of the second mode, and depletion, with more than 10% occupation in all but first 2 modes. Information regarding only the remaining atoms trapped in the well ($\lambda_{i}^{\mathrm{trap}}$) shows a large fraction of the BEC has escaped. Even though the many-body wavefunction would be considered over $30\%$ depleted by $t_{\mathrm{esc}}\approx 186$, the time at which the number of trapped atoms is $1/e$ of the initial value, our mean-field plots accurately bound the trend in many-body tunneling with a simple renormalization of the interaction parameter $g$ by about 25%. The effective mean-field simulation is able to qualitatively capture the overall trend in the number of trapped particles for tunneling, with less than $O(10^{-3})$ relative error in the trapped number of atoms, an error which would be indistinguishable within the error bars of many tunneling experiments including our own. While these results show that effective mean-field models can accurately reproduce low-order observables like trapped-atom density for depleted many-body wave-functions, we do not include loss, finite temperature, and other open-system effects in our many-body simulations.

In Section. III.2, we show how a full 3D mean-field treatment of the experiment captures the gross features of the dynamics, and in Section. IV we show how a JWKB model with an effective mean-field like parameter fits experiments runs well. These results are corroborated in this section, as we have shown how a renormalized mean-field interaction can adequately capture quantum tunneling of fragmented and depleted condensates, two effects that are likely present in an open quantum system such as our own, and compensated for with our JWKB parameter $a$. The applicability of mean-field seems to be larger than expected, and care must be taken to determine whether a given set of experimental data is purely mean-field. Future experiments with sufficient resolution may be able to distinguish between many-body and effective mean-field dynamics by reducing number fluctuations in the initial state to resolve effective interaction strengths, for example by post-selecting on atomic number in measurements or using atom interferometry Langen et al. (2015). Furthermore, coherence experiments on the escaped particles may further illuminate this distinction, as suggested by macroscopic occupation of more than one mode over the entire lattice in Fig. LABEL:sub@fig:MF_Plot:b.

VI Conclusions

We first provided a brief discussion about different MQT regimes according to four factors: statistical properties and quasiparticles, the role of interactions, type of trapping potential, and dimension of the system. In the discussion about atomic interactions, we suggested weak interactions can be described by mean-field-like theories, with fragmentation, entanglement and higher order quantum fluctuations increasing for stronger interactions and presenting experimentally untested regimes of MQT. We suggested different quasi-particle descriptions may result in distinct tunneling regimes. We outlined how the shape of potential wells define the modes which particles can occupy, before and after tunneling, leading to a strong distinction between e.g. Josephson physics and quasibound many-body dynamics, or quantum escape. We emphasize that the barrier is generally deformed by interactions and one must think at least in terms of an effective potential, as evidenced also by our experiment. Although tunneling is primarily a 1D effect, higher dimensions can affect tunneling by creating chaos in semiclassical paths behind the barrier, for instance. This brief survey set the tone for our own work and suggested future experiments and theoretical development needed, for instance, in the study of MQT of the unitary Fermi gas, perhaps in its holographic dual.

Then, we described the tunneling experiment and the resulting non-exponential decay of the trapped atoms. A mean-field description was validated by overall agreement between the experimental results and a 3D mean-field simulation, in decay curves and decay rates. The 3D mean-field simulation also calculated the kinetic term separately, which led to a correction of the kinetic term in the tight trap configuration. Our theoretical mean-field model reproduced the experimentally observed non-exponential decay, which was previously indicated to be the result of the participation of atomic interactions in the tunneling process. We then proceeded to explore the usefulness of a much simpler effective 1D model utilizing a modified version of the semiclassical approximation, or JWKB. A case study of the experiment with barrier height 190nK in the weak configuration demonstrated the effectiveness of our 1D model. We further divided the tunneling process into three sub-regions: an initial decay region corresponding to classical spilling over the two saddles appearing in our potential; a non-exponential quantum tunneling regime during which interactions heavily affected the tunneling rate; and a background-loss dominated region. This pattern appeared in all the experimental runs.

Subsequently, we described the modified JWKB model in detail. The transition point into JWKB is when the relative height of the saddle points equals the chemical potential. In the modified JWKB model we introduced an effective mean-field term which modified the barrier and produced non-exponential decay consistent with the experiment. Thus, this non-exponential decay can be described as generated by atomic interactions, which led to an effective dynamic barrier. Both 3D mean-field and 1D JWKB models confirmed an exponential relation between decay rate and chemical potential observed empirically in the experiments. JWKB used only two fitting parameters, an effective mean-field saddle height $V_{\mathrm{mf}}$ resulting from unitless mean-field parameter $a$, and saddle width $w$. The parameter $w$ grew wider with increasing barrier height, which followed the experimental trend. Increasing $V_{\mathrm{mf}}$ reflected the fact that higher peaks in the experiment have smaller maximum tunneling rates.

Finally, we showed how a renormalized mean-field theory is capable of capturing many-body quantum effects in low-order observables, by comparing discrete nonlinear Schrödinger and TEBD simulations, in analogy to optics contexts in which many modes create an effective scalar nonlinear-Schrödinger or Gross-Pitaevskii type description. Thus, experiments with large number fluctuations or error bars can have difficulty in discerning between mean-field and renormalized mean-field due actually to many-body effects. Future experiments which more precisely resolve atomic number and interaction strengths may be able to distinguish between the bare mean-field theories generally assumed to describe BEC dynamics and the effective or renormalized theories we have here suggested, in terms of a concrete observable, the number of atoms remaining in a many-body quasi-bound state as a function of time.

Acknowledgements.

The authors thank Rockson Chang, Gavriil Shchedrin, and David Spierings for helpful discussions and A. Stummer for technical support. Computations were performed on the gpc supercomputer at the SciNet HPC Consortium and Colorado School of Mines high performance computing resources as part of the Golden Energy Computing Organization. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund - Research Excellence, and the University of Toronto. This material is based in part upon work supported by NSERC, CIFAR, the US National Science Foundation under grant numbers PHY-1306638, PHY-1207881, and PHY-1520915 (LDC) , and the US Air Force Office of Scientific Research grant number FA9550-14-1-0287 (LDC), and China Scholarship Council (XXZ) and Northrop Grumman Aerospace Systems.