Effective Josephson dynamics in resonantly driven Bose-Einstein condensates

We show that the orbital Josephson effect appears in a wide range of driven atomic Bose-Einstein condensed systems, including quantum ratchets, double wells and box potentials. We use three separate numerical methods: Gross-Pitaevskii equation, exact diagonalization of the few-mode problem, and the Multi-Configurational Time-Dependent Hartree for Bosons algorithm. We establish the limits of mean-field and few-mode descriptions, demonstrating that they represent the full many-body dynamics to high accuracy in the weak driving limit. Among other quantum measures, we compute the instantaneous particle current and the occupation of natural orbitals. We explore four separate dynamical regimes, the Rabi limit, chaos, the critical point, and self-trapping; a favorable comparison is found even in the regimes of dynamical instabilities or macroscopic quantum self-trapping. Finally, we present an extension of the (t,t')-formalism to general time-periodic equations of motion, which permits a systematic description of the long-time dynamics of resonantly driven many-body systems, including those relevant to the orbital Josephson effect.


I. INTRODUCTION
In a bosonic Josephson junction (BJJ) two or a few single-particle states are coherently occupied by a macroscopic number of bosons [1]. These systems combine two qualities that make them interesting for experimental and theoretical studies. First, they are a quantum manybody system. Second, the underlying Hilbert space grows only linearly with the total particle number for the case of two single-particle states. Thus BJJs provide a simple framework to study the role of particle interactions in quantum gases. Due to the possibility of conveniently controlling the classical limit via the particle number, these systems are particularly well suited for the study of some aspects of the quantum-classical crossover [2][3][4].
The Josephson effect in Bose-Einstein condensates (BECs) can exist in a variety of qualitatively different forms. In an external Josephson junction [5], the modes involved refer to localized Wannier functions in a doublewell potential. A more robust realization is provided by the so-called internal Josephson effect [6,7], where the states are defined by the different electronic configurations of the gas atoms. Dissipation due to the exchange of non-condensed atoms leads to Ohmic damping of the Josephson oscillations [8]. Bosonic Josephson junctions can also be used for high precision measurement of physical quantities such as temperature [9], weak forces (gravitational [10] or electroweak [11]), or chemical potential differences [12]. Promising scenarios for the creation of macroscopic superposition states [13] or spin-squeezed states [14] are based on BJJs. Furthermore, in a recent experiment at NIST [15], a BJJ in an atom circuit was used to measure rotation.
Recently, a third kind of Josephson effect in BECs, not The resulting dynamics is that of a bosonic Josephson junction. In this schematic picture, the horizontal difference represents a chemical potential difference among the modes and the connecting line represents the coupling between them.
classifiable as external or internal, has been identified. In the orbital Josephson effect (OJE) the single-particle modes have identical internal structure and center-ofmass wave functions with strongly overlapping densities (see Fig. 1). The OJE can be realized when a timeperiodic driving potential induces a coupling among the unperturbed Floquet states of the system. Previous work on the OJE [16] has focused on the quantum ratchet on a ring potential because the OJE was originally identified in that context. This kind of system is realizable in current BEC setups. In fact an orbital Josephson junction can be obtained in a variety of static and driving potentials. The present work aims at presenting the OJE in a very general context. Starting with two very simple illustrative examples, a double well [5] and a box potential [17], we show that the OJE is a general concept that can be studied in a variety of existing experimental BEC systems. Then we explore different driving potentials in ring traps, inspired by recent experimental suc-cesses [15,[18][19][20][21][22][23][24].
For some selected cases, our theoretical model based on an effective few-mode Hamiltonian is compared to computationally demanding simulations from first principles performed with the Multi-Configurational Time-Dependent Hartree for Bosons (MCTDHB) algorithm [25][26][27]. For the present purposes, MCTDHB can be regarded as a method that enables a full manybody (FMB) study of the interacting boson problem. Thus we focus on the comparison between the truncated description of the orbital Josephson effect and a full many-body calculation that operates in a larger one-atom Hilbert space. Importantly, in both these pictures we work well-beyond the Gross-Pitaevskii (GP) approximation [28,29], which is also considered in some cases. We find that the MCTDHB results agree well with the truncated BJJ description in a variety of dynamical regimes. This includes regimes that show instabilities in the semiclassical limit such as chaos and unstable fixed points which are poorly described within the GP approximation. Interestingly, we find that the regime of macroscopic quantum self-trapping preserves its character in a non-truncated, full many-body description. This contrasts with the case of an undriven double-well potential (external Josephson effect), where self-trapping seems to be fragile within a multi-mode picture [30].
A major contribution of this paper is the generalization of the so-called (t, t ′ )-formalism [31,32] to an almost arbitrary type of equation of motion including the nonlinear Schrödinger equation. The (t, t ′ )-formalism has shown to be an efficient mathematical tool for the study of time-periodic [32,33] and aperiodic [34] linear systems. Here we employ the extended (t, t ′ )-formalism to analyze the dynamics of field operators in interacting many-body systems.
This paper is arranged as follows. Section II presents general two-mode orbital Josephson junctions, without specification of the trap geometry or the driving potential. It also includes a rather general presentation of the (t, t ′ )-formalism for field operators. Section III introduces two illustrative examples of the orbital Josephson effect: a rocked double well and a box potential perturbed with a modulated lattice. In Section IV the previously introduced concepts are used to explore in greater depth the case of a BEC in a ring trap subject to a ratchet potential. Up to three different driving potentials are discussed, each yielding a three-mode orbital Josephson system. Finally in Section V, we compare the predictions from the effective description in terms of an orbital Josephson junction with full many-body simulations using MCTDHB, focusing on one of the driving potentials introduced in Section IV. In Appendix A the extended (t, t ′ )-formalism is derived in great generality. Appendix B addresses technical details of the numerical MCTDHB study.

II. GENERAL CONCEPTS
In this section, we present a recipe for the realization of an orbital Josephson junction. We will keep our approach as general as possible. However, in Section IV, we will see that a variety of orbital Josephson junctions also exist which cannot be considered a specific realization of the type presented here.
We consider an initially static BEC, trapped in an arbitrary geometry, given by some static external potential V trap (r). At time t = 0, a time-periodic driving potential with periodicity T , is switched on. The dynamics of this many-body system is determined by the Heisenberg equation of motion for bosonic field operators, where time and spatial coordinates are made dimensionless. This can be realized by the choice of an appropriate length scale x 0 and then expressing all lengths, energies, frequencies, and times in units of x 0 , 2 /M a x 2 0 , /M a x 2 0 , and M a x 2 0 / , respectively, where M a is the atom mass. Furthermore, we set = 1, such that energies and frequencies have the same dimensions. The single-particle is a sum of kinetic energy and external potential. The initial field operator ψ(r) ≡ψ(r, 0) annihilates a bosonic particle at point r and satisfies the standard bosonic commutation relations [ψ(r),ψ † (r ′ )] = δ(r − r ′ ) and [ψ(r),ψ(r ′ )] = 0, which are preserved in time. The contact interaction strength λ is proportional to the s-wave scattering length of the gas atoms.
Our goal is to develop an effectively time-independent description that involves only two single-particle states. We do this by using an extended version of the (t, t ′ )formalism to describe the dynamics of the field operatorsψ(r, t) in second quantization [16]. The (t, t ′ )formalism has proven to be a powerful tool for solving the Schrödinger equation with time-dependent Hamiltonians [31,32]. The equation of motion for the field operators (2) reads in the extended (t, t ′ )-formalism where t ′ acts as an additional parameter which the field operators depend on. The single-particle part H(r, t ′ ) ≡ H 0 (r) + V (r, t ′ ) is the sum of kinetic energy, trapping potential, and a driving potential. Equation (3), together with the initial conditions, ψ(r, t ′ ; 0) ≡ψ(r) for all r and t ′ define a unique solution forψ(r, t ′ ; t), from which the physically relevant field operatorψ(r, t) can be obtained viaψ (r, t) =ψ(r, t ′ ; t)| t ′ =t .
Note that, the periodicity in t ′ , which is trivially imposed by the initial condition, Eq. (4), holds for all times, since the equation of motion (3) does not break this symmetry. Furthermore, under these initial conditions the initial commutation relations are also preserved, which read We assume the static single particle part H 0 (r) to be the predominant term in the equation of motion (2). This condition is met for systems with weak particle interactions λ, and a weak overall amplitude of the driving potential V (r, t). Furthermore, this assumptions implies, that the low energy eigenstates of the undriven (t < 0) system are condensates, with a negligible amount of depletion [35] and with the condensate orbitals given by eigenstates of H 0 (r).
Within the (t, t ′ )-formalism, this means that the predominant part is given by the unperturbed single-particle Floquet operator suggesting one ought to use a transformed representation of Eq. (3) that involves creation and annihilation operators with respect to the unperturbed Floquet states (eigenmodes of H 0 ).
A. Floquet representation in second-quantized form Therefore, we introduce the representation given bŷ where φ k (r) denotes an eigenmode of H 0 (r). The corresponding unperturbed single-particle quasienergy is given by E 0 km = ε 0 k − ωm, where ε 0 k is the energy of φ k (r). Quasienergies can be considered the time analog to quasi-momentum from Bloch theory. The inverse transformation to Eq. (8) readŝ where the sum over m runs over all integer numbers and k indicates the eigenstates H 0 (r). In this representation, the commutation relations (6) become In order to obtain an equation of motion forâ km (t) we enter these transformations into Eq. (3). For the rest of this section, we focus on driving of the form where the phase ϕ reflects the value of the modulation at the time point (t = 0) when the driving potential is switched on. Driving of this type is particularly relevant for experiments since it can be obtained naturally via inertial forces from shaking [36], or via intensity modulations [37]. To a first approximation, the latter possibility would cause an additional energy shift of the eigenmodes of H 0 (r). We note that the switching conditions can be tuned to generate the effective gauge factor e iϕ in Eq. (18). This generalizes a similar result derived in Ref. [38] within first quantization.
Note that, in principle, the concepts presented here can be generalized to driving potentials which include higher harmonics. This yields the equation of motion where the last sum is meant to run over all values of the indices {k 1 , k 2 , k 3 }, and {m 1 , m 2 , m 3 }. It involves the matrix elements of the driving potential, as well as the two-particle matrix elements for contact interaction We assume the system to be initially condensed in an eigenmode of the undriven single-particle part of the Hamiltonian H 0 (r), which will be indicated by the index k = 1. This is consistent with the previous assumption of weak particle interactions. Precisely, the condition λN ≪ ∆E has to be met in order to obtain a sufficiently condensed system; where ∆E is the energy difference between the initial mode and the neighboring modes.

B. Two-level description in a many-body framework
In order to obtain a dynamics that is governed by exactly two modes, we impose the following additional conditions on the driving field: First, there exists a further orbital, indexed as k = 2, to which the driving potential couples, i.e. V 1 2 = 0. Second, the driving frequency ω is near-resonant, i.e. close to the energy difference ω 0 ≡ ε 0 1 − ε 0 2 . This means that the driving frequency can be expressed as ω = ω 0 + ∆, where ∆ is a possible detuning from exact resonance. It should satisfy ∆ < 1 2 |ω 0 − ω 1 |, where ω 1 is the resonance closest to ω 0 . Given these two conditions, we say that the driving potential V (r, t) induces a resonant coupling between the modes |1 and |2 . As a final assumption, the driving shall not induce any resonant couplings between |1 or |2 and a third state. In practice, this just means any such couplings should be vanishingly small on the time scale of our simulations or the time scales for applying the driving in the case of BEC experiments.
For the case of vanishing particle interactions (λ = 0), it is known that these requirements justify a two-level description [39]. A typical approach to extract the longtime dynamics is given by the rotating wave approximation [39], which, within a single particle description, yields an equivalent, effectively time-independent, twolevel Hamiltonian. In this section, we perform the derivations within the (t, t ′ )-formalism instead of using the rotating-wave approximation. The reason is that in sections IV and V these concepts will be extended to driving potentials that contain a second harmonic in the time-like part and induce a coupling among three single-particle modes. Such systems are beyond the standard use of rotating-wave approximation, while the (t, t ′ )-formalism can be conveniently applied to more general types of driving.
The truncation to a two-level system can also be justified within the (t, t ′ )-formalism. When the driving frequency is tuned to exact resonance (or sufficiently close to it), the two modes |1, 0 and |2, 1 become degenerate or nearly degenerate with respect to the predominant Hamiltonian H 0 . Based on the assumptions above, we truncate the system of equations (13), in such a way that only the operatorsâ 10 ,â 21 , and their Hermitian adjoint participate. This yields an effective two-level Hamiltonian describing the many-body dynamics of a resonantly driven Bose system, and gives rise to the orbital Josephson effect.
In general, particle interactions induce a coupling to further modes and could invalidate the two-level description. Nevertheless, the single-particle part imposes a quasienergy difference between both modes |1 and |2 and further modes, which prevent the system from accessing these further modes. Therefore, it is reasonable to assume that the two-level description remains valid for weakly interacting particles.
Within this truncated space of resonant levels we will omit the second subindex m in the operatorsâ With a truncation of the set of modes to 1 and 2, the commutation relations (10) become those of standard bosonic creation and annihilation operators: In this truncated picture, Eq. (13) becomes wheren i =â † iâ i , and ϑ = ie iϕ V 1 2 . Furthermore, we have used the abbreviations w 1 ≡ W 11,11 , w 2 ≡ W 22,22 , and w 12 ≡ W 12,12 . These are the key interaction terms, corresponding to interactions within mode |1 , within mode |2 , and between modes |1 and |2 . For a contact interaction, the indices in the two-particle matrix elements are symmetric under a swap of the first two or the last two indices. For that reason, w 12 = w 21 has four equivalent permutations of indices, which gives rise to the prefactor 2. Terms like W 11,12 , W 22,21 , or other terms with permuted indices do not appear, because of the first Kronecker-delta in the interaction term in Eq. (13) makes these two-particle matrix elements vanish.
Since, within the truncated picture, standard commutation relations (16) are reobtained, the equations of motion (17) can be regarded as Heisenberg equations of motion coming from the truncated two-level-system Hamiltonian Schematic illustration of the discussed example systems: (a) driven two-mode system, (b) box-potential with driving. For illustrative purposes, the driven two-mode system (a) is here depicted as a double-well trap. The considerations in Section III A may also refer to an internal Josephson system. with andN is the total particle number operator. For systems with conserved particle number, the last term in Hamiltonian (18) has no effect on the dynamics and can be omitted. The constant term is only needed for relative comparisons of total energy, e.g. for systems with varying number of particles. When it is omitted, Hamiltonian (18) is that of a BJJ, with the peculiarity of a modedependent interaction strength. Note that this few-mode picture does not take into account initial depletion of the condensate. Therefore, we expect this effective description to work well in the limit of high particle numbers and small interaction strength λ, i.e., the mean-field limit λN = const., λ → 0, and N → ∞. If more than two modes are involved, it is still possible to incorporate the particle interactions in the same manner. However, in general, the particle interactions will include terms that induce a mixture among the modes. Furthermore, the mathematical expressions for the interaction term become lengthy when the participation of a third mode is incorporated in such a general manner. Under certain assumptions on the unperturbed Floquet states, it is still possible to obtain a similarly simple truncated description, cf. Section IV.

III. EXAMPLE SYSTEMS
In this section, we present two examples of specific setups that meet the conditions listed in the previous section and hence represent a realization of an OJE. A schematic illustration of the systems discussed here is shown in Fig. 2

.
A. Minimal example: driven two-mode system Our first explicit example of an OJE, is set up by a conventional (external or internal) BJJ. For convenience, we label the modes of this conventional BJJ by L and R. Our starting point is a driven bosonic two-mode Hamil-tonian (20) whereâ † i (â i ) creates (annihilates) a particle in one of the modes i = L, R and obey the usual bosonic commutation (20) can be considered a time-dependent driven version of the Lipkin-Meshkov-Glick model [40].
As for the general case discussed above, the driving is switched on at time t = 0. The eigenmodes of the undriven single-particle part of Hamiltonian (20) are given by the states, where the signs refer to positive (+) and negative (−) parity, if they are regarded as localized Wannier states. Their energy difference is ε − − ε + = J. Note that ε + < ε − . For weak driving strengths K, and when the system is initially condensed in one of the eigenmodes |± , we can apply the recipe presented in Section II and obtain an effectively time-independent description of the system. In this minimal example, the two-particle matrix elements have the simple form The remaining elements of W either can either be obtained by permutation of the indices or are zero. In any case, the above displayed elements are all we need to apply the prescription (19). Hence, we obtain an effective inversion of the interaction strength: u ± = − 1 2 . The matrix element of the driving is V + − = 1 2 K, from which follows the coupling element ϑ = i 2 Ke iϕ . With ∆ = ω − J, the effective two-level Hamiltonian (18) is fully specified.
In this minimal example, we already started with a two-mode system. However, the corresponding effective description as an orbital Josephson system, given by Hamiltonian (18), represents a truncated picture with respect to the unperturbed Floquet states. Figure 3 shows a comparison of the dynamics governed by the original many-body Hamiltonian (20) (dashed curve) and the time-independent effective Hamiltonian (18) (solid curve). The full dynamics features a fast wiggling onset whose frequency is similar to that of the external driving potential. However, the effective description does not reflect this wiggling. Rather, it refers to the long-term behavior of the system. Furthermore, there is also a discrepancy in the long-term behavior, which becomes more evident for stronger particle interactions. The precision of this effective description will be discussed in more detail in Section V.

B. Box potential
A second simple example is a BEC with a single spatial degree of freedom, trapped in a box potential with length 2π: Such a trap with a flat potential inside the confining walls has been realized in a recent experiment [17]. The eigenfunctions of the single-particle Hamiltonian are standing waves where k is a positive integer, and the related energy eigenvalue is ε k = k 2 /8. In this basis, the two-particle matrix elements for contact interaction are given by For the spatial part of the driving potential, we choose the lattice potential where κ > 0 is an integer valued parameter. The matrix elements, with respect to the standing waves (24) read This means that the potential V (x) induces a coupling among those modes that have momentum difference (k − l) of magnitude |κ|. A condensate that is initially prepared in a certain condensate state |k 0 , will couple either to the mode |k 0 + κ or to |k 0 − κ , if any. A simultaneous coupling to both modes is ruled out, because the driving frequency ω cannot match both energy spacings. At this point, it becomes clear why a harmonic trap is not suited for the realization of an OJE. For a harmonic trap, the spacing of neighboring energy levels is constant. Consequently, a resonant driving frequency would, in general, induce a subsequent coupling to all eigenmodes of the system. However, in actual experiments traps are only locally harmonic. Optical traps for example are built on Gaussians [41], and harmonic plus quartic traps have also been used [42]. So in general one does not require a box potential, and a series of other extant experimental systems are viable for the OJE. The example here can be adapted straightforwardly to those systems.
As in the previous example (Subsection III A), for the box potential, the interaction strength is effectively attractive, with u k = −1 independent of the specific modes k that form the orbital Josephson system.

IV. ORBITAL JOSEPHSON SYSTEMS WITH THREE MODES
In the following sections, we focus on a BEC in a ring trap, and consider three driving potentials that yield a bosonic Josephson system consisting of three angular momentum modes. Three illustrative examples are considered, suggesting that a variety of possible driving potentials exists which yield an OJE with this type of trapping potential. The ring trap [15,[18][19][20][21][22][23][24] shall be such that the motion of the atoms is effectively frozen along the radial degree of freedom. Consequently, the dynamics can be modeled as a one-dimensional system of length 2π with periodic boundary conditions. This means that we choose the radius R of the ring as our natural length scale (x 0 = R) and express all times, energies and frequencies accordingly -cf. Section II.
Initially, the potential energy along the single degree of freedom is flat. Therefore, the eigenstates of the singleparticle part of the Hamiltonian are plane waves (angular momentum eigenstates), characterized by a wave vector k and corresponding energy eigenvalue ε k = 1 2 k 2 . Initially, the Bose gas is assumed to be fully condensed in the k = 0 mode. At time t = 0, a weak driving potential is switched on, inducing coupling to modes with k = 0. Optimally, the applied potential provides the possibility to control the coupling strength to each of these modes separately. We consider three driving potentials, denoted as cases (a), (b) and (c), that combine these qualities. These are Case (a) has been realized experimentally on an extended optical lattice [37]. The second case (b) is the minimal driving potential that leads to the same first-and second-order processes as those of case (a). Finally, for case (c) the contribution of second-order processes can be neglected, which makes this driving potential particularly suitable to study the interplay between particle interactions and the truncated picture. For this case, the underlying processes are analogous to Bragg reflection.
In the adiabatic limit (ω → 0), each of the contributing waves could be regarded as a conveyor belt, and in classical systems they yield a Brownian surfer [43] if dissipation is added. A further important motivation of the third potential is that it can be straightforwardly implemented in experiments, because it results from the force of inertia that arises from an orbiting motion of the center of the ring along an ellipse. Common parameters of the potentials V a , V b , and V c are the wave vector κ of the first harmonic, which controls the angular momentum of the modes to which the driving couples. The driving frequency ω should be close to a resonant frequency according to the modes to which the driving couples. For V a and V b , the overall amplitude is controlled via the parameter K, and for V c we can define an overall amplitude as K = 1 2 (K + + K − ). All three of the driving potentials in Eq. 29 represent a realization of a ratchet potential. That is, they all break parity and time inversion symmetry [44,45]. For V a , both symmetries can be broken separately. The time-inversion symmetry is broken for β = 0 and ϑ = π/2, 3π/2 and parity is broken for α = 0 and ϕ = π/2, 3π/2. This potential was studied in Ref. [46] for the case ϕ = ϑ = 0. For V b , both inversion symmetries hold for γ = 0, while for γ = 0, time-inversion is broken with ϑ = π/2, 3π/2 and parity is broken for ϕ = π/2, 3π/2. Finally, V c simultaneously breaks both symmetries for K + − K − = 0. The phases ϕ ± have no control over the symmetry properties of V c .
All potentials have in common that they induce an OJE involving three modes that can be fully characterized via the angular momentum they carry. These modes are the initial center mode |0 , and two further modes with equal kinetic energy and opposed momenta, denoted as |± , cf. Fig. 4d. The coupling strengths for each of the constituent junctions can be tuned separately. Note that the systems presented in the following cannot be considered a specific realization of the type of orbital Josephson junctions considered in Section II. They differ in the number of modes that participate in the Josephson system. Furthermore, the driving potential is, in general, Fourier representation of three driving potentials (30) that yield an effective description given by Hamiltonian (31). The label of each panel indicates to which of the driving potentials listed in Eqs. (29) it refers to. The size of the points corresponds of the magnitude of the matrix element 00|V |km . The lower halves of panels (a-c) are omitted. Since the potentials (29) are real valued, all representation are point symmetric with respect to the origin. The circles highlight the principal participating modes, i.e. those that are intersected by the parabola of kinetic energy. The wave vector κ of the driving is κ = 1 in (a,b), and κ = 2 in (c). Panel (d) schematically shows the resulting three-level dynamics. not given by a product of a purely spatial and a purely time-like part; and in some cases (a,b), the contribution of a second harmonic in both parts become important. However, the driving potentials considered here are of course not the only possibilities to obtain an OJE with three modes.
As in previous sections, we perform our considerations in the (t, t ′ )-formalism. The matrix elements of the driving potentials (29) in the discrete representation (8) read (30c) Figure 4 shows a graphical representation of these matrix elements.
For noninteracting particles V a induces a coupling between the initial mode |k m = |0 0 , and the modes | ± 2κ 2 , where the effective coupling parameters can be derived with a perturbative calculation [46]. The second driving potential V b is equivalent to V a within a perturbational approximation that involves processes up to second order in K. Both perturbations V a and V b feature those matrix elements that are involved in the perturbational calculus presented in Ref. [46]. The potential V c consists of two counter-propagating sinusoidal waves. Each of these waves induces a coupling between the initial mode |0 0 , and | ± κ 1 . In contrast to the potentials V a and V b , these couplings can be directly obtained from a first order calculation.
Note that, since the calculations are performed within the (t, t ′ )-framework, each mode is characterized by two indices: the quantum number k indicating the angular momentum and the index m. As for the two-level system, within the truncated Hilbert space (0, ±) the second number m is determined by the angular momentum and will be omitted.
The dynamics resulting from each of the three drivings in Eq. (29) is governed by the effective three-level-system Hamiltonian [16] given that the total particle number is conserved. The index ν in the last sum takes values ± and 0. Here, Γ ± and µ are effective parameters that are determined via the parameters of the driving field. Their dependence for each driving potential is listed in Table I. In the case that either Γ + or Γ − vanishes, a two-level system is obtained and the resulting HamiltonianĤ 3LS becomes equivalent to Hamiltonian (18), but with a mode-independent interaction term. ∆   TABLE I. Effective parameters for the truncated picture in dependence of the driving field parameters. The left column indicates the driving potential V a/b/c . The listed functional dependences refer to the main resonance given by ω0 = 1 2 κ 2 for all potentials and ∆ denotes a possible detuning. Note that for V a and V b further resonant frequencies exist with a different functional dependence of Γ± and µ on the parameters [46].

V. NUMERICAL STUDY OF DYNAMICAL REGIMES
An effective description in terms of a BJJ is not guaranteed to always withstand a comparison with a more exact numerical simulation of the full system. For example, in the case of a BEC in a one-dimensional double well the related two-mode description, referring to localized Wannier functions, has been shown to be invalid near and in the regime where the two-mode model predicts self-trapping [30]. The reasons for that discrepancy are not clear since nominal criteria for the validity of the two-mode description were met in that study. In the following, we present a numerical check, where some results obtained with the truncated picture of an OJE are compared to the full many-body dynamics.
Because of its rich variety of dynamical regimes, we choose the ratchet system, presented in Section IV, for the numerical study. We focus on the dynamics governed by the driving potential V c -see Eq. (29c). This choice rules out possible effects of the intermediate modes that are involved in the higher-order perturbative calculation of the transition amplitudes Γ ± for the potentials V a and V b [46], and allows us to focus on the interplay between interaction and the truncated picture.
We have numerically solved the full many-body (FMB) dynamics involving a large number of modes, using MCTDHB [27] for various particle numbers and interaction strengths. The approximate three-level sys-tem (3LS), determined by Eq. (31), is solved numerically via exact diagonalization. We observe the occurrence of at least three qualitatively different dynamical regimes: (1) Rabi oscillations, between the initial mode |0 and |a ∝ Γ + |+ + Γ − |− , are expected to occur when the interaction term is weak compared to the driving strength, i.e., g/4π ≪ K, where g = λ(N − 1) is the mean-field interaction strength. The corresponding Rabi frequency is given by Ω R = 2 |Γ + | 2 + |Γ − | 2 + µ 2 /4. Its inverse T R = 2π/Ω R serves as a natural time scale in the truncated picture. (2) Chaotic dynamics are expected for intermediate interaction strength. (3) Self-trapping occurs when the particle interactions dominate over the driving strength g/4π ≫ K. The critical interaction strength for the occurrence of self-trapping is [16,47] g c = 8π max(|Γ + |, |Γ − |).
For our simulations, we choose κ = 1, which allows the use of a small single-particle basis for our FMB simulations (based on MCTDHB). The single-particle basis, used for the simulation presented here, is given by angular momentum modes, ranging from k = −4 to k = 4. Further details on the numerical method are presented in Appendix B. We choose a small overall driving amplitude K = 0.1, the regime for which the 3LS approximation should be a good one. The driving frequency is chosen to be exactly on resonance, i.e. ω = 0.5 for κ = 1 (∆ = 0). For the amplitudes of the constituent plane waves of the driving we choose K + = 0.07, K − = 0.03, and the phases ϕ ± = 0. We consider four different values of g, referring to qualitatively different dynamics. Those are (a) Rabi oscillations (g = 0.1), (b) chaos (g = 0.5), (c) dynamics near the critical point (g = 0.9), and (c) self-trapped dynamics (g = 1.5). The chaotic regime is identified by the calculation of the maximal Lyapunov exponent, which is zero for all other considered cases. The critical meanfield interaction strength for the self-trapping transition is, according to Eq. (32), g c ≃ 0.88. Particle numbers up to N = 40 are considered in the following analysis.
For a realistic experimental setup of a ring trap, loaded with 23 Na atoms (mass M Na = 22.99 amu = 3.8 × 10 −26 kg), with major radius R ≃ 10µm [21], the relevant time scale is t 0 ≃ 36ms. This means that the energy difference between the zero angular momentum mode and the first nonzero angular momentum mode corresponds to a frequency of about 900Hz. A driving potential that couples between the initial state and the first excited state would have a resonance at this frequency.
As was mentioned in Section II, in our 3LS calculations we do not include initial depletion of the condensate. Within FMB, we can analyze possible effects that originate from the initial depletion. The FMB simulations, presented in the following, take the ground state of the initial static Hamiltonian as an initial state, including particle interactions. Consequently, the initial state of the full system is depleted to some extent. This is a realistic choice for an initial state, since it can be experimentally obtained by the cooling of an interact- ing Bose gas. For N = 40 particles, the initial depletion for the considered interaction strengths stays below 1.4 × 10 −3 , which is reached for g = 1.5. We note that typical BECs in harmonic traps have depletion on this order or smaller [48]. Figures 5 and 6 show a comparison of FMB with 3LS over a time range of 50 driving cycles, corresponding to more than seven Rabi periods (see upper horizontal scale in Figs. 5 and 6.). Figure 5 shows the instantaneous mean current per particle I(t) ≡ −i dx ψ † ∂ xψ /N . Interestingly, for the unstable cases (g = 0.5 and g = 0.9), the discrepancies become quite large at certain times, but both curves revert to a good agreement at t ≃ 35T . For example, for g = 0.5 (Fig. 5b) the two descriptions deviate by about 40% after 15T . The largest discrepancies are observed in the self-trapped regime, where the relative difference amounts to 125% at the first maximum (t ≃ 1.2T ). Nevertheless, 3LS correctly reflects the drop of I(t) by at least one order of magnitude as compared to the other cases. Figure 6 shows the normalized occupation numbers of the natural orbitals, given by the eigenstates of the singleparticle density matrix (SPDM). The highest occupied orbital is the condensate orbital. For a normalized condensate occupation close to 1, the system can be well described by the GP equation, because this approximation can be derived by making the ansatz of a fully condensed system [28,29,49].
As can be seen in Fig. 6, for g = 0.1 and g = 1.5, the system remains condensed during the depicted time range. For the case of weak particle interactions, this is not surprising, since in the limit of zero particle interactions the system remains fully condensed for all times. For the values of g in between these extreme cases (Fig. 6bc), the dynamics is more complex and many-body effects become important. To be specific, a second orbital becomes macroscopically occupied by an amount of up to 30%. This indicates that in these cases, a mean-field treatment would fail and many-body considerations become indispensable. This observation is consistent with previous works, in which the validity of the GP equation for the case of chaotic dynamics was studied [16,[50][51][52].
One quantification of the relative disagreement between two curves f 1 (t) and f 2 (t) is given by the time- The disagreement,σ, of the condensate occupations, cf. Fig. 6, are 0.1%, 3%, 4%, and 1% for g = 0.1, 0.5, 0.9, and 1.5 respectively. Figure 7 shows the normalized occupation number of the outer modes i.e., the total occupation those modes that lie outside the three-level Hilbert space. Two different particle numbers (N = 20 and N = 40) and the mean-field approximation are considered for the same values of g as in Figs. 5 and 6. A general trend is that the participation of outer modes decreases for larger particle numbers. Furthermore, we see that the GP equation can also provide some information about the relative participation of outer modes. A clear trend of an increase of the participation of outer modes in time cannot be distinguished for the depicted cases. It seems to reach a saturated value after less than one Rabi period with a fluctuation onset, except for the Rabi regime, where the occupation of outer modes shows an oscillation. In the self-trapping regime, the GP equation clearly underestimates the contribution of outer modes to the dynamics. The participation n ′ (t) remains below 3 × 10 −3 for all cases depicted in Fig. 7.
We have performed additional FMB simulations, with a fully condensed initial state in the k = 0 mode, in order to see whether some of the discrepancies observed here between 3LS and the full dynamics can be assigned to the initial depletion, present in the initial state of the full system. For all interaction strengths considered, the discrepancies between the effective 3LS description and the full dynamics exceed the discrepancies observed here between the two versions of the FMB simulations (with and without initial depletion). Precisely, the relative disagreement between the condensate occupations of the initially depleted and the initially fully condensed system isσ = 3 × 10 −4 for g = 0.1 and is below 0.7% for the other three considered interaction strengths. These results suggests that initial depletion of the condensate plays a minor role in the comparison presented here. Figure 8 shows a comparison between 3LS and FMB with a weaker driving (halved values of K ± or K). The ratchet current is given in Fig. 8a, and Fig. 8b compares the occupation of the natural orbitals. The interaction strength is g = 0.25, such that the dynamics is completely analogous to the chaotic case shown in Figs 6b, but with a larger effective time scale. Note that the time range in Fig. 8 has doubled, compared to Fig. 6. It shows that both descriptions approach each other for smaller values of the driving amplitude. For example, the difference in the ratchet current between 3LS and FMB after three Rabi-periods amounts 0.1 for K = 0.1, while for K = 0.05 it has decreased to 0.05. Similarly, the disagreement, σ, of the normalized condensate occupation has halved to 1.5%. Thus the 3LS improves in accuracy for weaker driving amplitude.

VI. CONCLUSIONS
We have shown that the orbital Josephson effect is a general concept that can be realized in a variety of driven condensate setups. The OJE manifests itself when singleparticle states occupying the same region of space are coherently populated by a macroscopic number of resonantly driven bosons. It is distinct from the external or internal Josephson effect, which do not require external driving. We have listed the main characteristics that must be met by the trap and driving potentials to realize the OJE. We have discussed several trap geometries and for one of them we have considered three different driving potentials. Realistic experimental parameters and realizations were provided for these geometries and cases.
In two selected cases, the truncated Josephson description has been compared with the full many-body dynamics that encompasses a larger one-atom Hilbert space. In some cases, the effective, few-mode description is good approximation of the systems dynamics, and becomes increasingly better for weaker driving amplitudes. Even in cases where the discrepancies are large, the Josephson description still serves as an efficient way to capture the qualitatively different dynamical regimes and to predict in which parameter regimes they may be found. For some parameter regimes, many-body considerations (beyond mean-field) are necessary, even if depletion can be initially neglected. We have found that the effective (truncated) description can be correct in those cases as well. A few-mode description remains valid even near instabilities as long as the driving is sufficiently weak. Finally, we would like to remark that the regime of macroscopic quantum self-trapping is not an artifact of the effective description but is preserved within a full many-body calculation over simulation time scales.
We have extended the (t, t ′ )-formalism to a rather general case, which permits a convenient description of the coarse-grained, long-time dynamics of resonantly driven many-body systems. The extension can be applied, e.g., to the Heisenberg equation of motion of field operators in an arbitrary representation as well as to the timedependent nonlinear Schrödinger equation.

General case
Our goal is to map systems with an underlying timeperiodic equation of motion, with period T , to systems without explicit time dependence. The dynamics of a physical system is given by a trajectory in the phase space P of the system. For now, we do not need to make any further assumptions on P. It can have finite or infinite dimensions, and does not have to be a Hilbert space, i.e. no scalar product or norm needs to be defined on P. Further below, we will pay special attention to the case of square-integrable wave functions.
The state of the system is given by a vector v ∈ P. Its dynamics is fixed by the initial value problem: As mentioned above, we restrict our considerations to time-periodic systems: In any other sense, F is arbitrary. It can be nonlinear, discontinuous in v or in t, and of course it is allowed to have no explicit time dependence at all, in which case the period T can be chosen freely. Now we consider an arbitrary generalized loop in P, defined as v(t ′ ) ∈ P for all t ′ ∈ [0, T ] andv(0) =v(T ). (A3) We allowv(t ′ ) to be contracted, e.g., to a point (v(t ′ ) ≡ v 0 ) or to a line. Hence,v(t ′ ), strictly speaking, does not necessarily form a loop. This is why we assign the term generalized loop. This generalized loop may evolve in time, according to the initial value problem Note that F is a more general object than F , since it contains a partial derivative ofv with respect to t ′ . The dynamics of this generalized loop describes a surface with the structure of a generalized tube in the phase space of the considered system. Having our goal in mind, We are looking for an equation of motion that determines a parametrized family of curves through this surface. This family, with parameter τ , is given by Sincev(t ′ , t) and F are periodic in t ′ , it is sufficient to restrict ourselves to τ ∈ [0, T ]. An equation of motion for v ′ τ (t) can be obtained by deriving both sides of Eq. (A5a) with respect to t. The first argument ofv on the right hand side of the definition (A5) is itself a function of t. Therefore, we have to derive partially with respect to both arguments ofv and apply the chain rule for the case of the first argument. This yields where, t ′ τ is always meant to be a function of t, although it is not explicitly expressed. We have used Eq. (A4) to express ∂v(t ′ , t)/∂t as well as the identity ∂t ′ τ (t)/∂t = 1. The modulo operation within the argument of F (v, ·) does not have any effect due to the periodicity of F (v, ·). Hence, we obtain the equation of motion for v ′ τ (t), given by The initial conditions follow from Eq. (A4): v ′ τ (0) = v 0 . This means that for each τ ∈ [0, T ], v ′ τ (t) is the solution to the initial value problem In particular, the solution to the initial value problem, given in Eq. (A1), is obtained by Note that the t ′ -dependence ofv is originally restricted to the interval [0, T ], but as we imposed periodic boundary conditions in Eq. (A3), we can directly extendv(t ′ ) to be defined on the entire real axis.
The extended (t, t ′ )-formalism, derived here, is quite general. The underlying equation of motion could be a linear [31,32] or a nonlinear differential equation; it could be the Heisenberg equation of motion for field operators [16], the nonlinear Schrödinger equation (see Section A 2), or any other equation of motion from a very different context, including systems with dissipation.

Nonlinear Schrödinger equation
As mentioned above, these concepts can be applied to the dynamics of square integrable wave functions. In the case that the underlying equation of motion is linear, the standard (t, t ′ )-formalism [32] is obtained. Here, we will focus on the nonlinear Schrödinger equation, i∂ t ψ(r, t) = H 0 (r, t)ψ(r, t) + g|ψ(r, t)| 2 ψ(r, t), (A10) where H 0 (r, t) = V (r, t) − 1 2 ∇ 2 is the operator associated with the linear part of this equation. We consider it to be time periodic, with period T : H 0 (r, t + T ) = H 0 (r, t). The wave function ψ shall be normalized as dr |ψ(r)| 2 = 1, which can be shown to be conserved under equation of motion (A10).
In this case the phase space is the Hilbert space of square integrable complex-valued wave-functions ψ(x). Applying the recipe as outlined above in Eqs. A4-A1, we obtain the nonlinear Schrödinger equation within the (t, t ′ )-formalism: i∂ t ψ(r, t ′ ; t) = H 0 (r, t ′ ) − i∂ t ′ ψ(r, t ′ ; t) + g|ψ(r, t ′ ; t)| 2 ψ(r, t ′ ; t), (A11) with the normalization drdt |ψ(r, t ′ ; t)| 2 = T . This normalization is a direct consequence of the fact that the physically relevant wave function is obtained as ψ(r, t) = ψ(r, t, t). Equation A11 has the form of a nonlinear Schrödinger equation, with the linear part being static and the related wave function lives in an extended space given by the tensorial product of conventional Hilbert space and time-periodic functions, which was introduced by Sambe [53]. This implies that Eq. (A11) can formally be derived from the Hamilton functional (A12) Note that, due to the integral over over time (t ′ ), the functional H ′ GP has the unit of an action instead of energy. However, the equation of motion (A11) can be obtained from H ′ GP by applying Hamilton's equations. The minima of this Hamiltonian are stationary states of the equation of motion (A11). These stationary solutions yield states, whose associated physically relevant solutions (obtained by t ′ = t) are T -periodic up to a global phase factor e −iεT . They are the analog to Floquet states from the linear Schrödinger equation, and are named nonlinear Floquet states [3,54]. Accordingly, ε is the corresponding quasienergy.

Appendix B: Convergence study
Within MCTDHB the maximal number of natural orbitals that are taken into account is controlled via the parameter M , which determines the maximal number of nonzero eigenvalues of the SPDM. For M = 1 the GP approximation is obtained; M = 2 allows the description of a Bose gas occupying two natural orbitals, and so on. Accordingly, the value of M has to be chosen sufficiently high in order to obtain trustworthy results. Fig. 9a shows the time-evolution of the ratchet current for various values of M . We see that the curves undergo drastic changes when M is increased from 1 to 2, and also from 2 to 3. However, when M is increased further, these changes become minor. The difference in the value of the current for M = 4 and M = 5 stays below 0.06 for the considered time range. For the results presented in Section V, we chose M = 4. However, the relative difference σ(t) can approach values up to 200 near those time points where the current vanishes. Note that the relative difference between two curves diverges whenever their roots do not coincide exactly.
The influence of the single-particle basis has also been checked. For this study we have considered approximations up to M = 3 and used different numbers of modes. We observed that angular momentum modes beyond the k = 4 mode are only weakly occupied, which is why we dropped them for the more accurate calculations, and used a basis of nine angular momentum eigenmodes ranging from k = −4 to k = 4. Figure 9b shows the combined normalized occupation of those modes with the same magnitude of angular momentum, given by σ k ≡ n k +n −k /N.
Considered are those values of k that lie outside the truncated three-level description (k = 2, 3, 4). As can be seen, the contribution of outside modes drops roughly by factors of 100 for increasing value of k. This implies that for the precision of the results presented here, the contribution of the modes | ± 4 is sufficiently small. The normalized occupation of the angular momentum modes (k = ±4) is below 5 × 10 −6 during the considered time range. This upper bound for σ 4 during the first 50 driving cycles holds for all values of g. The employed MCTDHB code uses adaptive step-size integrators which are adjusted to a high accuracy, such that the results are trustworthy. We have checked the influence of integration accuracy in the beginning of the calculations and adjusted step size tolerance to a sensible value, which was used throughout our study in the following.