Collective excitations and tunneling dynamics in long bosonic Josephson junctions

TL;DR

This paper studies the low-energy excitations and tunneling dynamics of long bosonic Josephson junctions using advanced theoretical models, revealing complex oscillation behaviors beyond standard equations.

Contribution

It introduces a comprehensive analysis combining 2D Gross-Pitaevskii and Bogoliubov-de Gennes formalisms, along with a 1D hydrodynamic model, to better understand coupled Bose-Einstein condensates.

Findings

Multiple-frequency oscillations of atomic populations are explained by low-energy excitations.

The developed models accurately reproduce the excitation spectrum and population dynamics.

Generalization of the standard Josephson equations to include complex excitation effects.

Abstract

We investigate the low-energy dynamics of two coupled anisotropic Bose-Einstein condensates forming a long Josephson junction. The theoretical study is performed in the framework of the two-dimensional Gross-Pitaevskii equation and the Bogoliubov-de Gennes formalism. We analyze the excitation spectrum of the coupled Bose condensates and show how low-energy excitations of the condensates lead to multiple-frequency oscillations of the atomic populations in the two wells. This analysis generalizes the standard bosnic Josephson euqation approach. We also develop a one-dimensional hydrodynamic model of the coupled condensates, that is capable to reproduce the excitation spectrum and population dynamics of the system.

Equations

$$V(x,y,z)=\frac{m}{2}({\omega }_x^2{x}^2+{\omega }_y^2{y}^2+{\omega }_z^2{z}^2)+V_b{\mathrm{e}}^{-2y^2/{\lambda }^2}.$$

$$\mathrm{i}\hslash \text{partialderivative}t\Psi \left(r,t\right)=\left(-\frac{{\hslash }^2\text{laplacian}}{2m}+V\left(r\right)+g\text{absolutevalue}{\Psi }^2\right)\Psi \left(r,t\right),$$

$$g=\frac{4\pi {\hslash }^{2}a}{m}\sqrt{\frac{m{\omega }_z}{2\pi \hslash }},$$

$$Z(t)=\frac{N_1(t)-N_2(t)}{N},$$

$$N_j(t)={\int }_{A_j}\text{differential}r\text{absolutevalue}{\Psi \left(r,t\right)}^2,$$

$$\phi (t)={\theta }_2(t)-{\theta }_1(t),$$

$${\theta }_j(t)=\arg {\int }_{A_j}\text{differential}r\Psi \left(r,t\right).$$

$$\dot{Z}=-J\sqrt{1-Z^2}\sin (\phi ),$$

$$\dot{\phi }=\Lambda Z+J\frac{Z}{\sqrt{1-Z^2}}\cos (\phi ).$$

$${\omega }_p=\sqrt{J(\Lambda +J)},$$

$$\tilde{Z}({\omega }_n)=\sum _{m=0}^{M}Z(t_m){\mathrm{e}}^{-\mathrm{i}{\omega }_n{t}_m},{\omega }_n=2\pi n/t_M,$$

$$\Psi \left(r,t\right)={\mathrm{e}}^{-\frac{\mathrm{i}\mu t}{\hslash }}\left[{\psi }_0\left(r\right)+\delta \psi \left(r,t\right)\right].$$

$$\delta \psi \left(r,t\right)=\sum _k{c}_k\left[{\mathrm{e}}^{-\mathrm{i}{\omega }_{k}t}{u}_k\left(r\right)+{\mathrm{e}}^{\mathrm{i}{\omega }_{k}t}{v}_k^{\ast }\left(r\right)\right],$$

$$\hslash {\omega }_k{u}_k=({\hat{H}}_0+2g\text{absolutevalue}{{\psi }_0}^2-\mu )u_k+g{\psi }_0^2{v}_k,$$

$$-\hslash {\omega }_k{v}_k=({\hat{H}}_0+2g\text{absolutevalue}{{\psi }_0}^2-\mu )v_k+g{\psi }_0^2{u}_k,$$

$${\hat{H}}_0=-\frac{{\hslash }^2\text{laplacian}}{2m}+V\left(r\right).$$

$$\int d\mathbf{r}\left(|u_k{|}^2-|v_k{|}^2\right)=1.$$

$$Z=\frac{1}{N}⟨\Psi |\hat{Z}|\Psi ⟩,\hat{Z}\left(r\right)=\{\begin{array}{rl}1,& \mathbf{r}\in A_1(y<0),\\ -1,& \mathbf{r}\in A_2(y>0).\end{array}$$

$$Z(t)=\frac{2}{N}\sum _k{c}_k{D}_k\cos ({\omega }_{k}t).$$

$$D_k=⟨{\psi }_0|\hat{Z}|u_k+v_k⟩$$

$$\omega ={\omega }_x\sqrt{\frac{k(k+1)}{2}},$$

$$J_y(x,t)=\frac{\hslash }{2m\mathrm{i}}\int dy\left({\Psi }^{\ast }{\partial }_y\Psi -\Psi {\partial }_y{\Psi }^{\ast }\right).$$

$$J_y(x,t)=\sum _k{c}_k{j}_y^{(k)}(x)\sin ({\omega }_{k}t),$$

$$j_y^{(k)}(x)=\frac{\hslash }{m}\int dy\left[{\psi }_0{\partial }_y(u_k-v_k)-(u_k-v_k){\partial }_y{\psi }_0\right]$$

$$\Psi (x,y,t)={\Psi }_1(x,t){\chi }_1(y)+{\Psi }_2(x,t){\chi }_2(y).$$

$$\mathrm{i}\text{partialderivative}t{\Psi }_1(x,t)=\left(-\frac{\hslash }{2m}{\partial }_x^2+{\tilde{V}}_x+g_{1D}\text{absolutevalue}{{\Psi }_1}^2\right){\Psi }_1$$

$$-\left(K+F(\text{absolutevalue}{{\Psi }_1}^2+\text{absolutevalue}{{\Psi }_2}^2)\right){\Psi }_2-F({\Psi }_1^{\ast }{\Psi }_2+{\Psi }_2^{\ast }{\Psi }_1){\Psi }_1$$

$$+M\left(\text{absolutevalue}{{\Psi }_2}^2{\Psi }_1+({\Psi }_1^{\ast }{\Psi }_2+{\Psi }_2^{\ast }{\Psi }_1){\Psi }_2\right),$$

$$g_{1D}$$

$$K$$

$$F$$

$$M$$

$$\mathrm{i}\text{partialderivative}t{\Psi }_1(x,t)=\left(-\frac{\hslash }{2m}{\partial }_x^2+{\tilde{V}}_x+g_{1D}\text{absolutevalue}{{\Psi }_1}^2\right){\Psi }_1$$

$$-\left(K+F(\text{absolutevalue}{{\Psi }_1}^2+\text{absolutevalue}{{\Psi }_2}^2)\right){\Psi }_2,$$

Full text

I Introduction

Josephson effects in atomic Bose-Einstein condensates (BECs) receive considerable attention in experimental and theoretical studies as a prominent manifestation of quantum coherence on a macroscopic scale. These effects have been initially studied for coupled condensates in double-well traps Zapata1998 ; Smerzi1997 ; Raghavan1999 ; Albiez2005 ; Levy2007 ; LeBlanc2011 and coherently coupled spinor condensates PhysRevA.59.R31 ; PhysRevLett.105.204101 . In recent years Josephson effects were also observed and analyzed in more exotic systems such as fermionic superfluids Valtolina2015 , polariton condensates abbarchi2013macroscopic , spin-orbit coupled BECs Wang2018 ; PhysRevLett.120.120401 and condensates with attractive interaction Spagnolli2017 .

A common theoretical description of two coupled BECs is based on the bosonic Josephson equations. This model has been thoroughly studied theoretically Smerzi1997 ; Raghavan1999 . Various experiments confirm the predictions of this model, which include small-amplitude (plasma) oscillations, large-amplitude anharmonic oscillations and quantum self trapping Albiez2005 ; Levy2007 ; PhysRevLett.105.204101 . A formal derivation of this model relies on the two-mode description of the system, whereby each of the two condensates retains its density profile and a homogeneous phase. Such an approximation is valid if the coupling between the two BECs is much weaker than the energy required to create collective excitations inside each condensate. If this requirement is not fulfilled then internal collective excitations can be generated and influence the Josephson dynamics of the system. In order to describe the interference between internal and mutual collective motions of the condensates it is necessary to go beyond the simplified picture of the Josephson equations. Moreover the structure of collective excitations spectrum and consequently the results of such an interference may be considerably different depending on the geometric properties of the system. Existing theoretical studies address this question for certain specific types of collective excitations and trap geometries, such as phonon excitations in toroidal condensates Bidasyuk2016 , and higher modes of the harmonic trap PhysRevA.89.023614 . However, various aspects of collective excitations in coupled BECs and their relation to the tunneling transport have yet to be fully explored.

In the present work we analyze the collective excitations and the small-amplitude Josephson oscillations in a system of two highly anisotropic BECs. Two elongated condensates are placed in parallel to each other forming a so called long Josephson junction (LJJ). Such Josephson junctions were extensively studied for superconductors RevModPhys.51.101 ; PhysRevB.53.6622 ; PhysRev.164.538 and in atomic BECs mostly in the context of Josephson vortex dynamics PhysRevA.71.011601 ; Gil_Granados_2019 . Our interest in these systems is largely inspired by the experimental results of Ref. LeBlanc2011 where the populations of the wells were shown to oscillate at several distinct frequencies, indicating interference with internal collective excitations generated inside each condensate.

In order to describe the Josephson dynamics on the level of collective excitations in the system we need to relate the general physical picture of the Josephson equations with the Bogoliubov-de Gennes formalism. Following the developments in Refs. Danshita2005 ; Paraoanu2001 ; Burchianti2017 the Josephson plasma oscillations can be related to the lowest dipole mode in the two-well system. Along a similar line of arguments we show here how multiple Bogoliubov modes are involved in the low-energy Josephson dynamics of the anisotropic system.

We model the dynamics of the long bosonic Josephson junction in the mean-field framework of the two-dimensional Gross-Pitaevskii equation. The observed dynamical picture is complemented with an analysis of the collective excitation spectrum based on the Bogoliubov-de Gennes approach. We show how the specific structure of the Bogoliubov modes may lead to various physically relevant effects in the collective dynamics. In addition to enriched Josephson dynamics we identify for a certain range of barrier intensities unexpected quasi-degeneracy of the low-lying Bogoliubov modes with localization of the corresponding quasiparticles at the edges of the junction. Finally, we develop a simplified one-dimensional hydrodynamic model of the long bosonic Josephson junction, that is able to qualitatively reproduce all important features of the full simulations.

II Simulations of long Bose-Josephson junctions

First we briefly describe simulations of the dynamics of trapped condensed atoms in a bosonic LJJ. The trap potential is approximated by an anisotropic harmonic confinement in all spatial dimensions with a Gaussian barrier in $y$ direction,

[TABLE]

The barrier height $V_{b}$ is a tunable parameter to control the coupling between the two condensate clouds forming in this potential. This setup is patterned similarly to the experimental one described in Ref. LeBlanc2011 . We limit our considerations to values of $V_{b}$ larger than the chemical potential $\mu$ of the system such that classical hydrodynamic flow across the barrier is impossible, and the system remains in the tunneling regime.

The harmonic part of the potential is assumed to be cigar shaped ($\omega_{x}\ll\omega_{y,z}$). The barrier cuts the cigar into two parts along the $(x,z)$ plane parallel to the long axis of the trap. In such a configuration the trapping in $z$ direction remains harmonic. Disregarding the dynamics in this direction we may reduce the description of the system to two dimensions using a Gross-Pitaevskii equation (GPE)

[TABLE]

whith ${\bf\it r}=(x,y)$. The condensate wave function $\Psi\left({\bf\it r},t\right)$ is normalized to the total number of atoms in the system $\int\differential{{\bf\it r}}\absolutevalue{\Psi\left({\bf\it r},t\right)}^{2}=N$, which is a conserved quantity. The non-linear interaction strength in two dimensions is then given by Bao2003

[TABLE]

where $a$ is the s-wave scattering length and $m$ the mass of a condensate atom. We here consider ${}^{87}\mathrm{Rb}$ atoms with $a=5.819~{}\mathrm{nm}$ and $m=86.91~{}\mathrm{u}$ and choose $\omega_{x}=2\pi\cdot 1~{}\mathrm{Hz}$, $\omega_{y}=\omega_{z}=2\pi\cdot 50~{}\mathrm{Hz}$, a barrier $1/e^{2}$ half-width $\lambda=3.5~{}\mu\mathrm{m}$, and a total number of atoms $N=2.0\cdot 10^{4}$. However, most of our results should be qualitatively similar for any bosonic LJJ in the tunneling regime. For the calculations presented in this section we use a barrier amplitude $V_{b}/h=375~{}\mathrm{Hz}$. The ground state stationary solution of the GPE (2) shows the shape of two nearly one-dimensional parallel condensates (see Fig. 1). We find a chemical potential $\mu/h=273~{}\mathrm{Hz}$ and consequently $V_{b}/\mu=1.37$. Note that the chemical potential is still considerably larger than the trapping frequencies. Therefore, the mean-field treatment based on the GPE (2) is justified.

Josephson effects in bosonic systems lead to particle transfer between the wells. The relative population imbalance $Z(t)$ is the dynamical quantity which we use to quantify these effects

[TABLE]

where $N_{1}(t)$ and $N_{2}(t)$ denote the number of atoms in each well ($N_{1}(t)+N_{2}(t)=N$). These populations are related to the solutions $\Psi\left({\bf\it r},t\right)$ of the GPE. They are defined as integrals of the particle density over each well

[TABLE]

where $j=1,2$ and $A_{j}$ is the area subtended by well $j$.

The second dynamical quantity of interest is the phase difference between the condensates

[TABLE]

where each phase is defined by averaging over corresponding well

[TABLE]

One may show that this quantity is canonically conjugate to the population imbalance.

The time evolution of the population imbalance and relative phase of two BECs is often accurately described by the Josephson equations Smerzi1997 ; Raghavan1999

[TABLE]

Here, $\Lambda$ describes the intra-well interactions and $J$ characterizes the coupling between the two condensates. These quantities may be expressed in terms of the solutions of the GPE. A derivation of the Josephson equations with explicit definitions for $\Lambda$ and $J$ based on the two-mode model may be found in the Appendix. The ratio between $\Lambda$ and $J$ defines the physical region of validity of this model. It is expected to be valid in the region $1\ll\Lambda/J\ll N^{2}$ which is also called ‘Josephson regime’ PitaevskiiStringari2016 ; Leggett2001 . The value of this ratio is strongly dependent on the barrier height. For the range of configurations considered in the present study this ratio ranges from $\Lambda/J\approx 7$ for ($V_{b}/\mu=1$) to $\Lambda/J\approx 9800$ for ($V_{b}/\mu=2$) which is well inside the Josephson regime. The Josephson equations (8,9) have been successfully applied to describe various double-well experiments Albiez2005 ; Levy2007 ; PhysRevLett.105.204101 . Here we will use these equations as a basic reference and analyze their limitations when applied to anisotropic condensates.

In the case of a small initial population imbalance $\absolutevalue{Z_{0}}\ll 1$ or initial phase difference $\varphi_{0}\ll\pi/2$, the equations (8,9) reduce to a harmonic oscillator equation with the frequency

[TABLE]

known as plasma frequency. Thus the model predicts harmonic oscillations for small initial imbalances or phase differences. The value of the plasma frequency marks another important limitation of the model based on the system (8,9). In order for the Josephson oscillations to be decoupled from the internal oscillations inside each condensate it is necessary that $\omega_{p}\ll\omega_{x,y,z}$ PitaevskiiStringari2016 . This condition is not fulfilled for our system and therefore we can expect the effects of intra-well collective excitations to influence the Josephson dynamics.

For larger initial population imbalance or phase difference the Josephson equations predict the development of anharmonic large amplitude oscillations and the transition to the self-trapped regime at the critical population imbalance $Z_{cr}=2\sqrt{J/\Lambda(1-J/\Lambda)}$ Raghavan1999 .

In order to study the near-equilibrium dynamics of the junction we perform numerical simulations using the GPE (2). The initial state is prepared close to the ground state with an initial imbalance $Z_{0}=0.01$ inside of the plasma oscillation regime. In a real experiment such an initialization can be realized with an additional linear offset potential that is switched off at $t=0$. In our simulations the system evolves freely for $5~{}\mathrm{s}$ and the population imbalance $Z(t)$ is determined at equally spaced times $t_{m}$. The observed time series $Z(t)$ is shown in the upper panel of Fig. 2 revealing clearly population oscillations not determined by a single frequency.

We analyze the spectral properties of the oscillations by taking the discrete Fourier transform (DFT) of the time series $Z(t_{m})$

[TABLE]

where $M$ is the total number of time steps in the simulation. The Fourier-spectrum $\tilde{Z}(\omega)$ shows several well separated peaks (see Fig. 2, lower panel). We consider then $Z(t)$ to be a superposition of several harmonic oscillations. We extract the frequencies of these distinct oscillations simply by locating the maxima in the spectrum. It is worth noticing, that none of the detected frequencies can be related to the plasma frequency predicted by the two-mode model. In order to investigate this finding further, we vary the initial population imbalance and initial phase difference. The results plotted in Fig. 3 show that multiple peaks in the spectrum can be observed at arbitrarily small initial imbalance or phase difference. Their positions remain stable in the regions $Z_{0}\ll Z_{cr}$ and $\varphi_{0}\ll\pi/2$. Furthermore, the frequencies get smaller with increasing $Z_{0}$ or $\varphi_{0}$ as is similarly the case for the plasma frequency Raghavan1999 ; Martinez2018 . Finally, in the regime of large-amplitude oscillations distinct frequencies can hardly be distinguished and the DFT spectrum does not provide useful information. In the rest of the paper we will only consider the regime of small oscillations with $Z_{0}=0.01$ and $\varphi_{0}=0$, where the frequencies are well defined and stable.

III Bogoliubov-de-Gennes formalism

The results of the previous section show that several modes with different frequencies contribute to the population oscillations in the system. Furthermore, the number of these modes and their frequency is almost independent of the oscillation amplitude. This indicates that non-linear mode mixing effects are negligible in the considered regime. Nevertheless, such a multi-mode population oscillation spectrum is a clear indication that intra-well collective excitations influence the tunneling dynamics. In order to further analyze how different collective excitations influence the Josephson dynamics we employ the Bogoliubov-de Gennes formalism PitaevskiiStringari2016 ; Danshita2005 . We first write the condensate wave function in the form

[TABLE]

where $\psi_{0}$ is the stationary ground state, $\mu$ is the corresponding chemical potential and $\delta\psi\left({\bf\it r},t\right)$ is a perturbation of the form

[TABLE]

with the additional condition that the perturbation is small ($c_{k}\ll\sqrt{N}$). After inserting this ansatz into the GPE (2) and considering only the terms linear in $c_{k}$ we obtain the familiar system of Bogoliubov-de Gennes equations

[TABLE]

where

[TABLE]

The usual normalization condition for the Bogoliubov modes reads

[TABLE]

This linear system of equations can be solved numerically in order to obtain eigenmodes and their corresponding frequencies $\omega_{k}$. We solve the system (14) for different values of the barrier height $V_{b}$. The calculated spectrum as a function of $V_{b}$ is presented in Fig. 4. We also perform full GPE simulations with barrier heights in the same range to compare the frequencies extracted from the time series $Z(t)$ with the Bogolibov spectrum.

As one can see in Fig. 4 only some of the Bogoliubov modes correlate with the spectrum of population oscillations. In order to understand this fact let us see how the linear excitations of the form (13) enter into the dynamics of the population imbalance of the two-well system. To this end we introduce the following operator representation of the population imbalance

[TABLE]

which is valid for any wave function $\Psi$. We can now insert here the wave function of the form (12,13). We keep only terms up to first order in $c_{k}$ and assume (without loss of generality) that $\psi_{0}\left({\bf\it r}\right)$, $u_{k}\left({\bf\it r}\right)$ and $v_{k}\left({\bf\it r}\right)$ are purely real functions. Taking also into account that $\langle\psi_{0}|\hat{Z}|\psi_{0}\rangle=0$ the final expression yields

[TABLE]

Here

[TABLE]

is the excitation amplitude of the $k$’th Bogoliubov mode by an operator $\hat{Z}$. In other words the set of coefficients $D_{k}$ essentially characterizes the condensate response to the small perturbation proportional to $\hat{Z}$ PitaevskiiStringari2016 ; PhysRevA.56.587 .

The equation (16) shows that a populated Bogoliubov mode with index $k$ results in a transfer of population with the frequency $\omega_{k}$, but only if the coefficient $D_{k}$ is non-zero. To identify which Bogoliubov modes may correspond to non-zero values of $D_{k}$ we should consider symmetry properties of these modes. Both the potential (1) and the ground state are symmetric with respect to axes reflections in longitudinal ($x$) and transverse ($y$) directions. Therefore the solutions of (14) are expected to possess certain reflection symmetries in these directions as well. In Fig. 4 we distinguish three subsets of levels (plotted with different colors) based on their symmetry properties. These states are (A) transversely symmetric, (B) transversely antisymmetric but longitudinally symmetric and (C) antisymmetric in both directions. In the following we will refer to these states as A, B and C for shorter notation.

From the structure of Eq. (17) one can predict that $D_{k}$ may be non-zero only if $u_{k}+v_{k}$ is symmetric in $x$ and antisymmetric in $y$. Consequently the states of type B are the only subset of modes that can possibly contribute to the Josephson dynamics identified by $Z(t)$. This argument is confirmed by the results in Fig. 4 showing a nearly perfect match of the spectra extracted from GPE dynamics with the Bogoliubov modes of type B.

The above symmetry arguments provide necessary conditions for any of the collective excitations to be traceable in the population imbalance $Z(t)$. However, these conditions are by no means sufficient. The judgement on the number of non-zero amplitudes $D_{k}$ can be done only by calculating them numerically or measuring them. Such calculations show that there is always only a finite number of essentially non-zero values in this set. In Fig. 5 we show the values of $D_{k}$ and $\omega_{k}$ corresponding to the lowest type B excitations for the same barrier heights used in the previous section. Apparently only the five lowest excitations show amplitudes $D_{k}$ significantly larger than zero. One may also notice a similarity of Fig. 5 and the spectrum $\tilde{Z}(\omega)$ presented in the lower panel of Fig. 2. This shows that the spectrum of population oscillations can in general be predicted by analyzing the Bogoliubov modes of the system. This also means that the number of observable modes is a property of the system, and does not depend on the initial perturbation imposed in the numerical simulation.

IV General properties of the Bogoliubov spectrum

We now investigate further the complex structure of the calculated excitation spectrum. The first obvious observation is that the barrier affects differently the excitations of different symmetries. In particular, the excitations of type A which are symmetric in the transverse direction are practically unaffected by the changes of the barrier height. These solutions are node-less in $y$ and represent purely longitudinal excitations. They originate from the collective excitations of a single one-dimensional BEC and are well approximated by a simple analytical formula based on hydrodynamic approximation Menotti2002 (see Fig. 9 for direct comparison)

[TABLE]

The modes of types B and C are strongly affected by the barrier height but asymptotically they converge to the symmetric counterparts and are expected to match them exactly in the limit of two completely uncoupled condensates. In Fig. 6 we show the level spacings $\Delta\omega=\omega_{k}^{(B,C)}-\omega_{k}^{(A)}$ between the levels of type A and corresponding levels of types B and C. We see that in the high barrier limit all level spacings decay exponentially with growing $V_{b}/\mu$. Interestingly, all the decay rates are related to parameters defined in the two-mode model. The lowest pair of levels converge at the rate of plasma frequency $\omega_{p}$ but all the other pairs converge at the same rate as the coupling parameter $J$. In the limit of a very high barrier $V_{b}/\mu>1.9$ where $\omega_{p}<\omega_{x}$, higher modes disappear from $Z(t)$ dependence and the lowest oscillation frequency converges to the estimated value of the plasma frequency $\omega_{p}$ (see Fig. 4). This indicates that validity of the Josephson equations (8,9) is restored in this region.

Levels of types B and C also show an unexpected degeneracy in the region of intermediate barrier heights. The spacings between these levels are shown in Fig. 7. We show them in both linear and logarithmic scales to highlight that for all of the affected levels the minimal level spacing is observed at the same barrier height around $V_{b}/\mu\approx 1$. This quasi-degeneracy arises due to coupling between two condensates.

In order to understand the origin of this degeneracy we look into the distribution of the tunneling current initiated by these excitations. The transverse tunneling current density can be written in the following form

[TABLE]

Inserting here the wave function (12,13) we get

[TABLE]

where

[TABLE]

is the characteristic current density associated with $k$’th Bogoliubov mode. These current densities for the lowest Bogoliubov modes of types B and C are shown in Fig. 8. We see that in the region of quasi-degeneracy the tunneling current in these modes is localized at the edges of the junction and two edge-localized excitations basically decouple from each other. A similar effect was also described for the superconducting LJJs PhysRev.164.538 and attributed to the Meissner effect in the junction. One may also see from Fig. 4, that these edge-localized states are always involved in the dynamics of population imbalance $Z(t)$ (the crosses cover all regions of quasi-degeneracy).

V Effective one-dimensional model of elementary excitations

So far all of our descriptions for the bosonic LJJ have been in two spatial dimensions. However, the geometry of the system suggests that we can consider it as two coupled one-dimensional systems. For the symmetric modes the one-dimensional treatment of the system is quite simple and leads to the analytical expression (18). However, for the spectrum of antisymmetric modes to the best of our knowledge no analytical expressions exist. The model developed in Ref. Abad2013 describes dispersion relations of coupled infinite homogeneous 1D condensates, so while it captures some essential properties of the spectra, it can not reproduce the observed level degeneracy and edge-localized states. In the following we develop a simplified one dimensional hydrodynamic model that is still able to reproduce qualitatively the behavior of antisymmetric Bogoliubov modes. Our aim is to reduce the system to two coupled one-dimensional condensates. The reduction scheme that we propose is similar in spirit to the usual two-mode model. It relies on the additional assumption that the dimensions can be separated. So the long axis of the trap $x$ does not couple to the short axis $y$. In each well the wave functions can be written as a product state and the total wave function reads

[TABLE]

We assume the solutions in $y$ dimension $\chi_{1,2}(y)$ to be time independent. Additionally we define $\chi_{1}$ and $\chi_{2}$ such that they are real, orthogonal and normalized to unity. Our trapping potential is obviously separable $V(x,y)/\hbar=\tilde{V}_{x}(x)+\tilde{V}_{y}(y)$. Inserting this ansatz into the GPE (2) leads to a system of two coupled one-dimensional equations

[TABLE]

where for the second component the indices are interchanged. The coefficients that appear in these equations are defined as follows

[TABLE]

In the Josephson tunneling regime the overlap between functions $\chi_{1}$ and $\chi_{2}$ is small and we can safely assume $M\ll F\ll g_{1D}$. Then the equations simplify to

[TABLE]

These equations permit a one-dimensional description of the double-well system. They are similar to the commonly used model for two-component condensates with coherent coupling. However they contain additional terms proportional to $F$ which in general cannot be neglected. A more detailed study of the effects of these additional terms would require a separate study. We only mention here that $K$ is in general not sign-definite and neglecting the nonlinear terms proportional to $F$ may lead to unphysical results. This problem is usually rectified by taking $K$ as a phenomenological positive definite parameter or imposing additional approximations (see e.g. Ref. PhysRevA.81.025602 for a discussion and references).

In order to further simplify our model, we recall that for the parameters considered in the present work the 1D system is deeply in the Thomas-Fermi regime. This suggests that a hydrodynamic approximation may be valid to describe the system. We rewrite the equations (V) with the usual Madelung decomposition $\Psi_{1}(x,t)=\sqrt{\rho_{1}}\mathrm{e}^{\mathrm{i}\theta_{1}}$, where $\rho_{1}(x,t)$ is the density and $\theta_{1}(x,t)$ is the phase, which gives

[TABLE]

where we neglected the quantum pressure terms of the form $\partial_{x}(\sqrt{\rho})/\sqrt{\rho}$. The equations of the other component have interchanged indices.

To study the antisymmetric low energy excitations we consider a stationary state for each well and we add a perturbation with an opposite sign in each well:

[TABLE]

where $\tilde{\mu}=\mu/\hbar$. The Thomas-Fermi ground state density is

[TABLE]

in the region where $\tilde{\mu}+K>V_{x}$ and zero otherwise. With this ansatz we linearize the coupled equations,

[TABLE]

Combining the two equations we get one second-order equation for $\delta\rho$

[TABLE]

The term proportional to $(\partial_{x}\rho_{0})/\rho_{0}$ can be safely neglected in the Thomas-Fermi regime. We finally assume the solutions to be periodic in time $\delta\rho(x,t)=\cos(\omega t)\delta\rho(x)$, which leads to the following equation

[TABLE]

This equation provides a hydrodynamic approximation for the transversely antisymmetric collective excitations of the system. Using the explicit expression for the Thomas-Fermi ground state density (28) this equation can be easily recast in terms of the chemical potential and the trap potential, which makes it numerically tractable. In some limiting cases the equation (V) can be solved analytically. In particular, in the limit of uncoupled condensates ($K\rightarrow 0,F\rightarrow 0$) it is straightforward to show that the spectrum (18) is reproduced as is expected in this limit. Also if one considers the excitations of the same shape as the ground state ($\delta\rho(x,t)=Z(t)\rho_{0}$) then the result of the two-mode model is recovered and after integrating out the spatial dimension we get a single frequency $\omega=\omega_{p}$ identical to Eq. (10).

In figure 9 we compare numerical solutions of equation (V) with the Bogoliubov spectrum of the two-dimensional system. The parameters $g_{1D}$, $K$ and $F$ were estimated based on the ground state GPE solutions. The two approaches show reasonable agreement, which is better in the high barrier region. Specific features of the spectrum, such as formation of quasi-degenerate pairs of levels is also reproduced in the developed hydrodynamic approximation. This result justifies the general validity of this approach and applicability of the derived equation (V).

Using the ansatz (V) we can also write the population imbalance in the form analogous to (16)

[TABLE]

where the coefficients $D_{k}$ have a similar meaning to the excitation amplitudes discussed in the previous section and are defined by the solutions of Eq. (V) as follows

[TABLE]

The values of these coefficients are shown in Fig. 10. We see that they also qualitatively reproduce the results in Figs. 5 and 2. A number of additional zero-valued coefficients on the figure correspond to excitations of type C, which are also included in the spectrum of Eq. (V).

VI Conclusions

In the present work we investigated the low-energy dynamics of long bosonic Josephson junctions. In the simulations based on the time-dependent Gross-Pitaevskii equation we observe oscillations of the population imbalance consisting of multiple well-defined frequencies, that persist even for arbitrarily small initial perturbation.

In order to understand this behavior we analyze the spectrum of elementary excitations obtained from the Bogoliubov-de Gennes equations. By analyzing the condensate response to the imposed population imbalance we find the Bogoliubov modes that can contribute to the population transfer. This allows us to explain and predict the multi-mode spectrum of population oscillations observed in the dynamical simulations.

From the general structure of the excitations spectrum in the region of intermediate barrier heights we discover the development of quasi-degeneracy of low-lying levels. We connect this phenomenon with localization of corresponding modes at the edges of the junction. The geometry of LJJ allows such edge-localized excitations to decouple and form degenerate pairs. Quite interestingly, these edge-localized oscillations are also always observed in the dynamics of the population imbalance.

We also developed an effective one-dimensional model of two coupled condensates, that combines a hydrodynamic approximation in the longitudinal direction and quantum tunneling in the transverse direction of the barrier. This simplified model is shown to qualitatively reproduce all essential features of the excitation spectrum obtained using the Bogolubov-de-Gennes analyis: The excitation frequencies are in good agreement in a wide range of barrier heights. The model also correctly predicts multiple frequencies for the population dynamics of the condensate.

Appendix A The two-mode model

The two-mode model is a common approach to describe Josephson effects in a bosonic double-well system. We review it here to introduce the concepts and establish a notation. The central idea of the two-mode model is that the total wave function can be described with only two modes. Each one is (roughly) localized in one well and has one coherent phase

[TABLE]

$N_{1,2}$ is the occupation, $\theta_{1,2}$ is the phase and $\phi_{1,2}$ are the time-independent wave functions of each mode, which have unit norm and are orthogonal. One can obtain them with the standard ansatz of combining the ground state solution $\phi_{g}\left({\bf\it r}\right)$ with the first excited state that has a node in the barrier region $\phi_{ex}\left({\bf\it r}\right)$:

[TABLE]

The time-dependent variables of the two-mode model are the population imbalance $Z(t)$ and the relative phase $\varphi(t)$,

[TABLE]

where $N$ is the total number of atoms with $N_{1}+N_{2}=N$. Putting the ansatz (eqn. 33) into the GPE leads to the well-known two-mode model equations

[TABLE]

where the parameters are defined as

[TABLE]

Here $\Lambda$ is the on-site interaction, $J$ describes the coupling and $I$ is called interaction tunneling (see Ananikian2006 ). In case of small oscillations $\absolutevalue{Z}\ll 1$ and $\varphi\ll\pi/2$ the equations reduce to a harmonic oscillator equation with the frequency

[TABLE]

This is called plasma frequency and it characterizes the population transfer between the wells. In most configurations $\Lambda\gg J\gg I$, so $I$ can be neglected and one obtains $\omega_{p}=\sqrt{J(\Lambda+J)}$ as quoted in Eq. 10.

Results of experiments show that the two-mode model provides better predictions if the parameters are modified. Especially $\Lambda$ is usually too large with the definition 38. This is due to the approximation that the wave functions $\phi_{1}\left({\bf\it r}\right)$ and $\phi_{2}\left({\bf\it r}\right)$ are time-independent. Relaxing that condition by allowing the wave functions to depend on the population of the respective wells like this $\phi_{j}\left({\bf\it r}\right)\rightarrow\phi_{j}({\bf\it r},N_{j}(t))$ for $j=1,2$, provides a more accurate representation of the physics. This ansatz was implicitly introduced in Refs. Jezek2013 ; Nigro2017 and it was termed effective two-mode model. We provide a slightly different (and to our opinion more transparent) derivation. The ansatz for the total wave function is

[TABLE]

Here $N_{1},N_{2},\theta_{1}$ and $\theta_{2}$ are time-dependent. This can be seen as an adiabatic approximation because the wave function for one well is assumed to be in the ground state at the corresponding population. If we put this ansatz into the GPE (2) we obtain the equations

[TABLE]

where $J,\tilde{J},I,\Delta E$ and $\Lambda$ are functions of $Z$

[TABLE]

because $\phi_{1}({\bf\it r},Z)$ and $\phi_{2}({\bf\it r},Z)$ depend on it. For $Z=0$ the two-mode model parameters are recovered. The relation $\Lambda\gg J\gg I$ still holds, so $I$ can be neglected safely. Also the changes in $J$ and $\tilde{J}$ are expected to be small so $J(Z)\approx J(0)=J$ and $\tilde{J}(Z)\approx\tilde{J}(0)=0$. With that the equations of motion are

[TABLE]

The key of the effective model is to approximate the functions $\Delta E(Z)$ and $\Lambda(Z)$. However at this point it is not clear if and how they are related. For that let us introduce the local chemical potential of each well

[TABLE]

The difference between the local chemical potentials is indeed

[TABLE]

This relation clearly shows that $\Delta E$ and $\Lambda$ cannot be treated separately. A Taylor approximation at $Z=0$ gives

[TABLE]

The first term has to vanish due to the symmetry of the double well. The linear term can be changed to contain the appropriate populations $N_{1}$ and $N_{2}$ for each well

[TABLE]

Again the symmetry of the double-well causes the derivatives of $\mu_{1}$ and $\mu_{2}$ to be the same, so we finally get

[TABLE]

where

[TABLE]

This derivative can be obtained by calculating the ground state $\phi_{1}$ for several different populations $N_{1}$ and putting it in the functional $\mu_{1}$. We use this definition of the parameter $\Lambda$ for all two-mode model estimates presented in the present work.