On the low-energy description for tunnel-coupled one-dimensional Bose gases

We consider a model of two tunnel-coupled one-dimensional Bose gases with hard-wall boundary conditions. Bosonizing the model and retaining only the most relevant interactions leads to a decoupled theory consisting of a quantum sine-Gordon model and a free boson, describing respectively the antisymmetric and symmetric combinations of the phase fields. We go beyond this description by retaining the perturbation with the next smallest scaling dimension. This perturbation carries conformal spin and couples the two sectors. We carry out a detailed investigation of the effects of this coupling on the non-equilibrium dynamics of the model. We focus in particular on the role played by spatial inhomogeneities in the initial state in a quantum quench setup.


Introduction
The study of one-dimensional quantum many-body systems out of equilibrium has seen great progress in the past decades. Long-standing questions concerning the equilibration of observables, spreading of correlations and entanglement, and the emergence of statistical mechanics from microscopics have been successfully tackled using a range of innovative theoretical ideas [1][2][3][4][5][6][7][8][9], whilst spectacular advances in the ability to realize archetypical one-dimensional quantum many-body sytems using cold atoms [10][11][12][13][14] have made it possible to test many of these theoretical developments using tabletop experiments [15][16][17][18][19][20]. However, such experimental engineering of quantum many-body Hamiltonians relies on certain assumptions to make the experiments map onto a model of physical interest. These assumptions often include having a low energy density, at which an effective low-energy theory holds, and translational invariance, which can generally simplify the problem and specifically play an important role in the integrability of the low-energy theory. When studying non-equilibrium problems in finite quantum many-body systems, these two assumptions are sometimes brought into question.
Our interest lies in the effect of a finite tunnel barrier between the gases [14,[46][47][48]. This introduces a relevant perturbation and at sufficiently low energies leads to a decoupled theory of a Luttinger liquid describing the symmetric combination of Bose gas phases ("symmetric sector") and a sine-Gordon model [49] describing the relative phase ("antisymmetric sector"). The sine-Gordon model is of great theoretical importance as it is an exactly solvable, Lorentz invariant quantum field theory that exhibits a rich range of physical phenomena like dynamical mass generation and topological excitations and moreover has important applications to electronic degrees of freedom in solids [50]. Its behaviour out of equilibrium has received a lot of attention in the past decade. To be able to study dynamics, the very weakly interacting limit is amenable to a simple harmonic approximation [51][52][53], while the free fermion point can also be used to obtain exact results [51]. Integrability-based methods were used in Refs. [54][55][56][57] to study quenches from "integrable" initial states, whereas semiclassical methods [58,59] were applied to the study of the time-dependence of one and two-point functions as well as the probability distribution of the phase. The truncated conformal space approach [60] was employed in Ref. [61] to analyse the time evolution of two and four-point functions after a quantum quench. A first litmus test for the experimental realization of the sine-Gordon model using split Bose gas experiments was performed in an equilibrium situation: high order equilibirum correlation functions extracted from projective phase measurements in the classical limit have been found to agree well with classical field simulations [23]. A quantum many-body treatment of some such correlation functions is available as well, with possible generalizations to non-equilibrium situations [62]. For nonequilibrium initial conditions, however, experimental studies [24,63,64] have shown puzzling behaviour: when preparing two elongated Bose gases with an initial phase difference, applying a tunnel-coupling between them sets Josephson oscillations of density and phase in motion. These oscillations show a rapid damping, accompanied by a narrowing of the distribution function of the phase. To date, no satisfying theoretical explanation of this damping is known [65]. The damping seems incompatable with a description in terms of a translationally invariant sine-Gordon model, which fails to provide a mechanism for the observed strong and rapid damping in both a self-consistent harmonic treatment [66] and in a combination of truncated Wigner and truncated conformal space approaches [67].
In this work, we go beyond previous studies in two important ways: 1. We take into account the next most relevant perturbation at low energies. This perturbation induces an interaction between the symmetric and antisymmetric sectors.
2. We drop the assumption of translational invariance. To this end we place the model in a hardwall box geometry and consider inhomogeneous initial conditions.
Our strategy is to treat the resulting two-component model in the self-consistent time-dependent approximation (SCTDHA) as described in [66]. We consider the dynamics after initializing the system in a state in which the sectors are uncorrelated and observe how the new coupling term causes correlations between the two sectors to develop over time. In addition to this, energy starts to oscillate between the sectors. Depending on the initial density profile imprinted on the gas, Josephson oscillations of density and phase are affected by the presence of the additional term, showing modulations of the amplitude that differ from the ones observed in the SCTDHA treatment of isolated sine-Gordon dynamics [66]. This paper is organized as follows. In Sec. 2, we introduce the low-energy effective theory in a box geometry, the additional interaction term and the observable relevant for experiment. We also establish some notational conventions. In Sec. 3, we recapitulate the self-constistent time-dependent harmonic approximation as well as the framework to compute observables and some important distribution functions. In Sec. 4, we apply our formalism to an initial state which is commonly used in the literature, and present results on energy flow and growth of correlations between the sectors, along with the effect on Josephson oscillations, due to the additional interaction term. Sec. 5 summarizes our conclusions and discusses questions for further study.
2 Tunnel-coupled bose gases in a hard-wall box An appropriate model for the experiments carried out by the Vienna group is an interacting Bose gas confined in three-dimensional space by a tight harmonic potential in the z-direction, a double-well potential V ⊥ (y) in the y-direction and a shallow harmonic potential in the x-direction. We will refer to the x-direction as longitudinal, and to the remaining directions as transverse. To simplify the problem, we take the longitudinal potential to be an infinite square well Just like a shallow harmonic potential this breaks translational invariance in the longitudinal direction, but has the advantage to be considerably simpler to analyze. Our starting point is thus the following Hamiltonian where Ψ(x, y, z) are complex Bose fields obeying the usual bosonic commutation relations.

Low-energy effective theory
In situations where the transverse potentials are sufficiently tight, the dynamics in the y-and zdirections can be integrated out, in a way analogous to Ref. [70]. Details of this procedure will be reported elsewhere [71]. Projecting to the lowest two states of the transverse potential, and taking appropriate linear combinations of these, we obtain a Hamiltonian for two species of bosons, Ψ 1,2 , which are approximately localized in wells 1 and 2: Here the Bose fields Ψ(x) have commutation relations . The two Bose gases are coupled by a tunnelling term as well as contact interactions. The corresponding coupling constants Γ jklm follow from the details of the low-energy projection [71]. For our purposes, we will assume the diagonal elements to be equal to the usual Lieb-Liniger interaction constant, Γ jjjj = g ∀ j.
Hard-wall boundary conditions are imposed by restricting our problem to states |Φ where the density at the boundary has a vanishing eigenvalue: The one-dimensional model (3) gives an accurate description of the full theory H 3d at energies that are small compared to the energy E ⊥,2 of the second excited state of the transverse confining potential.
In the actual experiments this is a large energy scale. The physics of interest occurs at energies that are small compared to v/ξ E ⊥,2 , where ξ is the coherence length and v the speed of sound. This enables us to make a second low-energy projection by employing bosonization [42] This provides a low-energy description of (3) in terms of phase fields φ j and θ j with a cutoff length scale set by the coherence length of the gases, which for weak interactions is given by ξ = π/mv (the sound velocity v is defined below). The hard-wall condition is encoded in the boundary conditions of the θ-fields in a way that is described in Sec. 2.3. Let us first consider the case where interactions and tunnelling between the two gases are absent, meaning that both T ⊥ and the non-diagonal elements of Γ are zero. This leaves us with two Lieb-Liniger models in a hard-wall box, with interaction strength g. Under the mapping (5), the low-energy physics of this model maps to a pair of Luttinger liquids Here we have defined (anti)symmetric combinations of the phase fields by These fields are compact φ = φ + 2π, θ = θ + π and fulfil commutation relations This implies that the canonically conjugate fields to φ s/a are given by For weak interactions, the sound velocity v and Luttinger parameter K are related to the parameters in the Lieb-Liniger model in a simple way [72] Here γ = mg/ρ 0 is the dimensionless interaction strength and ρ 0 the average density of each of the two Bose gases.
In the next step we take into account the tunnelling term in (3) as well as "off-diagonal" interaction terms proportional to Γ ijkl with not all indices being equal. These introduce relevant perturbations (in the renormalization group sense) with respect to the critical Hamiltonian (6). Inserting the bosonization identity (5) and assuming Γ to be real, permutation symmetric and symmetric under 1 ↔ 2, we find that the perturbations with the lowest scaling dimensions can be written in the form where t ⊥ and σ depend on the microscopic parameters in (3). Importantly, the two terms in (11) get generated independently and we therefore will treat t ⊥ and σ as independent phenomenological parameters in the following. The Hamiltonian H s + H a + H ⊥ should be viewed as the result of integrating out high energy degrees of freedom in a renormalization group sense. As t ⊥ grows much faster than t ⊥ σ under the renormalization group it would be unphysical to consider very large values of σ. We have therefore restricted the numerical analyses reported below to the range 0 ≤ σ ≤ 2.
In addition to (11) there are other perturbations with higher scaling dimensions. Their systematic derivation as well as an analysis of their effects will be presented elsewhere [71]. In the case σ = 0 the full low-energy theory decouples into symmetric and antisymmetric sectors H = H s + H a , where H a is the Hamiltonian of a quantum sine-Gordon model [49] The non-equilibrium dynamics of this model was analyzed for the translationally invariant case in the framework of a SCTDHA in our recent work [66]. The additional σ-term in (11) couples the sine-Gordon model to the Luttinger liquid Hamiltonian H s . In the following we extend the analysis [66] to

Time-of-flight measurements
In the Vienna experiments [13, 14, 17-19, 23, 24, 73, 74] measurements are performed by turning off the trapping potential at some time t 0 , letting the gas expand freely and imaging the three-dimensional boson density after a time-of-flight t 1 . The outcome of each such "single-shot" measurement is determined by the eigenvalues e i 2 ϕa,s(x,t) of the bosonic vertex operators e i 2 φa,s(x,t 0 ) [26,28]. As shown in [26], the result of a single measurement of the boson density after a time-of-flight t 1 in the regime relevant for the Vienna experiments can be well approximated by Here d is the distance between the minima of the double well, x, x and r = (y, z) respectively denote longitudinal and transverse coordinates, and G(x, t) is the Green's function for a free particle The function f ( r, t) is an overall envelope whose precise from follows from the details of the trapping potential. By measuring tof , the system collapses to a simultaneous eigenstate of all e i 2 φa,s(x,t) . The outcome of such measurements can be simulated if one has access to distribution functions of the corresponding eigenvalues e i 2 ϕa,s(x,t) . Such distribution functions will be computed in Sec. 3.4. In principle, the observable (14) also contains small contributions from the density fields Π a,s (x) [26]. In order to treat these, the above description of a projective measurement has to be preceded by a diagonalization of the full observable, which now contains noncommuting fields. We do not pursue this further here because these effects are expected to be small in the regime where our low-energy approximation applies.
Experiments typically report results related to the quantity where we have defined

Mode expansions for the two-component Luttinger liquid
The free boson Hamiltonians H a,s are diagonalized by the mode expansions (see e.g. [72]) where The zero modes δN j have integer eigenvalues. The Hamiltonians then take the form Going back to Eq. (5), we see that the hard-wall condition (4) is guaranteed by choosing the c-number θ 0 such that It turns out to be useful in what follows to rewrite the mode-expansions in the form Here we have introduced a multi-index ν = (l, q) that runs over all positive momenta q ≥ 0 and the two sectors l = a, s and we have defined 3 Self-consistent time-dependent harmonic approxiation Our aim is to determine the non-equilibrium evolution after a quantum quench: the system is prepared in a density matrix ρ(0) that does not commute with the Hamiltonian (13). We moreover take the density matrix to be Gaussian for simplicity. The ensuing time evolution is described in the Schrödinger picture via the time evolving density matrix As our Hamiltonian of interest (13) is not solvable we resort to an analysis by means of a SCTDHA [41,66,[75][76][77]. Below we generalize the analysis of [66] to include the nonlinear interaction between the symmetric and antisymmetric sectors. The SCTDHA amounts to replacing the exact time evolution operator with where Here the functions g (1,2) (x, t) and h (1,2) (x, t) are determined self-consistently. In order to derive (29) we decompose the fields into their space and time dependent expectation values and their fluctuations Substituting this decomposition into the interaction part of the Hamiltonian (11) gives In the next step we expand the Hamiltonian to quadratic order in fluctuations following [66], which gives After re-expressing this in terms of the original fields φ a and Π s , we arrive at Eq. (29), where the functions h (j) (x, t) and g (j) (x, t) are determined self-consistently by Here we have defined two functions One subtlety associated with the SCTDHA concerns the zero mode φ a,0 . The spectrum of φ a,0 originally reflected the compact nature of the phase field φ a (x) = φ a (x) + 2π. The latter feature is lost in the SCTDHA, where fluctuations are assumed to be small but the fields themselves take arbitrary real values.

Gaussian initial states
In order to investigate the effects of the σ-term that couples the symmetric and antisymmetric sectors we want to start from a factorized state and study how correlations develop over time. An important requirement is related to our use of the SCTDHA: its accuracy strongly depends on the initial state obeying Wick's theorem. These two considerations lead us to consider the same class of initial states previously used in the literature [43][44][45] where ρ a (0) = |V, r, ϕ aa V, r, ϕ| is a Gaussian pure state It is useful to define new annihilation operators α a,k satisfying which are related to the b-operators via the canonical transformation In previous works it has been assumed that the symmetric sector is initialized in a thermal state [45]. We will follow this assumption, but in order to study the effects of spatial inhomogeneity we take our initial state to be given by a "displaced" thermal density matrix where the displacement operators are defined via This suggests the definition of displaced annihilation operators α s,k via a constant shift so that on the initial state. Since ρ s (0) satisfies Wick's theorem, it is then completely fixed by the vector R s,k along with connected two-point functions of the fields. Using the mode expansion of H s from Eq. (20) we simply find bosonic occupation numbers for q > 0, the anomalous expectation values b s,q b s,q c being zero. For the zero mode, the only expectation values on ρ s (0) that we will need are where the second identity implies ImR s,0 (0) = 0. As will become clear in the next section, expectation values involving the field φ s,0 will never be required for the computation of physical quantities.

Equations of motion
The SCTDHA allows for a closed-form expression of the equations of motion. We will work in the Heisenberg picture from here onwards. The SCTDHA guarantees that time evolving annihilation operators can always be written as where α µ are a set of bosonic creation and annihilation operators. We choose these to be given by where the α a,k are defined in (39) and the α s,k in (42). For (46) to be a canonical transformation we require The initial condition on R, S and T are given by T * ν,µ (0) = (sinh re iϕ ) pq if ν = (a, p), µ = (a, q), 0 else.
We note that the α µ 's satisfy Wick's theorem on the initial state, along with α µ = 0 for all µ.
The time evolution of any operator is then encoded in the time-dependence of the tensors R, S and T , which we will now determine. To this end, we write the SCTDHA Hamiltonian in the generic form The matrices A, B and vectors D, E depend on the self-consistency functions g (1,2) and h (1,2) , cf. Eqs. (34), and are given in Appendix A. Inserting the expansion (46) into the Heisenberg equation of motion yields a system of coupled, first order differential equations This system of ODE's is nonlinear : as a result of the self-consistency functions (34) on which the tensors A, B, D and E depend, these tensors are themselves functions of R, S and T , which therefore enter the system (52) in nonlinear combinations. To simplify some of the following equations we introduce linear combinations In terms of these functions mode expansions of the time evolved fields take the form The functions (35) can then be computed using Wick's theorem for the α-operators, based on the above expressions. This closes the system of ODE's (52). The zero mode in the symmetric sector φ s,0 reflects the compact nature of the phase field φ s and therefore needs to be treated separately from the finite momentum modes. We therefore define a field which time evolves as Importantly the zero mode φ s,0 does not get generated under Heisenberg time evolution of other fields. This is easily checked by inspection of the Hamiltonian (13) which is seen to not involve φ s,0 . This in turn implies that the zero mode cannot appear on the rhs of the Heisenberg equation of motion (51).
Since we can express the zero mode at t = 0 as we conclude that this linear combination of α-operators does not appear in the sums over modes in (54,55) except in the expansion for φ s (x, t), where it occurs in the term with ν = (s, 0). This directly leads to Re Q ν,(s,0) (t) = 0 ∀ ν = (s, 0), Im Q ν,(s,0) (t) = 0 ∀ ν.

One-point functions
As all relevant one-point functions of α ν and δN s are zero we have

Two-point functions
Comparing the definitions from Section 3.1 to the initial conditions (49), we find that for any ν, µ = (s, 0), If we define P (s) 0 to be the projector on the symmetric zero modes, along with its complement1 = 1 − P (s) 0 , we then find the following connected two-point functions In the above, indices on all matrices and vectors have been suppressed for conciseness. If we want to consider the field φ s instead of φ s , we need leave out the symmetric zero mode term. This leads, for instance, to and analogous modifications for φ s (x, t) φ s (y, t) c and φ s (x, t)φ a (y, t) c .

Full distribution functions
Individual measurement outcomes in interference experiments of interest [24] are fully determined by the eigenvalues ϕ a and ϕ s of the phase fields φ a and φ s [26], cf. Eq. (14). To model the outcomes of such measurements we therefore require the time-dependent distribution functions for ϕ a and ϕ s .
These can be determined in the framework of the SCTDHA [41,66]. For the case at hand, we first expand the eigenvalues of the phase fields as Fourier series, Here we have again used our multi-index notations µ = (j, q), where j = a, s labels the sector and q the momentum. Each measurement selects a particular set of Fourier coefficients and we denote the averages over many measurements by The mean values for the Fourier coefficients can be read off from the one-point functions calculated earlier, cf. Eqs. (60,61) The object of interest is then the time-dependent joint probability distribution P of Fourier coefficients {f µ }. Within the SCTDHA all cumulants of φ a,s other than the variance vanish, so that this probability distribution is Gaussian Here N is the total number of Fourier modes retained in (66). Noting that and comparing to Eq. (64), we can directly read off the covariance matrix as well: Having obtained a time-dependent probability distribution for the coefficients {f µ,t }, we can directly model experiments: we draw coefficients {f µ,t } from the distribution (69), reconstruct the corresponding eigenvalues (66), and insert these in the time-of-flight density (14) to compute the measured density profile. We note that in the notations used above the set {f µ } contains the non-physical Fourier coefficient f (s,0) . This quantity does not enter the observable (14), and can simply be discarded, whenever a set of coefficients is drawn from P ({f µ }, t). By repeating the above procedure for modelling a measurement many times over we can reconstruct the full distribution function of any observable that depends only on the phase fields φ a,s . In what follows, we will focus on the "interference term" in the spatially integrated density after time-of-flight R tof (x 0 , r, t 1 , t 0 ) defined in (16). The eigenvalues of this observable are proportional to where g ± (x) are defined in (17) and are related to the coefficients f µ via (66). Motivated by the experimental data analyses of Refs [21,22,65] we parametrize the interference term (72) as By drawing many sets {f µ } of coefficients from the distribution function P ({f µ }, t) and plotting the resulting values of Φ or C in a normalized histogram, we converge to probability distributions P Φ ,C for these quantities. These distribution functions can formally be written as 4 Results for experimentally relevant initial states

Choice of initial state
We now specialize to an initial state that is often used in the literature, see e.g. [43][44][45]. In these works, a quasi-classical argument is used to conjecture how the state of a pair of elongated Bose gases follows from the splitting process of a single gas. It is reasoned that when splitting a gas, each particle has an equal probability to end up in well 1 or in well 2. The relative particle number resulting from this poisson process is thus a stochastic variable with mean zero and variance proportional to the particle density. Assuming short-range correlations, one arrives at with η a phenomenological parameter which we will set to 1. Following [45], the delta function above is understood as a flat sum over plane waves running up to momentum π/ξ. To reproduce this initial two-point function, it suffices to use the initial state (37), with r a real and diagonal matrix and ϕ = 0. The resulting initial condition on Q, then leads to Π a (x, 0)Π a (y, 0) c = K L 2 e −2r 00 + j>0 qK πL cos (q j x) cos (q j y) e −2r jj .
Comparing Eqs. (76) and (49), we can thus read off for the antisymmetric sector. For the symmetric sector, we again follow Ref. [45]: the above quasiclassical splitting argument applies to the relative degrees of freedom, leaving the symmetric combinations of densities and phases unaltered. In [45], the symmetric sector is therefore taken to be in a finite temperature equilibrium state. We adhere to this conjecture here and use the thermal density matrix described in Section 3.1, thereby fixing the initial conditions for both T and S in conjunction with the above discussion. Finally, the initial conditions for R can be used to enforce various initial profiles on the density and phase fields in both sectors, which we will explore in Sec. 4.3 below.

Experimental parameters
We fix the parameters for our plots by following Ref. [21]: the one-dimensional density is taken to be ρ 0 = 45 µm −1 , the healing length is ξ = π/mv = π × 0.42 µm, the sound velocity is given by v ≈ 1.738 · 10 −3 m/s and the Luttinger parameter in our conventions is K ≈ 28. We take the onedimensional box size as large as we can achieve for a given value of the cutoff length scale, which amounts to L = 80 µm. This is comparable to the size reported in [21]. We work at a temperature of 5 nK throughout. In all figures, time is measured in units of the traversal time [4], t tr = L/2v, which is the time it takes for a light cone to reach the edge of the system from the centre of the box. We have chosen the value of the phenomenological tunnel coupling strength t ⊥ by considering a trade-off: we would like to maximize the Josephson frequency in order to follow as many density-phase oscillations as possible, whilst keeping the gap ∆ of the model's dispersion relation no larger than a small fraction of the energy cutoff in the Luttinger liquid. The latter is equal to c = vπ/ξ, with ξ the cutoff length scale. We have aimed for the ratio of the gap to the cutoff to be no larger than ∆/ c = 0.125, which we can guarantee by taking t ⊥ = 15 Hz for the above parameters. The only exception to the above is Fig. 2, where we take t ⊥ ≈ 1.17 Hz following Ref. [66], to enable a comparison with the case of periodic boundary conditions as presented there.

Time evolution
We now consider time evolution under the SCTDHA Hamiltonian (29), with the initial condition described in Sec. 4.1. Throughout, we choose R(0) such that The one-point functions φ s (x, 0) and Π s (x, 0) will be given different spatial profiles, to investigate the effects of broken translational invariance.
and σ = 0. This will serve as our benchmark, as it most closely resembles the translationally invariant scenario described in [66] in which the (anti)symmetric sectors remain uncorrelated. It is characterized by Josephson oscillations between density and phase, see Fig. 1(a), with a phase variance that initially grows, and then shows oscillating behavior, see Fig. 1(b).
To connect with our previous work [66] we include a comparison between results from that paper, where periodic boundary conditions were used, and the results derived for a box geometry in the present paper. Fig. 2 shows that the two geometries give extremely similar results in the centre of the trap for times below the traversal time, whereas deviations do occur after this time. It should also be noted that in [66] and Fig. 2, results are presented for smaller tunnel couplings (t ⊥ ≈ 1.17 Hz) than in the rest of this paper. The reason for choosing these values in [66] was that for a relatively shallow field potential, the inharmonicity of the cosine in the sine-Gordon model manifests itself more strongly, making deviations from the purely quadratic theory more apparent. For the purposes of this paper, however, it is more interesting to look at relatively large tunnel-couplings (t ⊥ = 15 Hz, see Sec. 4.2), as this enhances the coupling between the sectors in which we are interested.

Finite coupling between sectors (σ > 0) and homogeneous initial conditions
We next investigate different values of the coupling constant σ, and the resulting mixing between the sectors. Fig. 3 shows results for σ = 0, 1/2, 1, 3/2, 2, starting from completely flat profiles, as in Eqs.  (80), (81). When increasing σ, the phase oscillations remain essentially unchanged. A stronger effect is visible in the covariance between φ a and φ s , however. To quantify this, we define As can be seen in Fig. 3(b), the covariance C(x, t) increases to appreciable values as σ is increased. We also note that for larger values of σ, the variance of the relative phase increases somewhat for times below the traversal time, see Fig. 4.  Figure 4: Variance of the relative phase, for σ = 0 (blue) and σ = 2 (red). A slight increase in the variance is visible for the larger value of σ for times below t tr = L/2v It is also instructive to consider the energy flow between different terms in the Hamiltonian. To this end we define the following quantities We note that the total energy density, which is given by e sG (t) + e int (t) + H s /L, is independent of time, as required for a closed quantum system. Since we are interested in the time dependence of the various energy densities we subtract their values in the initial state and consider To quantify the effects of the σ-coupling on the flow of energy from and to the sine-Gordon model we show ∆e SG (t) in Fig. 5. To ascertain which fraction of the energy change is due to the kinetic and interaction parts of the sine-Gordon model we also show ∆e a (t) and ∆e ⊥,a (t) in Fig. 5(a). We observe that the change in ∆e SG (t) is very small, as significantly larger changes in ∆e a (t) and ∆e ⊥,a (t) largely compensate each other. In Fig. 5(b) we show how much of the energy from the sine-Gordon model ∆e SG (t) ends up in the new interaction term e int (t) and how much goes to H s (t) /L.

Finite coupling between sectors (σ > 0) and inhomogeneous initial conditions
As a next step, we investigate the effect of initial density profiles Π s (x) that are spatially inhomogeneous. These profiles will evolve in time as is shown in Fig. 6 (a,b). The profiles φ a (x) and Π a (x) are strongly affected by the strength of the σ-coupling to the inhomogeneous profile Π s (x) and develop inhomogeneities as a consequence. This is illustrated in Figs. 7(a,b) and has repercussions for the Josephson oscillations. The latter now displays spatial variations, which are caused by an effective Josephson frequency that has become σ-and position-dependent due to the presence of the space-dependent Π s (x)-field in the interaction term. This local and σ-dependent Josephson frequency is illustrated in Fig. 8. The spatial average of the phase, which is equal to the zero mode φ a0 , does not show any σ-dependence in its Josephson frequency, see Figs. 9,10. In this case, however, a σ-dependent modulation in the amplitudes is visible: the Josephson oscillations at different points in the box move out of phase due to the spatially varying Josephson frequency mentioned above. This leads to a decrease in the spatial average.
For an inhomogeneous profile of Π s (x, 0) , the covariance grows in time, in resemblance with the homogeneous case. This happens to an extent that is roughly proportional to σ. The same can be said of the energy flow between the (anti)symmetric sectors, as shown in Fig. 11. We see that the effects of the sector coupling term become stronger when we increase σ, but in the window of applicability of our bosonization based approach the effects remain small.

Distribution functions of the density after time of flight
As described in Sec. 3.4, our formalism allows the construction of distribution functions for the measured density after time-of-flight expansion. As a proof of principle we present such distribution functions in Fig. 12, for the observables Φ and C defined in Eq. (73).     : Distribution functions P Φ (α, t, t 1 ) (a) and P C (γ, t, t 1 ) for the observables Φ and C defined in Eq. (73). We choose a time of flight t 1 = 15 ms and integration length = 20 µm. The density profile used for these plots is homogeneous, with σ = 1.

Conclusion
We have extended the theory for non-equilibrium dynamics in pairs of elongated, tunnel-coupled Bose gases using a self-consistent time-dependent harmonic approximation (SCTDHA). In contrast to earlier works, we have studied the effect of a relevant perturbation which couples the (anti)symmetric sectors describing (anti)symmetric combinations of the two Bose gas phases. On top of this, we have dropped the assumtion of translational invariance by placing the system in a box and by imposing inhomogeneous initial density profiles. Starting from an initial state in which these sectors are uncorrelated, the coupling of the sectors under time evolution leads to a number of new but weak effects. First of all we observe the development of correlations between the sectors over time. This effect is present for all initial states we have considered, but the covariance between the sectors never reaches more than a few percent of the geometric mean of the variances. Second, the spreading of correlations is accompanied by a small transfer of energy between the sectors. And finally, the presence of the coupling term makes the dynamics in the antisymmetric sector susceptible to the breaking of translational invariance in the symmetric sector. The well-known Josephson oscillations of relative density and phase are modulated when taking an inhomogeneous initial density profile of the symmetric sector. This shows that the role of the trapping potential, which creates strong inhomogeneities, may play a more important role in experiment than was previously assumed. However, the model presented here does not capture the puzzling damping phenomenon observed recently [24,63,64]. This is not surprising given that our box potential is very different from the quadratic potential used in experiment. In future experiments, however, the box potential is likely to be used, which adds to the relevance of the calculations presented here.
We conclude that (i) the new term coupling the (anti)symmetric sectors leads to very weak effects. This means that the simulation of a sine-Gordon model using the setup described in this paper should not be severely hampered by the presence of this term. (ii) we have shown that it is possible to treat states with broken translational invariance in the SCTDHA formalism as presented in [66]. Combined with the sector coupling, we find that inhomogeneities in the density can have weak but nontrivial effects on the amplitude of Josephson oscillations. This means that the trapping potential is likely to have an effect on the dynamics probed in experiment. In a forthcoming paper, we will present a study of the projected Hamiltonian (3) in a microscopic analysis that takes a quadratic longitudinal potential into account. It would be interesting to combine such a microscopic approach with low-energy effective field theory calculations in the presence of such a quadratic trapping potential. However, the calculations using the box potential presented here may gain additional relevance when more experiments using a box potential, such as Refs. [68,69], are performed. VSB, Muller and Prins Bernhard Foundations.
The momentum indices q, k run within the blocks demarcated by solid lines, whereas the sector indices j = a, s change from one block to the other. The functions occurring above are defined via