Relaxation in an Extended Bosonic Josephson Junction

We present a detailed analysis of the relaxation dynamics in an extended bosonic Josephson junction. We show that stochastic classical field simulations using Gross-Pitaevskii equations in three spatial dimensions reproduce the main experimental findings of M. Pigneur et al., Phys. Rev. Lett. 120, 173601 (2018). We give an analytic solution describing the short time evolution through multimode dephasing. For longer times, the observed relaxation to a phase locked state is caused by nonlinear dynamics beyond the sine-Gordon model, induced by the longitudinal confinement potential and persisting even at zero temperature. Finally, we analyze different experimentally relevant trapping geometries to mitigate these effects. Our results provide the basis for future experimental implementations aiming to study nonlinear and quantum effects of the relaxation in extended bosonic Josephson junctions.


I. INTRODUCTION
The Josephson effect is a prominent example for the manifestation of macroscopic quantum effects. Originally formulated in the context of two weakly coupled superconductors [1] it shows a broad range of applications for systems featuring two coupled macroscopic quantum states. As such Josephson junctions were observed and analyzed in a variety of systems (see e.g. [2]). Increased interest over the last decades has been on its application to atomic systems, where two-body interactions enrich the dynamical behavior. This has led to a number of ongoing theoretical and experimental studies, from fermionic superfluids [3], macroscopic quantum self trapping [4,5], bosonic Josephson junctions [6][7][8][9][10][11], to different geometries [12,13] from small to extended junction arrays.
A recent experiment [14], studying an extended bosonic Josephson junction in an extended onedimensional (1D) geometry with two elongated (quasi-1D) 87 Rb superfluids in a double-well potential (see Fig. 1 and Fig. 2), makes the situation even less trivial. The observed Josephson oscillations (of the interwell atomnumber difference or of its conjugate variable, the global phase difference) were damped after only a few periods. More intriguing, this relaxation led to a phase locked state with strongly reduced fluctuations, observed experimentally through the almost straight interference fringes along the longitudinal axis of the trap (see Sec. II for details). This implies local damping of Josephson oscillations within each experimental realization and relaxation beyond simple dephasing dynamics. To date, this behavior could not be explained by various microscopic models [12,15,16] or its low energy effective description, the sine-Gordon model, not only in the semiclassical, but also in the quantum [17][18][19][20] treatment. * erne@atomchip.org Here we present a detailed numerical study of the quasicondensate dynamics, which explains the main results of the experiment [14]. For the experimentally relevant harmonic confinement we show that the system relaxes in two stages. The short time dynamics is fully described by multimode dephasing, already leading to a local damping of the oscillations. At longer times, we show that nonlinear dynamics beyond the sine-Gordon model causes the relaxation of the system to the observed phase locked state. We find this effect to persist even at zero temperature, which highlights the importance of understanding the relevance of classical nonlinear dynamics of thermally fluctuating fields when analyzing complex quantum many-body systems. Finally, we discuss different experimentally realizable trapping geometries to mitigate these effects.
The article is organized as follows. We start in Sec. II with an overview of our numerical simulations for the experimental system considered in [14] and give details of the numerical implementation and calculation of experimentally relevant observables. In Sec. III we present numerical results, reproducing the main findings of [14] over a wide range of initial conditions, and discuss the observed relaxation to a phase locked state. Lastly, in Sec. IV, we compare our numerical results to the effective one-dimensional description of the system, determining the microscopic origin of the observed relaxation and possibilities to mitigate the effect. We conclude our work in Sec. V.

II. NUMERICAL MODEL AND OBSERVABLES
Our numerical simulation follows closely the experiment in [14] where the extended bosonic Josephson junction was realized through an ultracold gas of 87 Rb atoms in an adjustable double-well potential on an atom chip. The gas was cooled well below degeneracy, such that the evolution of the system can be described by means of the  [14]. We consider two tunnel-coupled superfluids (red and blue ellipses) in a double-well (DW) potential (blue lines). By adjusting the barrier height the tunneling coupling J can be adjusted. The superfluids ψL,R are spatially separated and can be described in terms of density fluctuations δρL,R around a mean density profile ρ0 and a fluctuating phase φL,R (depicted in color). The tunneling-coupling J influences the antisymmetric phase φa and density fluctuations δρa. Considering only the global phase Φ and conjugate density difference n reduces the system to a two-mode model. The lower panel depicts schematically the experimental sequence (see also Fig. 2) from splitting a single condensate, introducing a phase difference (depicted in color) by applying a small tilt ∆E to the decoupled DW, to the Josephson oscillation regime where the finite tunneling-coupling J leads to coherent oscillations of particles between the two wells.

Gross-Pitaevskii equation (GPE)
The GPE describes the evolution for the order parameter Ψ(r, t) of a quasicondensate [21], capturing the contribution of nonlinear dynamics of classical fields. Here m and a s are the atom mass and s-wave scattering length of 87 Rb, respectively, and V is the external confinement potential modeled as The adiabatic radio-frequency potential V rf (r, u d ) on an atom chip [22] has a weak confinement along the longitudinal z-direction and can be continuously deformed from a single-well to a double-well (DW) potential (along the x-direction) for increasing values of u d (see Fig. 1). The control parameter u d , also known as the (normalized) dressing amplitude, determines the distance between the two wells, i.e. the height of the potential barrier in the DW. By means of the second control parameter u t it is possible to apply a small energy difference with respect to the transverse x-coordinate, leading to a tilted DW configuration used to initialize the Josephson oscillations in the experiment. Due to the tight radial confinement within each well, ν ⊥ ν , the system consists of two spatially separated elongated superfluids (see Fig. 1). If the typical energy scales of the gas are small compared to the energy of the first radially excited state, i.e. µ, k B T ω ⊥ , the system is in the quasi-1D regime, with the dynamics along the radial directions effectively frozen. Note that atomic repulsion may lead to a broadening of the radial wave function, which can be taken into account in the adiabatic limit [23]. We define the one-dimensional antisymmetric (relative) phase and conjugate 1D density difference with φ L,R and ρ L,R the longitudinal 1D phase and density profile of the left and right component, respectively (c.f. Fig. 1). In the following we drop the subscript 'a' when there is no risk of confusion. Note that, the potential in (2) is non-separable, which requires additional approximations when reducing the effective dimension of the system. Josephson effects in bosonic systems lead to a coherent oscillation of particles between the two sides of the DW via tunneling of atoms through the potential barrier (see Fig. 1). In the commonly used two-mode approximation the time evolution of the system is described by the Josephson equations (see e.g. [4,6,7]) where n(t) = (N L (t) − N R (t))/N and Φ(t) = Φ L (t) − Φ R (t) are the normalized atom number difference and conjugate relative phase difference between the two wells, respectively. Here µ is the on-site interaction energy and 2 J is the single particle tunneling coupling energy. In the following, we consider the Josephson regime 1 µ/(2 J) N 2 . For small initial amplitudes the Josephson equations describe oscillations of the phase and particle imbalance characterized by the plasma, or Josephson, frequency ω J ≈ 4Jµ/ .
For extended bosonic JJ the applicability of the twomode model strongly depends on the geometry of the system. In strongly elongated systems dynamics within each condensate may no longer be negligible, due to the high density of states along the extended direction. This naturally leads to a description of the system in the language of quantum field theory, where the Josephson equations (5) and (6) are replaced by the sine-Gordon model (see e.g. [24]) for the local fields φ(z, t) and δρ(z, t). A prominent example for the change in dynamical behavior is the breakdown of quantum self trapping (see e.g. [25]). Additionally, radially excited states may contribute to the dynamics, driving the system beyond the 1D regime. In order to reproduce the main findings of the original experiment we therefore decided to consider a full threedimensional classical fields simulation.
In the following, we first give an overview of the experimental sequence before discussing in the remainder of this section details of the numerical implementation and calculation of observables.

A. Overview of experimental sequence
We depict in Fig. 2 a typical evolution of the experimental system according to our numerical simulations. The experiment starts with a single quasi onedimensional condensate in thermal equilibrium. The initial trap at t = 0 has a weak harmonic confinement with frequency ν = 22 Hz along the longitudinal axis. The radial trapping frequencies are both given by ν ⊥ ≈ 3 kHz, such that the system is in the quasi-1D regime.
The trap begins to split at t > 0 and approaches at t = 21.5 ms its double-well configuration corresponding to the largest distance of the two halves of the quasicondensate. At this point the barrier of the DW is sufficiently high, such that the two 1D-BECs are completely decoupled. Unlike [26,27], splitting is assumed to be rather close to the adiabatic limit [28], and we end up in the situation when the quantum fluctuations are negligible to thermal ones. Therefore, the phase and density fluctuations after the splitting were found to be small in the experiment [14], signaling strong atom-number squeezing. This enables us to use finite-temperature classical statistical numerical methods, which do not include the quantum shot-noise, to describe the full experimental splitting process.
Next, a tilt along the transverse direction of the trap induces an energy difference between the two wells, which leads to the accumulation of a global phase difference Φ 0 (see Fig. 2 for t = 24.5 ms). The final phase difference Φ 0 can be adjusted by altering the maximum value of the tilt factor u t . Note that no particle number difference accumulates due to the complete decoupling of the two clouds. Experimentally the local relative phase φ can be extracted from the interference pattern resulting from the overlap of the wavefunctions after finite timeof-flight. Numerical results for the interference patterns after free expansion are depicted in the lower panel of Fig. 2. In accordance with the experiment, the straight interference fringes at t = 24.5 ms signal a well defined initial relative phase with negligible fluctuations along the whole condensate. This global relative phase Φ(t) can be extracted directly from the interference pictures by integrating along the longitudinal direction (white lines).
Lastly, within the period of t = 24.5 ms and t = 27.5 ms the distance of the two halves of the condensate is lowered to its final value (u d = u d and u t = 0), decreasing the barrier height. Tunneling through the barrier couples the two wells and the system starts to oscillate, realizing the extended bosonic Josephson junction in the relative phase φ and conjugate density difference δρ. The final potential for t ≥ 27.5 ms has a reduced harmonic confinement of ν = 12 Hz along the longitudinal direction and, within each well, an approximately harmonic confinement of ν ⊥ ≈ 1.5 kHz along the radial direction. Both frequencies ν and ν ⊥ are in nice agreement with the parameters of the original experiment [14].
In accordance with the experiment the global Josephson oscillation is rapidly damped and the system relaxes to a quasi-stationary state with a global relative phase Φ ≈ 0 and strongly reduced phase fluctuations. This so called phase-locked state is apparent in the stationary, almost straight interference fringes observed at long times (see last interference pattern in Fig. 2) and signals relaxation beyond local dephasing.

B. Numerical implementation
In order to prepare finite temperature initial conditions for Ψ(r, 0) we first compute ground state solutions of the GPE using imaginary time propagation [29]. In this context every wave function is normalized to a desired atom number N . The distribution of the total atom number N is assumed to be a normal distribution f (N ) = N (N , σ 2 N ) with the meanN = 3500 and standard deviation σ N = 0.16N . In the original experiment postselection further restricts the atom numbers to the range [N − δ N ,N + δ N ] with a given cut-off parameter δ N = 0.08N . In the numerical simulation the restricted atom number distribution is obtained by inversely sampling the cumulative distribution function F (N ) over the interval [F (N −δ N ), F (N +δ N )] using an equidistant distribution of n sr = 301 points representing the number of single runs.
Subsequently, the resulting zero temperature ground state solutions are propagated using the stochastic Gross-Pitaevskii equation (SGPE) (see e.g. [30] and references therein) until a new stationary thermal state is reached. Here, µ denotes the chemical potential of the eigenvalue problem at zero temperature and η is a complex random noise term. For the simple growth SGPE considered here, the positive constant γ can be freely tuned to improve the speed of convergence. We keep the atom number fixed within each SGPE realization by normalizing the wave function after each time step. This is done to preserve the exact atom number distribution already included in the above calculation of the ground state for each realization. Including the fluctuations of the SGPE would, however, only lead to a small broadening of the total atom number distribution.
The numerical propagation of the SGPE is based on a second-order accurate operator splitting and the spatial derivatives are approximated by means of the Fourier spectral collocation method. We note that the thermal noise term η is assumed to be constant for the duration of every time step. Depending on the desired temperature, several tens of thousands of time steps are necessary to reach a thermal state. In the simplest approximation, η denotes a complex, Gaussian, white noise process with correlations corresponding to a given temperature T . However, using the noise term in form of Eq. (8) results in unrealistic excitations of the quasicondensate along the tightly confined transverse directions of the trap. One obvious solution to this problem is to project the wave function onto a few of the lowest energy single particle eigenstates of the harmonic trap. This approach, which is known as the projected stochastic Gross-Pitaevskii equation [30], is prohibitively expensive in our three-dimensional setting. Due to the extremely strong transverse confinement of the initial trap k B T ω ⊥ , the main effect of radially excited single particle states is an increase in width of the Gaussian ground state wave function [23]. We therefore expect the desired thermal state to be an almost perfectly symmetric and smooth function with respect to the transverse directions x and y.
This assumption can be taken into account in the preparation of the thermal noise η(r, t) in Eq. (7). In this context, we first compute a complex field Ψ ⊥ (x, y, z) = Ψ(x, y, z)/ ρ(z) using the 1D density along the z-direction. For the noise term we finally employ the expression where λ(z) is one-dimensional Gaussian white noise with zero mean and variance given by Eq. (8) replacing δ(r − r ) with δ(z − z ). In the continuum limit an explicit cutoff for the noise is necessary due to the Raleigh-Jeans divergence, leading to a smeared out delta function. We checked independence of the cutoff for selected parameters. Convergence to the correct thermal state was verified by further evolving the system with the GPE after slightly disturbing the exact symmetry of the prepared states in the x, y-plane.
Once a set of thermal initial conditions Ψ(r, 0) has been prepared, the actual time evolution of all corresponding single runs is computed using the ordinary GPE (1). Similarly to the case of the SGPE we employ a second-order accurate operator splitting (Strang splitting) in combination with a Fourier spectral discretization of the spatial derivatives [31,32]. It is worth noting that the number of atoms N in every single run is conserved throughout the whole simulation as it is expected from the given numerical algorithm. We also checked that the total energy E is conserved (with high precision) as long as the two external parameters u t and u d remain unchanged (i.e. the Hamiltonian is time-independent).
The computational effort of a classical fields simulation in three spatial dimensions using several hundred initial states and on the order of 100 000 time steps is significant. However, the algorithm can be nicely implemented on a graphics processing unit (GPU) resulting in a dramatic speed-up in comparison with a CPU-based implementation on a multi-core CPU system.

C. Experimental observables
In the experiment, the evolution of the system is investigated by destructive measurements after time-of-flight (TOF) either looking at the interference pictures or the atom number difference between the two wells. We calculate the TOF expansion numerically, taking atom interactions into account for the first millisecond of the expansion. During this time, the system rapidly expands along the tightly confined radial direction, diluting the gas sufficiently such that the expansion becomes ballistic (i.e. non-interacting). Once the TOF-expansion has been computed, we integrate the density along the vertical direction y to obtain the desired interference pictures measured in the experiment. Additionally, to account for the finite imaging resolution in the experiment, we include a convolution of the numerical interference pictures with a Gaussian point-spread function The local relative phase φ(z, t) can be extracted by fitting a sinusoidal function to each pixel along the longitudinal z-direction. Finally, integrating along the longitudinal z-direction to obtain the profile n p (x), the contrast C tof and the global relative phase Φ tof are determined using where A, σ, C tof , k and Φ tof are found by solving a nonlinear least squares problem. The contrast C tof measures the visibility of the integrated interference fringes (c.f. Fig.2), i.e. it takes its maximum value C tof = 1 for negligible fluctuations of the relative phase along the whole condensate.
The computation of the global relative phase Φ tof by means of a TOF-simulation and formula (12) is a timeconsuming process. Alternatively, the global relative phase can be extracted from the in situ wave function (see e.g. [33,34] The results of both methods are practically indistinguishable and the weighting implicitly applied in Eq. (13) reflects the weighting involved in the procedure using the TOF-images reasonably well.
Consistently, the local phase profile φ(z, t) can be defined by omitting the integration over z in Eq. (13). Due to coherence along the tightly confined radial direction, the results are reasonably close to the line values Here, x 1 and x 2 are at the minimum of the radial potential in the left and right well (i.e. at the points of maximum density), respectively. Equivalently, the contrast C can be calculated from the in situ phase profiles via where we neglected the strongly suppressed density fluctuations [35]. We therefore consider in the following the in situ observables, omitting the time of flight expansion.

III. NUMERICAL RESULTS AND ANALYTIC ESTIMATES FOR THE EXPERIMENT
Numerical results of the experimental procedure outlined above are depicted in Fig. 3 for an initial global phase of Φ 0 = −1.25 rad. As mentioned earlier, the mean atom number amounts toN = 3500. Moreover, the thermal initial conditions correspond to a temperature of T = 20 nK.
The time-evolution of the local relative phase φ(z, t) of a selected (N = 3500) single run is shown in the upper panel of Fig. 3 (a). The corresponding time-evolution of the global relative phase Φ(t) is shown in the lower panel. In the real experiment it is impossible to observe the time-evolution of the relative phase in a selected single single run. Instead the experiment is repeated many times until meaningful statistical values can be extracted. We therefore depict in the upper and lower panel of Fig. 3 (b) the ensemble average of the local relative phase φ(z, t) and the global relative phase Φ(t) , respectively. We would like to recall that . . . denotes the ensemble average over n sr = 301 independent realizations, where we additionally take into account fluctuations of the total atom number N in accordance with the experiment. As the atom number difference is the canonical conjugate variable of the phase difference its time-evolution does not provide any new information and will therefore be omitted for brevity. The coherent splitting of a single condensate leads to the excitation of a common breathing mode, due to the halving of the atom number within each well. This breathing mode can be easily seen in Fig. 3 by the pair of black lines marking the region containing 75% of all particles. We note that the observed period T b ≈ 48 ms for ν = 12 Hz agrees perfectly with the theoretical prediction ν b = √ 3 ν for the breathing mode of a onedimensional condensate. For the timescales considered, we find the breathing mode to be sufficiently decoupled from the Josephson oscillation dynamics.

A. Local density approximation and dephasing
A first insight into the evolution of the system can be gained considering the local density approximation (LDA). In the Thomas-Fermi approximation the meanfield density profile is given by an inverted parabola ρ 0 (z) = n 0 (1 −z 2 ). Here n 0 is the peak density and we defined the dimensionless spatial coordinatez = z/R, where R is the Thomas-Fermi radius. Defining the local chemical potential µ(z) = gρ(z) leads to a spatially dependent Josephson frequency with ω J0 = 4Jn 0 / . The resulting phase profile describes an assembly of independent undamped Josephson junctions along the weakly confined longitudinal direction. The breathing of the condensate leads to a minor slow time dependence of ω J0 , which can be taken into account but does not significantly alter the results. We therefore in the following neglect the influence of the breathing for brevity. The spatial dependence of the Josephson frequency, decreasing towards the edges of the condensate, is clearly visible in Fig. 3 (a) through the bending of the local phase difference. Note that we are not in the self-trapped regime [36]. Therefore, different parts of the condensate exhibit local dephasing. This leads to a damping of the global phase difference even in the absence of dynamics along the longitudinal direction. Here H α (τ ) is the Struve function of the order α, predicting a power-law decay for the amplitude of the global Josephson oscillation [37].

B. Atom number fluctuations
The damping of the local and global relative phase in case of the ensemble averages is even stronger. Atom number fluctuations contribute to an additional dephasing between different realizations due to the dependence of the Josephson frequency ω J0 on the total particle number N . The most simple model Eq. (16) for the case of a harmonic longitudinal confinement yields ω J0 ∝ N κ with κ = 1/3. However, by fitting a whole series of 3D simulations covering a wide range of atom numbers, we find κ ≈ 0.43 (see Fig. 4), indicating an interwell tunneling that depends on the local density of the quasicon-  densates (c.f. e.g. [38]). Note that we chose the powerlaw dependence here based on the functional form predicted by Eq. (16) It is worth pointing out that fitting Eq. (19) to the ensemble average Φ(t) of our full numerical simulations, we find η to be reasonably close (within ≈ 10%) to its expected value η = κ δ N /N with κ ≈ 0.43.

C. Relaxation beyond dephasing
Despite being a good approximation for the evolution of Φ(t) , these effects alone do not correctly describe the observed relaxation of the system to a phase locked state. Firstly, the small variance of Φ(t) at late times (see Fig. 3(b)) shows that dephasing due to total atom number fluctuations is not the dominant effect for the damping. Equivalently, this implies damping of the Josephson oscillations within each realization. Secondly, while damping of Φ(t) is expected from Eq. (18), the observed local damping signals the breakdown of the LDA, which would predict an undamped oscillation for fixed position z (c.f. Fig. 3(a)). This difference can be clearly seen in the time evolution of the contrast C(t) depicted in Fig. 5, which increases at late times back close to its initial value. This, consistent with [14], is a clear indication for the relaxation to a phase locked state with small fluctuations of the relative phase along the whole condensate. Contrary, considering only local and atom number dephasing C(t) remains small despite the strongly damped global phase oscillations. Note that atom number fluctuations have only a minor influence on the evolution of the contrast C(t), as it is independent of the global phase Φ(t).

D. Comparison to the experiment
Finally, in a comprehensive simulation study we investigated the dependence of the damping time τ on the initial relative phase Φ 0 , the tunneling strength J and the mean atom numberN (see Fig. 6). In this context the damping time τ was determined by fitting the results of the numerical simulations using the model for a damped Josephson junction presented in [14], where a phenomenological damping term was added to Eq. (5) (see e.g. [39]). Due to the large number of different parameters it was necessary to reduce the number of single runs to n sr = 21. However, convergence of the extracted values was verified by (selected) computations using larger values of n sr . In all cases we found our results to be compatible with the experiment [14]. This includes a relatively weak dependence of τ on the initial global relative phase Φ 0 , a plateau-shaped dependence on the effective tunnel coupling strength J, and approximately τ ∼ N −0.5 . Note, however, that for the latter we find a slight additional dependence of the exponent on the single-particle tunneling coupling J.
Numerical results of a comprehensive simulation study illustrating the dependence of the damping time τ on the initial relative phase Φ0, tunneling strength J, and mean atom numberN . The final dressing amplitude u d determines the tunneling coupling (decreasing for larger u d ). If not explicitly depicted in the sub-figure, the remaining constant parameters areN = 3500, u d = 0.56, and Φ0 = −1.25, respectively. In all simulations the thermal initial conditions correspond to a temperature of T = 20 nK and the atom number distribution is characterized by the parameters σN = 0.16N and δN = 0.08N (see Sec. II B for details). Errorbars depict the standard deviation.

IV. MULTIMODE DEPHASING AND NONLINEAR RELAXATION
In order to study the fundamental mechanisms leading to the rapid local damping of the Josephson oscillations, we consider in the following a zero-temperature (T = 0) state with fixed particle number N = 3500. We further omit the preparation phase to mitigate the effect of a common breathing excitation and imprint the phase difference Φ 0 directly in the final trap configuration. These simplifications preserve the observed main results of rapid global and local damping of Josephson oscillations and enable us to compare results of the full 3D simulations with a tractable linearized, one-dimensional model (see e.g. [38,40]).

A. Effective one-dimensional description
We first proceed with the common dimensional reduction, reducing the system to two one-dimensional coupled quantum wires by integrating over the tightly confined transverse directions [23,41]. Regaining the dominant terms in the expansion, the system of two coupled 1D GPE equations can be described by the Hamiltonian The uncoupled evolution of the fields ψ L,R in the left (right) well of the double-well potential is described by where g 1D ≈ 2 a s ω ⊥ is the effective 1D interaction constant. The tunneling coupling is given by where we included the dominant nonlinear density dependence (see e.g. [38] for details). In general we have |J| |F max[ρ 0 ]|. Note, in particular, that also J can show an explicit density dependence due to radial swelling of the condensate [23]. We find the dominant effect of the nonlinear density dependence to be a moderate shift of ω J0 depending on the total atom number N (see Fig. 4). We therefore neglect in the following linearized model the nonlinear terms in Eq. (22), i.e. F = 0, and treat J as a free parameter to be determined from our 3D simulations.
Next, writing the fields in the Madelung representation ψ L,R = ρ 0 + δρ L,R e iφL,R , we expand in powers of the small density perturbations δρ and phase gradients |∂ z φ| (where | . . . | denotes the typical value) [42]. The Hamiltonian to quadratic order separates into a weakly coupled sum H ≈ H s + H a for the symmetric (s) and antisymmetric (a) degrees of freedom (DoF), defined as The Hamiltonians H s,a are given by where i = s, a and we suppressed the spatial and temporal dependence of the fields δρ(z, t), φ(z, t), and ρ 0 (z) for simplicity. Here ζ s(a) = 0.5 (1) for the symmetric (antisymmetric) DoF, δ ia is the Kronecker delta, and we included the quantum pressure and minor coupling correction (first and last term respectively) for completeness. If neglected, the symmetric and antisymmetric DoF are described by the Luttinger-Liquid (H LL ) and sine-Gordon (H sG ) model, respectively [24]. Note, in particular, that while second order in the small parameters δρ and |∂ z φ a | the sine-Gordon model is an interacting (quantum) field theory with H sG including terms with an arbitrary (even) number of fields φ a .

B. Analytic solutions for harmonically trapped systems
Within the framework of the sine-Gordon model the relaxation of Josephson oscillations in a homogeneous system, ρ 0 = const., was studied in [18][19][20]. This was, however, insufficient to describe the observed fast damping in the experiment [14]. Here, to investigate the influence of the longitudinal confinement, we present an analytic solution to the linearized equations of motion for a harmonically trapped system.
Under the assumption of (i) the Thomas-Fermi approximation (∂ z ρ 0 )/ρ 0 1, (ii) neglecting the nonlinear density dependence of the tunneling coupling, and (iii) small tunneling energy compared to the chemical potential J/µ 1 the linearized equations of motion for the relative phase field obey the eigenvalue equation where we assumed φ(t) ∼ e iωnt and defined the dimensionless eigenvalue λ 0n = √ 2ω n /ω and Josephson frequency γ = √ 2ω J0 /ω . We recall the definition of the dimensionless spatial coordinatez = z/R, where R is the Thomas-Fermi radius. The dominant correction beyond the LDA is given by the kinetic energy, first term in Eq. (26). If neglected, we recover the previous model of an assembly of undamped, uncoupled Josephson junctions with spatially dependent frequency. The kinetic energy term couples these oscillators, damping the rapid growth of the phase gradient caused by local dephasing.
Exact solutions to Eq. (26) are given by the angular oblate spheroidal wave functions S mn for m = 0 [37]. The eigenfrequencies ω n are depicted in the left panel of Fig. 7, featuring an increasing energy gap and twofold degeneracy of the spectrum for larger tunneling coupling J (cf. [38]). For vanishing tunneling coupling, γ = 0, Eq. (26) reduces to the Legendre differential equation, leading to the known spectrum for the excitations of the inhomogeneous Luttinger-Liquid model for a harmonically trapped quasi-condensate [43]. Imposing canonical commutation relations for δρ a and φ a determines the normalization of the mode functions and leads to the mode expansion of the phase and density quadratures φ a = n (2n + 1)µ 2R TF ω n ρ 0 S 0n (z) b n e iωnt +H.c. (27) δρ a = n (2n + 1) ω n ρ 0 2R TF µ S 0n (z) b n e iωnt +H.c. (28) In Fig. 7 (right panel) we show the occupation numbers |b n |, calculated for an initial constant phase difference φ(z) ≡ Φ 0 by projecting onto the quasiparticle basis, Eq. (27). In contrast to the spatially homogeneous system, where a constant global phase difference Φ 0 only has non-vanishing overlap with a single (k = 0) mode, here multiple modes are populated. Due to symmetry only even modes are occupied. Their amplitude decreases rapidly for n > 15, such that only the 10 lowest energy modes show significant population. We compare in Fig. 8(a) predictions of our analytic model to the evolution of the local relative phase φ(z, t) obtained from the full 3D GPE simulation. As previously mentioned, we consider a fixed particle number N = 3500, zero temperature T = 0, and omit the preparation phase to suppress the common breathing mode. While, e.g., the damping time τ depends on details of the initial state, these simplifications qualitatively preserve the main features of the relaxation (c.f. Fig. 3). In particular, the system still exhibits fast local damping of the Josephson oscillation and relaxes to a phase locked state on long timescales.
The time evolution in the harmonic approximation is given by the undamped oscillations between density and phase quadratures of the initially populated free quasiparticle modes (Eq. (27)). Each mode oscillates with its distinct eigenfrequency ω n = λ 0n ω / √ 2, where we treated the Josephson frequency ω J0 as a free parameter in Eq. (26), determined from the GPE simulations. The resulting multimode dephasing leads to a rapid initial damping of the global Josephson oscillation, similar to the previous local dephasing. Most notably, however, the dephasing of free quasiparticle modes already leads to a local damping of oscillations. At early times, we find excellent agreement between the full numerical simulation and the analytic predictions. For later times the linearized theory shows oscillations increasing again in amplitude, which propagate inwards from the boundaries. These are caused by (partial) rephasing dynamics as a result of the limited number of quasiparticle modes with significant population [44].
The absence of rephasing in the full 3D simulations signals the breakdown of the linearized model, either through nonlinear terms in the SG model or coupling of the symmetric and antisymmetric DoF. For the former the Hamiltonian of the system is still given by H ≈ H s + H a , such that the symmetric DoF can be neglected. Higher-order corrections couple the symmetric and antisymmetric DoF leading to a transfer of en- Note that, while the GPE obeys energy conservation, the depicted energies are calculated using Eq. (25) and hence are not strictly conserved. ergy between the two sectors, ultimately leading to complete equilibration of the system. This signals a definite breakdown for the effective description of the extended Josephson junction through the sine-Gordon model.
We quantify the coupling in Fig. 8(b), showing the time evolution of the energies E i (t) within each sector H s,a given by Eq. (25) at each instant of time. In addition, we display the energy of the zero-momentum (k = 0) mode of the antisymmetric DoF determining the spatially independent, global Josephson oscillation. Therefore, the initial energy E 0 is completely given by this zero-momentum mode. At early times the energies are approximately conserved within each sector, validating a description of the antisymmetric DoF in terms of the sine-Gordon model. For the parameters considered the linearized model (Eq. (26)) constitutes a good approximation in this regime, signaling that nonlinearities of the sine-Gordon model only lead to minor corrections. The rapid decline of the k = 0 mode is caused by the dephasing of quasiparticle modes, transferring energy to higher momentum states within the antisymmetric sector. Note that quasiparticle dephasing here leads to transport in momentum space since plane waves are not the eigenstates of the system. For 30 ms t 60 ms energy transfer from the antisymmetric to the symmetric DoF becomes dominant, revealing the breakdown of the SG model. At later times the system reaches equipartition of energy. This does however not imply complete thermalization of the system. We like to highlight that significant coupling already occurs at zero temperature due to the spatially inhomogeneous Josephson frequency. This is the dominant effect leading to the fast relaxation of the system to a phase locked state as observed in [14].

C. Flat-bottom potentials
The rapid equipartition of energy greatly impedes experimental studies of, e.g., the influence of quantum corrections or the long time evolution of extended bosonic Josephson junctions and the sine-Gordon model. Based on recent progress in shaping arbitrary trapping potentials [45,46] we now give an outlook on the possibilities for mitigating these effects in cold-atom experiments.
We show in Fig. 9(a) the time evolution of φ(z, t) and the normalized energies for a box-shaped potential, which is known to violate the system integrability, even if the potential is flat between the walls [47]. In accordance with typical experimental capabilities, we model the zdependence of V rf in Eq. (2) as with V 0 ≈ 15 kHz µ, length L = 100 µm, and finite wall width σ w = 2 µm. Spatial dependence of the Josephson frequency is limited to the edges of the condensate leading to disturbances emanating from the boundaries. Relaxation in flat-bottom potentials. Results of 3D-GPE simulations for parameters as in Fig. 8 for a boxshaped potential (a) and a ring-shaped potential (b) in the longitudinal z-direction. We adjusted the particle number N requiring equal Josephson frequencies ωJ (z = 0) for all trapping potentials. In neither case the system relaxes to a phase locked state since coupling of symmetric and antisymmetric degrees of freedom is highly suppressed. For (b) the Josephson oscillation completely decouples.
Once propagated inwards these disturbances lead to a rapid decline of the k = 0 mode caused by multimode dephasing. In contrast to the harmonic confinement, however, global Josephson oscillations prevail within the central region of the box (for t 20 ms). In particular, the amplitude of local phase oscillations remains high. Therefore, the contrast C(t) after its initial decrease remains small at longer times, i.e. the system does not relax to a phase locked state. Consistently, coupling between the symmetric and antsymmetric DoF is highly suppressed, with only ≈ 16% of the energy being transferred at t = 80 ms. Therefore, while the antisymmetric DoF shows multimode dynamics, decoupling from the symmetric DoF constitutes a good approximation.
Ring-shaped potentials further eliminate the influence of the boundaries, leading to undamped global Josephson oscillations at zero temperature (see Fig. 9(b)). In accordance, we find the total energy to remain in the k = 0 mode, with negligible coupling between the symmetric and antisymmetric DoF. Here, dephasing due to atom number fluctuations will dominate ensemble averages at late times, which can be mitigated through appropriate postselection. Naturally, this is the ideal setting to study the nonlinear dynamics and the influence of thermal and/or quantum fluctuations.

V. CONCLUSION
We gave a detailed discussion of the rich nonlinear dynamics in inhomogeneous extended bosonic Josephson junctions. We found our results for the full 3D-GPE simulations at finite temperature to reproduce the experimental findings [14] over a wide range of parameters.
A detailed analysis for the zero temperature case allowed us to distinguish two stages of the relaxation dynamics. The short time behavior was well described through the sine-Gordon model, i.e. the low-energy effective theory for the antisymmetric degrees of freedom of two tunnel-coupled one-dimensional superfluids. For the parameters considered, we found reasonable agreement with an analytic solution in the harmonic approximation. In contrast to the local density approximation, these solutions already predict the local damping of Josephson oscillations (in addition to the spatially dependent Josephson frequency). At later times, we explained the relaxation to a phase locked state, as observed in [14], through the breakdown of the sine-Gordon model description by coupling of the symmetric and antisymmetric degrees of freedom. Induced by the harmonic confinement, this coupling dominates the long time behavior of the system, even at zero temperature. Lastly, we showed that this coupling can be greatly reduced for box-or ring-shaped potentials along the longitudinal direction.
Our study is a crucial step when investigating the influence of thermal or quantum fluctuations on the relaxation dynamics in bosonic Josephson junctions. Our simulations provide, e.g., the relevant maximum timescale during which comparison to the sine-Gordon model is sensible. A detailed comparison to other microscopic mod-els, such as the self-consistent time-dependent Hartree-Fock approximation presented in [16], would be interesting but is hindered by the different regimes of applicability. Therein, results are presented in the small particle number regime N 200, wherein deviations of classical statistical simulations are expected to become relevant. Intriguingly, the behavior in [16] is similar to our harmonic 1D model, including local damping and recurrences of phase oscillations for harmonic confinements, and the absence of relaxation to a phase locked state. This, consistent with the considered small atom number limit, suggests that interactions play a less significant role in the regime of Ref. [16].
Our results highlight the importance of understanding the classical nonlinear dynamics, even at zero temperature, before comparing results to sophisticated quantum models. Our work outlines the experimental possibilities and steps necessary to study relaxation in extended bosonic Josephson junctions and the (quantum) sine-Gordon model out of equilibrium.