Shapiro effect in atomchip-based bosonic Josephson junctions

We analyze the emergence of Shapiro resonances in tunnel-coupled Bose-Einstein condensates, realizing a bosonic Josephson junction. Our analysis is based on an experimentally relevant implementation using magnetic double well potentials on an atomchip. In this configuration the potential bias (implementing the junction voltage) and the potential barrier (realizing the Josephson link) are intrinsically coupled. We show that the dynamically driven system exhibits significantly enhanced Shapiro resonances which will facilitate experimental observation. To describe the systems response to the dynamic drive we compare a single-mode Gross-Pitaevskii (GP) description, an improved two-mode (TM) model and the self-consistent multi-configurational time dependent Hartree for Bosons (MCTDHB) method. We show that in the case of significant atom-atom interactions or strong driving, the spatial dynamics of the involved modes have to be taken into account, and only the MCTDHB method allows reliable predictions.


Introduction
Coherent tunnelling dynamics of macroscopic many-body quantum states through a classically forbidden barrier is one of the most striking manifestations of quantum physics. In the presence of an external drive, such systems exhibit the Shapiro-effect, also known as photon-assisted tunneling [1,2]. Shapiro resonances have been studied extensively in the context of solid-state superconducting Josephson junctions where a resonant modulation of the energy bias leads to a DC current across the junction [3]. Such resonances can be used to exactly quantify an applied DC potential difference across the junction and are nowadays used to implement a voltage standard.
A similar effect occurs in tunnel-coupled Bose-Einstein condensates (BEC) in double well potentials, where a resonant modulation of the bias energy lets the condensates oscillate between both wells [2]. Such dynamics could be used for an accurate measurement of the chemical potential difference across the tunnel barrier [4]. Despite the close analogy to the superconducting system, the dynamics of the bosonic Josephson junction are strongly influenced by atom-atom interactions [5]. Additionally the presence of a finite atom reservoir inhibits the observation of a true DC atom flux.
In recent years much effort has been devoted to the theoretical and experimental study of the tunnelling dynamics of driven Bose-Einstein condensates [6]. In optical lattice experiments, arrays of bosonic Josephson junctions have been realized [7]. Shaking of the lattice allowed to control the effective tunnel coupling [8] and the superfluid to Mott-insulator transition [9]. This effect was also used for dynamical localization [10]. The direct analogy to the superconducting Josephson junction has been realized in an all-optical double well potential. Here the influence of atomatom interactions has been evidenced, leading to new dynamical modes [11] and subpoissonian number statistics [12].
In this work we analyze the Shapiro effect for magnetic double well potentials on atomchips in view of a recent experimental realization of a bosonic Josephson junction [13]. We focus on the experimentally relevant situation where a modulation of the energy bias between the left and right BEC is accompanied by a concurrent modulation of the tunnel coupling. As we show, this leads to a significant enhancement of the Shapiro resonances compared to the conventional driving [2], which will facilitate experimental observation.
Shapiro resonances have been theoretically analyzed in [14,2] within a two-mode (TM) model employing static condensate wave functions [15,16]. For the configuration investigated in this work it is important to take into account the spatial dynamics, as transverse excitations might occur due to the external drive, which have a strong influence on the tunneling dynamics. We therefore employ the multi-configurational time dependent Hartree equations for Bosons (MCTDHB) method [17,18], which is a self-consistent, essentially exact framework, taking into account the full dynamics of the system. In previous work, the MCTDHB method applied to bosonic Josephson junctions allowed to discover dynamics very different to that predicted by the Gross-Pitaevskii equation and the two-mode model, such as an inverse splitting regime [19], the decay of macroscopic quantum self-trapping [20], and violation of time-reversal symmetry when switching interactions from repulsive to attractive in bosonic Josephson junctions [21]. In this work we study another case, where for typical interaction strengths the Shapiro resonances predicted from MCTDHB appreciably differ from the simpler models.
The paper is organized as follows: We start with the physical system and its dynamic description in section 2. In section 2.1 we introduce the realization of the double well potential and how we implement the driving, and in section 2.2 the three models which will be used and compared in this work: the Gross-Pitaevski (GP) equation, the two-mode model, and MCTDHB. In section 3 the non-interacting driven system is discussed, and results which demonstrate enhanced Shapiro resonances are shown. It includes also a short discussion of how to choose experimental parameters on order to observe the effect. Results for the interacting case are presented in section 4. First we show how the resonances shift due to interactions in section 4.1. A discussion of the different dynamical models, as well as the damping of the number imbalance is finally discussed in section 4.2.

Double well potential
The following analysis is motivated by a recent experimental realization of a bosonic Josephson junction on an atomchip [13]. A symmetric double well potential is generated by a combination of static and oscillating magnetic near fields, making use of dressed adiabatic states [22,23,24] ‡. The system consists of two elongated traps with strong atomic confinement (ω ⊥ ∼ 2π · 2 kHz) in the x and y (transverse) directions, and a very weak confinement along the z (longitudinal) direction with ω x,z /ω z ∼ 100. Tunnelling dynamics takes place along the coupling direction x which connects the two potential minima of the double well. The system is assumed to remain in the (many-particle) ground state in the two orthogonal directions y and z which do not contribute to the dynamics. The distance d between the double well minima and the height of the tunnel barrier are adjusted by controlling the amplitude of the oscillating field component which is in the radio frequency (RF) range. This RF amplitude represents the dynamical control parameter λ(t) used to drive the Shapiro resonances.
To implement a potential difference ("voltage") between both wells, the double well can be tilted in space (making use of gravity) [22,25]. Further possibilities are a spatial inhomogeneity in the amplitude of the RF field or employing additional forces like DC electric near fields [26] or local optical dipole potentials [11,27]. The resulting potential is given by equation 10 of [23] and can be approximately described by a symmetric fourth order polynom V λ (x) for the double well with an additional linear gradient g that ‡ The considerations made in this work are equally valid for double well systems based on solely static magnetic fields or on optical dipole potentials. implements the tilt: (1) Note that the control parameter λ (RF amplitude) affects both the separation of the double well minima (d(λ)) and the steepness of the confinement (expressed by the coefficients c 1 d(λ) and c 2 d(λ)). By construction, a larger separation coincides with a higher potential barrier, both reducing coupling, making λ the most sensitive parameter in the system. The most straightforward implementation of a Shapiro experiment would be to dynamically modulate the potential gradient g. However in a purely magnetic implementation of the double well as in [13] this would also result in a significant spatial displacement of the trap minima, which might lead to excitations and uncontrolled dynamics. We therefore focus the analysis on the effects of a periodic modulation of the RF amplitude (splitting parameter λ(t)) keeping g fixed. This results in double well minima moving along a potential slope as indicated in figure 1 (upper right inset). As we will show, modulating both the tunnel coupling and the potential difference leads to a significant enhancement of the observed Shapiro resonances.

Model
2.2.1. Many-body Hamiltonian In this work we consider the dynamics in the splitting direction x, and neglect any dynamics in the other directions y and z [28]. The dynamics of the atoms in the double well is then governed by the many-body Hamiltonian [29,30] whereΨ(x) is the bosonic field operator, andĥ is the bare Hamiltonian. U 0 is the effective 1D interaction strength, obtained by integrating out the other spatial directions for the ground state [28]. The 1D approximation is reasonable since typically the dynamics in the other directions decouples to a good extent whenever the 3D-Potential factorizes, i.e., Direct solution of the dynamics due to this Hamiltonian is not possible for more than a few atoms, and thus we have to use approximation schemes. In this paper we will employ and compare several types of approximations which we will outline in the following. They can be considered as different ways of restricting the field operatorΨ(x) to a small number of modes.

Gross-Pitaevskii (GP) equation
The mean field dynamics [29] is obtained by restricting the field operator to a single mode function,Ψ(x) =â 0 φ 0 (x). The Heisenberg equation of motion forΨ(x), obtained from equation 2, yields the GP equation For a BEC in a double well, the GP wave function can be written as is localized in the left (right) well, and thus the many-body wave function reads The GP equation assumes that the number distribution between the left and right localized condensate has binomial number fluctuations √ N (perfectly "coherent state" [31]). Since for the initial state we consider in this work always a condensate localized in one well, the GP equation will well describe the short-time-dynamics. At longer times, however, a single-mode description will not be valid anymore.

Two-mode (TM) model
In a double well it is more natural to use a basis which comprises two modes instead of one: Here, the localized mode functions φ λ L (x) and φ λ R (x) are depend on the control parameter λ.â L (â † L ) annihilates (creates) an atom in the left well, and similar forâ R (â † R ). There exist several schemes for the choice of the modes φ λ L (x) and φ λ R (x). The simplest one uses superpositions of the two lowest lying eigenstates φ λ g (x) and φ λ e (x) of the single-particle Schrödinger equation of the symmetric potential V λ (x), which have gerade and ungerade symmetry: [15,16]. The Hamiltonian reads in terms of pseudo-spin operators: Hereby,Ĵ z = 1 2 (â † LâL −â † RâR ) measures the atom number difference between left and right well, andĴ x = 1 2 (â † LâR +â † RâL ) promotes an atom from the left to the right well and vice versa. Equation 5 depends on the generic parameters Ω(t), ∆E(t) and κ(t), which denote the tunnel coupling, energy bias, and nonlinear interaction energy, respectively. They are given as The time-dependence in these parameters is due to a time-dependent control λ(t).
Tunnel coupling and energy bias between left and right condensates are shown in figure 1 over a wide range of λ from the unsplit to the split case. § While the tunnel coupling decreases with increased splitting, the energy bias becomes larger. Typically, we choose a constant κ ≈ U 0 /2, which is a good approximation because this value depends less crucial on the shape of the modes as compared to Ω(t) and ∆E(t).
An improved two-mode model [32] can be obtained by using the first and second selfconsistent states of the GP equation φ GP,λ In the following we always use the improved two-mode model, unless stated otherwise.

2.2.4
. Two-mode model with self-consistent orbitals and occupations A dynamical description including also the spatial dynamics can be obtained by the multiconfigurational time dependent Hartree equations for Bosons (MCTDHB) [17], which represents a framework where the two modes are included self-consistently. Using two time-dependent orbitals has the crucial advantage that transverse spatial excitations can be included, in contrast to a model with two fixed orbitals. A full discussion of the working equations of MCTDHB [17,28] can be found elsewhere. We only sketch here the main ideas of the method. The ansatz for time- for the case of two modes, § We use units where = 1, mass of a 87 Rb atom m = 1, and atom energy and time is scaled by the confinement length a ho = /(mω ho ) = 1 µm and energy ω ho of a harmonic oscillator. It follows that the units of time and energy are, respectively, 1/ ω ho = 1.37 ms and ω ho = 2π · 116.26 Hz.
Other, more subtle corrections to the standard two-mode model (e.g., to κ), as suggested by Bergeman et al. [32], will not be used here since it does not lead to significant improvements. although in principle the number of orbitals can be chosen at will. ¶ MCTDHB then provides a way to determine the "best possible" shape of the orbitals at a given time. This is achieved by formulating an action integral based on the above ansatz, and then using a variational principle. For the orbitals one obtains nonlinear equations, and for the number part two-mode equations similar to the TM model. For the orbitals one obtains nonlinear equations, and for the number part two-mode equations similar to the TM model. Most importantly, the orbital part and the number part are coupled and thus have to be solved self-consistently.
An important quantity is the one-body reduced density Ψ † (x)Ψ(x) [35]. It can be diagonalized in terms of natural orbitals φ i (x) and natural occupations ρ i (i = 1, 2): Whenever one natural occupation dominates, we have a BEC [36] and one orbital suffices. MCTDHB then coincides with the GP equation. Whenever several eigenstates are finite, we have an m-fold fragmented BEC, and m-orbitals should be used. Whenever MCTDHB-calculations converge when increasing the number of orbitals, the results can be considered as an exact solution of the many-body Schrödinger equation.
In this work we measure the degree of fragmentation by the difference of the population of the natural orbitals. (ρ 11 − ρ 22 )/N = ±1 corresponds to a single BEC, whereas (ρ 11 − ρ 22 )/N = 0 corresponds to a fully fragmented system.

Observable
The Shapiro effect in a superconducting Josephson junction is related to a finite DC component of the tunnel current at the resonance frequencies. A similar effect is present in a bosonic Josephson junction [1]. However, current cannot be measured directly and furthermore the reservoirs, consisting of the atoms in the left or right well, are finite. Therefore, the 'current' changes its sign whenever one reservoir is empty, and in such a manner the atoms oscillate between both wells.
The initial configuration for the following investigations consists of all atoms in the lower well. We can characterize the Shapiro resonances by a time average of zero atom number difference between the wells. The time averaged imbalance is then, similar as in [2], A value of 0.5 then corresponds to no population transfer at all, whereas a value close to zero indicates a resonance.

Shapiro resonances in absence of interactions
We now discuss the emergence of Shapiro resonances in the atomchip geometry (i.e., by driving the double well separation in the presence of a fixed potential gradient) for the non-interacting case. This system can well be captured within the TM model.

Enhanced Shapiro effect
The time-dependence of the control parameter is chosen as with driving frequency ω. In the following we take λ 0 = 0.675, corresponding to a splitting distance of approximately 1 µm, and a gradient g = 1.5042, corresponding to ∆E 0 ≈ 2π · 280 Hz. Due to the linear relationship between λ and the tilt ∆E at λ 0 (see figure 1), we have to a very good approximation ∆E(t) = ∆E 0 + b · sin(ωt), where b is the driving amplitude of the bias energy. The tunnel coupling instead depends not linearly on the control at λ 0 , but rather in a polynomial fashion. Thus, the general form of the tunnel coupling is We now decompose the Hamiltonian, equation 2 for U 0 = 0, into H = H 0 (t) + H 1 (t), with Hereby, n is an integer corresponding to the order of the resonance. Next we transform into an interaction picture [37] with respect to H 0 (t): where we exploited thatĴ x couples only neighboring states |k and |k ± 1 . Then we insert the generating function of the ordinary Bessel functions J l (z) and obtain with an effective tunnel coupling Ω eff n , where A l = J l b ω . Thus, the contributions to the effective tunnel coupling Ω eff n add up in a vectorial fashion. Further, we haveĴ φn = cos (φ n )Ĵ x +sin (φ n )Ĵ y where the phase φ n has no influence on the considered Shapiro dynamics of J z (t) when initially all atoms are in the lower well.
We neglected fast oscillating terms and kept only the DC contribution, which is reasonable close to the resonance ω ≈ ∆E 0 /n, and whenever the driving frequency ω is much larger than the coupling Ω 0 [38,14].
Thus, the driven system can be approximately described as an undriven system, but with renormalized tunnel coupling [39,2,38]. At the resonance condition nω = ∆E, there is no energy bias between the left and the right condensate, and thus the atoms tunnel at rate Ω eff ("Shapiro current"). For example, the case n = 0 corresponds to the trivial resonance of an unbiased system (∆E 0 = 0) at ω = 0, with Lorentzian line shape characterized by Ω eff . The effect of the driving is then to modify the depth and the width of the Lorentzian [39]. A system with finite energy bias (∆E 0 = 0) has similar Lorentzian shaped resonances, but at discrete frequencies ω = ∆E 0 /n.
Most importantly is the oscillatory part of Ω(t), which gives additional crosscontributions to Ω eff , proportional to J n−m b ω . Since J n−1 (x) > J n (x) for small values of x, those additional contributions to Ω eff can be considerably larger than those due to the first term of equation 17, proportional to J n b ω , and these contributions can drastically enhance width and depth of the resonances.
A typical resonance structure is reported in figure 2 (solid line), showing distinct resonances up to order n = 5. We have chosen a driving amplitude which increases with the order of the resonance as λ 1 = 0.03 · ∆E 0 /ω, such that higher resonances become more distinct. The n = 1 resonance has not only contributions from J 1 b ω , but also from J 0 Hence, the combined driving with ∆E(t) and Ω(t) yields more than the sum of driving with only one of them, and we thus term these resonances enhanced Shapiro  resonances.

Optimal choice of parameters
We shortly discuss the regime of parameters which are best suited to find clear Shapiro resonances, taking into account limitations posed by an experiment realization.
3.2.1. Tunnel coupling Ω The tunnel coupling Ω, which is determined by the mean double well separation, controls how fast the atoms tunnel. Therefore it should be large enough in order that the averaging period for equation 10 is not too long (typically ΩT = 15). If Ω is too large, the resonances are shifted according to ∆E → √ ∆E 2 + Ω 2 [2]. Due to the larger amplitude Rabi oscillations, the resonances get wider and less clear. Another disadvantage of too large Ω is that for an interacting system, the initial state will not be localized [40]. This could be compensated by choosing a larger tilt, which has however other disadvantages as discussed below.

Asymmetry
A too small asymmetry (linear gradient g) is unfavorable in experimental realizations as it leads to tunnelling dynamics already during the system preparation. Furthermore the initial state is not well localized, and no clear resonances can be identified when averaging over time. A too large asymmetry on the contrary requires large driving amplitudes in order to induce a Josephson current.

Driving amplitude b
We have already seen in section 3 that a too small driving amplitude reduces the Shapiro current due to the Bessel function structure. A too large b, however, has the effect of broadening the resonances (∼ J n (b/ω)).

Resonance shifts
Shapiro resonances in a many-body system are shifted in frequency due to atomatom interactions. For weak interactions, as in [2], the shift is roughly given by the interaction energy U 0 N , with the exact resonance frequencies given in terms of elliptic functions [41,2]. However, for typical experimental interaction strengths U 0 N 1 as we consider here, the resonance frequency depends itself on the driving amplitude λ 1 [39,42].
The shift of the resonances due to interactions is demonstrated using MCTDHB calculations (solid lines in Figs. 3 (a-c)), where J z T is shown for several interaction strengths U 0 N = 1, U 0 N = 2, and U 0 N = 4, respectively. With increasing interaction strength, the resonances become shifted towards each other. For relatively strong interactions U 0 N = 4 ( figure 3 (c)), the initial state is not completely localized and the contrast of the resonances is reduced.
The dependence of the resonance shift on the amplitude λ 1 is shown in figure 3 (d) for U 0 N = 1, calculated from the TM model. The shift becomes smaller for larger λ 1 . The black line corresponds to λ 1 = 0.035 · ∆E 0 /ω, which was used in Figs. 3 (a-c) to map out the resonance most clearly.
For U 0 N = 1 and for the TM model, we compare also to the case of artificially driving only with Ω(t) or ∆E(t), similar as in figure 2. We find that the enhanced Shapiro resonances are much more distinct than the usual Shapiro resonances, and the difference is much larger than in the non-interacting case. We also find that the resonance shifts of the usual Shapiro resonances (i.e., at constant Ω) are different from those of the enhanced Shapiro resonances and the case of constant ∆E. This is similar to the findings in [42], where it has been shown that an oscillating Ω(t) leads to a shift in the effective interaction strength.

Spatial dynamics and damping
Finally we discuss the differences between several models for calculating the resonances. The best description is provided by the MCTDHB method, such that the MCTDHB results in Figs. 3 (a-c) (solid lines) serve as a reference. The main question we want to answer here are the limitations of simpler models in describing the driven system.
We compare in Figs For U 0 N = 2 the GP fails to be a good description, especially for the first resonance. For U 0 N = 4 finally, the TM model does neither correctly reproduce depth and location of the resonances. The GP equation gives slightly better results than the TM-model, but still appreciably misses the correct results. From these simulations we see that, whereas for U 0 N = 1 and U 0 N = 2 the coincidence between TM and MCTDHB is very good for all times, for U 0 N = 4 the coincidence is lost at short times. Thus, whereas for weaker interactions Ω(t) and ∆E(t) as deduced from static orbitals are sufficient, the full self-consistency of MCTDHB is needed for stronger interactions (U 0 N = 4).
The spatial dynamics of the condensates as shown in figure 4 (a) is adiabatic to some extent, since the driving frequency is much smaller than the transverse trap frequency ∆E 0 ω x,y . However, as we have seen, for strong interactions (U 0 N = 4) the exact details of the condensate oscillations become important and have to be taken into account.
We show in Figs  of the TM model differs completely from the MCTDHB one after 10 ms. In order to point out the importance of the renormalization to the tunnel coupling according to the last term in equation 7, we also indicate results where it is not included (bright dashed line figure 4 (c)). The GP equation is in principle not a good description of the Shapiro dynamics, as it does not take into account damping of the Shapiro oscillations on a time scale of tens of ms. We show in the upper panels of Figs. 4 (b-c) the width of the number distribution ∆J z (t) for the TM calculations by the bright band around the mean value J z (t). Whereas for short times the number distribution is very narrow, it starts to become much broader after a time which depends on interactions. At the same time the GP results start to deviate from the exact ones. Finally, the distribution stays very broad and the mean value does not really change anymore, while the GP equation predicts oscillations of J z (t) with approximately constant amplitude.
Thus, the time averaged imbalance J z T does not fully characterize the system dynamics, and it is more appropriate to study the time evolution of the imbalance. However, it has been expressed in the literature [2] that for weak interactions and the "standard" Shapiro effect, the averaged imbalance J z T of the TM model and the GP equation should be the same, although the dynamics are very different in each case. For the configuration we discuss here we find this to be justified only for U 0 N = 1 (see figure 3 (a)). In this case the damping is slower than the period of the oscillations, and J z T corresponds approximately to the average of the extremal values, see the middle panel of figure 4 (b). However, for stronger interactions, the imbalance becomes stuck at an early point of the dynamics and thus the results for the averaged imbalance are significantly different from those predicted by the GP equation.
In order to have another perspective on the dynamics of damping, we plot the population difference of the natural orbitals φ i (x) (i = 1, 2). As can be seen in the lower panels of figure 4, during the Shapiro dynamics the BECs become fragmented into two incoherent ( â † 1â 2 = 0) condensates due to the nonlinear interactions. We find that indeed the GP description starts to fail whenever the system starts to fragment.
In addition we plot in figure 4 (d) also results for N = 1000 atoms (thin line with symbols), in order to demonstrate that fragmentation is clearly important also for larger numbers of atoms. Thus, convergence of MCTDHB results to GP results happens only at very large N . +

Conclusion
In conclusion, we have studied Shapiro resonances in a configuration where not only the bias potential, but also the tunnel coupling is driven dynamically. This is typical for double well potentials realized on atomchips and thus our findings are directly relevant for future experiments. We show that this configurations has favorable properties as it leads to enhanced Shapiro resonances. Due to a spatial deformation of the potential induced by the driving, the question of transverse excitations is of great interest. We find that, at least for significant interactions, the realistic MCTDHB method has to be used instead of simpler models in order to properly capture spatial dynamics of the involved modes. supported in part by NAWI GASS and the ESF Euroscores programm: EuroQuaser project QuDeGPM. We also acknowledge support by the Austrian Science Fund within projects P21080-N16, P22590-N16 and F41, by the EU project MIDAS, and by the Viennese Fund for Science and Technology (WWTF) projects MA45 and MA07-7. J.G. acknowledges support by the Humboldt-foundation.