Orbital Josephson effect and interactions in driven atom condensates on a ring

In a system of ac-driven condensed bosons we study a new type of Josephson effect occurring between states sharing the same region of space and the same internal atom structure. We first develop a technique to calculate the long time dynamics of a driven interacting many-body system. For resonant frequencies, this dynamics can be shown to derive from an effective time-independent Hamiltonian which is expressed in terms of standard creation and annihilation operators. Within the subspace of resonant states, and if the undriven states are plane waves, a locally repulsive interaction between bosons translates into an effective attraction. We apply the method to study the effect of interactions on the coherent ratchet current of an asymmetrically driven boson system. We find a wealth of dynamical regimes which includes Rabi oscillations, self-trapping, and chaotic behavior. In the latter case, a full many-body calculation deviates from the mean-field results by predicting large quantum fluctuations of the relative particle number.


Introduction
The Josephson effect (JE) is a fundamental quantum phenomenon which reflects the coherent occupation of two or a few single-particle states by a macroscopic number of bosons. It was first predicted [1] and observed [2] for superconductors, and seen later in superfluids [3]. In Bose-Einstein condensates (BECs) one may refer to the internal or external JE, depending on whether the two coherently connected states occupy the same region of space with different internal atomic states, or two different regions with the same spin state. Both the internal [4,5] and the external [6] JE have been observed. The internal JE holds the promise of generating highly entangled quantum many-body systems [7]. Here we report on a new type of JE which involves the macroscopic coherent occupation of Floquet states in ac-driven BECs. This orbital Josephson effect is properly neither external nor internal, since the connected Floquet states occupy the same region of space with the same internal atom structure, their sole difference residing in their orbital state.
The Josephson effect requires interactions in order to display truly collective behavior associated with the macroscopic occupation of two states [8,9]. In the absence of interactions, the resulting Rabi dynamics is merely an amplification of the dynamics undergone by a single atom [4,9]. The role of interactions in ac-driven many-body systems is difficult to treat, and frequently simplifications are made such as the two-mode approximation [10], mean field theory [11,12,13], and its first [14,15] and second [16,17] order corrections, or the use of an effective description [18] in terms of a static manybody Hamiltonian with renormalized parameters. This latter approach is valid when the driving frequency is the dominant energy scale. Alternatively such systems can be numerically studied by the exact simulation of small clusters [18], or by recently developed techniques such as the time-dependent density-matrix renormalization group (t-DMRG) method [19], or multiconfigurational time-dependent Hartree for bosons (MCTDHB) [20], both of which are able to treat larger systems. Here we develop a description of the coarse-grained dynamics of the quantum field operator in an ac-driven many-boson system, thus going beyond a mean-field treatment and its first corrections. We find that, in the case of resonantly connected Floquet states, the long-time evolution can be described by conventional quantum dynamics, by which we mean one deriving from an effective interacting Hamiltonian. If the driving provides the only external potential, we find that a repulsive interaction in real space translates into an attractive interaction in the truncated Hilbert space of a discrete number of resonantly connected Floquet states. We apply these findings to understand in depth the role of interactions in asymmetrically driven BECs which exhibit the quantum ratchet effect [11,12,21,22], together with a wealth of dynamic regimes.
-+ -+ Figure 1. Schematics of the orbital Josephson effect. (a) Resonant driving induces coupling between the zero-momentum and two opposite momentum eigenmodes. The resulting Josephson link between these orbitals, with Rabi frequencies 2Γ ± , has a negative effective interaction energy. (b) A possible departure from exact resonance lifts the degeneracy of the quasi-energies.

(t, t ′ ) formalism in second quantization
Our starting point is the equation of motion for the field operator of a system of interacting bosons, where H(x, t) = − 1 2 ∂ xx + V (x, t) includes a time-dependent potential, and λ is the effective strength of the contact interaction in a quasi-one-dimensional ring of radius R. We set = 1 and measure all energies and frequencies in units of 2 /MR 2 , with M the atomic mass. Furthermore, our unit of length is R, so that the circumference of the ring is 2π. Such rings, as schematically depicted in Fig. 1a, have been proposed and experimentally realized -see e.g. Refs. [23,24], among many other ring experiments. A possible particle flow can be observed via standard time-of-flight measurement techniques [23]. The field operatorψ(x, t) obeys periodic boundary conditions in x, and its algebraic structure is given by the standard equal-time bosonic commutation relations: In order to study the long-term dynamics of the many-body system, we generalize the (t, t ′ ) formalism to treat second-quantized operators. This formalism was originally developed for single particles [25], and later adapted to treat helium [26] always within first quantization. We do this extension by introducing an additional parameter t ′ which the field operators depend on, and impose periodic boundary conditions in t ′ with period Tψ We defineψ(x, t ′ ; t) as the solution of the equation of motion together with the initial condition The periodic boundary conditions in t ′ , which are imposed initially via Eq. (4), are guaranteed to be preserved over time t due to the periodicity of H(x, t ′ ). We note that the field operatorψ(x, t ′ ; t) acts on the same Fock space as the Heisenberg field operatorsψ(x, t), for all t ′ and all times t.
We may choose t ′ to be restricted to be a function of t given by t ′ τ (t) = t + τ , with τ ∈ [0, T ]. Then the resulting operatorψ τ (x, t) ≡ψ(x, t ′ τ ; t) =ψ(x, t + τ ; t) is the solution of the Heisenberg equation (1), but with a shifted switching H τ (x, t) = H(x, t + τ ) and initial conditionψ τ (x, 0) =ψ(x, 0). This can be shown by calculating the total time derivative along the line t ′ = t ′ τ (t): where we have used Eq. (3) and the identity dt ′ τ /dt = 1. In other words,ψ(x, t ′ τ ; t) = ψ τ (x, t) is but a solution of the Heisenberg equation of motion for a Hamiltonian shifted in time by τ with the initial conditionψ τ (x, 0) =ψ(x, 0). In particular t ′ = t (τ = 0) gives the solution to Eq. (1). This implies the following commutation relations for the solutions of Eqs.
The advantage of the (t, t ′ ) approach is that it provides a natural separation of timescales. The t ′ coordinate describes the behavior of the system on timescales shorter than the driving period T , while the long-time dynamics is described by t -see Eqs. (7)-(8) below. Furthermore, the originally time-dependent problem [see Eq. (1)] is mapped to a formally time-independent one [see Eq. (3)], which opens the possibility of using resolution methods developed for such type of problems [21,25,26].
The field-operator extension of the (t, t ′ )-formalism which we have presented is rather general. It is neither restricted to bosonic systems, nor to the ring geometry we study here. It can be straightforwardly translated to the case of fermionic or mixed systems, as well as to other trap-geometries.
Due to the periodic boundary conditions in t ′ , it is natural to express the t ′dependence ofψ(x, t ′ ; t) via a Fourier decomposition, and we may approximate the dynamics by incorporating only a few modes.
To include the effect of a weak time-periodic perturbation, V (x, t) = V (x, t + T ), it is convenient to work in the representation of unperturbed (V = 0) Floquet states. The first-quantized version of (3) was studied in Ref. [21] in the absence of interactions. The unperturbed stationary states evolve as where ε 0 ℓm = 1 2 ℓ 2 − mω and ℓ, m are integers labeling the Fourier modes which are the unperturbed Floquet states with quasi-energy ε 0 ℓm , with ω = 2π/T . The operators in this representation, satisfy the equation of motion and the commutation relations (6) translate to If one were to view m as just an additional standard (orbital or spin) quantum number, these commutation relations would appear unconventional (very much as Eqs. (6) would look unconventional if t ′ were regarded as an extra space variable). However, Eqs. (11) follow naturally from the Fourier transformation of Eqs. (6).
Here V lm,l ′ m ′ is the matrix element of the driving operator between two unperturbed states, as given in Eq. (8). When the driving frequency is such that the system is at or close to resonance, only a few states (all with the same or similar value of ε 0 ℓm ) are relevant [11,12,21]. Importantly, in that subspace the index m is uniquely determined by ℓ. Thus within that truncated space we can drop the index m. As a result, Eqs. (11) become equivalent to the standard bosonic commutation relations [â ℓ (t),â † ℓ ′ (t)] = δ ℓℓ ′ and [â ℓ (t),â ℓ ′ (t)] = 0.

Resonant driving
If we calculate the effective matrix elements connecting the various resonant unperturbed Floquet states, we are left with a conventional few-mode boson problem whose dynamics, given by Eq. (10), can be studied using established techniques. The resulting dynamics between Floquet states reflects the coarse-grained, long-time dynamics of the true quantum state evolution.
We remark that the derivation of an effective Hamiltonian does not require ℓ to describe a plane wave. In such a case, the interaction term in Eq. (12) would contain a sum over four orbital indices, and the matrix elements would depend on all four involved orbitals.
The conservation of total particle number ℓn ℓ = N, together with the identity leads to the interesting result that, up to a constant (2N 2 − N)λ/4π, the Hamiltonian (12) is equivalent tô Thus a repulsive interaction in real space translates into an attractive interaction in (angular) momentum space. The resulting energy gain in the macroscopic occupation of a single state in gases with repulsive interactions underlies the stability of Bose-Einstein condensation, as discussed in Ref. [27].
To derive (14), we have used the selection rules (stemming from the simultaneous conservation of momentum and quasi-energy) which are implicit in Eqs. (10)- (12). Interestingly, in those dynamical regimes where the truncated resonant space is effectively reduced to two modes, Eq. (14) applies for orbitals of arbitrary shape (not necessarily plane waves), but with an interaction strength that is state dependent.
When a few modes are macroscopically occupied, Eq. (12) describes a Josephsontype link between states ℓℓ ′ with Rabi frequency 2Γ ℓℓ ′ . Remarkably, this Josephson link takes place between atom states with the same internal state and occupying the same region of space. Thus it is not appropriate to describe it as an internal or external JE. The Hamiltonian (14) describes what can be termed the orbital Josephson effect, which here occurs between resonant Floquet states in ac-driven Bose condensates.

Application to the quantum ratchet
As a particular application of the method developed here to treat the many-body problem in driven systems, we consider a BEC subject to the asymmetric driving which has been numerically studied [12] and experimentally implemented by Salger et al. [28] on an extended optical lattice. Since the driving in Ref. [28] conserves quasi momentum, the here considered periodic boundary conditions can be related to the experimentally more convenient lattice system. For ϕ = 1 2 π, 3 2 π it yields a coherent quantum ratchet provided that both α and β are nonzero. We consider a system initially prepared in |ℓ = |0 . It was analytically shown in Ref. [21] that, if we drive it with resonant frequency ω = 1 and small amplitude K, the subsequent evolution will mix the initial state only with the states |ℓ = | ± 2 . These three states all satisfy the resonance condition ε 0 ℓ = 0, and we label them with the indices ±, 0. The resulting three-level system Hamiltonian iŝ where ν takes values ±, 0, and ∆ = ω − 1 accounts for a possible small detuning that shifts the quasi-energies in Eq. (14) to ε 0 ℓ = −(ℓ 2 /2)∆. The resulting level structure is schematically depicted in Fig. 1b. Near resonance the tunneling parameters Γ ± can be calculated analytically [21]. For ∆ ≪ 1, one finds This result demonstrates that the quantum ratchet current (which requires |Γ + | = |Γ − |) originates in the interference between first-and second-order processes in the driving strength [21]. In Eq. (17) we have neglected the possible effect of interactions on the second-order tunneling terms. Without interactions, the initial state |0 couples only to a state |a which is an asymmetric combination of |± , yielding an average current which is half that carried by |a [21]. It has been noticed that interactions destroy the coherent current [12]. Below we study the role of interactions in depth and show that they yield a rich variety of dynamical regimes.

Numerical results
We employ three different approximations to investigate the many-body problem: full numerical resolution of the time-dependent Gross-Pitaevskii equation (FGP) using a fourth-order Runge-Kutta method with a fixed time-step; study of the GP equation in the truncated space of the three resonant states 0, ± (3GP); and resolution of the many-body problem in the three-level system (3LS), via exact diagonalization of the effective interacting Hamiltonian (16) with up to 40 atoms. Importantly we note that 3GP can be obtained both as a truncation of FGP or a mean-field version of 3LS. We always assume the condensate to be initially in state 0. Figure 2 shows the expectation value of the ratchet current I(t) -given by the mean momentum per particle -over the first 200 driving cycles. The regimes of weak, moderate, and strong interaction [as characterized by g ≡ λ(N − 1)], are all calculated with the three methods described above. For weak interactions, the two truncatedspace calculations (3GP and 3LS) yield similar Rabi oscillations, both differing from FGP in that they do not display high-frequency dynamics, as revealed both in the Fourier spectrum and in the time dependence of the current. The fast dynamics of FGP disappears if the stroboscopic current is plotted, as shown in the shaded region. For strong interactions, the system tends to remain trapped in the initial state ‡, as may be expected since Eq. (16) has the structure of a Josephson link with negative interaction energy [29]. The three numerical methods predict differing behaviors of the resulting small current, but agree on the position and strength of the main peak in the ‡ In the absence of decoherence the initial state 0 would mix with resonant states ± after a very long time. spectrum. For intermediate interaction, discrepancies between the three methods soon appear, which suggests chaotic behavior. This is confirmed and studied further below.

Time-dependent current
Within the accuracy limited by the ω-sampling, the positions of the main peaks coincide exactly in the Rabi regime, and vary by about 10% in the self-trapping case. The peak heights vary by about 6% in the Rabi case, and about 9% between FGP and 3GP in the the self-trapped regime. However for 3LS, the discrepancy with FGP (and 3GP) is worse (differing by a factor of 3), which coincides with the observation that for the self-trapped dynamics one has to consider particle numbers far beyond N = 40 in order to get an agreement between the full many-body description and the mean-field result [30].

Time-averaged current
after 400 cycles, computed within 3GP, as a function of g and K at exact resonance. We assume |Γ + | > |Γ − |; see Eq. (17). For small K, the perturbative scheme of [21] applies, and at weak interaction a nonzero ratchet current results from Rabi oscillations between |0 and |a . For stronger interactions the ratchet current disappears, because the system remains in the initial state |0 due to the large negative energy associated with the initial macroscopic occupation of |0 [see Eq. (14)]. This can be viewed as a form of macroscopic quantum self-trapping [31], albeit with a negative effective interaction energy. We find that just below the critical interaction strength the system oscillates between |0 and |+ . This is consistent with the fact that the threshold interaction value obtained from that assumption, g c = 8π|Γ + |, agrees well with the numerical simulation (see solid curve in Fig 3a). Finally, we notice the existence of a region of suppressed or weakly reversed current for stronger driving and intermediate interaction. Figure 4a shows the effect of departing from exact resonance. In general, the effect of ∆ > 0 is that of shifting the behavior of current towards higher interactions. In particular we note that, starting from the white-blue region of suppressed or reversed current, interactions tend to restore the positive ratchet current. This effect can be understood if one notes that, for positive detuning, the effective degeneracy between the three states ±, 0 disappears because |0 acquires a higher energy. Degeneracy is restored due to the massive initial occupation of |0 , which lowers its energy via the attractive mean-field interaction in Eq. (14). It thus becomes clear that the effect of a positive detuning is counteracted by an increase of the effective attraction. In particular,  the starting degeneracy recovered with the help of interactions permits the onset of Rabi oscillations known to be essential for the emergence of a ratchet current [21].

Phase dependence and interaction symmetry
In Fig. 5 we show how the the ratchet current depends on the choice of initial driving phase. A change ϕ = 0 → π corresponds to K → −K and thus to a reversal of the current. On the other hand, we note that for ϕ = 1 2 π, 3 2 π the ratchet current is always zero, because parity is restored for those values of ϕ. Figure 5 also shows that the long-time averaged current (but not the instantaneous one) is an even function of the interaction strength.  Figure 3b shows the population n 3 of the third eigenstate of the time-averaged oneparticle density matrix (calculated in 3GP)

Chaotic dynamics
where A µ (t) (µ = 0, ±) are the time-dependent expansion coefficients of the condensate orbital. Comparison with Fig. 3a shows a clear correlation between the occupation of a third state and the suppression or weak reversal of the ratchet current. The dynamics of a macroscopic condensate in a three-level system is that of two coupled non-rigid pendula, which beyond the harmonic regime can be expected to be chaotic. This behavior is confirmed from the inspection of Figs. 4b-f. In Fig. 4b we show a plot ofĪ as a function of ∆ for g = 0.1, obtained with the three different calculation methods described before. For the values of ∆ marked with labels c-f, we show in Figs. 4c-f the corresponding Poincaré cross sections. The one point which shows chaotic behavior is e, which falls in the white-blue region of Fig. 4a. We have just noted (see Fig. 3a) that such a region is correlated with the occupation of a third eigenstate and thus with the possible emergence of chaotic behavior on the coarse-grained time scale. This underlines the connection between three-level occupation, chaotic behavior, and current reduction or reversal. A calculation of the maximal Lyapunov exponent λ M (given by black dots in Fig. 4b) further confirms the chaotic character of this parameter region, as λ M acquires positive finite values there. For more details, see Appendix A. Figure 4b also indicates that in the chaotic regime the convergence in time is slower than in the regular regime, as can be seen by the discrepancy between 3GP calculations performed over 400 or 10 4 driving periods. We emphasize that the current reversal feature in the chaotic region is not just a transient but reflects a lasting behavior; we have checked it up to 10 5 periods in some cases. Finally, we note that all calculations performed in the FGP approximation indicate that, for large K and g close to the self-trapping transition g c (beyond which the system essentially remains in |0 ), the behavior undergoes a significant influence from states beyond the three-state basis 0, ± (not shown). is shown in (e) -same for the maximum of C 11 (red triangles) and the minimum of C 01 (black circles). In (a-c) parameters are the same as in Fig. 2c-d. In (d-e), the parameters are also those of Fig. 2c-d except for the total particle number N , which is here plotted as a variable.

Two-particle correlation functions
It was noted in Ref. [32] and hinted in a previous work [14] that chaotic motion within GP indicates the end of validity of this approximation. This was recently reexamined in Refs. [33,34]. Within the 3LS, we are able to test the validity of GP directly for relatively high particle numbers. We do this by analyzing the relative occupation numbers of the natural orbitals, which are those obtained by diagonalizing the reduced one-particle density matrix, defined in a many-body context as For the two limiting cases of Rabi oscillations and self-trapping, the relative occupation of the condensate orbital (defined here as the most occupied eigenstate of ρ (1) ) stays above 98% during the first 400 driving cycles, which justifies a mean-field treatment. In contrast, for the chaotic regime, we find that the occupation of the condensate orbital decreases abruptly after some 120 cycles (about 2 Rabi periods) and a second natural orbital gets macroscopically occupied -see Fig. 6a. We find a low sensitivity of this behavior to the particle number.
In this chaotic regime, the macroscopic occupation of more than one orbital does not occur via a simple fragmentation of the condensate. This can be seen in the particle number correlation functions C ij ≡ n i n j − n i n j , where i, j refer to the natural orbitals. In a fragmented condensate, C ij ≡ 4C ij /N 2 would be zero (for N → ∞), while for a noon (cat) type state such as e.g.
its absolute value approaches unity, reflecting the creation of macroscopic particle number correlations. Figures 6b-c show the time-evolution of C ij for the same parameter set as in Figs. 2c-d. We see that, with the decrease of the condensate fraction, some particle number correlations increase in a non-negligible way. In order to check the scaling of this behavior with the particle number, we focus on the local maximum of the variance of the condensate occupation (C 00 ) for different particle numbers. Figure 6d shows that this maximum gets shifted to later times for higher N. The value of this maximum stays above 0.2 for all considered particle numbers and increases monotonically for N > 25, as shown in Fig. 6e. This indicates a scaling of the correlations C ij with N 2 , which in turn reveals noon-like (21) behavior, i.e., a dynamics dominated by a few many-body configurations differing by macroscopically large relative particle numbers. Similar observations have been made by Weiss and Teichmann [35], investigating a driven double-well system. Thus we find that in this chaotic regime more than one orbital gets macroscopically occupied. At the same time the particle number correlation functions scale with the square of the total particle number. This indicates a complex dynamics in which the macroscopically large particle number correlations become a substantial feature.

Conclusions
In summary, we have introduced a formalism to account for interactions in ac-driven many-body systems. In the case where only resonant states effectively intervene, we derive a Hamiltonian involving conventional operators which describes the long time dynamics of the driven system. Although we have focused on the case where the undriven states are plane waves, the method we have developed applies to an arbitrary set of states that become resonant in the presence of ac-driving. We have applied the method to the calculation of the coherent ratchet current carried by an asymmetrically driven atomic condensate, and compared it with continuum and truncated mean-field descriptions. We have found a rich dynamical behavior with crossovers from self-trapping to chaotic behavior to regular oscillations. In the latter case, we find a strong departure from the mean-field picture and in particular an appreciable increase of particle number correlations.
where |ψ 0 (t) and |ψ 1 (t) are both solutions of the GP equation. Because of the Hamiltonian structure of the nonlinear Schrödinger equation, it is convenient to allow only perturbations that yield states with the same value of the Hamiltonian [33]. Note that in our effective description the Hamiltonian reflects the system's quasi-energy instead of its energy. Within 3GP, the system's dynamics is given by the 3-tuple A(t) = (A + (t), A 0 (t), A − (t)), reflecting the expansion coefficients of the condensate orbital with respect to the three modes ±, 0. In order to extract the Lyapunov exponent from the system's dynamics, we consider the time evolution of the two initial states A δ+ = N δ (iδ, 1, 0) and A δ− = N δ (0, 1, iδ), where N δ = 1/ √ 1 + δ 2 is the normalization factor. The parameter δ parameterizes the initial distance, √ 2δ, between the trajectories A δ± (t). The Hamiltonian of both the states A δ± has the same value, and for δ = 0, we recover the initial state considered by us. We denote their time-dependent separation by d δ (t) ≡ d 2 (A δ+ (t), A δ− (t)), and use this quantity to estimate the Lyapunov exponent. The limit ψ 1 → ψ 0 in the definition (A.2) translates into the limit δ → 0. Figure A1 shows the time evolution of d δ (t) on a logarithmic scale for chaotic (a), and regular dynamics (b). Two different values of δ are considered. For chaotic dynamics the distance between the two trajectories d δ (t) grows over several orders of magnitude, until it saturates at a value around unity. This saturation occurs because the distance d δ (t) is bounded from above by the value 2, since both states A δ± are normalized to unity. The value of λ M is estimated by fitting an exponential to the rising slope of d δ (t) before saturation is reached. For the case shown we obtain λ M ≈ 0.61 (0.51) for δ = 10 −6 (10 −5 ). The estimated values thus have a relative difference of about 18%, which implies a high relative uncertainty. Nevertheless, the obtained value for λ M is clearly positive, indicating chaotic dynamics.
In the regular regime d δ (t) deviates by less than 10% of its initial value during the depicted time range, for both values of δ. Chaos can be precluded when the global maximum of d δ (t) can be decreased by choosing smaller values for δ.
In all cases, the accuracy of the simulation was checked by propagating both the states A δ± backward in time after t = 3000 T . The simulation is reliable, when d δ (t) comes back to its starting value after t = 6000 T . This check is necessary to ensure that the observed separation of trajectories is not due merely to numerical errors.