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

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.


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 con-finement 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 constant 1 . 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 Josephsonvortex dispersion, a new unstable branch of solitary waves arises at finite cross-component nonlinear interactions where both dark solitons and Josephson vortices are stable.  1) and (6). (a) A scalar Bose gas confined in an elongated double-well potential. A tight binding approximation leads to Eq. (1), where J describes tunneling through the potential barrier and the cross-interaction vanishes (g c = 0, γ = 1). (b) An atomic Bose gas with two accessible (hyperfine) spin components in a cigar-shaped quasi-one-dimensional trap. The two components are slightly off-shited for clarity. Linear coupling of the spin-components with constant J is achieved by coherent driving with a radio-frequency field shown in blue. In this case, the cross-component interaction g c is typically of similar magnitude to g (g c ≈ g, |γ| 1) but can be tuned to different values by means of a Feshbach resonance for certain atomic species [41,46]. 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.

The model
We are considering a quasi-one-dimensional two-component Bose gas with identical boson mass m for the two components and linear coupling J in the Gross-Pitaevskii approximation with the energy functional where Ψ j (x, t) is the order parameter of component j ∈ {1, 2}, µ is the common chemical potential and for simplicity, the intra-component interaction constant g [47] is also taken as common for the two components, but the cross-component interaction constant g c may be different. We assume that J > 0, which is the physical case when the coupling is achieved by tunnelling through a potential barrier [48]. In the general case, the phase of J can be absorbed into the definition of Ψ 2 . 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 mg/ 2ñ 0 1, i.e. the Lieb-Linger parameter is small [49], whereñ 0 = (µ + J)/(g + g c ) 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 Ψ j = δW [Ψ * j , Ψ j ]/δΨ * j : We are interested in solitary wave solutions that translate with a constant velocity V s . With the aim of adimensionalising the equations of motion we make the ansatz where are dimensionless variables. Further introducing the dimensionless linear and nonlinear coupling constants we obtain the dimensionless Gross-Pitaevskii equation for uniformly translating solutions We are looking for solutions to these equations that asymptotically tend to the stable constant background solution for z → ±∞ [48]. The boundary conditions can thus be written as 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 bears physical significance.

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 A key quantity is the excitation energy associated with the solitary wave where E 0 = −2L(1 + ν) 2 is the (dimensionless) energy of the constant background, and E sol is the energy of the solitary wave solution. For a localised solitary wave solution that heals to the constant background, E s is independent of the box size 2L for sufficiently large L. It will depend on the soliton velocity v s but there may be multiple solutions for each v s . Another useful observable is the momentum, which is scaled by µ/(g + g c ). The background solution has zero momentum and we introduce the following dimensionless observables: where P s is the physical momentum of the solitary wave with boundary conditions (7). The momentum difference ∆P 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 P cf 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 P c is the momentum that the solitary wave excitation has with periodic boundary conditions but it is also significant for the open boundary conditions (7) due to the relation The canonical momentum provides a convenient way of parameterising the solitary-wave solutions and the relation of E s vs. P c 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 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 ω z by where m P 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 m P = mN d | vs=0 2 . The (missing) particle number of the solitary wave N d is obtained by integrating the density and subtracting the background where n k (z) = |ψ k (z)| 2 are the dimensionless particle densities in the two components and n 0 = 1 + ν =ñ 0 (g + g c )/µ 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 µ/ √ µm(g + g c ).

Analytically known solitary-wave solutions
Several exact solutions of (6) are known. The lowest-energy constant solution is ψ 1 = ψ 2 = √ 1 + ν, which we shall refer to as the background. We remark that γ = 0 separates the miscible (γ > 0) and immiscible (γ < 0) phases of the system. In the miscible regime, ν = 0 corresponds to a degenerate mean-field ground state with undefined spin polarisation, with any ν = 0 leading to a background solution polarised along the x-direction. For this work we consider the miscible regime with γ ≥ 0 and a polarised ground state along the +x-direction obtained with ν > 0, while analogous results hold for ν < 0, where polarisation along the −x-direction is obtained.

Dark solitons
The coupled BECs system supports dark soliton solutions which satisfy ψ = ψ 1 = ψ 2 and are given by [1] ψ 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, v B = √ 1 + ν. Note that the dark soliton solutions are independent of the cross-interaction parameter γ. The zero-velocity case is visualised in Fig. 2 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 takes a maximum value at zero velocity and vanishes at v s = ±v B . The phase difference is π for stationary solitons and reaches the extremal values 0 and 2π at the limiting velocities ±v B . The missing particle number of the dark soliton evaluates to  The canonical momentum of the dark soliton is varying in the interval P c ∈ [0, 4π(1 + ν)] = [0, 2πn 0 ]. The inertial mass evaluates to m I = −8 √ 1 + ν.

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 and only exists for 0 < ν < 1 3 . The parameter value ν = 1 3 marks a bifurcation point where the Josephson vortex solution becomes identical to the dark soliton solution (16). For ν < 1 3 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 z = 0. Indeed, tracing the phase along the +z direction in ψ 1 and along −z in ψ 2 amounts to a total phase change of 2π if the origin z = 0 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 The energy and momentum difference for the Josephson vortex at v s = 0 are where the −(+) stands for the Josephson vortex ψ (anti-vortex ψ * ).

Manakov solitons
When the cross-and intra-component nonlinearities are equally strong (i.e. γ = 0), 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 χ 1,2 = 1 √ 2 (ψ 2 ± ψ 1 ), we re-write equations (6) for the new variables: where the two different signs in front of ν are to be taken with the two different indices, k = 1, 2. A trial solution of the form [57] χ 1 = αi + β tanh(ηz), is found to satisfy (23) if the parameters are given by In fact, χ 2 may be multiplied by an arbitrary phase factor, e iθ , and the resulting solution still satisfies the differential equations. Transforming back to the ψ-fields gives Notice that for the parameters in (25) to be real (and the solution to be non-trivial) we need ν < 1/3 and v 2 s < 4ν. The phase angle θ remains a free parameter, indicating the large degeneracy of these solutions. For v s = 0 and θ = ∓π/2 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 γ = 0 analytically: the excitation energy, angular momentum, phase difference, canonical momentum, and missing particle number are The inertial mass evaluates to m I = −2 5ν+1 √ ν and is independent of velocity. Note that the limits of P c are the same as for dark solitons. The Manakov solitons are illustrated in Fig. 2 (e) & (f).

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 z ∈ [−L, L] as described in Sec. 2 3 . As boundary conditions, we require zero first derivatives at ±L for both fields and choose L large enough for the solutions to settle in to the constant background.
After a solution is obtained, we check that the densities n k at ±L are within 0.01 of the background density, and that the phases of the two fields at ±L are within 0.01 of each other (φ k (±L) = φ 3−k (±L)). If either condition is not fulfilled, L 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 v s , ν, γ in small steps.
We found that out of the entire θ-spectrum of analytic Manakov solutions at γ = 0, only the θ = 0, ±π and θ = ±π/2 solutions extend to positive, finite γ. When γ = 0, the stationary Manakov solution is identical to the zero-velocity Josephson vortex solution if θ = −π/2. Indeed, following the θ = −π/2 solution from γ = 0 to γ > 0 yields the Josephson vortex branch obtained by following Josephson vortices from γ = 1 to γ < 1. On the other hand, following the θ = 0 solution from γ = 0 to γ > 0 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)).

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. In all cases for γ = 1, ν ≥ 0.15 the solutions were obtained by starting from the known zero-velocity Josephson vortices (21) and increasing velocity at a fixed ν. Note that physically, at P c = 2π(1 + ν), v s = 0, 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. Note that, as shown in section 7, Fig. 5 (a), as ν goes to zero, the Josephson vortex dispersion relation is "split in half" as E s (P c = 2π(1+ν)) drops to zero. At ν = 0 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 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      As is clear from Fig. 5 (a), when γ = 1 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 0 < γ < 1 the staggered soliton branch is unstable and lies between the stable dark soliton & Josephson vortex branches.
In 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 ν ≈ 0.14 − 0.15 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 ∆P (and are therefore identified as vortices), while dark and staggered solitons have ∆P = 0, and are hence classified as solitons. ν ≤ 0.33 and ν = 0.005, 0 ≤ γ ≤ 1. 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 γ = 1, Kaurov and Kuklov [29] found that zero-velocity Josephson vortex solutions only exist for ν < 1/3, 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 γ = 1, 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 γ < 1 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 γ < 1 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 P c -γ 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 γ = 0 (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 0, 4π(1 + ν), the limits of P c . These are regions where staggered solitons and Josephson vortices do not exist and dark solitons are stable (as in the regime γ = 1). The development of these slivers of stability is seen in Fig. 9 (b) as a dip of the critical velocity, starting at about γ = 0.86.
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 ν = 0.005 as in Fig. 10, the result is γ = 0.975, 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 γ = 1 when ν = 0, smoothly decreases and reaches γ = 0 at ν = 1/3, so that outside of 0 < ν < 1/3, neither Josephson vortices nor staggered solitons exist.
Thus, there is a region of bistability for γ < 1 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 ν = 0.005, v s = 0, P c = 2π(1 + ν) 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 (9) when |ψ 1 | 2 = |ψ 2 | 2 . 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.

Inertial mass and missing particle number
In this section we focus on the first parameter range (γ = 1) 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 1/m I curve (red solid line in Fig. 12) crosses zero. This happens at ν = 0.1413. 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 experimentallyaccessible quantity is m I /N d , 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 N d at the extrema of the Josephson vortex dispersion relation as a function of ν, and Fig. 14 shows the ratio N d /m I 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 m I /N d → ±∞ on either side of the critical point. This implies that essentially, the Josephson vortices become infinitely heavy.

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 v s = 0, P c = 2π(1 + ν), i.e. in the immediate vicinity of the known analytical solution (21). We take the variational ansatz a form general enough to capture dark solitons, zero-velocity Josephson vortices and Manakov solitons. One then has to evaluate L = E s − v s P c for this variational guess and take away L for the background state, resulting in the difference, ∆L . Differentiating ∆L with respect to all five variational parameters (A, α, B 1 , B 2 , ε) 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 v s : we set where the zeroth order parameters are chosen to correspond with the solution (21): 1+ν , ε 0 = α 0 = 0. The zeroth-order terms in the linearized equations thus cancel, and it remains to set the first order terms (in v s ) to zero. IntroducingB ± =B 1 ±B 2 , we replace the equations resulting from d∆L /dB 1,2 = 0 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:Ã =B + = 0, and α = √ ν 216ν 2 (π 2 − 8) + 6γν 168 − 888ν + 4π 2 (19ν − 5) + π 4 (1 + ν) Linearising the variational equations in v s is an approximation that is of the same order as keeping terms up to O(v 2 s ) in E s (the excitation energy) and O(v s ) in P c (the total momentum). Making such an expansion we can calculate the inertial mass m I = 2 dEs which is plotted in Fig. 12 alongside the numerical results. Using the zero-velocity solution (21), we can compute the missing particle number at v s = 0 as N d = −8 √ ν. The ratio N d /m I from this calculation is shown in Fig. 14 as the magenta dash-dotted line. Note that Ref. [45] predicted N d /m I = (5ν − 1)/(1 + ν), which is also displayed in Fig. 14 for comparison.

The sine-Gordon equation
The second parameter regime that we have investigated (ν = 0.005) 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 where The Hamiltonian density can be obtained in the usual way: , where P φ is the canonical conjugate coordinate to φ a . The Euler-Lagrange equation yields the sine-Gordon equation: Rewriting in dimensionless form (see (5)) and in a frame moving at v s , the sine-Gordon equation becomes The solution is given by The Hamiltonian density is and the excitation energy is Next, using we get the canonical momentum as We can eliminate v s to get the dispersion relation: or if we choose to write (in analogy to a relativistic particle) then we identify as the"mass" and "speed of light" of the sine-Gordon soliton, respectively.

Relativistic behavior
We have seen that at γ = 1 and small ν, the coupled-BECs Josephson vortex dispersion relation develops a dip about P c = 2π(1 + ν) (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 P c = 2πn 0 . Next, we would like to compare the sine-Gordon parameters m SG and c SG 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 c JV from using data about P c = 2π(1 + ν), and then obtain m JV as The "relativistic mass" m JV calculated this way (for ν ≤ 0.14) is indistinguishable from m I 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 (m SG ) in Fig. 12, it appears that the Josephson vortex mass m JV indeed approaches the sine-Gordon result as ν → 0. 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", c JV , 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 ν = 0 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 where k is a dimensionless wavenumber. If for some sufficiently small k the frequency ω becomes imaginary, the vaccum state is unstable. Thus, the vacuum can become unstable if γ < 0. The speeds associated with each branch are the speed of sound, c B = √ 1 + ν, and c RP = γ(1 + ν) + 2ν, which can be interpreted as a "speed of light". Both the elementary speeds are shown in Figs. 16 and 17 for comparison with sine-Gordon and Josephson vortex results. Notice that c RP is never equal to (the variational) c JV for γ > 0. Figure 17 finally explores the regime of finite cross-nonlinearity, γ < 1. Here we compare the sine-Gordon "speed of light" c SG to its equivalent from the Gross-Pitaevskii model c JV (showing both a numerical calculation and a variational approximation), and to the elementary speeds c B and c RP . We can see that the difference between c SG and c JV remains constant as a function of γ (it only depends on ν) and that both the Josephson vortex and sine-Gordon "speeds of light" exhibit a square-root dependence on γ (recall that c SG = γ(1 + ν)) while c B is independent of γ. Thus, by decreasing γ at a small ν we can decouple two fundamental speeds in the Gross-Pitaevskii model dispersion relation.
We can also use the results of section 10 to obtain an approximate analytical expression for the "speed of light". The excitation energy of the stationary vortex (21) which together with m I of equation ( In the absence of cross-nonlinearity (γ = 1), there is a critical point at ν ≈ 0.1413 where the Josephson vortex dispersion relation at P c = 2π(1 + ν) changed concavity. This corresponds to the inertial mass changing sign, going through ±∞. Thus, very "heavy" Josephson vortices can be created by tuning the coupling strength in this range. The heavy solitonic vortices observed experimentally in [9,10] are closely related, but there it is not possible to change the sign of m I by tuning a parameter.
The coupled BECs Josephson vortex dispersion relation at small ν and γ = 1 can be compared to the dispersion relation of the integrable sine-Gordon equation. The sine-Gordon parameters may be expressed through the Gross-Pitaevskii model parameters by deriving the sine-Gordon model from the Gross-Pitaevskii model in the small ν limit. We found that the Josephson vortex dispersion relation about P c = 2π(1 + ν) became equivalent to the sine-Gordon one asymptotically for ν → 0 but that the characteristic velocity scales of the Josephson vortices and the gapped linear excitations differ for finite ν. This challenges the widely-used approximation, or at least suggests some caution in its application. However, by working in the small ν regime, Josephson vortices may open the possibility for experimental study of "relativistic particles" (to a good approximation) using collective excitations of ultracold atoms.
When γ < 1, there exists a γ-and ν-dependent region where dark solitons and Josephson vortices are both stable, separated (in energy) by the unstable staggered solitons. For γ = 1, dark solitons are always unstable and Josephson vortices are stable. Therefore, observing dark solitons in such a coupled BECs system could be difficult because they would quickly decay into two opposite-circulation Josephson vortices. If one worked in the bistable region, however, since dark solitons are dynamically stable they would not decay. This could potentially enable one to observe dynamics and interaction of Josephson vortices with dark solitons experimentally. Note, however, that in the bistable regime dark solitons are still thermodynamically unstable, and thus the presence of a thermal cloud in a cold-atom experimental setting would lead to the eventual decay of the dark soliton, due to its negative inertial mass [1,59]. Josephson vortices around the minimum of the dispersion relation, however, should be stabilized at finite temperature.

A Stability calculation
This appendix gives details of how the stability of a solution to equations (6) is determined. We start from the Gross-Pitaevskii equations in dimensionless form, allowing for additional time dependence, other than mere translation at v s : To find out whether a solution is stable we must add a variation to the wavefunction: The right-hand side of (53) is substituted into (52); zero-order terms in δψ k give the unperturbed equations (52), terms of second order in δψ k and higher are discarded, and the first order terms give two linear equations for δψ k : We then make the ansatz Substituting (55) into (54) and separating out terms proportional to e −iλτ from those proportional to e iλ * τ (in light of orthogonality), we obtain four equations: where When these equations are written in matrix form (in the basis a 1 , b 1 , a 2 , b 2 ), it becomes clear that solving for the λ's reduces to diagonalizing the following matrix: (58) M is a matrix of operators, each of which must also be represented by a matrix. Let us consider these constituent operators first. These operate on the spatial dimension, discretized in steps of h. If the interval [−L, L] is discretized in to N grid points, then ν appearing in M is in fact ν multiplied by the N × N identity matrix. The wavefunctions, in turn, are represented by N × N matrices with ψ on the main diagonal. Products of ψ's are achieved by multiplying the appropriate ψ matrices together.
To construct ∂ z and ∂ zz we use a five-point stencil. In particular, if f (x) is some function and x is discretized in steps of h, then the first and second derivatives are approximated as  ) for the second-derivatives. In order to avoid boundary effects, on the second and pre-last rows we use a three point stencil: On the first and last rows, we also use the three point stencil but with additional assumptions. For the first derivative, we are forced to take a one-sided derivative, and for the second derivative, assume that f (x + h) = f (x − h). This is because only one of x ± h is part of the discrete grid when x is the first or the last point.
To find out whether a solution is stable or not, we need to know if there are any complex eigenvalues. The accuracy of the calculation is limited by h, and in our case, h = L/100 where 2L is the size of the system. h is usually 0.01, but for the largest systems can get up to 0.05 or 0.06. Note that the coupled Gross-Pitaevskii equations in this discrete representation are satisfied to order h 2 : the norm of the residuals is of order 10 −4 . In light of this, the cut-off for deciding whether the complex part of an eigenvalue is spurious or real is set to 0.01. Then, for each complex eigenvalue, the mod-squared eigenvector is inspected. If it is peaked in [−L/2, L/2], it is assumed to be an actual unstable mode. If it peaks outside this range, the complex eigenvalue is assumed to be spurious.
In the γ < 1 parameter regime, some extra care has to be taken when computing stability. For dark solitons, spurious unstable modes sometimes satisfy our conditions for true instability defined in the paragraph above. To distinguish them from real unstable modes, we required the eigenvector mod-squared at ±L to have decayed to one hundredth of the maximum value or more. The spurious modes have undamped oscillations beyond the region where the dark soliton is localized and are therefore ruled out by this extra condition. The next issue occurs for both dark solitons and staggered solitons: when the eigenvalue of a true unstable mode goes to zero as a function of some parameter, at some point it inevitably crosses our threshold of 0.01 (set in the paragraph above). This was suspected to occur in the high velocity limits. Therefore we checked that the pure imaginary eigenvalue belonging to the only potentially unstable eigenvector decayed smoothly as a function of velocity to zero. This confirmed that the mode in question was indeed unstable, even though its imaginary eigenvalue was less than 0.01.

A.1 Analytical stability for dark solitons
We are able to analytically determine the boundary between the stable and unstable regions in parameter space for the known dark soliton solutions. This calculation is not completely general, as in order for it to work, we are forced to assume that the dark solitons are stationary, thus fixing one of the parameters; ν and γ remain arbitrary, though.
Numerically we observe that the unstable eigenvector for dark solitons always has zero real part, and therefore, when dark solitons change stability (i.e. when the imaginary part of the eigenvalue goes through zero), the entire eigenvalue is zero. We are thus interested in solving M = 0. Defining the change of basis matrix D −D * + γψ 2 − γψ * 2D +D * − γψ 2 − γψ * 2 D +D * + γψ 2 + γψ * 2D −D * − γψ 2 + γψ * 2 , and we will denote U = [ã,b] T . The choice v s = 0 guarantees that ψ 2 = ψ * 2 andD =D * , and hence the diagonal elements of (64) vanish. The resulting equations read These equations have the same form as the (time-independent) Schrödinger equation, i.e. the eigen-problem for the Hamiltonian. In addition to the usual kinetic term we have a sech 2 potential -known as the Rosen-Morse potential after the authors who first solved this problem analytically [60], and a constant term which can be interpreted as the eigenvalue. The energy spectrum consists of a few discrete bound states (localized and square-integrable), followed by a continuum of higher-energy, unbound states (delocalized). When the parameters are just right for the bound energy eigenvalues of the Hamiltonians to match the eigenvalue terms in the equations, the two equations (65) are satisfied with localized solutions. In other words, for such parameter values a zero eigenvalue of (61) exists and dark solitons switch stability. Reference [60] derives the following results: the equation has discrete, bound eigenvalues where n = 0 or n ∈ N, n ≤ κ + 1 4 − 1 2 . For direct comparison of (65) with this result, we must rewrite the potential terms through sech 2 and change to the scaled position coordinatez = √ 1 + νz. This procedure yields 4 [ν + γ(1 + ν)] 1 + νã = ∂zz + 2(1 + 2γ) sech 2 (z) ã, 4ν 1 + νb = ∂zz + 2 sech 2 (z) b .
Examining the equation forb and comparing to the Rosen-Morse results, n can only be 0 or 1. Moreover, we easily compute 0 = 1 and 1 = 0. Next we set each n equal to the eigenvalue 4ν 1+ν and see what conditions this imposes on our parameters. Doing this for 1 leads to ν = 0 and for 0 leads to ν = 1/3. These are well-known points at which dark solitons do change stability: at ν = 0 Josephson vortices appear and dark solitons change from stable to unstable while the reverse process occurs at ν = 1/3. Now let us compare the equation forã to the Rosen-Morse results: n can be 0, 1 or 2, the latter only if γ ≥ 5/8. Setting n equal to the eigenvalue of theã equation gives We can use this condition to check our numerical results. Setting ν = 0.005, and taking n = 0, 1, 2 in turn, we plot the left-and right-hand sides of (69) as a function of γ, looking for the intersection point. For n = 0 (69) is satisfied at γ ≈ 0.975, for n = 1 the lines do not cross and for n = 2 they cross at γ ≈ 0.1565 < 5/8, so n = 2 is not actually possible at this point in parameter space. Thus, this calculation predicts that stationary dark solitons at ν = 0.005 will change stability once, at γ ≈ 0.975. This point is added to Fig. 10 (red square) and fits perfectly on the numerical curve (blue circles).

B Derivation of the sine-Gordon equation
In this appendix we show how one can obtain the sine-Gordon equation from the Gross-Pitaevskii model of section 2. The Lagrangian density of the coupled BECs system is given by where the energy density (also see (9)) is and The Gross-Pitaevskii equations (2) can be recovered from the Euler-Lagrange equations for the fields Ψ k and Ψ * k . To proceed, we take the following ansatz for the wavefunctions: In terms of the new fields, (70) becomes (∂ x φ s ) 2 + (∂ x φ a ) 2 + u 2 2 cos(2Θ)∂ x φ s ∂ x φ a +µu 2 + Ju 2 sin(2Θ) cos(φ a ) − g 2 u 4 cos 4 (Θ) + sin 4 (Θ) −g c u 4 cos 2 (Θ) sin 2 (Θ).
We now assume that the densities of the two wavefuncitons are almost the same, i.e., we take Θ(x, t) = π 4 + y(x, t), where y is a field of small magnitude. We expand L to second order in y: Expanding out all the brackets in (76), we keep only the 2 nd , 6 th , 9 th , 12 th , and 14 th terms. This selection is based upon whether or not the term is needed in the reduced Lagrange density in order for it to yield the sine-Gordon equation. The reduced Lagrangian reads We write down the Euler-Lagrange equations for y and φ a , make the approximation that u is a constant, eliminate y between the two equations and get where u was set to the background value, Equation (38) is identical to (78).