Beating dark-dark solitons in Bose-Einstein condensates

Beating dark-dark solitons in Bose-Einstein condensates D. Yan,1J.J. Chang,2C. Hamner,2M. Hoefer,3P.G. Kevrekidis,1P. Engels,2V. Achilleos,4D.J. Frantzeskakis,4and J. Cuevas5 1Department of Mathematics and Statistics, University of Massachusetts, Amherst, Massachusetts 01003-4515, USA 2Washington State University, Department of Physics & Astronomy, Pullman, Washington 99164, USA 3North Carolina State University, Department of Mathematics and Statistics, Raleigh, North Carolina 27695, USA 4Department of Physics, University of Athens, Panepistimiopolis, Zografos, Athens 157 84, Greece 5Nonlinear Physics Group. Escuela Polit´ ecnica Superior. Departamento de F´ ısica Aplicada I. Universidad de Sevilla. C/ Virgen de ´ Africa, 7. 41011-Sevilla (Spain) Motivated by recent experimental results, we study beating dark-dark solitons as a prototypical coherent structure that emerges in two-component Bose-Einstein condensates. We showcase their connection to darkbright solitons via SO(2) rotation, and infer from it both their intrinsic beating frequency and their frequency of oscillation inside a parabolic trap. We identify them as exact periodic orbits in the Manakov limit of equal inter- and intra-species nonlinearity strengths with and without the trap and showcase the persistence of such states upon weak deviations from this limit. We also consider large deviations from the Manakov limit illustrating that this breathing state may be broken apart into dark-antidark soliton states. Finally, we consider the dynamics and interactions of two beating dark-dark solitons in the absence and in the presence of the trap, inferring their typically repulsive interaction. I. INTRODUCTION One of the principal themes of study in the emerging field of atomic Bose-Einstein condensates (BECs) is the examination of the coherent structures that arise in them [1–4]. When such explorations started over a decade ago [5–9], they were considerably hindered by either geometric or thermal effects, which were detrimental towards the lifetime of dark solitons and vortices that can be formed in repulsive BECs. Yet, the newer generations of experiments have enabled considerable strides towards the observation of dynamics and interactions of such nonlinear waveforms [10–15]. In addition to the above context of single-component BECs, soliton and vortex states may also arise in multi-component condensates, such as the two-component pseudo-spinor BECs, or the three and higher component spinor BECs [1, 2, 4]. One DB DB DD DB DB DD DD Density [a.u.] 0x[microns] 100 200 DB DB DD DB 2,2 1,1 2,2 1,1 2,2 - 1,1 - DB DB DD DB DB DD DD Density [a.u.] 0x[microns] 100 200 DB DB DD DB 2,2 1,1 2,2 1,1 2,2 - 1,1 - FIG. 1: (Color online) Prototypical experimental images of dark-bright and dark-dark solitons in a two-component BEC. The two components are vertically offset for separate imaging. All dynamics occur with vertically overlapped components before the imaging procedure. Clear examples of dark-bright and dark-dark solitons are marked as DB and DD respectively. In the fourth panel, the red (thick) line shows a radially integrated cross section of the upper component in the boxed region of the third panel, while the black (thin) line shows the cross section of the lower component. The |F, mF⟩hyperfine states used for these images are given to the right of each component. 2 of the prototypical examples of a soliton state in these settings is the so-called dark-bright (DB) soliton [16, 17]. Experimental images of DB solitons in a two-component BEC are presented in Fig. 1. The BEC in this figure is comprised of two different hyperfine states of 87Rb, and the solitons are generated by subjecting the BEC to inter-component counterflow; details of this technique are described in Refs. [18, 19]. In each panel, the atom clouds of the two components are vertically offset for imaging only, while all the dynamics leading to the soliton formation occurs in overlapped clouds. Clear examples of darkbright solitons are marked as “DB” in the figure and they consist of a dark soliton in one component that is coupled to a bright soliton in the second component. These structures can be thought of as “symbiotic” (or even parasitic) states because their bright component cannot be supported alone in the case of repulsive interactions [3]; in fact, the bright soliton is only sustained because of the presence of its dark-counterpart, which operates as an external trapping potential. Although dark-bright solitons (and even a prototypical interaction thereof) were first observed some time ago in the context of nonlinear optics [20, 21], their observation in recent atomic BECs experiments [10] triggered a sizeable burst of research activity centered around them. Topics of study included (but were not limited to) multi-DB soliton solutions from the viewpoint of integrable systems [22], numerical study of DB soliton interactions [23], discrete DB solitons [24], experimental realizations of DB soliton trains [18], DB soliton oscillations and interactions [25, 26], as well as interaction of DB solitons with localized impurities [27]. Recently, a “cousin” of these DB solitons, namely the dark-dark (DD) soliton — which involves two dark solitons but with potentially a breathing oscillation between their densities was also experimentally observed [19]. Pertinent examples are marked as “DD” in Fig. 1. These solitons show interesting dynamics in which they periodically change their form, from the one shown in the first panel to the one shown in the second panel, and back (note the order of the hump/notch features in each of the DD’s component; see also Fig. 5 below). Such “beating dark-dark solitons” are expected to emerge in the integrable two-component (so-called Manakov) limit of the relevant mean-field theoretic models [28] and were, in fact, earlier observed in numerical experiments involving the dragging of defects through the binary condensates [29]. The current experimental advances, such as the ones leading to the soliton images of Fig. 1, motivate the present theoretical study, in which we revisit DD soliton states at the integrable Manakov limit and extract information from their connection to the DB solitons (section II). These results are corroborated by the identification of such single DD soliton states, as genuine periodic orbits of the Manakov case (with and without a trap) and the study of their stability, internal modes and associated near-equilibrium dynamics (section III). In addition, we examine the dynamics of individual such solitons, upon departure from the integrable limit (section IV). Experimentally it has also become possible to generate several solitons, and even solitons of different types, in a single BEC – see, e.g., the third panel of Fig. 1 which demonstrates the coexistence of dark-bright and dark-dark solitons. Although the experimentally exploited counterflow between the two components is beyond the scope of our current analysis, these experimental findings motivate our investigation of the interactions between two dark-dark solitons (section V). Finally, conclusions of our study, as well as a number of interesting perspectives for future work, are also presented (section VI). II. DARK-BRIGHT AND DARK-DARK SOLITONS: THEORETICAL BACKGROUND We consider a two-component elongated (along the x-direction) repulsive BEC, composed of two different hyperfine states of the same alkali isotope. In the case of a highly anisotropic trap (i.e., if the longitudinal and transverse trapping frequencies are such that ωx≪ω⊥), this system can be described by two coupled Gross-Pitaevskii equations (GPEs) of the form [1]: i~∂tψj=(−~2 2m∂2 xψj+V(x)−µj+ 2 ∑ k=1 gjk|ψk|2)ψj.(1) Here, ψj(x, t)(j= 1,2) denote the mean-field wave functions of the two components (normalized to the numbers of atoms Nj=∫+∞ −∞ |ψj|2dx), mis the atomic mass, and µjare the chemical potentials; furthermore, gjk = 2~ω⊥ajk are the effective one-dimensional (1D) coupling constants, ajk denote the three s-wave scattering lengths (note that a12 =a21) which account for collisions between atoms belonging to the same (ajj) or different (ajk, j =k) species, while V(x) = (1/2)mω2 xx2is the external trapping potential. Let us now assume that the two-component BEC under consideration consists of two different hyperfine states of 87Rb, such as the states |1,−1⟩and |2,1⟩used in the experiment of Ref. [30], or the states |1,−1⟩and |2,−2⟩used in the experiments of Refs. [18, 19, 25]. In the first case the scattering lengths take the values a11 = 100.4a0,a12 = 97.66a0and a22 = 95.00a0, while in the second case the respective values are a11 = 100.4a0,a12 = 98.98a0and a22 = 98.98a0(where a0is the Bohr radius). In either case, it is clear that the scattering lengths take approximately the same values, say aij ≈a. This way, measuring the densities |ψj|2, length, time and energy in units of 2a,a⊥=√~/ω⊥,ω−1 ⊥and ~ω⊥, respectively, we may cast Eqs. (1) into 3 the following dimensionless form, i∂tu1=−1 2∂2 xu1+V(x)u1+ (|u1|2+|u2|2−µ)u1,(2) i∂tu2=−1 2∂2 xu2+V(x)u2+ (|u1|2+|u2|2−µ)u2,(3) where we have also assumed that the chemical potentials characterizing each component are equal. Note that the potential in Eqs. (2)-(3) is now given by V(x) = (1/2)Ω2x2, where Ω = ωx/ω⊥is a natural small parameter of the system. The above system of Eqs. (2)-(3) is invariant under SU(2) rotations [28]. In particular, let us first recall that a general matrix element of SU(2) takes the form U=(α−β∗ β α∗), where αand βare complex constants such that |α|2+|β|2= 1. Then, it can be shown that if (u1, u2)Tare solutions of Eqs. (2)-(3), then, (u′ 1 u′ 2)≡U(u1 u2)=(αu1−β∗u2 βu1+α∗u2), are also solutions of Eqs. (2)-(3). This suggests that we may start from the exact dark-bright (DB) soliton solution (which exists in the absence of the potential) and, derive the beating dark-dark (DD) soliton solution. More specifically, in the absence of the external potential (V(x) = 0), and for the boundary conditions |u1|2→µand |u2|2→0as |x| → ∞, Eqs. (2)-(3) possess an exact analytical single DB soliton solution of the following form: u1(x, t) = √µ{cos ϕtanh ξ+isin ϕ},(4) u2(x, t) = ηsech ξexp{ikx +iθ(t)},(5) where ξ=D(x−x0(t)),ϕis the dark soliton’s phase angle, cos ϕand ηrepresent the amplitude of the dark and bright solitons, and Dand x0(t)are associated with the inverse width and the center position of the DB soliton. Furthermore, k=Dtan ϕand θ(t)are the (constant) wavenumber and phase of the bright soliton, respectively. The above parameters of the DB soliton are connected through the following equations: D2=µcos2ϕ−η2,(6) ˙x0=k=Dtan ϕ, (7) ˙ θ=1 2(D2−k2),(8) with ˙x0and ˙ θdenoting the DB soliton velocity and angular frequency, respectively (overdots denote time derivatives). Thus, the DB soliton (4), (5) is characterized by three free parameters (seven parameters µ,ϕ,η,k,D,˙x0,˙ θand four constraints (6)-(8)). Notice that the amplitude ηof the bright soliton, the chemical potential µof the dark soliton, as well as the (inverse) width parameter Dof the DB soliton are connected to the number of atoms NBof the bright soliton by means of the following equation: NB≡∫R|u2|2dx=2√µη2 D.(9) According to the above arguments, one may start from the DB soliton and construct SU(2) rotated solutions, in the following form: u1(x, t) = α√µ{cos ϕtanh ξ+isin ϕ}−β∗ηsech ξexp{ikx +iθ(t)},(10) u2(x, t) = β√µ{cos ϕtanh ξ+isin ϕ}+α∗ηsech ξexp{ikx +iθ(t)}.(11) With the additional four parameters α,β∈Cand the constraint |α|2+|β|2= 1, the solution (10)-(11) is characterized by six free parameters. Introducing a new parameter c, the velocity of the background fluid, another solution can be constructed from Eqs. (10)-(11) via a Galilean boost: exp[i(cx −c2t/2)]u1,2(x−ct, t). Thus, in the most general case, this DD soliton solution is characterized by seven free parameters. One natural set of parameters can be found from the far-field, |x| → ∞ behavior consisting of two densities, an overall fluid velocity, and four phases. 4 Due to Galilean invariance and phase invariance, u′ j(x, t) = eiφjuj(x, t), we will assume, without loss of generality, that the background is at rest (c= 0) and focus, more specifically, on the case of the SO(2) rotated DB soliton. In this case, the corresponding orthogonal matrix is given by: U=(cos(χ)−sin(χ) sin(χ) cos(χ)),(12) where χis an arbitrary angle. This way, the relevant SO(2) rotated soliton solution takes the form: u1(x, t) = cos(χ)√µ{cos ϕtanh(D(x−x0(t))) + isin ϕ}−sin(χ)ηsech(D(x−x0(t))) exp{ikx +iθ(t)},(13) u2(x, t) = sin(χ)√µ{cos ϕtanh(D(x−x0(t))) + isin ϕ}+ cos(χ)ηsech(D(x−x0(t))) exp{ikx +iθ(t)},(14) The solution (13)-(14) is a DD soliton solution characterized by 4 free parameters. The asymptotics of this solution are |u1|2→ µcos2(χ)and |u2|2→µsin2(χ)as |x| → ∞. The densities of the above dark solitons read: n1≡ |u1|2=µcos2(χ)−(µcos2(χ) cos2ϕ−η2sin2(χ))sech2ξ −√µη sin(2χ){sin ϕsin[kx +θ(t)] + cos ϕcos[kx +θ(t)] tanh ξ}sechξ, (15) n2≡ |u2|2=µsin2(χ)−(µsin2(χ) cos2ϕ−η2cos2(χ))sech2ξ +√µη sin(2χ){sin ϕsin[kx +θ(t)] + cos ϕcos[kx +θ(t)] tanh ξ}sechξ, (16) while the total density ntot of the DD soliton is given by: ntot =n1+n2=µ−D2sech2ξ. (17) Notice that the total density of the DD soliton is time-independent and has the form of a dark soliton density of depth D2on top of a background density µ. The above density is, in fact, identical to the density of the DB soliton; this is due to the fact that under SO(2) rotation the total density, as well as all other conserved quantities of the system, remain unchanged. This will be particularly important when considering the motion of the DD soliton in a trap — see below. On the other hand, one may consider the individual dark soliton densities, n1and n2, across the trajectory of the DD soliton, i.e., for ξ= 0: in such a case, x=x0(t) = kt and the densities read: n1(ξ= 0) = µcos2(χ) sin2ϕ+η2sin2(χ) −√µη sin(2χ) sin ϕsin [1 2(k2+D2)t],(18) n2(ξ= 0) = µsin2(χ) sin2ϕ+η2cos2(χ) +√µη sin(2χ) sin ϕsin [1 2(k2+D2)t].(19) It is clear that n1,2(ξ= 0) are periodic functions of time; the relevant angular frequency (which constitutes the internal beating frequency of the DD soliton) is given by: ω0=1 2(k2+D2) = 1 2(µ−η2sec2ϕ),(20) where we have also used Eq. (6). The frequency ω0is bounded by two limiting values. First, in the case η→0, the beating DD soliton becomes a regular DD soliton, characterized by a width D=√µcos ϕand a velocity k=√µsin ϕ; in this case, ω0→(1/2)µ. Second, in the limiting case D→0, the beating DD soliton is reduced to a plane wave; in this case, ω0→(1/2)k2. In other words, the intrinsic oscillation frequency take values in the range: 1 2k2< ω0<1 2µ. (21) III. DARK-DARK SOLITONS AS PERIODIC ORBITS IN THE MANAKOV MODEL In this section, we analyze the existence, stability and dynamics of single beating DD solitons in a trap of the form V(x) = 1 2Ω2x2, considering them as periodic orbits. In the presence of the trap, the dynamics of the center of mass x0(t)of the beating DD soliton is still described by the dynamics of the original (unrotated) DB soliton center x0. This is due to the fact that the GPEs (2)-(3) are invariant under SO(2) rotations even in the presence of V(x), and so are all conserved quantities of the system, 5 −20 −15 −10 −5 0 5 10 15 20 −1 −0.5 0 0.5 1 x u1(x), u2(x) −20 −15 −10 −5 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 x |u1(x)|2, |u2(x)|2 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Im(λ) Re(λ) FIG. 2: Left panel: Profiles and densities of the beating dark-dark soliton solution with Ω=0,ω0= 0.5and µ= 1.5at t= 0. Right panel: Floquet multipliers spectrum for the dark-dark soliton displayed in the left panel. such as the total energy. Since the derivation of the equation of motion for the DB soliton center x0in Ref. [16] was relying on the change of energy (due to the presence of the trap), it is clear that the evolution of the beating DD soliton center follows the same dynamics: it performs a harmonic oscillation in the trap according to the equation ¨x0+ω2 oscx0= 0, where the oscillation frequency ωosc is given by [16]: ω2 osc = Ω2(1 2−r 8√1 + (r 4)2)(22) where r=NB √µis a measure of the ratio of the number of atoms in the bright and dark soliton component. In order to compute the soliton profile and determine its stability, we consider the solution of Eqs. (2)-(3), with g11 =g22 =g12 = 1, as a Fourier series expansion of period ω0[32], namely, u1(x, t) = ∞ ∑ k=−∞ zk(x)eikω0t, u2(x, t) = ∞ ∑ k=−∞ yk(x)eikω0t,(23) with {zk},{yk} ∈ R. Then, the dynamical equations are reduced to a set of coupled equations: [µ−kω0−V(x)]zk+1 2∂2 xzk=∑ p∑ q (zpz∗ q+ypy∗ q)zk−p+q(24) [µ−kω0−V(x)]yk+1 2∂2 xyk=∑ p∑ q (zpz∗ q+ypy∗ q)yk−p+q(25) where we have used the notation zk≡zk(x),yk≡yk(x). If the trap is absent, it is straightforward to see that z0(x) = √µ 2tanh(√2ω0x) = y0(x),(26) z1(x) = −√µ 2−ω0sech(√2ω0x) = −y1(x),(27) zj(x) = yj(x)=0,|j|>1or j =−1,(28) is actually the solution (13)-(14) for χ=π/4,ϕ=k= 0, and ω0=D2/2. In order to numerically find a DD soliton solution in the system with the trap, the previous solution (with the dark component {zk}multiplied by the Thomas-Fermi cloud with uT F 1=√max(µ−V(x),0)) is introduced as a seed for a fixed-point method in the system of Eqs. (24)-(25). Throughout this section, we have considered —for convenience— a trap strength Ω=0.2in order to consume less time in the numerical calculations, as will be explained below. Figures 2 and 3 show the periodic orbit for t= 0 without and with a trap potential, respectively. It is worth remarking that solutions in the trap exist for µ > 2ω0, as predicted in the end of section II. Once a periodic solution is found, its (linear) orbital stability can be analyzed by means of Floquet analysis. To this end, the time evolution of a small perturbation ξ1(x, t), ξ2(x, t)to a periodic solution {u1,0(x, t), u2,0(x, t)}must be traced. For the double indices of ui,j,irepresents the component index, i.e., i= 1,2is the first and second component of the dark-dark soliton 6 −20 −15 −10 −5 0 5 10 15 20 −1 −0.5 0 0.5 1 x u1(x), u2(x) −20 −15 −10 −5 0 5 10 15 20 0 0.2 0.4 0.6 0.8 1 x |u1(x)|2, |u2(x)|2 −1 −0.5 0 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Im(λ) Re(λ) FIG. 3: Same as Fig. 2 but in the trapped case with Ω = 0.2. solution respectively. The index j= 0 denotes that this is the (numerically) exact periodic solutions for u1and u2, around which we linearize in our Floquet analysis. The perturbations are introduced in the dynamical equations (2)-(3) as: u1(x, t) = [u1,0(x, t) + δξ1(x, t)] , u2(x, t) = [u2,0(x, t) + δξ2(x, t)] ,(29) and the resulting equation at order O(δ)reads: i∂tξ1= [−1 2∂2 x+ 2|u1,0|2+|u2,0|2−µ+V(x)]ξ1+u2 1,0ξ∗ 1+u∗ 2,0u1,0ξ2+u2,0u1,0ξ∗ 2,(30) i∂tξ2= [−1 2∂2 x+|u1,0|2+ 2|u2,0|2−µ+V(x)]ξ2+u2 2,0ξ∗ 2+u∗ 1,0u2,0ξ1+u1,0u2,0ξ∗ 1.(31) Then, in the framework of Floquet analysis, the stability properties of periodic orbits are resolved by diagonalizing the monodromy matrix M, which is defined as:    Re(ξ1(x, T)) Im(ξ1(x, T)) Re(ξ2(x, T)) Im(ξ2(x, T))   =M   Re(ξ1(x, 0)) Im(ξ1(x, 0)) Re(ξ2(x, 0)) Im(ξ2(x, 0))   .(32) with T= 2π/ω0. As the system is symplectic and Hamiltonian, the linear stability of the solutions requires that the monodromy eigenvalues, λ(also called Floquet multipliers) must lie on the unit circle (see, e.g., [35, 36] for details). The Floquet multipliers can also be written as λ= exp(iΘ), with Θbeing denoted as the Floquet exponent. An internal mode of the soliton corresponds to a spatially localized solution of Eqs. (30)- (31), with its oscillation frequency related to the Floquet exponents as ωm= Θω0/(2π). Figures 2 and 3 show a typical Floquet multiplier spectra, indicating stability of the periodic orbits. All the analyzed solutions (i.e. with Ω = 0 and Ω = 0.2) are stable. The choice of a trap strength Ω = 0.2for studying the stability of periodic orbits instead of the value Ω = 0.01 as in the rest of the paper is twofold. On the one hand, as indicated by Eq. (22), the oscillation period scales with Ω−1; consequently, decreasing 20 times the trap strength implies an integration time 20 times larger; on the other hand, the Thomas-Fermi radius (RT F =√2µ/Ω), which measures the dark soliton size, would also increase 20 times, so the number of equations to integrate also increases at this rate. With Ω=0.2, obtaining the monodromy matrix takes around 30 minutes. Consequently, the computation of the monodromy matrix for Ω = 0.01 would take about 200 hours. Some interesting results can be extracted by the analysis of the internal modes of the periodic orbits. Figure 4(left) shows the dependence of three internal modes of the Floquet spectrum with respect to µfor ω0= 0.5. The blue line is close to the frequency predicted by Eq. (22) [depicted as a dashed red line]. Indeed, perturbing the beating DD soliton with the corresponding eigenmode, we have confirmed that this perturbation leads to an oscillation of the soliton in the trap with a frequency equal to that of the eigenmode (cf. left panel of Fig. 5). It can be observed that the agreement between the numerical eigenfrequency and that predicted by Eq. (22) improves when µincreases, as expected. This is because the assumption that the solitary wave is a particle inside a Thomas-Fermi cloud is one of increasing validity the deeper one is within the Thomas-Fermi limit of large µ. The right panel of Figure 4 shows the dependence of the frequency of the internal mode corresponding to the oscillation of the trap with respect to ω0for fixed µ= 5 and compares it with the frequency predicted by Eq. (22). We note here, as an aside in the case Ω = 0, that the internal soliton modes are neutral modes located at (1,0) on the unit circle. In particular, the mode associated with the oscillation of the DD soliton in the trap becomes in this case a neutral mode associated with the translation of the soliton. The algebraic multiplicity of the multiplier at (1,0) in the case of Ω = 0 is 8, while in the trapped case (due to the lifting of translational invariance) it is 6. 7 1.5 2 2.5 3 3.5 4 4.5 5 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 µ ωm 0.5 1 1.5 2 2.5 0.04 0.06 0.08 0.1 0.12 0.14 0.16 ωo ωm FIG. 4: (Left panel) Dependence of the eigenfrequencies of some internal modes (solid lines) with respect to µwith fixed ω0= 0.5(The black, blue, green and red lines are the upper solid, middle solid, lower solid and the dashed respectively). The right panel shows that dependence with respect to ω0with fixed µ= 5 for the mode in blue in the left panel. In both panels, the red dashed line corresponds to the oscillation frequency (22) in a trap with Ω = 0.2. x t 0 50 100 150 −10 −8 −6 −4 −2 0 2 4 6 8 10 x t 0 50 100 150 −10 −8 −6 −4 −2 0 2 4 6 8 10 x t 0 50 100 150 −10 −8 −6 −4 −2 0 2 4 6 8 10 x t 0 50 100 150 −10 −8 −6 −4 −2 0 2 4 6 8 10 x t 0 50 100 150 −10 −8 −6 −4 −2 0 2 4 6 8 10 x t 0 50 100 150 −10 −8 −6 −4 −2 0 2 4 6 8 10 FIG. 5: Contour plots showing the evolution of the densities of the DD soliton components (first and second rows depict respective components) when perturbed by three different eigenmodes: the left panel corresponds to a soliton perturbed along the blue mode and leads the soliton to harmonically oscillate near the center of the trap; in the central panel, the perturbation (along the black mode) leads to a breathing of the soliton width, whereas in the right panel (perturbation along the green mode), the outcome corresponds to an oscillation of the whole condensate. In both cases, µ= 3,ω0= 0.5and Ω = 0.2. In order to observe the properties of other internal modes, we have perturbed the beating DD soliton with the corresponding eigenmodes. In particular, a perturbation along the direction of the localized mode depicted in black in the left panel of Fig. 4, leads to a breathing in the width of the soliton — see central panel of Fig. 5. On the other hand, a perturbation along the direction of the mode depicted in green in the left panel of Fig. 4, leads to an oscillation of the condensate along the trap, leaving the beating DD soliton unaffected (i.e., the soliton stays at the trap center) — see right panel of Fig. 5. For progressively weaker traps, the modes of the background condensate and of the dark-bright solitary wave essentially decouple and in fact two of the frequencies shown in Fig. 4 (green and blue) tend to 0, as the corresponding motions (of the solitary wave through the background or of the background through the solitary wave) become neutral. Finally, we make a remark about the way we have calculated the value NBthat must be introduced in Eq. (22). The procedure consists in performing an SO(2) rotation with χ=−π/4to the periodic DD soliton at t= 0. This solution is shown in the left panel of Fig. 6, whereas the rotated solution is depicted in the right panel of the same figure. Thus, NBis the norm of the bright component of the rotated solution. It can also be inferred from the Fourier coefficients of the periodic orbit: NB=∫R|u2|2dx=1 2∫R [|∑ k zk|2+|∑ k yk|2−2Re((∑ k z∗ k)(∑ k yk))]dx. (33) 8 −20 −15 −10 −5 0 5 10 15 20 −2 −1 0 1 2 x u1(x), u2(x) −20 −15 −10 −5 0 5 10 15 20 0 0.5 1 1.5 2 2.5 x |u1(x)|2, |u2(x)|2 −20 −15 −10 −5 0 5 10 15 20 −2 −1 0 1 2 x u1(x), u2(x) −20 −15 −10 −5 0 5 10 15 20 0 1 2 3 x |u1(x)|2, |u2(x)|2 FIG. 6: (Left panel) Profiles and densities of a periodic orbit at t= 0 with µ= 3,ω0= 0.5and Ω = 0.2. The right panel shows the dark-bright soliton arising by rotating with χ=−π/4the dark-dark soliton of the left panel. x g11:g12:g22=1:1:1 100 200 300 400 −5 0 5 t x 100 200 300 400 −5 0 5 g11:g12:g22=1.03:1:0.97 100 200 300 400 −5 0 5 t 100 200 300 400 −5 0 5 0.2 0.4 0.6 0.2 0.4 0.6 0.2 0.4 0.6 0.2 0.4 0.6 FIG. 7: The comparison between the integrable case g11 :g12 :g22 =1:1:1(left column) and the non-integrable case g11 :g12 :g22 = 1.03 : 1 : 0.97 (right column), is demonstrated. The upper panels show the densities of the first dark soliton component while the lower ones show the second dark component. Here η= 0.6, χ =π/4, θ = 0, k = 0, µ = 1. Based on the similarity of the relevant dynamics, we will focus on the case of unit nonlinear coefficients. IV. SINGLE BEATING DARK-DARK SOLITONS NEAR AND FAR FROM THE MANAKOV LIMIT We now turn to a numerical study of the properties of the beating DD soliton states. Firstly, in the absence of a trap, we are going to compare the integrable case with equal scattering lengths g11 :g12 :g22 = 1 : 1 : 1 to the non-integrable case g11 :g12 :g22 = 1.03 : 1 : 0.97 (see Ref. [30]). From Fig. 7, we observe that both of the dark components are oscillating with fixed frequencies and these two cases are very similar [38] All of the runs reported below for one of these parameter sets have been repeated with the second one and in all cases we have observed a close similarity between the dynamical phenomenology of these two cases. To highlight the fact that substantial variations of the scattering length–which can be imposed by virtue of a Feshbach resonance–may have a significant impact on the robustness of the beating DD solitons, we consider scattering lengths in the set with ratios g11 :g12 :g22 =g: 1 : 1. In particular, we take g= 1.1,1.2,1.6in Fig. 8. When gis not so large, i.e., g= 1.1,1.2, the beating DD soliton oscillates and, as tincreases, the change in the scattering length results in mobility of the coherent structure. However, more dramatic events can arise when gis relatively large, e.g., for g= 1.6. There, we can see that the soliton is finally split into two fragments (upon growth of the intrinsic beating oscillation which eventually induces the splitting) and results in two states that resemble dark-antidark solitons [31] (see also Ref. [29]). In particular, each of the components acquires a dark soliton coupled to an anti-dark soliton, i.e., a density hump (instead of a dip) on top of a finite background, in the second component. In Fig. 9, we show a particular example of the DD soliton in the trap, which oscillates around its center; the parameter values are µ= 1,η= 0.6, initial soliton position x0(t=0)=2.5, and trap strength Ω = 0.05. Note that for these runs, the initial profile of the beating DD soliton in the trap is approximated by the numerically found (in trap) ground state — i.e., the Thomas- 9 g=1.1 200 400 600 800 −5 0 5 200 400 600 800 −5 0 5 g=1.2 200 400 600 800 −5 0 5 200 400 600 800 −5 0 5 t g=1.6 200 400 600 800 −20 0 20 t 200 400 600 800 −20 0 20 0.2 0.4 0.6 0.2 0.4 0.6 0.2 0.4 0.6 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 1 FIG. 8: The dynamics of the densities of a DD soliton in the absence of a trap but for g= 1.1, g = 1.2, g = 1.6respectively (g≡g11); upper and lower panels show the densities of the first and second dark soliton component, respectively. When g= 1.1or g= 1.2, the soliton is set into translational motion. However, for g= 1.6the coherent structure executes a growing oscillation which eventually results in its splitting into a pair of dark-antidark solitons (i.e., a dark soliton in the one component coupled to a lump in the other). The parameters used here are the same as in Fig. 7. t x 200 400 600 800 −20 0 20 t x 200 400 600 800 −20 0 20 0 500 −2 0 2 t x0 2 4 6 8 10 12 14 0 0.02 0.04 r ωosc 0 0.5 0 0.5 FIG. 9: Motion of a dark-dark soliton in a trap of strength Ω = 0.05. Other parameters are the same as in Fig. 7. The upper two panels show the oscillation of the soliton initially centered at x0= 2.5(the chemical potential is µ= 1). The lower left panel demonstrates the center of mass of the DD in the upper panels. The analytical oscillation frequency, given by Eq. (22), is 0.03123, while the numerical frequency, calculated by Fourier transform, is 0.03238. The lower right panel yields the comparison between the analytically calculated frequencies (red line) versus the numerical obtained ones (the blue triangles), as rvaries between 0.1 and 14. Fermi cloud — multiplied by the beating DD solution (without a trap) of Eqs. (13)-(14). Then via a time-stepping algorithm (a fourth-order Runge-Kutta scheme), we obtain the time evolution of the densities of the oscillating solitons in the upper two panels. Moreover, the left lower panel shows the center of mass of the beating DD soliton in the trap. Using Fourier analysis, we can infer the numerical frequency of in-trap oscillation, which can, in turn, be compared to the analytical one, cf. Eq. (22). As shown in the bottom right panel of the figure, there is very good agreement between the two. Next, we consider the in-trap dynamics of a single beating DD soliton but for the non-integrable cases. Again, when g11 : g12 :g22 = 1.03 : 1 : 0.97, we observe a nearly identical phenomenology to that of unit gij’s. For the more significant deviations from that case of the form g11 :g12 :g22 =g: 1 : 1 where g= 1.1,1.2,1.6, the results are reported in Fig. 10. For lower values of g= 1.1,1.2, the behavior of the DD is similar to the case with g= 1, however, we progressively observe more significant radiative emissions which also affect the oscillation frequency. However, once again the modifications of the phenomenology are most dramatic in the case of g= 1.6of the bottom panels. There, the radiation emission is accompanied by growing intrinsic oscillations which eventually result in the breakup and formation of a single dark-antidark solitary wave.