# How an interacting many-body system tunnels through a potential barrier to open space

See allHide authors and affiliations

Edited by Eric Johnson Heller, Harvard University, Cambridge, MA, and approved July 3, 2012 (received for review January 24, 2012)

## Abstract

The tunneling process in a many-body system is a phenomenon which lies at the very heart of quantum mechanics. It appears in nature in the form of α-decay, fusion and fission in nuclear physics, and photoassociation and photodissociation in biology and chemistry. A detailed theoretical description of the decay process in these systems is a very cumbersome problem, either because of very complicated or even unknown interparticle interactions or due to a large number of constituent particles. In this work, we theoretically study the phenomenon of quantum many-body tunneling in a transparent and controllable physical system, an ultracold atomic gas. We analyze a full, numerically exact many-body solution of the Schrödinger equation of a one-dimensional system with repulsive interactions tunneling to open space. We show how the emitted particles dissociate or fragment from the trapped and coherent source of bosons: The overall many-particle decay process is a quantum interference of single-particle tunneling processes emerging from sources with different particle numbers taking place simultaneously. The close relation to atom lasers and ionization processes allows us to unveil the great relevance of many-body correlations between the emitted and trapped fractions of the wave function in the respective processes.

The tunneling process has been a matter of discussion (1⇓–3) since the advent of quantum mechanics. In principle, it takes place in all systems whose potential exhibits classically forbidden but energetically allowed regions. See, for example, the overview in ref. 4 and Fig. 1. When the potential is unbound in one direction, the quantum nature of the systems allows them to overcome potential barriers for which they classically would not have sufficient energy and, as a result, a fraction of the many-particle system is emitted to open space. For example, in fusion, fission, photoassociation, and photodissociation processes, the energetics or life times are of primary interest (5⇓⇓⇓–9). The physical analysis was made under the assumption that the correlation between decay products (i.e., between the remaining and emitted fractions of particles) can be neglected. However, it has to be stressed first that, at any finite decay time, the remaining and emitted particles still constitute one total many-body wave function and, therefore, can be correlated. Second, in contrast to the tunneling of an isolated single particle into open space, which has been amply studied and understood (4), nearly nothing is known about the tunneling of a many-body system. In the present study, we demonstrate that the fundamental many-body aspects of quantum tunneling can be studied by monitoring the correlations and coherence in ultracold atomic gases. For this purpose, the initial state of an ultracold atomic gas of bosons is prepared coherently in a parabolic trapping potential which is subsequently transformed to an open shape allowing for tunneling (Fig. 1, *Upper*). In this tunneling system the correlation between the remaining and emitted particles can be monitored by measuring deviations from the initial coherence of the wave function. The close relation to atom lasers and ionization processes allows us to predict coherence properties of atom lasers and to propose the study of ionization processes with tunneling ultracold bosons.

In the past decade, Bose–Einstein condensates (BECs) (10, 11) have become a toolbox to study theoretical predictions and phenomena experimentally with a high level of precision and control (12). BECs are experimentally tunable, the interparticle interactions (13), trap potentials (14), number of particles (15), their statistics (16) and even dimensionality (17⇓–19) are under experimental control. From the theoretical point of view, the evolution of an ultracold atomic cloud is governed by the time-dependent many-particle Schrödinger equation (TDSE) (20) with a known Hamiltonian (*Methods*). To study the decay scenario in Fig. 1, we solve the TDSE numerically for long propagation times. For this purpose, we use the multiconfigurational time-dependent Hartree method for bosons (MCTDHB) (for details and literature see *Methods*). Our protocol to study the tunneling process of an initially parabolically trapped system into open space is schematically depicted in Fig. 1, *Upper*. We restrict our study here to the one-dimensional case, which can be achieved experimentally by adjusting the transverse confinement appropriately. As a first step, we consider *N* = 2,4, or 101 weakly repulsive ^{87}Rb atoms in the ground state of a parabolic trap (see *Methods* for a detailed description of the considered experimental parameters).

The initial one-particle density and trap profile are depicted in Fig. 1, *Upper*. Next, the potential is abruptly switched to the open form *V*(*x*,*t*≥0) indicated in this figure, allowing the many-boson system to escape the trap by tunneling through the barrier formed. Eventually, all the bosons escape by tunneling through the barrier, because the potential supports no bound states.

The initial system, i.e., the source of the emitted bosons, is an almost totally coherent state (21). The final state decays entirely to open space to the right of the barrier, where the bosons populate many many-body states, related to Lieb–Liniger states (22⇓–24), which are generally not coherent. It is instructive to ask the following guiding questions: What happens in between these two extremes of complete coherence and complete incoherence? And how does the correlation (coherence) between the emitted particles and the source evolve? By finding the answers to these questions we will gain a deeper theoretical understanding of many-body tunneling, which is of high relevance for future technologies and applied sciences. In particular, this knowledge will allow us to determine whether the studied ultracold atomic clouds qualify as candidates for atomic lasers (11, 25⇓–27) or as a toolbox for the study of ionization or decay processes (5⇓⇓–8).

Because the exact many-body wavefunctions are available at any time in our numerical treatment, we can quantify and monitor the evolution of the coherence and correlations of the whole system as well as between the constituting parts of the evolving wave packets. In the further analysis we use the one-particle density in real ρ(*x*,*t*) and momentum ρ(*k*,*t*) spaces, their natural occupations and correlation functions *g*^{(1)}(*k*^{′}|*k*; *t*) (28⇓–30). We defer the details on these quantities to *Methods*. To study the correlation between the source and emitted bosons we decompose the one-dimensional space into the internal IN and external OUT regions with respect to the top of the barrier, as illustrated by the red line and arrows in Fig. 1, *Lower*. This decomposition of the one-dimensional Hilbert space into subspaces allows us to quantify the tunneling process by measuring the amount of particles remaining in the internal region in real space as a function of time. In Fig. 2 we depict the corresponding quantities for *N* = 2 and *N* = 101 by the green dotted curve. A first main observation is that the tunneling of bosonic systems to open space resembles an exponential decay process.

The key features of the dynamics of quantum mechanical systems manifest themselves very often in characteristic momenta. Therefore, it is worthwhile to compute and compare evolutions of the momentum distributions ρ(*k*,*t*) of our interacting bosonic systems. Fig. 3 depicts ρ(*k*,*t*) for *N* = 2,4,101 bosons. At *t* = 0 all the initial real space densities have Gaussian-shaped profiles resting in the internal region (Fig. 1, *Upper*). Therefore, their distributions in momentum space are also Gaussian-shaped and centered around *k* = 0. With time the bosons start to tunnel out of the trap. This process manifests itself in the appearance of a pronounced peak structure on top of the Gaussian-shaped background, see Fig. 3, *Black Framed Upper*. The peak structure is very narrow (similar to a laser or an ionization process), the bosons seem to be emitted with a very well-defined momentum. For longer propagation times a larger fraction of bosons is emitted and more intensity is transferred to the peak structure from the Gaussian background. Thus, we can relate the growing peak structures in the momentum distributions to the emitted bosons and the Gaussian background to the bosons in the source. We decompose each momentum distribution into a Gaussian background and a peak structure to check the above relation (*Methods*). The integrals over the Gaussian momentum background, , are depicted in Fig. 2 as a function of time. The close similarity of the , characterizing the amount of particles remaining in the internal region in real space, and confirms our association of the Gaussian-shaped background with source bosons and the peak structure in the momentum distribution with emitted ones (Fig. 3, *Black Framed Upper*).

Now we are equipped to look into the mechanism of the many-body physics in the tunneling process with a simplistic model. When the first boson is emitted to open space one can estimate its total energy as an energy difference, *E*^{N} - *E*^{N-1}, of source systems made of *N* and *N* - 1 particles. Assuming negligible coupling between the trapped and emitted bosons the energy of the outgoing boson reads: The quantity μ_{1} is simply the chemical potential of the parabolically trapped source system. For noninteracting particles its value would be independent of the number of bosons inside the well. Because the boson tunnels to open space where the trapping potential is zero, it can be considered for the time being as a free particle with all its available energy, μ_{1}, converted to kinetic energy . We have to expect the first emitted boson to have the momentum . The value of the estimated momentum agrees excellently with the position of the peak in the computed exact momentum distributions, see the arrows marked *k*_{1} in Fig. 3, *Lower* (*i–iv*)). This agreement allows us to interpret the peak structures in ρ(*k*,*t*) as the momenta of the emitted bosons. As a striking feature, also other peaks with smaller *k* appear in these spectra at later tunneling times (see Fig. 3, *Lower* (*i*, *ii*, and *iv*)). We can use an analogous argumentation and associate the second peak with the emission of the second boson and its chemical potential, μ_{2}. Its kinetic energy and the corresponding momentum, *k*_{2}, can be estimated as the energy difference μ_{2} = *E*^{N-1} - *E*^{N-2} between the source subsystems made of *N* - 1 and *N* - 2 bosons, under the assumption of zero interaction between the emitted particles and the source. The correctness of the applied logic can be verified from Fig. 3, *Lower* (*i* and *ii*), where for the *N* = 2 (*N* = 4) particles the positions of the estimated momenta *k*_{2} (*k*_{3}, *k*_{4}) of the second (third, fourth) emitted boson fit well with the position of the second (third, fourth) peak in the computed spectra. The momenta and details on the calculation are collected in the *SI Text*. The momentum spectrum for *N* = 101 bosons shows a similar behavior—the multipeak structures gradually develop with time starting from a single-peak to two-peaks and so on, see Fig. 3, *Lower* (*iv*). However, from this figure we see that the positions of the peaks’ maxima, and with them the momenta of the emitted bosons, change with time. On the one hand we see that the considered tunneling bosonic systems can not be utilized as an atomic laser: The initially coherent bosonic source emits particles with different, weakly time-dependent momenta. In optics such a source would be called polychromatic. On the other hand we can associate the peaks with different channels of an ionization process and their time-dependency with the channels’ coupling. Thus, we conclude that it is possible to model and investigate ionization processes with tunneling ultracold bosonic systems. We mention that such processes can occur in the multielectron photoionization dynamics of molecules in a laser field. Namely, several electrons can be ionized by tunneling through a barrier, see the recent experiment in ref. 31.

Let us now investigate the coherence of the tunneling process itself. In the above analysis of the momentum spectra we relied on the exact numerical solutions of the TDSE for *N* = 2,4,101 bosons. We remind the reader that in the context of ultracold atoms the Gross–Pitaevskii (GP) theory is a popular and widely used mean-field approximation describing systems under the assumption that they stay fully coherent for all times. In our case the GP approximation assumes that the ultracold atomic cloud coherently emits the bosons to open space and keeps the source and emitted bosons coherent all the time. To learn about the coherence properties of the ongoing dynamics it is instructive to compare exact many-body solutions of the TDSE with the idealized GP results, see Fig. 3, *Lower* (*iii*). The strengths of the interboson repulsion have deliberately been chosen such that the GP gives identical dynamics for all *N* studied. It is clearly seen that for short initial propagation times the dynamics are indeed coherent. The respective momentum spectra obtained at the many-body and GP levels are very similar, see Fig. 3, *Lower* (*i–iv*) for ρ(*k*,*t*_{1} = 100). At longer propagation times (*t* > *t*_{1}), however, the spectra become considerably different. This difference means that with time, the process of emission of bosons becomes less coherent.

Next, to quantify the coherence and correlations between the source and emitted bosons we compute and plot the momentum correlation functions |*g*^{(1)}(*k*^{′},*k*|*t*)|^{2} in Fig. 4 (*Left*) for *N* = 101 (for *N* = 2,4 they look almost the same). Let us stress here that the proper correlation properties cannot be accounted for by approximate methods. For example the GP solution of the problem gives |*g*^{(1)}|^{2} = 1, i.e., full coherence for all times. For the exact solution we also obtain that at *t* = 0 the system is fully coherent, and thus |*g*^{(1)}(*k*^{′}|*k*; *t* = 0)|^{2} = 1. Hence, Fig. 4, *Upper Left* is also a plot for the GP time evolution. However, during the tunneling process the many-body evolution of the system becomes incoherent, i.e., |*g*^{(1)}|^{2} → 0. The coherence is lost only in the momentum-space domain where the momentum distributions are peaked, the *k*-region associated with the emitted bosons (Fig. 4, *Left*). In the remainder of *k*-space the wave function stays coherent for all times. We conclude that the trapped bosons within the source remain coherent. The emitted bosons become incoherent with their source and among each other. Therefore, the coherence between the source and the emitted bosons is lost. A complementary argumentation with the normalized real-space correlation functions is deferred to *SI Text*.

In the spirit of the seminal work of Penrose and Onsager on reduced density matrices (29), we tackle the following question: how strong is the loss of coherence in many-body systems? The natural occupation numbers, , obtained by diagonalizing the reduced one-body density (*Methods*) define how much the system can be described by single, two, or more quantum mechanical one-particle states. The system is condensed and coherent when only one natural occupation is macroscopic and it is fragmented when several of the are macroscopic. In Fig. 4 (*Right*) we plot the evolution of the natural occupation numbers for the studied systems as a function of propagation time. The initial system is totally coherent—all the bosons reside in one natural orbital. However, when some fraction of the bosons is emitted, a second natural orbital gradually becomes occupied. The decaying systems lose their coherence and become twofold fragmented. For longer propagation times more natural orbitals start to be populated, indicating that the decaying systems become even more fragmented, i.e., less coherent. In the few-boson cases, *N* = 2 and *N* = 4, one can observe several stages of the development of fragmentation (32, 33)—beginning with a single condensate and evolving toward the limiting form of *N* entangled fragments in the end which means a full fermionization of the emitted particles. In this fermionization-like case each particle will propagate with its own momentum. For larger *N* the number of fragments also increases with tunneling time. The details of the evolution depend on the strength of the interparticle interaction and number of particles. By resolving the peak structure in the momentum spectra of the tunneling systems at different times we can directly detect and quantify the evolution of coherence, correlations, and fragmentation.

Finally, we tackle the intricate question whether the bosons are emitted one by one or several at a time? By comparing the momentum spectra ρ(*k*,*t*) depicted in Fig. 3, *Lower* (*i*, *ii*, and *iv*) at different times it becomes evident that the respective peaks appear in the spectra sequentially with time, starting from the most energetical one. If multiboson (two or more boson) tunneling processes would participate in the dynamics, they would give spectral features with higher momenta which are not observed in the computed spectra depicted in Fig. 3. A detailed discussion and a model are given in *SI Text*. This model suggests that the bosons tunnel out one by one. However, the fact that the peaks’ heights and positions evolve with time indicates that the individual tunneling processes interfere, i.e., they are not independent. The origin behind this interference is the interaction between the bosons. We conclude that the overall decay by tunneling process is of a many-body nature and is formed by the interference of different single-particle tunneling processes taking place simultaneously.

We arrive at the following physical picture of the tunneling to open space of an interacting, initially coherent bosonic cloud. The emission from the bosonic source is a continuous, polychromatic many-body process accompanied by a loss of coherence, i.e., fragmentation. The dynamics can be considered as a superposition of individual single-particle tunneling processes of source systems with different particle numbers. On the one hand, ultracold weakly interacting bosonic clouds tunneling to open space can serve as an atomic laser, i.e., emit bosons coherently, but only for a short time. For longer tunneling times the emitted particles become incoherent. They lose their coherence with the source and among each other. On the other hand, we have shown the usage of tunneling ultracold atoms to study the dynamics of ionization processes. Each peak in the momentum spectrum is associated with the single particle decay of a bosonic source made of *N*, *N* - 1, *N* - 2, etc. particles—in close analogy to sequential single ionization processes. These *N* discrete momenta comprise a total spectrum—in close analogy to total ionization spectra.

We have focused here on the many-body tunneling process of an initially coherent bosonic cloud. Nevertheless, it is interesting to inquire whether the many-body tunneling mechanism would change when the interaction between the bosons is increased, e.g., by employing a Feshbach resonance. To this end, we have repeated our studies for *N* = 2,4,101 bosons with sevenfold stronger interactions for which the initial state is still close to condensed, as well as for *N* = 2,4 200-fold stronger interacting bosons, for which the initial state is already fermionized. The results are shown in *SI Text*. The tunneling mechanism is not altered by the stronger interactions and our model predicts the positions of the momentum peaks of the escaping bosons well.

As an experimental protocol for a straightforward detection of the kinetic energy of the emitted particles one can use the presently available single-atom detection techniques on atom chips (34) or the idea of mass spectrometry (see discussion in *SI Text*). Summarizing, in many-particle systems decaying by tunneling to open space the correlation dynamics between the source and emitted parts lead to clearly observable spectral features which are of great physical relevance.

## Methods

### Hamiltonian and Units.

The one-dimensional *N*-boson Hamiltonian reads as follows: Here λ_{0} > 0 is the repulsive interparticle interaction strength proportional to the *s*-wave scattering length *a*_{s} of the bosons and *x*_{i} is the coordinate of the *i*-th boson. Throughout this work λ = λ_{0}(*N* - 1) = 0.3 for all considered *N* is used (further results for *N* = 2,4,101 bosons with sevenfold stronger interaction, λ = 2.1, and for *N* = 2,4 bosons with 200-fold stronger interaction, λ = 60, are discussed at the end of the main text and *SI Text*). For convenience we work with the dimensionless quantities defined by dividing the dimensional Hamiltonian by , where is Planck’s constant, *m* is the mass of a boson, and *L* is a chosen length scale. At *t* < 0 trap is parabolic , the analytic form of *V*(*x*,*t*) after the opening is given in ref. 35.

In this work we consider ^{87}Rb atoms for which *m* = 1.44316·10^{-25} kg and *a*_{s} = 90.4*a*_{0} without tuning by a Feshbach resonance, where *a*_{0} = 0.0529·10^{-9} m is Bohr’s radius. We emulate a quasi one-dimensional cigar-shaped trap in which the transverse confinement is *w*_{⊥} = 2,291.25 Hz, which is amenable to current experimental setups. Following ref. 36, the transverse confinement renormalizes the interaction strength. Combining all the above, the length scale is given by , and the timescale by . The relation between the (dimensionless) interaction parameter λ_{0} and the (dimension-full) scattering length *a*_{s} is hence given by .

### The Multiconfigurational Time-Dependent Hartree for Bosons Method.

The time-dependent many-boson wave function Ψ(*t*) solving the many-boson Schrödinger equation is obtained by the MCTDHB, see refs. 37 and 38. Applications include unique intriguing many-boson physics such as the death of attractive soliton trains (39), formation of fragmented many-body states (40), and numerically exact double-well dynamics (41⇓–43). Recent optimizations of the MCTDHB, see, e.g., refs. 44 and 45, allow now for the application of the algorithm to open systems with very large grids (here 2^{16} = 65,536 basis functions), a particle number of up to *N* = 101, and an arbitrary number of natural orbitals (here up to 14). We would like to stress that even nowadays such kind of time-dependent computations are very challenging.

The mean-field wave function is obtained by solving the time-dependent Gross–Pitaevskii equation, which is contained as a special single-orbital case in the MCTDHB equations of motion, see refs. 37 and 38. To ensure that the tunneling wave packets do not reach the box borders for all presented propagation times the simulations were done in a box [-5; 7465]. In the dimensional units we thus solve a quantum mechanical problem numerically exactly in a spatial domain extending over 8.29 mm.

### Many-Body Analysis of the Wave Function.

With the many-boson wave function Ψ(*t*) at hand the various quantities of interest are computed and utilized to analyze the evolution in time of the Bose system. The reduced one-body density matrix of the system is given by , where is the usual bosonic field operator creating a boson at position *x*. Diagonalizing ρ^{(1)}(*x*|*x*^{′}; *t*) one gets the natural orbitals (eigenfunctions), , and natural occupation numbers (eigenvalues) from the expression . The latter determine the extent to which the system is condensed (one macroscopic eigenvalue) or fragmented (two or more macroscopic eigenvalues) (32, 33, 46). The diagonal part of the reduced one-body density matrix ρ(*x*,*t*) ≡ ρ^{(1)}(*x*|*x*^{′}; *t*) is the system’s density. The first-order correlation function in coordinate space quantifies the degree of spatial coherence of the interacting system (28, 30). The respective quantities in momentum space, such as the momentum distribution ρ(*k*,*t*) and the first-order correlation function in momentum space , are derived from ρ^{(1)}(*x*|*x*^{′}; *t*) via an application of a Fourier transform on its eigenfunctions.

In real space the density-related nonescape probability is given by (see ref. 35 for the non-Hermitian results). The momentum-density related nonescape probability is obtained by least-squares fitting a Gaussian function ρ^{Gauss}(*k*,*t*) = *Ae*^{-(Bx)2} to ρ(*k*,*t*) in the *k*-space domain [-∞,0]. *A* and *B* are the fit parameters. We then define the momentum-density-related nonescape probability as

## Acknowledgments

We thank Shachar Klaiman and Julian Grond for a careful reading of the manuscript and comments as well as Lincoln Carr for discussions. Computation time on the bwGRiD and the Cray XE6 cluster Hermit at the High Performance Computing Center Stuttgart (HLRS), and financial support by the Heidelberg Graduate School of Mathematical and Computational Methods for the Sciences (HGS MathComp) and the Deutsche Forschungsgemeinschaft (DFG) also within the framework of the Enable fund of the excellence initiative at Heidelberg university are greatly acknowledged.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: axel.lode{at}pci.uni-heidelberg.de.

Author contributions: A.U.J.L., A.I.S., K.S., O.E.A., and L.S.C. designed research, performed research, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1201345109/-/DCSupplemental.

## References

- ↵
- ↵
- ↵
- Kramers HA

- ↵
- Razavy M

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Ullrich J,
- Shevelko VP

- ↵
- Pitaevskii LP,
- Stringari S

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Glauber RJ

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Streltsov AI,
- Sakmann K,
- Lode AUJ,
- Alon OE,
- Cederbaum LS

- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics