Non-linear mixing of Bogoliubov modes in a bosonic Josephson junction

We revisit the dynamics of a Bose-Einstein condensate in a double-well potential, from the regime of Josephson plasma oscillations to the self-trapping regime, by means of the Bogoliubov quasiparticle projection method. For very small imbalance between left and right wells only the lowest Bogoliubov mode is significantly occupied. In this regime the system performs plasma oscillations at the corresponding frequency, and the evolution of the condensate is characterized by a periodic transfer of population between the ground and the first excited state. As the initial imbalance is increased, more excited modes -- though initially not macroscopically occupied -- get coupled during the evolution of the system. Since their population also varies with time, the frequency spectrum of the imbalance turns out to be still peaked around a single frequency, which is continuously shifted towards lower values. The nonlinear mixing between Bogoliubov modes eventually drives the system into the the self-trapping regime, when the population of the ground state can be transferred completely to the excited states at some time during the evolution. For simplicity, here we consider a one-dimensional setup, but the results are expected to hold also in higher dimensions.


I. INTRODUCTION
Two weakly coupled Bose-Einstein condensates (BECs) in a double-well potential constitute a paradigmatic system for investigating the physics of bosonic Josephson junctions [1][2][3][4]. Owing to the nonlinear character of the interactions, this system exhibits different dynamical behaviors, ranging from Josephson plasma oscillations (in the limit of very small imbalance between the population of the two wells) [5], to macroscopic selftrapping where -above a critical value of the imbalance -the population of the two wells is almost locked to the initial value [3,4,6]. Due to the conceptual importance of these phenomena, BECs in double-well potentials and arrays of coupled boson Josephson junctions have been extensively investigated in the last two decades both theoretically [2,3, and experimentally [4,6,[30][31][32][33][34][35][36][37][38], as well as their counterparts with fermionic superfluid atomic samples [39][40][41][42].
The physics of these systems is well captured by a two-mode approximation of the Gross-Pitaevskii equation, each mode being localized in one of the two wells, which allows for an effective description in terms of only two parameters, namely the population imbalance z(t) and the phase difference φ(t) between the left and right components. Here we provide a complementary description by means of the quasiparticle projection method of Ref. [43], extending the Bogoliubov treatment of Ref. [29] to the case of arbitrary initial imbalance. For the sake of simplicity, we shall restrict the analysis to the case of a (quasi) one-dimensional condensate [44] We find that in the regime of a small initial imbalance, where only one Bogoliubov mode is significantly occupied and the system performs plasma oscillations at the corresponding frequency [29], the evolution of the condensate is characterized by a periodic transfer of popula-tion between the ground state and the first excited state. As the initial imbalance is increased, more Bogoliubov modes get coupled during the evolution of the system, and their population also varies with time, contrarily to what happens in a linear system. As a consequence, the frequency spectrum of the imbalance turns out to be still peaked around a single frequency which is shifted towards lower values, rather than getting relevant contributions at higher frequencies, where Bogoliubov modes are located. By further increasing the initial imbalance, the population of the ground state can be completely transferred to the excited states at some time during the evolution, driving the system into the macroscopic selftrapping regime.
The paper is organized as follow. In Sec. II we introduce the formalism, reviewing the definition of the two-mode approach (II A) and of the quasiparticle Bogoliubov expansion (II B). Then, in Sec. III we present the results by discussing the behavior of the system in the regime of Josephson plasma oscillations (III A), the selftrapping regime (III C), and that intermediate between the former two (III B), highlighting the role of non-linear mixing (III D). Final considerations are drawn in the conclusions.

II. MODEL
Let us consider the following (dimensionless) Gross-Pitaevskii equation [29] i∂ t ψ(x, t) = − 1 2 The prediction of the TM model [3] in Eq. (9) is also shown as a reference (dotted line).
and dx|ψ(x)| 2 = 1, describing a (quasi) one dimensional condensate trapped in a double-well potential. The latter is composed by a harmonic potential term, plus a barrier of intensity V 0 and width w, with δx providing a relative shift between the two. Here we are interested in describing the dynamics triggered by an initial population imbalance between the two wells. This can be obtained by preparing the system in the ground state ψ g (x) = ψ(x, 0) of the above potential with δx = 0, and then suddenly switching δx = 0 at t = 0 (only the harmonic potential is shifted, the barrier does not move). The ground state ψ g (x) is obtained from with µ being the condensate chemical potential. As for the parameters, here we choose w = 0.3 and V 0 = 50, that correspond to a double-well configuration within reach of current experiments (see e.g. Ref. [42]), whereas the interaction strength u 0 and the initial shift δx are taken as free parameters, and will be varied for exploring different regimes (see later on). In particular, δx, is chosen in order to produce the desired initial imbalance z 0 . As it is known, in the limit of very small initial imbalance the system performs Josephson plasma oscillations [1], and eventually enters a self-trapping (ST) regime at a critical imbalance [3] whose specific value depends on the strength u 0 of the nonlinear term (see Fig. 1). The dynamics of the system will be analyzed by means of an expansion over the Bogoliubov modes, by comparing with the exact evolution and the two-mode (TM) approach.
A. Two-mode model Usually, the dynamics of a condensate in a double-well potential is treated by means of the two-mode approach, which consists in writing the condensate wave function as (see e.g. [29,45] and references therein) where the functions ψ L,R (x) are localized in the left and right well, have unit norm, and are orthogonal to each other, ψ L |ψ R = 0. Though somewhat approximate,and not entirely justified from the formal point of view [29] -the two-mode model provides an effective description of the double-well system in several respects, and will be used in the rest of the discussion as a reference. Here we construct the two-modes ψ L,R (x) from the ground state ψ g (x) (symmetric) and the first-excited solution ψ 1 (x) (antisymmetric) of the stationary GP equation [46]. Namely, we take the following linear combination, , corresponding to the most common approach in the literature [9,13,18,33,45,[47][48][49]. Then, by defining (α = L, R) and one gets the following equations for the phase difference In the following, this set of equations will be referred to as the full two-mode (fTM) model. When the terms Λ 1 and Λ 2 can be neglected, it reduces to the well-known two-mode (TM) model by Smerzi et al. [3]. We recall that the TM model predicts that the system enters the ST regime when the parameter Λ exceeds a critical value Λ c . For φ 0 = 0, it takes the following value: Here, we shall use as independent parameter u 0 rather than Λ (which will depend on u 0 ). Then, the previous equation can be easily inverted, yielding As shown in Fig. 1, this formula provides a good estimate for the actual critical imbalance extracted from the GP equation, in the whole range considered (u 0 ∈ [1, 200]).

B. Bogoliubov approach
As a complementary description, here we employ the quasiparticle projection method introduced by Morgan et al. in Ref. [43]. It amounts to a Bogoliubov expansion [50,51] where the condensate and quasiparticle populations are allowed to vary with time, namely The functionsũ i (x) andṽ i (x) are the Bogoliubov eigenmodes, with the tilde indicating that they are chosen to be orthogonal to ψ g (x) [52]. They are solutions of (from now on we fix ψ * g = ψ g without loss of generality) [53] L gψ 2 with The solutions of Eq. (12) satisfy the following orthogonality relations [54] dx The coefficients b g (t) and b i (t) are given by [43] b In the linear regime, when the modes remain decoupled during the whole evolution, the coefficients b i (t) are solutions of iḃ i (t) = ω i b i (t), namely [43,50] where the coefficients b i0 ≡ b i (0) (which do not depend on time) are fixed by the initial conditions, see Eq.
Imbalance. To construct the population imbalance between the right and left wells we start by integrating the particle density over the positive and negative x semi-axis. By taking into account the symmetries of the problem we have with Remarkably, only the Bogoliubov excitations with odd i contribute to the imbalance, owing to the symmetries of the system. In the linear regime we have (using also the fact that in our caseũ i ,ṽ i can be chosen real without loss of generality)

III. RESULTS AND DISCUSSION
Here we shall discuss the evolution of the imbalance for different values of its initial value z 0 ≡ z(0) (throughout this work we set φ 0 ≡ φ(0) = 0), discussing the behavior of the system in terms of the quasiparticle projection method [43] introduced in the previous section. The general behavior of z(t) has already been extensively studied, at least in the framework of the TM model, see e.g. the seminal Ref. [3]. In the rest of this paper we fix the ratio µ/V 0 ≡ 0.25, a value that characterizes a typical Josephson regime (with the chemical potential much lower than the barrier height [45]). The explicit behavior of the imbalance evolution is shown in Fig. 2 for u 0 = 4 and z 0 = 0.1, 0.3, 0.5, 0.7 (empty squares in Fig. 1), ranging from the regime of Josephson plasma oscillations (z 0 0.1), to the ST regime (z 0 0.62). A detailed description of the different dynamical behaviors and of the various lines plotted in the figure is given in the following.

A. Josephson plasma oscillations
In the limit of very small imbalance, the system performs Josephson plasma oscillations characterized by a frequency ω J . This frequency corresponds to the energy of the lowest Bogoliubov mode [29]. In fact, in this limit only one Bogoliubov mode is occupied, the system is in the linear regime, and z(t) is well reproduced by Eqs. (23), (24) with the only contribution of B 01 , namely [29] z(t) = 4B 01 cos(ω 1 t).
This is shown in Fig. 2(a), where the GP prediction (solid purple line) is perfectly reproduced by that of Eq.  Fig. 2(a). This difference may increase further by increasing u 0 [29].
In this regime we also have , where n g (t) ≡ |1 + b g (t)| 2 represents the (relative) population of the ground state, and n e1 (t) ≡ |b 10 | 2 dx |ũ 1 | 2 + |ṽ 1 | 2 + 2b 2 10 cos(2ω 1 t) dxũ 1ṽ1 that of the first Bogoliubov excitation, the other excited modes being essentially irrelevant. The evolution of n g and n e1 is plotted in Fig. 3(a) for z 0 = 0.1 (the other three panels will be discussed later on). A sinusoidal oscillation -with frequency 2ω 1 , see Eq. (28) -is clearly visible in Fig. 3(a). It corresponds to a (small) periodic transfer of population between the ground state and the first excited state, contrarily to what happens in a truly linear system, where the occupation number of each energy level is constant.

B. Intermediate regime
In general, when one increases the initial imbalance z 0 , the form of z(t) changes, and its frequency (the inverse of the period) changes as well [3], see Fig. 2(b) and 2(c). In particular, before entering the ST regime, the imbalance is still characterized by periodic oscillations, but with a frequency ω that is shifted with respect to the plasma value ω J as an effect of the nonlinearity, see Fig. 4(a). These changes are reflected in the change of the Fourier spectrum, in Fig. 4(b). Remarkably, in this regime the spectrum is still peaked around a single frequency, that is shifted continuously towards lower frequency values with respect to ω J ≡ ω 1 , contrarily to the naive expectation of having more Bogoliubov modes macroscopically occupied (the first Bogoliubov frequencies is indicated by the solid (red) point on the horizontal axis of Fig. 4(b), higher modes lye far outside the present range). In fact, we find that the system exits the linear regime, namely Eq. (24) fails in reproducing the actual behavior of z(t) -see Fig.  2(b) and 2(c) -even if higher Bogoliubov modes have an initial population that is still below 1% that of the lowest mode. In other words, the linear approach fails not because some of the other excited modes are initially macroscopically occupied (as it would be the case for a truly linear system), but because of the nonlinear mixing during the evolution of the system. In this regime, the analog of the decomposition (28) becomes more complicated as the contribution of all the excited modes to the total density is now meaning that it is not possible to write the total density as the sum of separate contributions of each Bogoliubov mode. The evolution of n g (t) and n e (t) is shown in Fig.  3 for increasing values of the initial imbalance z 0 . This figure shows that the transfer of population between the ground and the excited states increases by increasing z 0 . Initially, when the system exits the linear regime but z 0 is not too large [e.g. z 0 = 0.3, Fig. 3(b)], n g (t) and n e (t) are still characterized by sinusoidal oscillations. For larger values of z 0 , the oscillations in the population deviates from this behavior, as it does the corresponding imbalance [see e.g. Fig. 3(c) and Fig. 2(c)]. In any case, the oscillations of n g (t) are always in phase with those of |z(t)| (that is, the maximal imbalance is obtained when the population of the ground state is maximal).

C. Self-trapping regime
By further increasing the initial imbalance z 0 , the period of z(t) gets larger and larger (see also [3]), and eventually diverges at the critical value z 0c where ω ∝ 1/T → 0, see Fig. 4(a) (z 0c ≃ 0.62 in the present case). Notice that the value of z 0c obtained from the solution of the GP equation is reproduced with great accuracy by the fTM model (we have verified that this holds true even for values of u 0 larger than that considered in the present paper). Remarkably, the onset of ST corresponds to a situation in which the population of the ground state can be transferred completely to the excited states, namely when n g (t) = 0 at some t during the evolution, see Fig.  3(d). This feature is indeed a distinctive characteristic of the ST regime. In this regime the imbalance is stuck on the positive side (or the negative one, depending on the initial conditions), still oscillating, but with an irregular pattern [3]. The latter reflects in the shape of the frequency spectrum, that significantly broadens and acquires a relevant contribution from the low frequency region, ω ≃ 0 [dotted-dashed orange line in Fig. 4(b)].

D. Non-linear mixing
Owing to the coupling between the different Bogoliubov modes, see Eq. (29), we argue that |b i (t)| 2 cannot be identified with the occupation number of the i-th quasiparticle level, contrarily to the the interpretation given in Ref. [43]. However, since the coefficients b i (t) represent the quasiparticle amplitudes in the sense of the expansion (11), in the following we shall consider their modulus squared |b i (t)| 2 as a measure of the weight of each mode in the system dynamics. Their evolution (for i = 1, 2, 3, 4) is shown in Fig. 5(left), along with the corresponding time-averaged values for the same values of the initial imbalance as in the previous figures, namely z 0 = 0.1, 0.3, 0.5, 0.7. In Fig.  5(right) we show the corresponding region of the complex plane spanned by the real and imaginary parts of b i (t) (here normalized to b 1 (0), for easiness of visualization) during the evolution of the system. In Fig. 5(e) we also show the trajectory of the lowest Bogoliubov mode (i = 1) for z 0 = 0.005, indicating that in the limit z 0 → 0 the expected behavior is recovered: in this case the coefficient b 1 (t) is constant in modulus as dictated by Eq. (18) for the linear regime, and the contribution of all higher excited modes is negligible [29]. As z 0 is increased, the dynamics in the complex plane becomes chaotic, each mode spanning a larger portion of the plane. Notice also the change in the orbit shape, from circular (in the limit z 0 → 0) to elliptical (for z 0 0.1). Both left and right panels evidence a mixing between different modes, especially those with i = 2 and i = 3. Remarkably, the mode  i = 2 -which, being even, does not contribute directly to the imbalance, see Eqs. (22) and (23) -indeed affects it through the mixing with other modes during the evolution of the system.

IV. CONCLUSIONS
We have analyzed the dynamics of a (quasi) onedimensional Bose-Einstein condensate in a double-well potential, from the regime of Josephson plasma oscillations to the self-trapping regime, by means of the Bogoliubov quasiparticle projection method [43]. In the limit of very small initial imbalance, the system performs Josephson plasma oscillations characterized by the frequency of the lowest Bogoliubov mode (the only Bogoliubov mode being significantly occupied) [29]. In this regime, the evolution of the system is characterized by a periodic transfer of population between the ground state and the first excited state. As the initial imbalance is increased, the system still performs periodic oscillations between the left and right wells, but with a frequency that is continuously shifted towards values lower than the plasma frequency. This occurs because of the nonlinear mixing of the Bo-goliubov modes during the evolution of the system, and not because some of the excited modes (besides the lowest one) are initially macroscopically occupied, contrarily to what happens in a linear system. The frequency spectrum of the imbalance is therefore still peaked around a single frequency, and the corresponding period diverges when the system enters the self-trapping regime. This corresponds to a situation in which the population of the ground state can be transferred completely to the excited states at some time during the evolution. This feature is indeed a distinctive characteristic of the ST regime. The present picture is expected to hold also in higher dimensions.