Josephson oscillations in split one-dimensional Bose gases

We consider the non-equilibrium dynamics of a weakly interacting Bose gas tightly confined to a highly elongated double well potential. We use a self-consistent time-dependent Hartree--Fock approximation in combination with a projection of the full three-dimensional theory to several coupled one-dimensional channels. This allows us to model the time-dependent splitting and phase imprinting of a gas initially confined to a single quasi one-dimensional potential well and obtain a microscopic description of the ensuing damped Josephson oscillations.

Given that the sine-Gordon description only applies in an appropriate scaling limit a crucial question is how close the experiments are to this regime. In equilibrium correlation functions obtained from time-of-flight measurements of the boson density were found to be in good agreement with classical field simulations of the sine-Gordon model [11]. In non-equilibrium situations like the ones studied in Refs. [10,14,15] the situation is much less clear. In these experiments two elongated Bose gases are prepared in a quantum state characterized by a phase difference between the two gases. A tunnel coupling between the gases is then applied, which induces Josephson-like oscillations of density and phase. These oscillations quickly damp out and the distribution function of the phase is seen to narrow. Various studies based on the sine-Gordon model have so far failed to account for these observations [62,65,66]. In particular, taking into account Gaussian fluctuations on top of the solution of the classical field equations in a self-consistent way produces only very weakly damped Josephsonlike oscillations [62,66]. Given this state of affairs it is natural to question whether the experiments are in the right regime for a sine-Gordon based description to apply. In the experiments one deals with three dimensional bosons in a time-dependent confining potential. An obvious question is how good the low-energy projection to two one-dimensional Bose gases is in the experimentally relevant parameter regime. Another important issue is that the initial state that is prepared after splitting the gas and imprinting a phase difference is in fact not known, as the splitting process has so far only been modelled in a qualitative phenomenological way [67,68], or via methods that rely on a two-mode approximation [69], a classical field approximation [70] or a restriction to the transverse direction only [69]. In order to start addressing these questions we return to the drawing board and consider a gas of weakly interacting bosons subject to a tight harmonic potential in the z-direction, a time-dependent double well potential V ⊥ (y, t) in the y-direction and a shallow harmonic potential in the x-direction. This leads to the Hamiltonian H 3D (t) = d 3 rΨ † ( r ) − ∇ 2 2m + mω 2 x 2 x 2 + V ⊥ (y, t) + mω 2 z 2 z 2 Ψ ( r ) where r = (x, y, z) is the 3D coordinate andÛ ( r ) is the effective interaction potential We will always consider elongated gases with ω x ω z and refer to the x-direction as the longitudinal, and the remaining coordinates r ≡ (y, z) as the transverse directions. In order to make contact with the experiments of Ref. [14] we use with the values c 1 = 2π · 2.52 kHz, c 2 = 2.17 µm and I c = 0.4. For I(t) = I c , V ⊥ is a quartic potential with a flat bottom and for I > I c , it develops a double well structure. The term F (t)y is used to imprint a phase difference between the gases in the two wells (a precise description of what we mean by this is given below). The idea is to use (1) to describe the splitting of the gas, the phase imprinting and finally the subsequent non-equilibrium dynamics, but to take advantage of the fact that (i) interactions are weak; (ii) the confinement is tight in the y-and z-directions. The combination of these two allows us to project the full three dimensional theory to a small number of one-dimensional channels The resulting Hamiltonian is time-dependent but retardation effects are negligible. Some comments on the projection procedure are provided in Appendix A. We stress, that as a result of working with the instantaneous basis of single-particle eigenstates of −∂ 2 y /2m + V ⊥ (y, t) our effective one-dimensional field operatorsψ † a (x, t) have an explicit time dependence. After working out how to obtain H proj (t) from (2) we proceed as follows: 1. We treat the interactions in a time-dependent self-consistent Hartree-Fock approximation (SCHFA).
As we are dealing with an effective one-dimensional system we do not allow for the formation of long-range order. The resulting approximation is quite different from Gross-Pitaevskii theory.
The main attraction of the SCHFA is that it can be implemented straightforwardly, while its main limitation is that it treats interaction effects in a rather crude way. However, it is nonetheless expected to provide a good description as long as the interaction strength is sufficiently weak and the energy density in the system is not too low. We first identify a corresponding parameter regime and then model the Josephson oscillations experiments in this regime.
2. We start with a confining potential that forms a single elongated well and initialize our system in a thermal low-temperature state.
3. We evolve the state under a time-dependent transverse potential V ⊥ (y, t) that models the splitting and phase imprinting protocols used in the experiments. This provides us with a characterization of the "initial state" used in the Josephson-like oscillation experiments.
4. Finally we consider the non-equilibrium evolution of the split, phase-imprinted state. We observe damped Josephson-like oscillations.
This paper is organized as follows. In Sec. 2, we describe the low-energy projection used to arrive at the Hamiltonian (4), and the rationale for considering one-dimensional field operators carrying explicit time-dependence. In Sec. 3, we introduce the observables relevant to experiment, and connect them to the Green's functions of one-dimensional field operators that are computed in this paper. Sec. 4 introduces the self-consistent time-dependent Hartree-Fock approximation, and the resulting nonlinear partial differential equations that govern the time-evolution of the experimentally relevant Green's functions. Sec. 5 describes how the initial state of the system is modelled by preparing a gas in a thermal state of a single well and splitting it by a deformation of the trapping potential. The time-dependent definition of the one-dimensional field operators is shown to be an important tool in enabling this model for the preparation sequence. In Sec. 6, numerical results are presented for the time-evolution of experimentally relevant observables after the preparation stage. Density-phase oscillations are observed to be strongly damped over timescales that are comparable to those seen in the experiment.
2 Time-dependent projection to one-dimensional channels We start from the 3D Hamiltonian (1) with δ-interactions (2), where we have definedD u = −∂ 2 u /2m + mω 2 u u 2 /2 for u = x, z, and The 3D Bose fieldΨ( r ) satisfies the usual bosonic commutation relations. We use a double well potential of the form (3) for V ⊥ (y, t) throughout this paper, where the phase imprinting is implemented by the imbalance potential F (t)y. To arrive at an effective 1D model, we expand the 3D field operator in an instantaneous basis of single-particle eigenstates of the quadratic part of the Hamiltonian Here the single-particle eigenstates fulfil and we have defined canonical 1D Bose field operators bŷ At this stage we are dealing with an infinite number of Bose fields. We now exploit the fact that the single-particle eigenvalues ω z (c + 1/2) and b (t) constitute very large energy scales for c ≥ 1 and b ≥ā, and that interactions are weak. This allows us to truncate the expansion of the 3D Bose fields (7) to a small, finite number of channels. The rationale for working with explicitly time-dependent single-particle states rather than working in a fixed basis is that a subset chosen to span the low-energy subspace of the free part of the Hamiltonian (5) at t = 0 will in general only span the low-energy subspace at later times if we include a large number of channels. This would make the truncation much less efficient. Let us now give the details of the truncation procedure outlined above. If ω z ω x , we expect the dynamics to be frozen into the lowest single-particle eigenstates in the z-direction. We can then project to the corresponding low-energy subspace by truncating the expansion (7) to the c = 0 term. In the y-direction, the double well V ⊥ (y, t) gives rise to more states in the low-energy sector than just the ground state. We therefore need to retain multiple single-particle eigenstates Φ b (y, t). These wave functions are explicitly time-dependent eigenstates of the double well operatorD y (t) from Eq. (6), with eigenvalues a (t). If these eigenvalues show a gap above energy a−1 (t) that is large compared to all other energy scales in the problem for all times, the expansion (7) can be truncated at a = a. The resulting projection of the 3D field operator to the low-energy sector then readŝ where we have definedψ for all times. When starting from the 3D Hamiltonian (5), inserting the projected operator (10) and integrating in the y, z-directions leads to a model for a species of bosons, with coupling constants that are given by overlap tensors Corrections to (12) will be negligible as long as the following conditions hold: • Interactions are small. This holds by construction.
• The initial occupation numbers where ρ(0) is the density matrix at time t = 0, are very small for c > 0 and b ≥ā. We ensure that this is the case by working with an initial thermal density matrix at a sufficiently low energy density compared to ā (t). Experimentally this condition could be fulfilled by making V ⊥ (y, t) sufficiently tight.
• The transverse potential is changed slowly enough so that Tr[ρ(t) b † a,b,c (t)b a,b,c (t)] remain very small for c > 0 and b ≥ā. This provides a (rather obvious) restriction on the experimental protocol.
In equilibrium it is straightforward to evaluate the corrections to (12) and we outline the necessary steps in Appendix A. Perturbatively integrating out the high-energy channels above some cutoff generates all two, four and six boson interactions between the low-energy channels allowed by particle conservation. The interactions are very slightly retarded and non-local in space (the corresponding scales are set by the cut-off energy) but are negligible compared to the terms retained in (12). An analogous analysis can be in principle be carried out in the time-dependent situation of interest here, but as we don't require the corrections we do not follow this line of enquiry further.

Connection to previous literature
For time-independent double-well potentials with a very high tunnel-barrier, the lowest two singleparticle eigenstates Φ 0,1 (y) are approximately given by symmetric and anti-symmetric combinations of wave packets g L,R (y) that are localized in the left and right wells We can then define left-and right-localized one-dimensional Bose operatorsψ L,R ≈ (ψ 0 ±ψ 1 )/ √ 2. Inserting these definitions into Eq. (4) with a = 2 leads to the model of two Lieb-Liniger Hamiltonians connected by a tunnel-coupling term. The coupling constant of the Lieb-Liniger model is given by All other overlap tensors involving four combinations of g L,R (x) vanish if the two wells are separated by a high tunnel barrier, so that the only coupling between the left and right gases is given by the tunneling term proportional to ( 1 − 0 )/2. Eq. (15) is the Hamiltonian that is studied in most of the literature, following Ref. [45]. In this paper, we will instead focus on the more general Hamiltonian (12).

Three channel model
So far we have kept the numberā of one-dimensional channels in our theory arbitrary. In practice it turns out to be sufficient to work withā = 3 in order to accommodate the experimental situation realized in the Vienna group. The trapping geometries in these experiments are chosen so as to strongly suppress the occupation of the second excited level (a = 2). By including this suppressed level in our simulations and staying close to the experimental energy scales and trapping frequencies we can therefore ensure that we are in a regime where the occupation of the (time-dependent) third single-particle excited level can be safely neglected. The energy scales relevant to this reasoning are displayed in Fig. 1. This figure shows that at t = 0, the single-particle energy of the third level, which is neglected in our simulations, differs from the single-particle ground state energy by 3 Measurements and Green's functions

Measured operator in time of flight
The Bose gases in the double well can be probed through matter-wave interferometry [31,71,72]. After a tunable time t 0 spent in the double well, the Bose gases are released by turning off the trapping potential. This causes them to expand and overlap in three-dimensional space, and eventually their ( combined density is measured by absorption imagining after a "time of flight" t 1 . The theoretical description of this measurement process in the framework of a low-energy theory description is described in detail in Ref. [38] in the simpler case when two transverse single-particle states, given by Gaussian wave packets in the left and right wells respectively, are kept in the projection to an effectively onedimensional model. We will briefly recapitulate this construction, before expanding it to the case of more than two levels that are not perfectly localized in the wells. The absorption imagining can be thought of as a von-Neumann measurement of the boson density at time t 0 + t 1 The density operator is diagonal in the position eigenbasis {|r 1 , . . . , r N } which implies that the measurement outcomes are particle positions N j=1 δ(r−r j ) and the associated probability distribution is where (t 0 + t 1 ) is the density matrix of the system at time t 0 + t 1 . The moments of this probability distribution are The density matrix at time t is given by In the Heisenberg picture we have where the Heisenberg-picture field operators are given bŷ The approach of Ref. [72] is to relate the quantum state of the system after time-of-flight to the state at the time of trap release by assuming that interactions are negligible during the time of flight. This is a reasonable assumption since the gas, which is no longer constrained and expands in 3D, very quickly becomes highly dilute. The free, transverse expansion is then effectuated by the evolution operator This allows us to relate the Heisenberg picture field operators at times t 0 + t 1 and t 0 where G 0 (r, t) is the Green's function of the non-interacting boson Hamiltonian describing the free expansion. Now we exploit the fact that the initial density matrix ρ(0) involves only the low-energy sector, i.e. states in which only very few transverse modes are occupied. This allows us to project the field operatorsΨ (H) (r, t 0 ) and concomitantly the operatorsÔ (H) (r 1 , . . . , r n , t 0 + t 1 ) in the expression (22) for the moments M n to the low-energy description where S is some set of indices labeling single-particle states g a (y, t) in the transverse direction. The equations of motion of the Heisenberg picture operatorsψ a (x, t) are derived below in section 4. In [38], this set S = {L, R} refers to single-particle states with no explicit time-dependence that are localized in the left and right wells, respectively. In the following we will focus on the average over many absorption images Carrying out the convolutions we obtain Here we have defined and an analogous expression is obtained for Ξ 0 (z). The higher moments M n>1 can be related to expectation values in the low-energy description in the same way. In many works [44,72] it is assumed that the longitudinal expansion has little effect (even though it can be straightforwardly taken into account in a low-energy field theory framework in some cases [38]). This assumption is based on the state at the time of trap release: since the gas is spatially very constrained in the transverse directions, its momentum distribution in these directions is much broader than in the longitudinal direction. As a result, the time scale for expansion in the longitudinal direction far exceeds that for the transverse directions. If the time of flight t 1 is short, this suggests the approximation of neglecting longitudinal expansion altogether, replacing the free evolution operator (24) by This results in a simplified expression for the operatorρ We will use this approximate expression in much of the remainder of this work.

Green's functions of interest
In what follows, we will derive equations of motion for the Green's functions of the 1D Bose fields, defined as Solving these numerically gives us access to the expectation value ofρ (H) tof (r, t 0 + t 1 ) as well as averages of Fourier-transformed quantities like (41). In order to connect to the experimental works [10,14,15] we have to account for the fact that the data extracted from absorption imaging has been analyzed in terms of the number/phase representation for an effective two-channel model. Denoting the corresponding Bose field operators byψ L,R we can define averages of the relative density and phase via In order to connect to these quantities we need to expressψ L,R in terms of operators in our threechannel model. As interactions are weak this transformation can be taken to be linear. To be specific let us work with an effective three-channel model, i.e.ā = 3. We then carry out a change of basis such thatψ We have introduced a third, "excited" boson speciesψ e to be able to span the full space of 3 transverse levels. The set of labels used in Eq.
The transformation matrices c (α) j (t) are chosen with orthonormal rows and columns, so that they translate between the basis of single-particle eigenstates Φ 0,1,2 (y, t) of the transverse operatorD y (t) and another basis that contains left-and right-localized wave functions g L,R (y, t) as well es a third wave function, g e (y, t).
In [38], the wave functions g L,R (y, t) were simply given by (anti)symmetric combinations of Φ 0 and Φ 1 . However, the presence of the third wave function g e (y, t) now creates ambiguity, meaning that the c (α) j (t) can be defined in multiple ways. We will give two options here.
Choice 2: Since the double well is centered around y = 0, we find the vector c j (t)| 2 = 1. This fixes g L (y, t) as the single-particle wave function in the space spanned by Φ 0,1,2 (y, t) with the smallest possible probability for the particle to be found at y > 0, i.e. in the right well. Mutatis mutandis for c   j (t) and for I ≥ 0.55, the values are practically indistinguishable. We will therefore present results for the much simpler Choice 1, and comment on the changes that occur for Choice 2 wherever they are relevant.
Using Choice 1 and Eq. (33), the average relative density and phase (33) are given by Another quantity of experimental interest is the mean interference contrast, which we define as

Experimental data analysis and its relation to Green's functions
As discussed above the average over many absorption images gives access to in the three-channel model. We refer to Eq. (28) for the case when longitudinal expansion is taken into account. We will now show that the phase ϕ(x, t 0 ) of interest in Eq. (37) can be extracted from M 1 (r) by taking a suitable partial Fourier transform In the simpler case of gases whose transverse single-particle states are given by Gaussian wave packets in the left and right wells respectively the choice q = md/t 1 , where d is the distance between the wells' minima, gives access to the relative phase, see e.g. Ref. [38]. It is important to check how this situation is affected by the presence of the third channel and by the fact that the localized single-particle states are not given by perfect Gaussian wave packets. Studying the amplitudes A ij numerically for a given double well potential, we can establish which terms in (40) contribute at the wave vector Q = md/t 1 for a realistic potential in the three channel model. As shown in Fig. 2, A LR has a marked peak in Fourier space around q = md/t 1 . The diagonal terms ∼ A ii only contribute around q ≈ 0. The terms ∼ A Le and ∼ A Re do contribute at higher wave vectors, but their Fourier transforms both become very small around |q| = md/t 1 for all values of the double well (3) we consider. Moreover, the occupation of the "excited" transverse wave function g e (y, t) is much smaller than that of the wave functions g L,R (y, t). For these reasons, the Fourier transform (41) at q = md/t 1 is well approximated by (in the sense of taking expectation values of powers of the operator in states that belong to the low-energy subspace) whose argument then provides the relative phase of interest If on the other hand one works with a trapping potential where A Le and A Re do not have small 2.00 q ( ) | q( g * L g R) | | q( g * L g e) | | q( g * L g L + g * R g R) | q = md/t 1 Figure 2: Fourier transformed products of single-particle wave functions after time-of-flight g L,R (y, t 1 ) occurring in Eq. (31), for the parameters given below Eq. 3 with I = 0.5. The cross term g * L (y, t 1 )g R (y, t 1 ) (green) shows a peak around q = md/t 1 , whereas g * L (y, t 1 )g e (y, t 1 ) (cyan) becomes small there. The same can be said about the other cross terms involving g e . This allows to extract ϕ a (x) using Eq. (43).
Fourier components at q = md/t 1 and if the occupation of g e (y, t) is not small, the identification (42) might fail. Another likely scenario is that the value of |q| = md/t 1 cannot be established to sufficient precision. In such cases, Eq. (40) shows how different boson bilinears contribute to the measured density after time-of-flight, using numerical evaluations of the amplitudes A ij (y, t 0 , t 1 ).
Our way of extracting the relative phase from the average over many absorption images should be contrasted to the way in which the experiments [10,14,15] have been analyzed. In these works a value for a relative phase φ(x, t 0 ) is extracted for each (typical) absorption image by fitting the observed density profile to an expression of the form where f (y, z, t 1 ) is a Gaussian envelope. The data is then analyzed in terms of the average φ(x, t 0 ) over many shots. An interesting open question is to establish the precise relation between φ(x, t 0 ) and ϕ(x, t 0 ) extracted from the average over many images.

Hartree-Fock time evolution
Having established how Green's functions are related to averages over experimental measurements, we now consider their time evolution. We do so in the Heisenberg picture, indicated with a superscript (H), and consider the equations of motion for the 1D field operators, Here U (t) is the time-evolution operator and the additional, explicit time-derivative is nonzero due to the time-dependent definition ofψ a (x, t), via the corresponding eigenstates Φ a (y, t) of the transverse potential V ⊥ (y, t). In order to work out the last term on the right hand side of (45) we revert to the expansion of the Bose field into channels without projection (7) Using the orthonormality of the single-particle wave functions we obtain Using our assumption that the transverse potential is varying sufficiently slowly in time we can project these equations to our model withā transverse channels Physically, this term in the equation of motion (45) keeps track of transitions a → b to different levels due to time-dependence in V ⊥ (y, t). In what follows, we will drop the superscript (H) and fix a = 3. We now make the Hartree-Fock approximation for the interaction term, Using the symmetry of Γ abcd (t) the truncated Heisenberg equation (45) then yields the self-consistent equations describing the time evolution of the Green's functions of interest, namely C ab (x, x , t) with a, b = 0, . . . , a − 1. The HF approximation is equivalent to neglecting all higher connected n-point functions other than these Green's functions. The self-consistency of the HF scheme is implemented by the effective potentials with B(t) given by Eq. (49).
The system of Eqs. (51) can be solved numerically. In our implementation, we use a mixed implicitexplicit method for the propagation in time, employing a Crank-Nicholson scheme for the terms linear in Green's functions and a first order forward Euler method for the nonlinear terms. We work on a 2D square spatial grid of linear size 250 µm, using 1000 × 1000 grid points and approximating spatial derivatives by fourth order finite differences. We have checked convergence with respect to the grid size as well as the time step, which is 0.015 ms in the figures presented below. At each time step during the preparation sequence, the matrix B(t) given by Eq. (49) is computed for the lowest a eigenfunctions corresponding to the potential V ⊥ (y, t).

Quality of the SCHF approximation in equilibrium
An important question is how well we expect the HF approximation to work. It is well known [76] that at sufficiently low temperatures, 1D Bose gases form a quasi-condensate which is not well captured in the HF approximation. Specifically, the 1D boson density develops a central density peak which is underestimated by HF calculations. To make this precise we consider the simpler case of the Lieb-Liniger model in a harmonic trap V (x), where we can compare finite-temperature HF computations to results using Yang-Yang thermodynamics combined with the Local Density Approximation (YY+LDA). The LDA method is expected to provide highly accurate results in the appropriate parameter regime and its application to the Lieb-Liniger model have been described in detail in [75]. It has been successfully tested in experimental settings [77] and we will use it to compute the quantities The expectation values · HF are computed by the methods of Sec. 5.2 and using Wick's theorem. The expectation values · YY+LDA , on the other hand, are computed by numerically solving the thermodynamic Bethe Ansatz equations at finite temperature [78], using a chemical potential that is slowly varying in space µ(x) = µ 0 − V (x). For ∆ 2 , the Hellman-Feynman theorem must be used in addition [75]. The criterion for LDA to be applicable [75] can be checked a posteriori, and is found to be satisfied everywhere away from the boundaries of the gas for our parameters. A comparison between HF and YY+LDA for density profiles ρ 0 = ψ † (x)ψ(x) of a single gas is presented in Fig. 3. We see that while the HF approximation works quite well overall, it does underestimate the central peak. This failure occurs above a certain particle number, and the number where this cross-over occurs decreases with temperature. We will therefore work at a relatively high temperature of T = 60 nK in what follows. To make sure our particle number does not exceed the crossover where HF fails, we have plotted ∆ 1,2 for a range of particle numbers and longitudinal trapping frequencies in Fig. 4. This allows to monitor the quality of HF in the initial state for the parameters of our simulation. In particular, ∆ 2 shows whether the connected 4-point function, which is zero in HF, has a significant value in the initial state. For T = 60 nK and N 200, Fig. 4(b) shows it to be small. We note that our self-consistent Hartree-Fock results can in principle be improved upon for weak interactions and small particle numbers using only the self-consistently determined Green's function ψ † (x)ψ(x) HF , combined with perturbation theory. Rather than simply using Wick's theorem to compute the expectation values · HF occurring Eq. (53) in terms of ψ † (x)ψ(x) HF , one could include perturbative corrections to this results, working in powers of the interaction strength and performing contractions using the self-consistently determined Green's function ψ † (x)ψ(x) HF . We have checked the first order term and observed that it brings the result closer to the YY+LDA result for small particle numbers (N 300), whereas the correction starts to diverge for larger particle numbers. to ω = 2π · 7.5 Hz (green), ω = 2π · 10 Hz (cyan), ω = 2π · 12.5 Hz (blue) and ω = 2π · 15 Hz (red).
5 Initial state and gas splitting

Preparation sequence
We now have an equation of motion at hand for the relevant Green's functions that enter observables. Starting from an appropriate initial state, we can thus simulate the effect of the gas splitting, phase imprinting and free evolution performed in the experiments [10,14,15]. We implement these manipulations through the functions I(t) and F (t) which are present in the definition (3) of V ⊥ (y, t). We distinguish a number of stages: 1. A single gas is prepared in a thermal state. The transverse confining potential is a single well with a flat bottom, given by (3) with I = I c and F = 0.
2. We raise the double well barrier over some time t r by increasing I linearly from I c to I max . At t = t r we are left with a split gas and a high tunnel barrier.
3. We raise one of the wells over a time t imp by increasing F (t) linearly from 0 to F max > 0. Physically, this imprints a phase difference between the wells.
4. We remove the imbalance between the wells by tuning F (t) back down to zero in time t imp .
5. Finally we lower the tunnel barrier somewhat to enable tunneling on the relevant time scales, by decreasing I from I max to I f in a time t low .

Numerical determination of the initial state
At stage 1, the system is initialized in a thermal state of the Hamiltonian (12), subject to the HF condition (50). This state is determined as follows. We expand the field operators aŝ where χ α are real eigenfunctions of the harmonic oscillator potential in the x-direction, and we keep n m + 1 such modes. The Hamiltonian (12) subject to (50) can then be written as with the tensors Reshaping h aα,bβ and diagonalizing the resulting matrix numerically yields a canonical transformation The new creation and annihilation operators diagonalize the Hamiltonian H which combined with (56) forms a self-consistent system of equations. We proceed by iteration: starting from an initial guess b † cγb dδ 0 , which we take to be thermal with respect to the non-interacting Hamiltonian, we diagonalize h aα,bβ and compute (58) with the resulting P and E. Reinserting into (56) leads to the next iteration, and we repeat until convergence is reached.
A major hurdle in the above procedure is presented by the overlap tensor Γ αβγδ . As we use n m = 1000 modes, this tensor is too large to store numerically. However, using known identities for Hermite polynomials [74], we can write (56) as where the tensors B pm αβ are 0 if α + β − 2m − p is odd and/or negative, and otherwise given by The considerably smaller tensors A p αβ can now be separately contracted with other terms in (56), leading to a great memory gain. Even so, evaluating and storing the tensors (61) is still a very slow process for n m = 1000. We therefore make a simplifying assumption: we set for some Λ, which we choose to be 40 in our numerics. To see how this is justified, we note that the Hamiltonian (55)-(56) implies the relation on the canonical transformations P for all a, a, α, α. The assumption (62) is therefore valid if the P aα,bβ become very small whenever |α − β| Λ. This is reasonable since the weak interactions are not expected to couple harmonic oscillator modes that have widely different numbers of nodes. We check a posteriory that this assumption is consistent and well within the range set by Λ. We have also checked the assumption explicitly for the case of n m = 400. Finally, we have verified that the Green's functions resulting from the above procedure remain time-independent when they are propagated in time under (51) with a time-independent potential V ⊥ (y, 0).
The above procedure yields a set of Green's functions C ij (x, y) which characterize the state of the system at t = 0. In the central region of the trap, with |x| < 3 µm, we find exponential decay of the Green's functions C ii (x, −x) for the parameters presented in Sec. 6.1. The associated correlation length is roughly 0.5 µm.

Josephson oscillations
We are now in a position to model the full experimental sequence. To do so, we first fix the values for various constants and parameters.

Experimental parameters
The transverse potential V ⊥ (y, t) is described by Eq. (3) and its time evolution follows Sec. 5 with I max = 5.8, I f = 0.5, t r = 4 ms, t imp = 2 ms and t low = 2 ms. This means that after a time t r + 2t imp + t low = 11 ms, the confining potential becomes time-independent, and the 1D field operators lose their explicit time-dependence as a result. We consider a temperature of 60 nK and take the transverse confining potential in the z-direction to be harmonic with ω z = 2π · 1.7 kHz. The s-wave scattering length and atomic mass for the experimental system of 87 Rb atoms [15] are a s ≈ 5.2 nm and m ≈ 1.4 · 10 −25 kg, respectively. This fixes all parameters in the problem.

Assessment of time-dependent truncation errors in a toy model
In our full model, the initial thermal state contains three different transverse levels which mutually interact. An example of the resulting initial density profiles is given in Fig. 5(a), with occupation of the higher levels being suppressed as expected thanks to their larger energy cost. For t > 0, the occupations can change in a way that is both due to interactions and to the non-adiabaticity of the deformation of V ⊥ (y, t). The latter is modelled by the additional term (49) in the equations of motion (45), which are truncated at a = a = 3. To assess the error made in this truncation, we briefly consider the quantum mechanical problem of bosons in a double well V ⊥ (y, t). We discard the x-direction and set interactions to zero, so that the problem is given by Eq. (51) in the absence of x-dependence and with Γ ijkl = 0. This problem can be integrated numerically for any value of the truncation index a. Results for a = 3 (as we use in the full model) and a = 15 are compared in Fig. 5(b). The lines remain close, showing that the truncation error has a very small effect on transitions induced by the time-dependence of V ⊥ (y, t). ii (x, x,

Characterization of the quantum state after the preparation sequence
In our Hartree-Fock approximation, the state of the system at time t is fully determined by the Green's functions C ij (x, y, t), with i, j = 0, . . . , a − 1. Having these at hand thus allows us to give a full, quantitative description of the state of the system after the splitting and phase imprinting procedure, at the level of Hartree-Fock. This is a major improvement beyond existing, more phenomenological [67,68] or approximate [69,70] methods. As an illustration of the ability of our method to provide the full Green's functions, we plot C 00 (x, y, t) and |C 01 (x, y, t)| at time t = t r + 2t imp + t low = 11 ms, that is, after the preparation sequence (see Fig. 6). One sees that the Green's functions are strongly peaked around the main diagonal. To further illustrate their behavior, it is therefore instructive to plot the diagonal (Fig. 7) and anti-diagonal (Fig. 8) of the Green's functions of interest.
An equivalent way to express the same Green's functions is through the occupation numbers whereb † aα creates a particle in instantaneous eigenstate ξ α (x)Φ a (y, t), as defined via Eq. (54). As the occupation numbers are strongly suppressed away from the diagonal, we display this diagonal (Fig. 10) and the anti-diagonal pertaining to α + β = 20 (Fig. 9) for all occupation numbers of interest.
In a nutshell, the preparation sequence described above provides us with an initial state characterized by very short-ranged correlations. In the centre of the trap the correlation length is roughly 0.5µm, in line with the correlation length at t = 0.

Damping of density-phase oscillations
By monitoring the observables from Sec. 3, we can follow the relative density and phase between the gases. As soon as the barrier is lowered (step 5. in Sec. 5), oscillations in the relative density and phase can be observed, cf. Fig. 11(a), with an offset of a quarter period between the two. Importantly, the amplitude shows an initial period of damping, for all particle numbers we have considered. The mean interference contrast C(x, t), on the other hand, shows only very limited time-dependence. We C 00 (x, −x) Figure 8: The spatial anti-diagonal of Green's functions C ij (x, y) of interest, with the same parameters as in Fig. 6.  α,β , defined in Eq. (64). The chosen parameters are the same as in Fig. 6.
have fitted the density-phase oscillations at the center of the trap between t = 11 ms and t = 35 ms to and extracted the damping time τ and frequency ω. We stress that this is by no means a full description of the phase oscillations but merely a phenomenological formula to quantify the time scale τ of the damping observed in the early oscillation stage of the HF simulation. The dependence of this damping time τ on N is displayed in Fig. 11(b), whereas the dependence of the frequency ω on N is displayed in Fig. .12. There is a range of values of N for which the damping time as a function of N is in (a) qualitative agreement with the power-law dependence reported in [10,14]. For N ∼ 300, the behavior suddenly changes. This transition coincides with the breakdown of HF in the initial state: around this particle number, the errors ∆ 1,2 between HF and YY+LDA from Eq. (53) start to increase to significant values. This is displayed in the inset to Fig. 11(b). We thus conjecture that the deviation of τ (N ) from a power law for N 300 is due to a breakdown of HF in that regime. A number of additional observations can be made. The frequency of density-phase oscillations is highest at the center x c of the trap in the x-direction. Away from this point, the frequency is smaller, as displayed in Fig. 13(a). This figure also shows that the damping during the first few periods is somewhat weaker at points away from the trap center, where the gas density is smaller. Second, the gas as a whole shows a breathing motion. This can be shown by studying the squared longitudinal size of the left and right gas profiles, Fig. 13(b) shows that the squared longitudinal sizes of the left and right gases oscillate out of phase with one another. On top of this, there is an overall breathing motion of the gas with a frequency that depends monotonously on ω x . This breathing gets damped over a timescale that is large compared to the breathing period of the separate left and right gases. It is instructive to investigate the effect on the damping that various aspects of our set-up might have. First, there are two possible definitions of left-and right-localized bosonsψ L,R , as described in Sec. 3. As mentioned there, we stick to Choice 1 (cf. (36)) by default. Do our results, and the observed damping in particular, change if we switch to Choice 2? Fig. 11(a) shows results for Choice 2 in red. The curve is shown to lie very close to the blue curve, which was computed with Choice 1. This behavior occurred for all performed simulations, showing that the choice between Choices 1 and 2 does not significantly affect our results.
Second, we can investigate the effect of the second excited level, by turning off the corresponding couplings (13), setting Γ 2jkl = 0 for all permutations of indices. This completely shields the lowest two levels 0 and 1, and hence the relative density and phase (37), from any effects which level 2 might have. The resulting curves for ϕ fall on top of the curves for nonzero interaction with the second excited level, as exemplified by Fig. 14(a). We conclude that the effect of the additional boson species on the damping is negligible.
Third, we can study the effect of the longitudinal potential on the damping. This effect turns out to be very significant. In Fig.11(b), we see that the τ (N )-curves are shifted upwards as the strength of the potential is decreased. A weaker potential thus leads to a decrease in the damping effect. This suggests that in a box potential, the damping effect might be completely absent (within the SCHF approximation). We have therefore performed the same simulations in a box potential, by imposing hard wall boundary conditions at x = x c ± L/2 on the PDE (51). Fig. 14(b) shows a representative result, with parameters that are chosen to closely match those of Fig. 11(a). In particular, the bulk density is chosen to match the peak density from the initial condition of Fig. 11(a). The result is striking: in the box, no damping is visible at all. In fact, a very slight increase in the amplitude of the density-phase oscillations is observed.  Figure 14: (a) the same curve as the phase ϕ from Fig. 11(a), presented alongside the same quantity, but computed with Γ 2jkl = 0 for all permutations of indices. (b) oscillations of relative density n and phase ϕ for the same parameters as Fig. 11(a) but in a hard-wall box potential of size L = 80 µm. The bulk density is chosen to match the peak density from the initial condition of Fig. 11(a) 7 Beyond self-consistent Hartree-Fock The main attraction of the self-consistent Hartree-Fock approximation is its simplicity. However, it is not expected to provide a quantitatively accurate account of the non-equilibrium dynamics and it would be interesting to improve on it. A good way forward would be to implement the second Born approximation [79] following Refs. [80,81]. The two significant complications compared to these works are the absence of translational invariance and the explicit time-dependence of the Hamiltonian during the splitting and phase imprinting sequence. As a first step we consider the non-equilibrium evolution after the phase imprinting, which is described by a time-independent Hamiltonian (12) We now expand in harmonic oscillator modes notation and substitute this back into the expression for the Hamiltonian. Introducing a multi-index we can rewrite the Hamiltonian in a very compact form Here we have defined where Γ abcd andΓ ijkl are given by (13) and (56) respectively and k 1 = (a, i), k 2 = (b, j), k 3 = (c, k) and k 4 = (d, l). The second Born approximation for the single-particle Green's function can then be derived by generalizing the steps given in [81,82] to the case at hand. This results in the following set of equations of motion Y (k, p; q 1 , . . . , q 4 )e itE(q 1 ,...,q 4 ) G(q 1 , q 3 , 0)G(q 2 , q 4 , 0) L(k, p; q 1 , . . . , q 6 |t − s) G(q 1 , q 2 , s)G(q 3 , q 4 , s)G(q 5 , q 6 , s) , where E(k 1 , . . . , k 4 ) = ε(k 1 ) + ε(k 2 ) − ε(k 3 ) − ε(k 4 ) and the integral kernels are given by L(k, p; q 1 , . . . , q 6 |t) = 8 p X(k, p; q 1 , q 3 , q 6 , p; p, q 5 , q 2 , q 4 |t) + 16 p X(k, p; q 1 , q 3 , q 2 , p; p, q 5 , q 4 , q 6 |t) , K(k, p; q 1 , q 2 , q 3 , q 4 |t) = 8 k 1 ,k 2 X(k, p; k 1 , k 2 , q 2 , q 4 ; q 1 , q 3 , k 1 , k 2 |t) , The set of integro-differential equations (73) is clearly much more difficult to solve numerically than the self-consistent Hartree-Fock equations. The time integration is crucial at short times, while for sufficiently late times (73) ought to be reducible to a matrix quantum Boltzmann equation [81].
Integrating (73) is beyond the scope of this paper.

Conclusions
In this work we have developed a microscopic theory for the non-equilibrium evolution of bosons confined by a time-dependent quasi-one-dimensional trapping potential. Using that the transverse confinement is tight we have projected the full three-dimensional theory to a finite number of coupled, one-dimensional channels. By employing a time-dependent projection the number of channels that need to be retained in experimentally relevant parameter regimes is very small: three channels suffice. We then analyzed the resulting theory by means of a self-consistent time-dependent Hartree-Fock approximation and showed how the resulting Green's functions are related to averages of experimentally measured quantities. The Hartree-Fock approximation is expected to apply only for sufficiently weak interactions and sufficiently high energy densities. We have tried to identify a corresponding parameter regime by comparing the SCHF approximation to results obtained by combining the exact solution of the Lieb-Liniger model with a local density approximation in the trapping potential. On the basis of these considerations we restricted our initial states to temperatures of at least 60 nK and to particle numbers below ∼ 200. In this parameter regime we expect the HF method to work well at least at short times, when the neglected higher connected n-point functions have not had time to grow substantially. Our method has a number of attractive features. First, it allows to include the effects of various longitudinal potentials. Second, it can account for higher excited levels of the transverse confining potential which are normally neglected. Finally, it allows us to model the gas splitting and phase imprinting in a fully microscopic way. To our knowledge, such a model has not been presented before, and one of our main results is a characterization of the quantum state of the system after gas splitting and phase imprinting in terms of single-particle Green's functions of the one-dimensional channels. The second main result of our work is the description of the density-phase oscillations that ensue after the splitting and phase imprinting. In particular we find that these are damped over a few oscillation periods. These damped oscillations agree with recent measurements [10,14,15] in multiple ways. First, the damping time is inversely related to the number of particles, following a curve compatible with [10]. Second, the oscillation frequency decreases away from the center of the trap, as observed in [15]. We have shown that the coupling to the second excited level has very little effect on the damping. On the other hand, the longitudinal trapping potential is seen to play a very important role: the weaker the longitudinal trapping frequency, the weaker the damping. In a hard wall box, no damping is observed at all within our Hartree-Fock approximation. This suggests that damping effects are suppressed in this geometry for weak interactions. It therefore would be very interesting to repeat the experiments [10,14,15] in a hard-wall box potential. Such potentials are indeed under development [12,16] and our model can serve as a direct theoretical prediction for such setups.
The main limitation of our method is the way interactions are treated. In order to access the parameter regime of the experiments [10,14,15], in which the particle number was significantly higher than in our simulations, it is necessary to go beyond the SCHF approximation used here. A significant improvement would be provided by the second Born approximation discussed in section 7, but this is much harder to implement numerically. Ideally one would want to employ a controlled approximation scheme like [83,84] for our fully time-dependent problem.
Our work has several implications for attempts to describe Josephson oscillations in tunnel-coupled one-dimensional Bose gases based on the sine-Gordon model. A key finding of our work is the strong effect the longitudinal confining potential has on the damping of Josephson oscillations in the param-eter regime studied here. This suggests that the low-energy field their calculations based on the sine-Gordon model [62,[64][65][66] should be extended to account for the longitudinal confinement. Secondly, our characterization of the "initial state" after splitting and phase imprinting provides very useful information on what initial states to consider in the sine-Gordon framework. In the first instance one should consider Gaussian states with very short correlation lengths that reproduce the single-particle Green's functions reported here. Our microscopic modelling of the splitting process enables us to provide the same kind of information also in previously studied cases without phase-imprinting. where The situation we are interested in is where the eigenvalues k of the transverse confining potential constitute a large energy scale and the transverse level spacings |ε k − ε j | between highly excited transverse states are large too. We can then "integrate out" the transverse degrees of freedom above some cutoff Λ. Let us denote the first eigenvalue above Λ by ā , and rewrite the action as where S < is the part of the action that only involve the fields Ψ k , Ψ † k with 0 ≤ k <ā, S > the quadratic part of the action that involves only fields with k ≥ā, and S int are the remaining quartic terms that mix channels below and above the cutoff and describe interactions between channels above the cutoff. Defining we can eliminate the degrees of freedom above the cutoff using that we are dealing with weak interactions. Up to second order in S int we have the following expression for the low-energy part of the action The first order term generates hopping between the low-energy channels where W k 1 ,k 2 = ∞ n=ā 4V n,k 1 ,n,k 2 ∞ j=0 |χ j (x)| 2 e β ε k +ω(j+1/2) − 1 We see that this is small compared to the interaction strength c because the Bose occupation factors are by construction negligible. The second order term in S int contains all possible quadratic, quartic and sextic interactions involving ψ k (x, τ ) and ψ * k (x, τ ) compatible with particle number conservation, e.g.ā −1 where U k 1 ,k 2 ,k 3 ,k 4 (τ − τ , x, x ) = −8 ∞ n 2 ,n 2 =ā V n 1 ,k 1 ,n 2 ,k 2 V n 2 ,k 3 ,n 1 ,k 4 × G n 1 (τ − τ, x , x)G n 2 (τ − τ , x, x ) , As ε k > Λ the Matsubara Green's function of the high energy channels is very short-ranged in both imaginary time and space, so that retardation effects can be neglected and working with a purely local interaction between the low-energy channels remains justified. Hence the quartic terms generate only a very small renormalizations of the interaction terms already present between the low-energy channels.