A Stability calculation

# Quasiparticles of widely tuneable inertial mass: The dispersion relation of atomic Josephson vortices and related solitary waves

Quasiparticles of widely tuneable inertial mass: The dispersion relation of atomic Josephson vortices and related solitary waves

S.S. Shamailov1, J. Brand1*

1 Dodd-Walls Centre for Photonic and Quantum Technologies, Centre for Theoretical Chemistry and Physics and New Zealand Institute for Advanced Study, Massey University, Private Bag 102904 NSMC, Auckland 0745, New Zealand

* J.Brand@Massey.ac.nz

March 6, 2018

## Abstract

Superconducting Josephson vortices have direct analogues in ultracold-atom physics as solitary-wave excitations of two-component superfluid Bose gases with linear coupling. Here we numerically extend the zero-velocity Josephson vortex solutions of the coupled Gross-Pitaevskii equations to non-zero velocities, thus obtaining the full dispersion relation. The inertial mass of the Josephson vortex obtained from the dispersion relation depends on the strength of linear coupling and has a simple pole divergence at a critical value where it changes sign while assuming large absolute values. Additional low-velocity quasiparticles with negative inertial mass emerge at finite momentum that are reminiscent of a dark soliton in one component with counter-flow in the other. In the limit of small linear coupling we compare the Josephson vortex solutions to sine-Gordon solitons and show that the correspondence between them is asymptotic, but significant differences appear at finite values of the coupling constant. Finally, for unequal and non-zero self- and cross-component nonlinearities, we find a new solitary-wave excitation branch. In its presence, both dark solitons and Josephson vortices are dynamically stable while the new excitations are unstable.

## 1 Introduction

The concept of inertial (or effective) mass [1] is commonly used in condensed matter physics: it captures the response of a quasiparticle in an interacting system to an applied force, encapsulating the emergent Newton’s equations of quasiparticle dynamics. Atomic Bose-Einstein condensates (BECs) [2, 3] provide a platform for the emulation of other quantum-many-body systems under highly controllable conditions [4, 5, 6]. The possibility of adjusting the inertial mass of localized excitations in BECs by tuning experimental parameters could potentially open the way to interesting applications. Solitary waves with large inertial mass, identified as solitonic vortices [7, 8], have recently been observed in superfluid Fermi gases [9, 10] and BECs [11]. Since the inertial mass of a solitonic vortex depends on the mass density of the superfluid and the length scale of the transverse confinement [10, 12, 13] it is tuneable through the trap geometry within certain limits, but it cannot change sign. In this work we study the properties of solitary waves in a one-dimensional two-component atomic superfluid where an adjustable linear coupling between the two components provides a convenient control parameter that can be used to tune the inverse inertial mass through zero, and thus achieve large positive and negative inertial masses.

Quasi-one-dimensional BECs have been prepared experimentally almost two decades ago [14], and more recently, two coherently-coupled one-dimensional BECs have been demonstrated [15, 16, 17, 18]. Dark solitons [19] are localized density depletions with a phase drop across them that propagate at constant speed, preserving their shape. They have been observed in single-component BECs [20, 14, 21, 22] and require quasi-one-dimensional confinement to be stable and long-lived [23, 7, 24]. Josephson vortices are known to occur as quantised magnetic flux lines corresponding to a vortex of the superconducting order parameter inside a long Josephson junction between two bulk superconductors [25, 26], and their theoretical description is often reduced to a sine-Gordon equation [27]. Josephson vortices in atomic superfluids were first discussed as domain walls of the relative phase in two-component BECs with a weak coherent coupling between the components [28] and later as vortices that have entered the extended barrier region of a single-component BEC in a double-well geometry [29, 30] (see Fig. 1). While observation of regular vortices is by now common place in BECs [31, 32] and superfluid Fermi gases [33], proposals for the observation of Josephson vortices in BEC were made in [34, 35, 36]. Very recently, spontaneously created Josephson vortices were identified through interference patterns in linearly coupled BECs [37].

The model of linearly-coupled two-component BECs describes two different physical realizations: a spinor BEC with two spin components and coherent coupling achieved through radio-frequency or microwave radiation driving a hyperfine transition [28] or a single-component BEC in a double-well potential [38, 39, 40, 6, 36]. Both options are illiustrated in Fig. 1. Theoretical considerations ranged from testing the Kibble-Zurek mechanism in a ring geometry [36], to modeling the decay of an unstable vacuum to a universe with structure [6, 41], to metastable domain walls [28, 38], to the dynamical response to periodic modulation [39], and to tunneling quenches leading to breather modes forming out of quantum fluctuations [40]. Out of these studies, Refs. [28, 39, 40, 6] reduced the model to the integrable sine-Gordon case, assumed to be applicable in the small tunneling limit, in order to obtain their main results. The experiment [37] has assumed likewise. Testing the validity of this approximation is part of the work presented in the current paper.

After the work presented in this paper was completed we became aware of the related recent work by Qu et al. [42], which considers Josephson vortices (called magnetic solitons in [42]) in a regime of weak tunneling and almost equal self- and cross-nonlinearities, where the particle number density is approximately constant1. They find analytical and numerical results for solitary wave solutions, their dispersion relations and their dynamics under harmonic trapping. Their findings are consistent with our own results, which cover wider and more general parameter regimes. Furthermore, we have independently obtained similar predictions for the dynamics of the Josephson vortex in a trapped condensate [44].

In this work we consider several families of solitary nonlinear waves in quasi-one-dimensional linearly-coupled two-component BECs. Son and Stephanov [28] discussed the existence of a domain-wall of the relative phase in this system, which was later called a Bose Josephson vortex by Kaurov and Kuklov [29, 30]. The latter authors found exact stationary solutions to the coupled Gross-Pitaevskii equations and discussed the bifurcation of the dark soliton solution of this model, which is stable above a critical strength of the linear coupling, into stable Josephson vortices which coexist at smaller coupling strength with the now unstable dark soliton. Qadir et al[45] added a detailed stability analysis and approximated the properties of moving Josephson vortices for small velocities. Exact solutions and the full dispersion relation for moving Josephson vortices have so far not been available and approximate dispersion relations have only been found for the regime of weak coherent coupling and nearly-equal nonlinear interactions [42]. Here we present numerical solutions for the complete dispersion relation of solitary wave solutions including moving Josephson vortices (already used in Ref. [43] to simulate collisions of Josephson vortices). In addition to the stable Josephson-vortex dispersion, a new unstable branch of solitary waves arises at finite cross-component nonlinear interactions where both dark solitons and Josephson vortices are stable.

We show that there is a critical value of the linear coupling at which the Josephson vortex dispersion relation changes from having a single maximum to having three (local) extrema: a maximum, minimum and another maximum. At this point the inertial mass of the Josephson vortex at the center of the dispersion relation changes sign and diverges to on either side of the bifurcation. The two maxima appearing at smaller coupling strength correspond to small-velocity solitary waves with negative inertial mass that resemble a dark soliton in one component with a back-flow current in the other.

In the small coupling limit, we test whether the central part of the Gross-Pitaevskii Josephson vortex dispersion relation approaches the sine-Gordon dispersion relation. The latter can be described by only two parameters – the “mass” and the “speed of light”. The inertial mass of the Josephson vortex approaches the sine-Gordon “mass” parameter quite rapidly as the tunneling is decreased, but the “speed of light” does not – the approach to the common value at zero tunneling occurs with different slopes.

The paper is structured as follows. Section 2 introduces the model equations, 3 defines useful observables for characterizing the solutions, and section 4 summarizes known analytical solutions and their properties. Next, section 5 explains how we numerically obtain solutions, which are visualized in section 6. Section 7 discusses the dispersion relations, while 8 discusses parameter regimes and a stability analysis of the solutions. Section 9 addresses the inertial mass and missing particle number of Josephson vortices. Section 10 presents a variational approximation to slow-moving Josephson vortices, yielding analytical expressions for the key numerical results. Section 11 summarizes some central results regarding the sine-Gordon equation, and 12 examines the validity of approximating the Gross-Pitaevskii equations with the sine-Gordon model. Discussion and conclusions are given in 13. Appendix A presents details of how we perform the stability calculation and appendix B derives the sine-Gordon equation from the Gross-Pitaevskii model, thus enabling a direct comparison of the two.

## 2 The model

We are considering a quasi-one-dimensional two-component Bose gas with identical boson mass for the two components and linear coupling in the Gross-Pitaevskii approximation with the energy functional

 W=∫{∑j=1,2[ℏ22m∣∣∂xΨj∣∣2+g2|Ψj|4−μ|Ψj|2]+gc|Ψ1|2|Ψ2|2−J(Ψ∗1Ψ2+Ψ∗2Ψ1)}dx, (1)

where is the order parameter of component , is the common chemical potential and for simplicity, the intra-component interaction constant [47] is also taken as common for the two components, but the cross-component interaction constant may be different. We assume that , which is the physical case when the coupling is achieved by tunnelling through a potential barrier [48]. In the general case, the phase of can be absorbed into the definition of . Two possible realisations of this model are illustrated in Fig. 1. Even though the homogeneous one-dimensional Bose gas never completely fully condenses, the Gross-Pitaevskii (mean-field) approximation is justified when , i.e. the Lieb-Linger parameter is small [49], where is the background particle density in each component.

The time evolution of the order parameter is described by the Gross-Pitaevskii equation, which can be formally obtained from :

 iℏ∂tΨ1 = −ℏ22m∂xxΨ1−μΨ1+g|Ψ1|2Ψ1+gc|Ψ2|2Ψ1−JΨ2, iℏ∂tΨ2 = −ℏ22m∂xxΨ2−μΨ2+g|Ψ2|2Ψ2+gc|Ψ1|2Ψ2−JΨ1. (2)

We are interested in solitary wave solutions that translate with a constant velocity . With the aim of adimensionalising the equations of motion we make the ansatz

 Ψj(x,t)=√μg+gcψ(z), (3)

where

 z =ξ−vsτ, ξ =√mμℏx, τ =μℏt, vs =√mμVs, (4)

are dimensionless variables. Further introducing the dimensionless linear and nonlinear coupling constants

 ν =Jμ, γ =g−gcg+gc, (5)

we obtain the dimensionless Gross-Pitaevskii equation for uniformly translating solutions

 −ivs∂zψ1 = −12∂zzψ1−ψ1+12(1+γ)|ψ1|2ψ1+12(1−γ)|ψ2|2ψ1−νψ2, −ivs∂zψ2 = −12∂zzψ2−ψ2+12(1+γ)|ψ2|2ψ2+12(1−γ)|ψ1|2ψ2−νψ1. (6)

We are looking for solutions to these equations that asymptotically tend to the stable constant background solution for [48]. The boundary conditions can thus be written as

 limz→±∞ψ1= limz→±∞ψ2, limz→±∞|ψ1|2= limz→±∞|ψ1|2=1+ν. (7)

Note that this leaves a complex phase undetermined at each end, and while Eqs. (6) are invariant under the change of an overall phase, the phase difference

 Δϕ=argψj(z→−∞)ψj(z→∞), (8)

bears physical significance.

## 3 Physical observables

Let us now define several useful quantities that shall be evaluated later on for the numerical solutions. The energy functional in dimensionless units is evaluated as

 E= √μm(g+gc)ℏμ2W= = L∫−Ldz ∑k=1,2{12|∂zψk|2−|ψk|2−νψ∗kψ3−k+14(1+γ)|ψk|4}+12(1−γ)|ψ1|2|ψ2|2. (9)

A key quantity is the excitation energy associated with the solitary wave

 Es=Esol−E0, (10)

where is the (dimensionless) energy of the constant background, and is the energy of the solitary wave solution.

For a localised solitary wave solution that heals to the constant background, is independent of the box size for sufficiently large . It will depend on the soliton velocity but there may be multiple solutions for each .

Another useful observable is the momentum, which is scaled by . The background solution has zero momentum and we introduce the following dimensionless observables:

 Ps = L∫−L(p1+p2) dz, ΔP = L∫−L(p1−p2) dz, pk = −i2[ψ∗kdψkdz−ψkdψ∗kdz], Pcf = 2(1+ν)Δϕ, Pc = Ps+Pcf, (11)

where is the physical momentum of the solitary wave with boundary conditions (2). The momentum difference between the two components indicates a degree of symmetry breaking. In the scenario where the two components are spatially separated it has the significance of an orbital angular momentum (see Fig. 1). The quantity is the momentum of the counter flow that has to be added in periodic boundary conditions (ring geometry) in order to compensate for the phase step . The canonical momentum is the momentum that the solitary wave excitation has with periodic boundary conditions but it is also significant for the open boundary conditions (2) due to the relation

 dEsdPc=vs. (12)

The canonical momentum provides a convenient way of parameterising the solitary-wave solutions and the relation of vs.  is known as the dispersion relation. Examples of dispersion relations that summarise the results of this work are presented in in section 7, Fig. 5, panels (a), (c), (e). In the framework of Landau’s quasiparticle picture, the dispersion relation determines the dynamics of the solitary waves in a slowly-changing environment as long as the nature of the solitary wave changes adiabatically such that the solitary wave solutions are well approximated by the stationary solutions of Eqs. (6) at any one time and the energy stored in the solitary wave is conserved [50, 12].

Of particular interest is the inertial mass of the quasiparticle, which is given by

 mI=dPcdvs=2dEsd(v2s) = (d2EsdP2c)−1. (13)

It is directly related to the measurable oscillation frequency of the solitary wave under the influence of a weak harmonic trap in the longitudinal direction with frequency by

 ω2zΩ2=mImP, (14)

where is the physical mass [51, 52, 12, 53]. At zero velocity, the physical mass is proportional to the particle number depletion of the solitary wave by 2. The (missing) particle number of the solitary wave is obtained by integrating the density and subtracting the background

 Nd=L∫−L[n1(z)+n2(z)−2n0] dz, (15)

where are the dimensionless particle densities in the two components and is the dimensionless background density. Note that Eq. (15) yields the particle number in terms of the reduced dimensionless quantities and is scaled by a factor .

## 4 Analytically known solitary-wave solutions

Several exact solutions of (6) are known. The lowest-energy constant solution is , which we shall refer to as the background. We remark that separates the miscible () and immiscible () phases of the system. In the miscible regime, corresponds to a degenerate mean-field ground state with undefined spin polarisation, with any leading to a background solution polarised along the -direction. For this work we consider the miscible regime with and a polarised ground state along the -direction obtained with , while analogous results hold for , where polarisation along the -direction is obtained.

### 4.1 Dark solitons

The coupled BECs system supports dark soliton solutions which satisfy and are given by [1]

 ψ=√1+ν−v2stanh[√1+ν−v2sz]+ivs. (16)

This corresponds to identical dark soliton solutions in each component with . The maximal velocity at which a dark soliton can travel is the Bogoliubov speed of sound of the system, . Note that the dark soliton solutions are independent of the cross-interaction parameter . The zero-velocity case is visualised in Fig. 2(a) & (b).

The soliton’s properties can be calculated by direct integration from the analytical solution and are well known [1]. In our dimensionless units they explicitly depend on the coupling parameter . The excitation energy

 Es=83(1+ν−v2s)3/2, (17)

takes a maximum value at zero velocity and vanishes at . The phase difference

 Δϕ=π−2tan−1[vs√1+ν−v2s] (18)

is for stationary solitons and reaches the extremal values and at the limiting velocities . The missing particle number of the dark soliton evaluates to

 Nd=−4√1+ν−v2s, (19)

and the momentum difference vanishes (). The velocity dependence of the dark soliton’s properties is shown in section 7, Figs. 68 as green lines.

The canonical momentum of the dark soliton is

 Pc=2π(1+ν)−4vs√1+ν−v2s−4(1+ν)tan−1(vs√1+ν−v2s), (20)

varying in the interval . The inertial mass evaluates to .

### 4.2 Stationary Josephson vortex

The stationary Josephson vortex is found as a complex solution of Eq. (6), which breaks the symmetry between the two components [28, 29]. It is given by , with

 ψ=√1+νtanh(2√νz)+i√1−3ν sech(2√νz), (21)

and only exists for . The parameter value marks a bifurcation point where the Josephson vortex solution becomes identical to the dark soliton solution (16). For two degenerate solutions are obtained from and , which can be interpreted as vortices of opposite circulation. The vortex nature is most clearly seen in the double-well potential scenario of Fig. 1(a), where a phase singularity sits in the middle of the double-well barrier at . Indeed, tracing the phase along the direction in and along in amounts to a total phase change of if the origin is included; note that both components have equal phase far away from the Josphson vortex. This can be seen following the phase profiles shown in Fig. 2 (c) & (d).

The energy and momentum difference for the Josephson vortex at are

 Es =83√ν(3−ν), ΔP =∓2π√1+ν√1−3ν, (22)

where the stands for the Josephson vortex (anti-vortex ).

### 4.3 Manakov solitons

When the cross- and intra-component nonlinearities are equally strong (i.e. ), the coupled equations (6) can be mapped on to the integrable vector non-linear Schrödinger equation known as the Manakov system [55, 56]. In this limit, a whole family of solutions can be found analytically [57]. Defining , we re-write equations (6) for the new variables:

 −ivs∂zχk=−12∂zzχk−(1±ν)χk+12(|χ1|2+|χ2|2)χk, (23)

where the two different signs in front of are to be taken with the two different indices, . A trial solution of the form [57]

 χ1 = αi+βtanh(ηz), χ2 = δsech(ηz)eiεz (24)

is found to satisfy (23) if the parameters are given by

 α = √1+ν2νvs, β = √(4ν−v2s)(1+ν)2ν, η = √4ν−v2s, δ = √(4ν−v2s)(1−3ν)2ν, ε = vs. (25)

In fact, may be multiplied by an arbitrary phase factor, , and the resulting solution still satisfies the differential equations. Transforming back to the -fields gives

 ψ1,2=1√2[αi+βtanh(ηz)±eiθδsech(ηz)eiεz]. (26)

Notice that for the parameters in (25) to be real (and the solution to be non-trivial) we need and . The phase angle remains a free parameter, indicating the large degeneracy of these solutions. For and the solution (26) reduces to the stationary Josephson vortex (anti-vortex). We will refer to the family of solutions (26) as the Manakov solutions, even though the presence of the linear coupling provides a point of difference to the solutions of the original Manakov system.

As for the dark soliton, it is possible to calculate all the quantities of interest for the Manakov solutions at analytically: the excitation energy, angular momentum, phase difference, canonical momentum, and missing particle number are

 Es = 4√4ν−v2s[23(4ν−v2s)−(3ν−1)], ΔP = 2π√1+ν√1−3ν sech(πvs2√4ν−v2s)sinθ, Δϕ = π−2tan−1[vs√4ν−v2s], Pc = 2π(1+ν)−4{vs√4ν−v2s+(1+ν)tan−1[vs√4ν−v2s]}, Nd = −4√4ν−v2s. (27)

The inertial mass evaluates to and is independent of velocity. Note that the limits of are the same as for dark solitons. The Manakov solitons are illustrated in Fig. 2 (e) & (f).

## 5 Numerical methods

In order to extend the analytical solutions into unknown parameter regimes we numerically solve the boundary value problem with open boundary conditions and as described in Sec. 23. As boundary conditions, we require zero first derivatives at for both fields and choose large enough for the solutions to settle in to the constant background.

After a solution is obtained, we check that the densities at are within 0.01 of the background density, and that the phases of the two fields at are within 0.01 of each other (). If either condition is not fulfilled, is increased and the solver is called again.

For the numerical procedure, an appropriate guess for the wave-function has to be provided. The initial guess is obtained from one of the analytically known solutions and is then followed in one of the parameters. All parts of the dispersion relation could be conveniently accessed by changing either of the controlling parameters , , in small steps.

We found that out of the entire -spectrum of analytic Manakov solutions at , only the and solutions extend to positive, finite . When , the stationary Manakov solution is identical to the zero-velocity Josephson vortex solution if . Indeed, following the solution from to yields the Josephson vortex branch obtained by following Josephson vortices from to . On the other hand, following the solution from to gives an entirely new branch, which we shall refer to as staggered solitons, due to the fact that the centers of the density dips are shifted with respect to each other (see Fig. 4 (c)-(f)).

## 6 Visualizing the solutions

In order to visualize the solutions, we show surface plots where the width of the two cylinders is related to the density of the two fields and the phase is encoded as a color map. In addition, we provide one-dimensional plots of the density and phase profiles to better resolve the finer details. We choose representative examples that illustrate the different solutions in all distinct regions of parameter space.

Figure 3 (a) & (b) show a moving Josephson vortex for . In all cases for the solutions were obtained by starting from the known zero-velocity Josephson vortices (21) and increasing velocity at a fixed . Note that physically, at , the Josephson vortex is centered exactly half way between the two parallel BEC lines. Its distinctive features are an equal dip in the density and an equal-but-opposite phase step in each condensate.

Figure 3 (c) & (d) show a stationary Josephson vortex at the maximum of the dispersion relation for and panels (e) & (f) show a moving Josephson vortex for the same and . The solutions for were obtained by starting from the previously-calculated wavefunctions at , and at each velocity gradually decreasing .

Note that, as shown in section 7, Fig. 5 (a), as goes to zero, the Josephson vortex dispersion relation is “split in half” as drops to zero. At each “wing” of the dispersion relation corresponds to a dark soliton in one of the two BEC lines. This can be seen clearly in Fig. 3 (c) & (d) where the density of one condensate is practically flat and the other has a strong dip. We therefore refer to the quasi-particles around the maxima of the Josephson vortex dispersion relation as “Josephson vortex maxima”, and interpret them as single-strand dark solitons. At the maxima of the dispersion relation, the vortex is exactly crossing one of the BEC strands as it moves out (perpendicularly to the BECs) from in between the two strands. Conversely, the dark soliton dispersion relation consists of a dark soliton in each of the BEC strands and the Josephson vortex dispersion relation merges with it as .

Figure 4 (a) & (b) show an example of a moving Josephson vortex for . The solutions for were obtained by starting from the previously-calculated wavefunctions at , and at each velocity gradually decreasing . Once part of the dispersion relation was available at each value, if necessary, we could complete it by following in .

Figure 4 (c) & (d) show a stationary staggered soliton for and panels (e) & (f) show a moving staggered soliton for the same and . These solutions were obtained by starting from the analytical Manakov wavefunctions at , and at each velocity gradually increasing . This gave us the central part of the dispersion relation at all values, which we then extended in at each constant .

## 7 Dispersion relation and other observables

Figure 5 panels (a), (c), (e) show the dispersion relations of dark solitons, Josephson vortices and of the staggered solitons, the latter only for . From panel (a) it is clear that for , the Josephson vortex dispersion relation changes concavity at at around . The same process is observed in reverse as with , as we move from panel (c) to (e). At , the equations reduce to the Manakov case, which is solved analytically in section 4.3 and indeed the Manakov solitons have a dispersion relation with a single central maximum.

As is clear from Fig. 5 (a), when the Josephson vortex dispersion relation bifurcates from the dark soliton one at some critical velocity that depends on , examined later in Fig. 9 (a). In panel (b) for the staggered soliton branch is unstable and lies between the stable dark soliton & Josephson vortex branches.

In Fig. 5 panels (b), (d), (f) we compare the energy of dark solitons, Josephson vortices, Josephson vortex maxima and staggered solitons at the extrema of the dispersion relations (which necessarily implies zero velocity) as a function of . In (b), for , the Josephson vortex and Josephson vortex maximum lines merge at around . It may be expected that this bifurcation point depends on , and this is indeed found to be the case. Panel (d) shows that at , the bifurcation point has now moved from to around . Notice that the staggered solitons bifurcate from the dark solitons, which accounts for their similar properties. Finally, at in panel (f), only the dark soliton-Josephson vortex bifurcation remains: Josephson vortices join the dark soliton line at , which is independent of .

Note that in (d), both the Josephson vortex and the Josephson vortex maximum solutions are stable, but the Josephson vortex maxima have , unlike all other solutions shown. The staggered soliton solutions only exist below about where they are unstable, while the higher energy dark solitons are stable. For , staggered solitons disappear and dark solitons become unstable.

The energy, missing particle number, angular momentum and phase difference are plotted as a function of velocity in Figs. 6-8 for the three parameter sets that were used in Figs. 3 & 4. The color code is identical to that used in Fig. 5. The change in concavity of the Josephson vortex branch as goes down through in Fig. 5 (a) is seen as the development of a loop in velocity-energy plots (compare panels (a) of Figs. 6 and 7). In Fig. 8, we see that staggered solitons and Josephson vortices merge and terminate at common end points. Only Josephson vortices have non-zero (and are therefore identified as vortices), while dark and staggered solitons have , and are hence classified as solitons.

## 8 Parameter regimes, types of excitations and their stability

Dark soliton solutions are analytically known for all parameter values. We have numerically obtained all translating Josephson vortex solutions in two parameter regimes: , and , . Staggered solitons were found in the second regime; this branch always has zero angular momentum and energy higher than Josephson vortices but lower than dark solitons. It is understood to be a transitory state through which Josephson vortices are able to reverse their circulation. Wherever the staggered soliton branch does not exist, dark solitons perform the role of the transitory state.

When , Kaurov and Kuklov [29] found that zero-velocity Josephson vortex solutions only exist for , at which point Josephson vortices merge into dark solitons. In fact, a bifurcation exists for all velocities, but the critical value of depends on velocity. A more natural point of view for us will be to say that for any given value of the tunnelling strength , the Josephson vortex and dark soliton dispersion relations merge smoothly at some critical momentum (associated with some critical velocity), and for larger momenta, the Josephson vortex branch does not exist. This is illustrated in Fig. 9 (a) where we plot the maximal velocity reached by the Josephson vortex branch (the critical velocity) as a function of . We found that, with , whenever Josephson vortices and dark solitons coexist, Josephson vortices are stable and dark solitons are unstable and when Josephson vortices cease to exist, dark solitons become stable. This fact was exploited in Ref. [45] where the authors present a similar plot to Fig. 9 (a) based on a stability calculation for dark solitons. An outline of the stability calculation is presented in Appendix A.

When we see that once again there exists a critical momentum beyond which the Josephson vortex solutions do not exist, but the Josephson vortex dispersion relation now terminates by touching the dark soliton dispersion relation non-tangentially (i.e. the slopes of the curves are different). The critical velocity is plotted as a function of in Fig. 9 (b). The staggered soliton branch terminates at the exact same critical momentum and velocity as the Josephson vortex branch.

In the regime, Josephson vortices are again always stable, but the situation for dark solitons is quite different. Figure 10 shows a numerically-determined boundary line (plotted in blue circles) in the - plane such that above this curve, dark solitons are unstable and below it they are stable. As soon as dark solitons become stable, staggered solitons appear. These are always unstable except exactly at (the entire Manakov family of solutions is always stable). There exist small regions of stability in Fig. 10, bounded by the almost vertical sections of the stability-flip curve and , the limits of . These are regions where staggered solitons and Josephson vortices do not exist and dark solitons are stable (as in the regime ). The development of these slivers of stability is seen in Fig. 9 (b) as a dip of the critical velocity, starting at about .

For zero velocity dark solitons, we can analytically compute the points in parameter space where stability changes – this is done in Appendix A.1. For as in Fig. 10, the result is , in agreement with numerical calculations (this point has been added to Fig. 10 as a red square). In fact, the analytical calculation also allows one to see that this stability-flip point starts at when , smoothly decreases and reaches at , so that outside of , neither Josephson vortices nor staggered solitons exist.

Thus, there is a region of bistability for where Josephson vortices (lowest energy) and dark solitons (highest energy) are both stable with the unstable staggered soliton branch (intermediate energy) between them. An illustration is given in Fig. 11 where we fix and plot the energy as a function of . The energies of dark solitons and Josephson vortices are constant since the solutions (16) and (21) are independent of , as is the energy functional (3) when . Overall, this has the familiar shape of a bistability bifurcation diagram with a fold. The unusual features are that the upper branch continues to the right past the fold and that the three lines do not make a single, smooth curve.

## 9 Inertial mass and missing particle number

In this section we focus on the first parameter range () and examine some overall properties of the dispersion relations. To start with, we can calculate the inertial mass of Josephson vortices and Josephson vortex maxima (evaluating the derivatives in (13) at the minimum and maximum of the dispersion relation, respectively) as a function of , which yields Fig. 12. The blue and red solid curves were obtained from the numerical Josephson vortex solutions. We define the bifurcation point at which the central part of the Josephson vortex dispersion relation changes concavity by the value at which the curve (red solid line in Fig. 12) crosses zero. This happens at . The magenta dash-dotted line shows the variational approximation for Josephson vortices (see section 10). The black dashed line will be described in section 12.

The inertial mass is a useful characteristic of an excitation, but the experimentally-accessible quantity is , the ratio of the inertial mass to the number of particles in the excitation, as it relates to the experimentally-measurable frequency ratio of small amplitude oscillations in the presence of weak harmonic trapping [50, 9, 10, 51]. With this in mind, Fig. 13 shows at the extrema of the Josephson vortex dispersion relation as a function of , and Fig. 14 shows the ratio obtained by combining the data from Figs. 12 and 13. The magenta dash-dotted line shows the variational approximation for Josephson vortices, described in section 10. It is clear that the red curve certainly crosses zero, which means that on either side of the critical point. This implies that essentially, the Josephson vortices become infinitely heavy.

## 10 Variational calculation for Josephson vortices

In light of the results of the previous section, we endevour to find a variational approximation for Josephson vortices near , i.e. in the immediate vicinity of the known analytical solution (21). We take the variational ansatz

 ψ1,2=√1+ν{isin(α)+cos(α)tanh(Az)±iB1,2sech(Az)eizε}, (28)

a form general enough to capture dark solitons, zero-velocity Josephson vortices and Manakov solitons. One then has to evaluate for this variational guess and take away for the background state, resulting in the difference, . Differentiating with respect to all five variational parameters () and setting the resulting expressions to zero, we obtain a system of five coupled non-linear equations. These are quite complicated, and a direct solution is impractical. Instead, we linearize the equations in : we set , where the zeroth order parameters are chosen to correspond with the solution (21): . The zeroth-order terms in the linearized equations thus cancel, and it remains to set the first order terms (in ) to zero. Introducing , we replace the equations resulting from by the sum and difference of these two equations. The five equations we must now solve decouple into two sets: two- and three-coupled equations. The solutions are: , and

 Ω = −48{−γ2+2(γ−2)γν+ν2[24+γ(44+3γ)]}[3ν+γ(6ν−2)] (29) − 4π2[2ν+γ(3ν−1)][3γν(7−29ν)−54ν2+5γ2(1+ν)(3ν−1)] + 3γπ4(1+ν)(γ−2ν−3γν)2, ~α = √ν{216ν2(π2−8)+6γν[168−888ν+4π2(19ν−5)+π4(1+ν)] (30) + γ2(3ν−1)[−96(1+13ν)+4π2(13ν−5)+3(1+ν)π4]}/Ω, ~ε = 72(1+2γ)ν2[6ν(π2−8)+γ(3ν−1)(3π2−32)]/Ω, (31) ~B− = 144γ(1+2γ)ν3/2π(3ν−1)/Ω. (32)

Linearising the variational equations in is an approximation that is of the same order as keeping terms up to in (the excitation energy) and in (the total momentum). Making such an expansion we can calculate the inertial mass , to obtain

 mI = 8√ν{48[γ+2ν(3+γ)+ν2(30+49γ)][3ν+γ(6ν−2)]−3γ(1+ν)2π4(2ν+γπ4(3ν−1)) (33) − 4π2[27ν2(1+5ν)+3γν(ν(137ν−14)−7)+γ2(5+ν(ν(309ν−133)−5))]}/Ω,

which is plotted in Fig. 12 alongside the numerical results. Using the zero-velocity solution (21), we can compute the missing particle number at as . The ratio from this calculation is shown in Fig. 14 as the magenta dash-dotted line. Note that Ref. [45] predicted , which is also displayed in Fig. 14 for comparison.

## 11 The sine-Gordon equation

The second parameter regime that we have investigated () is particularly interesting in terms of how it compares to the analytically solvable sine-Gordon model. In order to carry out such a comparison, we first give a brief review of the sine-Gordon equation.

In Appendix B we derive the sine-Gordon equation from the model of section 2 by assuming that the densities of the two fields are practically equal to each other and are almost constant. In addition, we isolate the terms from the Lagrangian density that contribute to the relative phase sector, which asymptotically decouples from the total phase sector in the limit of vanishing tunneling. While the total phase sector supports gapless elementary excitations, it is the relative phase sector that is captured by the sine-Gordon model and is relevant for the Josephson vortices. This selection of terms is partly justified a posteriori by the success of the analysis we perform in section 12.

The derivation of Appendix B allows one to express the sine-Gordon parameters through the Gross-Pitaevskii model parameters, thus enabling a direct comparison of the two models. In this section we will present some analytical results for the sine-Gordon equation [58], written with parameters determined by the procedure in Appendix B.

The Lagrangian density of the sine-Gordon model is

 L=ℏ24(g−gc)(∂tϕa)2−ℏ24mμ+Jg+gc(∂xϕa)2+2Jμ+Jg+gccos(ϕa), (34)

where

 ϕa=ϕ1−ϕ2. (35)

The Hamiltonian density can be obtained in the usual way:

 Pϕ = ∂L∂(∂tϕa), H = Pϕ(∂tϕa)−L, (36)

where is the canonical conjugate coordinate to . The Euler-Lagrange equation

 ∂L∂ϕa−∂x∂L∂(∂xϕa)−∂t∂L∂(∂tϕa)=0 (37)

yields the sine-Gordon equation:

 ∂ttϕa−γm(μ+J)∂xxϕa=−4Jγ(μ+J)ℏ2sin(ϕa). (38)

Rewriting in dimensionless form (see (5)) and in a frame moving at , the sine-Gordon equation becomes

 [v2s−γ(1+ν)]∂zzϕa+4νγ(1+ν)sin(ϕa)=0. (39)

The solution is given by

 ζ = √4νγ(1+ν)γ(1+ν)−v2s, ϕa = 4tan−1(eζz). (40)

The Hamiltonian density is

 H=14[v2sγ+1+ν](∂zϕa)2−2ν(1+ν)cos(ϕa), (41)

and the excitation energy is

 Es=8ν(1+ν)ζ+2ζ(1+ν+v2sγ). (42)

Next, using

 Pc(vs)=vs∫0d¯vs1¯vsdEsd¯vs, (43)

we get the canonical momentum as

 Pc=4vsγζ. (44)

We can eliminate to get the dispersion relation:

 E2s=(1+ν)[γP2c+64ν(1+ν)], (45)

or if we choose to write (in analogy to a relativistic particle)

 E2s=m2SGc4SG+c2SGP2c, (46)

then we identify

 mSG = 8√νγ, cSG = √γ(1+ν), (47)

as the“mass” and “speed of light” of the sine-Gordon soliton, respectively.

## 12 Relativistic behavior

We have seen that at and small , the coupled-BECs Josephson vortex dispersion relation develops a dip about (see Fig. 5), similar in shape to the central part of the dispersion relation of the sine-Gordon equation. The equivalence of the two models in this regime has been suggested before [30], and now that we have the sine-Gordon dispersion relation expressed through the Gross-Pitaevskii model parameters, we are in a position to check this statement.

First, we can compare the dispersion relations visually. This is shown in Fig. 15, and the Josephson vortex dispersion relation indeed seems to be very close to the sine-Gordon curve near the zero-velocity point . Next, we would like to compare the sine-Gordon parameters and to their equivalents in the coupled BECs model as a function of . A sensible way of extracting these parameters from the Josephson vortex dispersion relation is to first obtain from

 cJV=√max(dE2sdP2c), (48)

using data about , and then obtain as

 mJV= ⎷E2s(Pc=2π(1+ν))c4JV. (49)

The “relativistic mass” calculated this way (for ) is indistinguishable from obtained as a derivative using equation (13) plotted in Fig. 12 as a red solid line. Comparing the red line to the black dashed line () in Fig. 12, it appears that the Josephson vortex mass indeed approaches the sine-Gordon result as . Note that we are unable to compute numerical Josephson vortex solutions at smaller because the excitation length-scale becomes unmanageable.

As for the “speed of light”, , Fig. 16 shows that the functional dependence on is completely different for the coupled BECs and sine-Gordon models, and it is clear that the two only become equal at but the slopes remain different. We therefore conclude that the Gross-Pitaevskii model approaches the sine-Gordon model only asymptotically.

There are two fundamental speeds in the coupled-BECs model, which can be found by computing linearized excitations about the vacuum state, as was done in [6]. The authors find two elementary excitation branches: gapless Bogoliubov phonons (subscript “B”) and a gapped relative-phase excitations (subscript “RP”). A standard Bogoliubov calculation (such as the one in Appendix A) leads to the dimensionless oscillation frequencies

 ωB = √1+ν ⎷12k2(12(1+ν)k2+2), (50) ωR