Avalanches of Bose-Einstein Condensates in Leaking Optical Lattices

One of the most fascinating experimental achievements of the last decade was the realization of Bose-Einstein Condensation (BEC) of ultra-cold atoms in optical lattices (OL's). The extraordinary level of control over these structures allows us to investigate complex solid state phenomena and the emerging field of ``atomtronics'' promises a new generation of nanoscale devices. It is therefore of fundamental and technological importance to understand their dynamical properties. Here we study the outgoing atomic flux of BECs loaded in an one-dimensional OL with leaking edges, using a mean field description provided by the Discrete Non-Linear Schrodinger Equation (DNLSE). We demonstrate that the atom population inside the OL decays in avalanches of size $J$. For intermediate values of the interatomic interaction strength their distribution ${\cal P}(J)$ follows a power law i.e. ${\cal P}(J)\sim1/J^{\alpha}$ characterizing systems at phase transition. This scale free behaviour of ${\cal P}(J)$ reflects the complexity and the hierarchical structure of the underlying classical mixed phase space. Our results are relevant in a variety of contexts (whenever DNLSE is adequate), most prominently the light emmitance from coupled non-linear optics waveguides.


Introduction
One of the most fascinating experimental achievements of the last decade was the realization of Bose-Einstein condensation (BEC) of ultra-cold atoms in optical lattices (OLs) [1]- [3]. The extraordinary level of control over these structures allows us to investigate complex solid state phenomena [2]- [8], and the emerging field of 'atomtronics' promises a new generation of nanoscale devices. It is therefore both of fundamental and technological importance to understand the dynamics and transport properties of BECs in OLs.
A fundamental source of physical information is time resolved decay measurements of quantum mechanical systems which are coupled to a continuum via leads or absorbing boundaries. While radioactive decay is a prominent paradigm, more recent examples include atoms in optically generated lattices and billiards [9,10], the ionization of molecular Rydberg states [11], photoluminescence spectroscopy of excitation relaxation in semiconductor quantum dots and wires [12], and pulsed electromagnetic propagation studies [13].
In this paper, using a mean field description provided by the discrete nonlinear Schrödinger equation (DNLSE), we study the statistical properties of the outgoing atomic flux of BECs loaded in OLs with leaking edges (for experimental realizations see [14,15]). We find three types of dynamical behavior: for very weak interatomic interaction strength , the population decay is smooth. When is increased the atom population N (τ ) inside the OL decays in bursts of size J ≡ δ N , which are triggered by collisions between discrete breathers (DBs) and other excitations supported by the OL. For large values of the distribution P(J ) is exponential. But for a whole range of intermediate 's, the bursts resemble scale-free avalanches and P(J ) follows a power law, i.e. P(J ) ∼ 1/J α , characterizing systems at a phase transition. This scalefree behavior of P(J ) reflects the complexity and the hierarchical structure of the underlying classical mixed phase space. We mark that avalanches of BECs have been discussed in different contexts in [16], while the possibility to control the phase of BEC via localized dissipative defects in order to generate different types of coherent structures has been studied recently in [17]. Other research propositions of dissipation as a mechanism able to support or induce nonlinear phenomena can be found in [18].
Finally, we point out that although our focus is given to atomic BECs, our results are also relevant in a large variety of contexts (whenever the DNLSE is adequate). The most prominent among them is light emittance from arrays of coupled nonlinear waveguides [19] where boundary dissipation can be realized with suitable mirrors.
The structure of the paper is as follows: in section 2, we present the mathematical model that describes the dynamics of an atomic BEC in a leaking OL and commend on its limitation to describe BECs in one-dimensional OLs. In section 3, we introduce the notion of DBs as solutions of the DNLSE and analyze the temporal behavior of the participation ratio as a measure of the atomic occupation inside the leaking OL. The analysis of the avalanche distribution is done in section 4. In section 5, we present a heuristic model that explains the distribution of avalanches based on the phase-space structure of a three-site lattice. In section 6, we present our conclusions.

Mathematical modeling of leaking OLs
The Bose-Hubbard Hamiltonian is the simplest model that captures the dynamics of a dilute gas of bosonic atoms in a deep OL with chemical potential small compared to the vibrational level spacing. In the case of weak interatomic interactions (superfluid limit) and/or a large number of atoms per well (so that the total number of atoms N ∼ O(10 4 -10 5 ) is much bigger than the number of wells M), a further simplification is available since the BECs dynamics allow a semiclassical mean field description [20,21]. The resulting Hamiltonian is where n is the well index, |ψ n (t)| 2 ≡ N n (t) is the mean number of bosons at well n, U = 4πh 2 a s V eff /m describes the interaction between two atoms on a single site (V eff is the effective mode volume of each site, m is the atomic mass and a s is the s-wave atomic scattering length), µ n is the well chemical potential and T is the tunneling amplitude. Above, we have assumed that only the fundamental energy band is involved in the dynamics [22].
The 'wavefunction amplitudes' ψ n (t) ≡ √ N n (t) exp(−iφ n (t)) can be used as conjugate variables with respect to the Hamiltonian iH leading to a set of canonical equations i∂ t ψ n = ∂H/∂ψ * n ; i∂ t ψ * n = −∂H/∂ψ n . (2) Substituting (1) into the canonical equations of motion, we get the DNLSE.
To simulate the leaking process at the two edges, we supplement the standard DNLSE with a local dissipation at the two edges of the lattice [14,23,24]. The resulting leaking DNLSE is given by where N n = |ψ n | 2 is the mean atomic population at site n. The time τ , the interatomic interaction χ = 2U/T , the well chemical potential µ n and the atom emission probability γ describing atomic losses from the boundary of the OL [14] are measured in units of the tunneling rate. The results reported below correspond to µ n = 0 and dissipation rate γ = 0.2 (for some estimations on γ see [14]). Nevertheless, we have checked numerically that the same qualitative behavior results for other values of γ . An OL with leaking edges can be realized experimentally, by the action of two separate continuous microwave or Raman lasers to locally spin-flip atoms inside the BEC to a nontrapped state [14,15]. Spatially localized microwave fields focused below the wavelength can be obtained at the tip of tapered waveguides. The spin-flipped atoms do not experience the magnetic trapping potential and are released through gravity in two atomic beams at the ends of the OL.
Below, we study the decay and the statistical properties of the total atom population inside the OL, as a function of the initial effective interaction per well, i.e. = χρ, where ρ = N (τ = 0)/M is the initial average density of atoms in the OL. In our numerical experiments we have used initial conditions with on-site populations N n (τ = 0) having only small random fluctuations across the OL while the phases φ n of the 'wave function amplitudes' ψ n = √ N n exp(−iφ n ) have been drawn from a binary distribution with values φ = 0, π taken with equal probability. We mark that in the superfluid regime the phases between neighboring lattice sites can be adjusted by applying a magnetic field gradient (Bloch oscillation) [26]. It is worth noting that the 4 results reported below remain unchanged if the initial phases between neighboring sites alternate between 0 and π or (in the other extreme) if we assume that they are random variables uniformly distributed between [0, 2π ]. Such preparations have to be contrasted with the situation where all sites have the same phase (i.e. the initial preparation is randomized only via the on-site population). In the latter case, we were not able to observe breathers and our analysis below is inapplicable. Further randomization of the initial state is achieved, by imposing a conservative (i.e. γ = 0) transient period, typically τ = 500. After this transient is completed, the dissipation at the lattice boundaries is switched on, leading to a progressive loss of atoms. In all our calculations, we have normalized the wave functions such that N (τ = 0) = 1.
Finally, we would like to comment on the validity of the DNLSE (2, 3) to describe actual experiments of BECs in OLs. Experimental realizations have been performed for values of χ = 2U/T in the range 10 −5 -10 −3 , while the number of atoms is typically N ∼ O(10 4 -10 5 ). These estimations lead to experimentally feasible parameters 1 for which the DNLSE is a good approximation. For example, the experiment of [3] shows that the DNLSE, describes very well the BEC dynamics in an OL with parameters N ∼ 2 × 10 5 , is the recoil energy and k L is the laser mode which traps the atoms), and M ≈ 200 wells, correspond to a value of the control parameter ≈ 0.5. Below, for completeness of the theoretical analysis, we present results for a larger parametric regime of -values. Nevertheless, we point out that the interesting -regime for which avalanches of BEC population are observed within the DNLSE approximation, is within the limits of applicability of the DNLSE to describe BEC dynamics in one-dimensional OLs.

DBs and leaking OLs
An exciting result appearing in the frame of nonlinear lattices is the existence of spatially localized, time-periodic, stable solutions, termed DBs, which emerge due to the nonlinearity and discreteness of the system. DBs were observed in various experimental setups [27]- [32] while their existence and stability were studied thoroughly during the last decade [33], [35]- [38]. Their importance was recognized in [39], where it was shown that they act as virtual bottlenecks which slow down the relaxation processes in generic nonlinear lattices [38]- [40]. Further works [14,40] established the fact that absorbing boundaries can take generic initial conditions towards self-trapped DBs.
Since we are interested in the effects of DBs on the relaxation process of the DNLSE (3), we introduce a localization parameter PR which is a measure of the relative number of sites that are still occupied by the remaining atoms in the leaking OL. It is defined as where · · · indicates an average over different initial conditions. For γ = 0, equation (5) gives the standard participation ratio. Accordingly, the more evenly the atoms spread over the lattice, the closer PR is to a constant of O(1). PR approaches two limiting values: (a) for U = = 0 (linear regime), the norms are exponentially distributed, leading to PR = 1/2 and (b) for → ∞ (or T → 0), the wells are un-coupled and the norms are uniformly distributed, leading to PR = 5/9, which might be viewed as the formation of O(M) DBs (multiple breather regime). For γ = 0, the transition between these limiting cases is smooth (see figure 1). In the open system (γ > 0) the quantity PR = PR(τ ) is of course time dependent. After a transient, however, it rapidly approaches a constant value PR S indicating the formation of a quasi-steady state (see inset of figure 1). In the following, we study PR S as a function of the interaction strength . Instead of a smooth transition between the two extremes (as in a closed system), we observe a sharp drop of PR S at a critical interaction strength of b ≈ 0.15 (our numerics indicate that this drop indeed becomes a step function in the limit M → ∞). This is associated to the self-trap transition phenomena analyzed in the frame of one-dimensional DNLSE lattices (see for example [41,42] 5 ). PR S drops down to its lowest possible value (∝ 1/M) corresponding to a single occupied site, i.e. the final state consists of a single selftrapped DB. We mark that for < b , the atomic population N (τ ) decays smoothly to zero, following the same behavior as for the = 0 even though the DNLSE allows for breather solutions in this range as well 6 . For a representative evolution of the atomic population for < b , see figure 2(a). In the opposite limit of b (see figure 2(c)), we observed the formation of O(M) DBs (multiple breather regime), as in the case of closed systems.
The most interesting dynamics emerges for interatomic interaction strengths in a critical range of L < < U (see figure 2(b)), where we find multiple breather states with a scale-free distribution of single site norms N n . Figure 3 shows P N n (x) for two different system sizes. It clearly displays an inverse power law with a value β = 1.9 ± 0.05 given by a best least square fit (the cutoff for small N n seen in figure 3 can be shifted to arbitrarily small values for larger OLs). In our numerics (with system sizes up to M = 4096), we found U ≈ 2 and L ≈ 0.5, with a strong indication that in the limit M → ∞ this lower bound approaches b . We note again that the DNLSE can describe satisfactory the BEC dynamics in an OL for 1 (see discussion in section 2).

Leaking dynamics and avalanches at the critical -range
At the critical -regime, we have found that the atom population N (τ ) leaks out of the lattice in sudden bursts (see the lower left inset of figure 4) of magnitude δ N = J . We are interested in x -1.9 x P(x) Figure 3. Power law distribution of norms P(x = |ψ n | 2 ) ∼ x −β for = 1. The best least square fit (dashed blue line) indicates that β = 1.9 ± 0.05 ≈ α.
their distribution P(J ) for an ensemble of initial conditions. To this end, we have numerically defined the burst by a threshold D * in the derivative D(τ ) = |dN (τ )/dτ |: whenever D(τ ) rises above D * we register a burst until it drops again below the threshold. In most of our calculations we have used D * ≈ 10 −5 . We have additionally checked that our numerical results for the distribution P(J ) are stable for other choices of the threshold and for different sampling intervals of the N (τ ) time-series. In all cases studied, we had at least 10 4 trajectories for statistical processing. We have found that the avalanches follow a scale-free distribution, as demonstrated in figure 4 for the case of = 1. In fact, the size J of the observed avalanches is proportional to the site-norms N n of the critical states (see next paragraphs). Therefore, one expects that P N n (x) and P(J ) follow the same decay law. Indeed, the best least square fit to our numerical data give α ≈ 1.86 ≈ β. For smaller , the distribution P(J ) shows clear deviations from the power law. By increasing M, however, these deviations become smaller (see figure 4 for = 0.5), indicating L → b in the M → ∞ limit. For L > L U the distribution P N n (x) shows a very nice exponential behavior (see upper right inset of figure 4). Let us study in more detail the dynamics that lead to the creation of an avalanche in the critical regime. One such event is depicted in figure 5. An excitation, coming from the bulk of the lattice, collides with the outermost self-trapped DB. As a result the breather is shifted inwards by one site while at the same time a particle density proportional to the density of the excitation tunnels through the DB. This particle density eventually reaches the leaking edge of the OL and decays in a form of an avalanche. Details of the migration process are still under investigation [25,44].
To investigate theoretically the collision process that involves the self-trapped DB, we will make use of a general concept called the local ansatz [33,34]. The DB solutions in the DNLSE are very well localized and the most basic and important breathers occupy only three sites. Within the local ansatz, only the central site and the two neighboring sites of a breather-like object are considered, allowing us to turn a high dimensional dynamical problem involving M sites, into a reduced problem with three degrees of freedom (nonlinear trimer corresponding to equation (3) with M = 3 sites and γ = 0). A detailed analysis of the reduced problem [33,34], shows that (i) the breather-like object corresponds to a trajectory in the phase space of the full system which is practically embedded on a two-dimensional torus manifold, thus being quasiperiodic in time and (ii) the breather-like object can be reproduced within a reduced system. In figure 6 (inset) we show a representative Poincaré section of the nonlinear trimer for ≈ 1. We make use of the variables (N 2 , φ 3 − φ 2 ) where N 2 is the norm at site n = 2 of the M = 3 DNLSE (3), (the central site of the DB), and φ 3 − φ 2 is the phase mismatch between the second and the third site. The Poincaré section correspond to the plane φ 1 = φ 3 anḋ φ 1 >φ 2 of the energy surface. The phase space consists of two components involving islands of regular motion embedded in a chaotic sea. Trajectories inside the islands correspond to DBs, provided that their frequency is outside the linear spectrum. In contrast, chaotic trajectories have continuous Fourier spectra, parts of which overlap with the linear spectrum of the infinite lattice [33]. As long as the DB is stable, it acts as a barrier which prevents atoms from reaching the leaking boundary. A necessary condition for an avalanche event is thus the destabilization of the DB. This is caused by a lattice excitation with particle density N prt 1 greater than a threshold, which can push the regular (quasiperiodic) orbit out of the island towards the chaotic sea. One needs a perturbation that is of the order of the linear size S of the island (which represents a DB). At the same time a portion J ≈ N max 3 ∝ N pert 1 ∼ S is transmitted through the DB (see figure 6).

Distribution of island sizes in mixed phase space
The above analysis (in the frame of the local ansatz), allows us to turn the problem of the analysis of P(J ) into the study of the distribution of island sizes P(S) of the reduced system, with M = 3 in the -regime where the phase space is mixed. Let us consider a simple heuristic model that mimics the hierarchical ('island over island') structure of a typical mixed phase space, where at each hierarchy level k > 0 there exists a main island of size S k−1 with a number n k of subislands. At the kth hierarchy level, the fraction of sizes of the main island to a sub-island is f k = S k−1 S k . Choosing S 0 = 1 and making the additional simplifications f k = f and n k = n for each hierarchy level, we get that S k = f −k while for the total number of islands p k we get p k = k j=1 n j = n k . Hence, the distribution of island sizes reads P(S) = n k(S) = S −α where α = ln n/ ln f. Due to the self-similar structure of the phase space, we have to require that the total number of islands diverges. This condition, together with the fact that the total phase-space volume is finite, provides two bounds for the power-law exponent 1 < α < 3, (9) which is in excellent agreement with our numerical findings (see equation (7)). The power-law distribution, equations (8,9), is supposed to be a generic feature of systems with a hierarchical mixed phase space in 2D, of which the Poincaré section of the trimer is an example. This seems a likely supposition, however, due to the complicated phase-space structure of mixed chaotic systems it is a formidable task to directly verify it numerically.
To get further confidence in our heuristic modeling, we sought an indirect verification of the power-law distribution of island sizes particular suited to our line of questioning. We did this for a paradigmatic model of mixed phase-space dynamics, the kicked rotor [43]. The classical Hamiltonian has the form of a pendulum in a gravitational field which is given as a periodic sequence of kicks:  (12)) for the kicked rotor with kicking strength K = 3.5, corresponding to a mixed phase space (inset). The island sizes were evaluated by studying the co-evolution of close-by initial conditions as a function of their separation for successively longer times τ .
In the limit of τ → ∞ the distribution converges to a power law.
Above p is the angular momentum, θ is the angular displacement, I is the moment of inertia (below we will assume I = 1) and k is the perturbation strength. Such a time-dependent potential is very convenient both for analytical and numerical calculations and allows one to use the following mapping of angular momentum and displacement: where P n = p n τ , K = kτ and both (P, θ) are taken modulo 2π (torus geometry). The investigation of this map is computationally much more convenient than the nonlinear trimer. For K = 3.5 map (11) possesses a typical mixed phase space (see inset of figure 7). Using the kicked rotor, we have investigated the escape of nearby trajectories in the surroundings of an island: we evolved a random ensemble of pairs of trajectories, that were initially (n = 0) separated by a distance D in phase space. We then calculated the probability density P in,out (D) of one trajectory staying in a given region surrounding an island, up to a large cut-off time n = T (i.e. it starts in or very near an island) while the other trajectory leaves the region before n = T (thus being a chaotic trajectory). Provided that equations (8,9) apply, one finds that this probability (in the limit of T → ∞) takes the form where the constant term C 0 reflects the existence of an upper cut-off for the island sizes. We present in figure 7 the outcome of our numerical calculations for increasing values of T .
Our numerical data confirms that in the asymptotic limit T → ∞ the distribution P in,out (D) approaches the form equation (12). This also confirms the power law scaling equation (8) of the island sizes over more than three orders of magnitude.

Conclusions
In this paper, we have investigated the statistical properties of the population decay N (τ ) of an atomic BEC in an OL with leaking boundaries. We have found that for interatomic interactions greater than a critical value b , the atomic population N (τ ) decays in sudden bursts (avalanches) of size J = δ N (τ ) which are followed by time-intervals for which N (τ ) remains almost constant. This behavior is related to the formation of DB. Specifically, the DBs closest to the leaking boundaries act as barriers that trap the atoms inside the bulk of the OL. Collisions of these DBs with excitations created in the bulk of the OL lead to their destabilization. During this process, a portion of the atoms hitting the DB from the inside tunnel through it and reach the boundary to create an avalanche.
For very large values of the interatomic interactions > U , the distribution of avalanches P(J ) is exponential, while for moderate -values, L < < U we find a power law avalanche distribution P(J ) ∼ 1/J α ; the latter indicating the existence of a phase-transition. Using the local ansatz approach, we were able to relate the appearance of the scale-free avalanches, with the hierarchical structure of the underlying mixed phase space. Specifically, we have argued that the avalanche size is proportional to the size of the phase-space islands, associated with the out-most DB. A heuristic model was constructed which explained the power law distribution of island sizes, being a generic property of dynamical systems with mixed phase space. The outcome of the heuristic modeling has been further confirmed against numerical calculations with the standard map (kicked rotor) in the mixed phase-space regime.
We expect that our results will shed new light on the problem of relaxation in nonlinear lattices. In the specific case of an atomic BEC loaded within a leaking OL, the time derivative of N (τ ) is equal to the outgoing atomic current i.e. J (t) = − ∂ N (τ ) ∂τ , which is an experimentally accessible quantity. A measurement of the emitted atoms can be used to infer information about the interatomic interactions that dictate the dynamics inside the OL.