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


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.

Figure 1: Possible realisations of a linearly-coupled two-component Bose gas as described by Eqs. (1) and (6). (a) A scalar Bose gas confined in an elongated double-well potential. A tight binding approximation leads to Eq. (1), where describes tunneling through the potential barrier and the cross-interaction vanishes (, ). (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 is achieved by coherent driving with a radio-frequency field shown in blue. In this case, the cross-component interaction is typically of similar magnitude to (, ) 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.

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


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 :


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




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 [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.

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


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


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:


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


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


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


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


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]


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).

Figure 2: Stationary solitary-wave solutions () in a coupled two-component BEC: (a) & (b) Dark soliton, (c) & (d) Josephson vortex, and (e) & (f) Manakov soliton at . Left column: The density is encoded in the width of the tube and the phase in the color. The upper tube corresponds to component 1 and the lower one to component 2. Specifically we define and plot isosurfaces at the value . Right column: density (top panels) and phase (bottom panels) profiles are shown directly as a blue dashed line (component 1) and a red dashed line (component 2). Other parameters are for all plots.

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 . The phase difference


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


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


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


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


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:


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


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


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


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


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: Numerical solutions from the Josephson vortex family corresponding to various points of interest on the dispersion relations of Fig. 5. Left column: The density is encoded in the width of the tube and the phase in the color. The upper tube corresponds to component 1 and the lower one to component 2. Specifically we define and plot isosurfaces at the value . Right column: density (top panels) and phase (bottom panels) profiles are shown directly as a blue dashed line (component 1) and a red dashed line (component 2). (a) & (b) Moving Josephson vortex: , (c) & (d) stationary Josephson vortex maximum: , (e) & (f) moving Josephson vortex: .
Figure 4: Numerical solutions from the Josephson vortex family corresponding to various points of interest on the dispersion relations of Fig. 5. Left column: The density is encoded in the width of the tube and the phase in the color. The upper tube corresponds to component 1 and the lower one to component 2. Specifically we define and plot isosurfaces at the value . Right column: density (top panels) and phase (bottom panels) profiles are shown directly as a blue dashed line (component 1) and a red dashed line (component 2). (a) & (b) moving Josephson vortex: , (c) & (d) stationary staggered soliton: , (e) & (f) moving staggered soliton: .

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.

Figure 5: Families of solitary-wave solutions. Left column: dispersion relations with excitation energy vs. canonical momentum for [or as indicated in panel (a)]. Right column: vs. the coupling parameter for stationary solitary waves (, marked by dots in the left column). (a), (b): (no cross-interaction, ). The highest energy branch corresponds to dark soliton solutions (green, labeled ‘DS’); the Josephson vortex (red, labeled ‘JV’) corresponds to the lowest available energy state for given . Different values of are shown in panel (a). For three local extrema coexist, each corresponding to a stationary solitary wave solution. There are two degenerate maxima (blue, labeled ‘JV(M)’), and one minimum corresponding to the stationary Josephson vortex solution. (c), (d): . Staggered soliton solutions (unstable) appear at intermediate energies for (black, labeled ‘SS’). (e), (f): . Josephson vortices and staggered solitons merge into the highly degenerate branch of “Manakov” solutions (stable). Throughout: Full (dashed) lines indicate stable (dynamically unstable) solutions; square markers indicate bifurcation points.

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.

Figure 6: No cross-interaction, intermediate coupling regime (). Energy (a), missing particle number (b), angular momentum (c) and phase difference (d) as a function of velocity. Properties of Josephson vortices (labeled ‘JV’) are plotted in red and those of dark solitons (labeled ‘DS’) in green. Solid lines indicate stable solutions, while dashed lines depict unstable ones.
Figure 7: No cross-interaction, weak coupling regime (, ). Energy (a), missing particle number (b), angular momentum (c) and phase difference (d) as a function of velocity. Properties of Josephson vortices (labeled ‘JV’) are plotted in red and those of dark solitons (labeled ‘DS’) in green. Solid lines indicate stable solutions, while dashed lines depict unstable ones.
Figure 8: Intermediate cross-interaction, weak coupling regime (). Energy (a), missing particle number (b), angular momentum (c) and phase difference (d) as a function of velocity. Properties of Josephson vortices (labeled ‘JV’) are plotted in red, those of dark solitons (labeled ‘DS’) in green, and quantities associated with staggered solitons (labeled ‘SS’), in black. Solid lines indicate stable solutions, while dashed lines depict unstable ones.

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.

Figure 9: Critical velocity (maximum allowed speed) for Josephson vortices as a function of with (a) and as a function of with (b). The red square is the analytical Manakov result and the full (blue) lines are numerical results.

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.

Figure 10: Stability diagram for dark (grey) solitons at . Dark solitons are unstable in the shaded region. The boundary points are calculated analytically (red square, see Appendix A.1) and numerically (blue dots) and the black line is a guide to the eye (spline fit to the data points). In the stable region below the black line unstable staggered solitons coexist with dark solitons. In the stability region at small and large all dispersion branches have merged and only stable grey solitons exist.
Figure 11: Excitation energy of stationary (unstable) staggered solitons as a function of at (black dashed line). The energy of the stable dark soliton (top green line) and Josephson vortex (lower red line) are also shown.

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.

Figure 12: Inverse of the inertial mass for Josephson vortices (upper red full line, labeled ‘JV’) and Josephson vortex maxima (lower blue full line, labeled ‘JV(M)’) at as a function of tunneling strength with obtained from numerical solutions. The magenta dash-dotted line is an approximate result obtained from a variational calculation for Josephson vortices (labeled ‘var JV’), equation (33). The black dashed line shows from (47), discussed is 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.

Figure 13: Missing particle number for Josephson vortices (red upper curve, labeled ‘JV’) and Josephson vortex maxima (blue lower curve, labeled ‘JV(M)’) as a function of tunneling strength with , evaluated at the extrema of the dispersion relation (i.e. ).
Figure 14: Missing particle number over inertial mass for Josephson vortices (red lower solid line, labeled ‘JV’) and Josephson vortex maxima (blue upper solid line, labeled ‘JV(M)’) as a function of tunneling strength with , evaluated at the extrema of the dispersion relation (i.e. ). The magenta dash-dotted line is an approximate result obtained from a variational calculation for Josephson vortices (labeled ‘var JV’), and the black dashed line is a prediction from Ref. [45].

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


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


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


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




The Hamiltonian density can be obtained in the usual way:


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


yields the sine-Gordon equation:


Rewriting in dimensionless form (see (5)) and in a frame moving at , 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 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.

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


using data about , and then obtain as


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