Coherent impurity transport in an attractive binary Bose–Einstein condensate

We study the dynamics of a soliton-impurity system modeled in terms of a binary Bose–Einstein condensate. This is achieved by ‘switching off’ one of the two self-interaction scattering lengths, giving a two component system where the second component is trapped entirely by the presence of the first component. It is shown that this system possesses rich dynamics, including the identification of unusual ‘weak’ dimers that appear close to the zero inter-component scattering length. It is further found that this system supports quasi-stable trimers in regimes where the equivalent single-component gas does not, which is attributed to the presence of the impurity atoms which can dynamically tunnel between the solitons, and maintain the required phase differences that support the trimer state.


Introduction
Multi-component matter plays host to a plethora of novel phenomena, at both the classical and quantum mechanical level. The coexistence of several coupled, interacting degrees of freedom can facilitate different phases of matter, such as the miscible-immiscible phase-separation of binary fluids, arising from energetic competition between the differing components of the fluid [1].
Quantum fluids-systems of interacting particles comprised of Fermions or Bosons cooled below their respective degeneracy temperature, can now be used to give direct insight into many analogous systems due to their high degree of experimental controllability. In particular, it is now feasible to engineer the dimensionality [2, 3], particle interactions [4] and potential landscape [5] of these macroscopic systems. Complementary to this, the optical manipulation of these systems has reached maturity-oppurtunities now exist to emulate complex phases of matter in the prescence of gauge fields [6,7], which form a key ingredient for many condensed matter effects of interest.
Solitary waves have been produced experimentally in both single and multi-component condensate systems. In the former case, quasi-stable soliton states have been generated, comprising single [8] as well as trains of bright solitons [9,10]. Further work demonstrated bright solitons sensitivity to surface physics in the form of both repulsive [11] and attractive potentials [12]. Understanding the observed stability of these fragile systems has revealed the important role the complex phase of the matter-wave plays in these systems [13,14]. Matterwave solitons have been touted for applications in metrology, where these state's inherent coherence advocates them as strong candidates for engineering matter-wave interferometry [15][16][17][18]. This in particular has led to the realisation of a matter-wave bright soliton Mach-Zehnder interferometer with a 85 Rb condensate [19], as well as proposals to controllably split solitons [20], and very recently schemes to realise bright soliton states with minimal noise have appeared [21]. The purity of cold atom systems has also been exploited to gain insight into the role disorder plays for the dynamics of bright solitonic states in cold atomic gases [22,23].
There have also been experimental realizations of solitary wave structures in multi-component systems. Early theoretical work studied the properties of dark-bright and bright-bright solitons [24,25] the first of which was subsequently realized individually [26] and also in the form of trains [27]. As well as this, studies have focussed on the role of potential barriers in the dynamics of vector solitons [28]. Theoretical work has predicted that the single component focussing nonlinear Schrödinger equation can possess chaotic solutions in the presence of an axial harmonic potential [29,30], as well as the observed interaction induced frequency shift of pairs of trapped bright solitons [31] in the experiment of [13]. Complementary to this, theoretical work has focussed on solitary waves in higher spin systems, revealing the existance of integrable points in the full parameter space of the spin-1 condensate, in the form of so-called 'polar' bright solitons [32,33]. Although solitons are usually studied as the solutions to one-dimensional nonlinear models, there have also been predictions of stable two-dimensional solitary wave solutions in dipolar Bose-Einstein condensates [34,35], where the additional nonlocal nonlinearity provides the stabilizing mechanism for these solitons. Very recently the Jones-Roberts soliton was realized experimentally, a true two-dimensional solitary wave structure [36].
The realization of artificial electromagnetism with cold gases, and in particular spin-orbit coupling for Bose-Einstein condensates opens a novel route towards studying nonlinear wave structures. Here, the coupling of the condensates momentum to a quasi-spin leads to stripe-like soliton phases, related to the underlying immiscible phase of these systems [37,38]. Spin-orbit coupling forms a key ingredient in simulating more exotic scenarios, such as Dirac-like equations, where confined solutions have been predicted [39] that resemble their bright soliton cousins in single component condensates.
Atomic condensates benefit from being exceptionally pure systems-this in turn allows one to investigate the effects of disorder and defects with an unprecedented level of control. The presence of impurities in ensembles of ultracold matter has led to predictions of impurity-molecules and lattices at the mean-field level [40], as well as the role of many-body correlations for a single impurity out-of-equilibrium [41]. Experimental work has studied the role that spin impurities have in the strongly correlated Tonks-Girardeau limit [42] and also magnetic spin models [43], which have also been the focus of subsequent theoretical investigations [44][45][46]. Complementary to this, recent experimental advances have led to the realization of trapping one matter-wave inside another, where a degenerate Fermi gas of 6 Li atoms was confined inside a 133 Cs Bose-Einstein condensate [47].
The ability to both prepare and control ultracold gas experiments gives access to physical regimes that mimic and go beyond those associated with conventional condensed matter physics. Impurities play a central role in condensed matter, since most materials will contain some imperfections. One important example drawn from this field is the polaron, a quasi-particle that consists of an electron and the distortion caused by the passage of the electron through the ionic lattice. Impurities in the form of polarons can act as a sensitive probe within manyparticle systems, and can be used to explore the correlations of these systems. Additionally it should be noted that polarons are not necessarily dependent on the prescence of impurities in a material, they can also appear in ideal crystals. Over the last few years, ultracold gas experiments have succeeded in simulating the physics of polarons, including the pioneering experimental realisation of polarons of both bosonic [48,49] and fermionic [50] gases. The physics of polarons has also formed an ongoing focus of theoretical investigations. Optical lattices yield access to many models of interest in condensed matter physics, however as they are constructed from the interference of two counter-propagating laser modes, they do not naturally yield lattice vibrations (phonons), a key ingredient for polaron physics. This important question was investigated in [51], which proposed a methodology to overcome this drawback. Further work investigated the effect of dimensionality on the self-trapping of impurities, revealing regions where stable polarons can exist [52]. Very recently, a theoretical investigation has revealed the universal behavior of the bosonic polarons energy and its dependence on the Efimov parameter [53].
In this publication we will outline the collisional dynamics of multicomponent soliton-impurity systems, and how binary or triplet collisions might be exploited to perform deterministic population transfer operations on the impurity, providing a toolkit for future applications to metrology and quantum computation. The soliton-impurity system at the heart of our work is shown schematically in figure 1, where two soliton isosurfaces are shown with the delocalized impurity component. The paper is organized as follows. In section 2 we examine the stability of this system using a full three dimensional variational approach in order to understand the regimes where stable dynamics can be realized. Then in section 3, we state the model for the two component system in terms of coupled mean-field Gross-Pitaevksii equations for the dynamics. After this in section 3.1 we explore the ground states of the binary system, following which in section 3.2 we undertake a scattering analysis of a single soliton molecule carrying an impurity with an 'empty' soliton. We then proceed to show how soliton molecule complexes can be built using three solitons in section 4, and study the resulting nonlinear dynamics of the solitons and impurity as a function of the relative phase and inter-component scattering length, revealing the coherent nature of the impurities dynamics. We also discuss the conditions under which this state is stable to thermal fluctuations, before demonstrating that the impurity undergoes a novel localization transition. We conclude with a summary of our findings in section 5.

Soliton-impurity stability
The majority of experiments with atomic condensates are realized with repulsive inter-particle interactions confined by harmonic potentials. Under these conditions, the condensate is unconditionally stable. The introduction of attractive interactions can lead to a collapsed state, originating in the dispersive kinetic energy of the gas being overwhelmed by the attractive interactions between particles. We consider a two component model, where the second component of the system can be modeled as an 'impurity', since the mass (number of atoms) of either component can be independently varied [54,55]. We consider a two component (binary) system forming a Bose-Einstein condensate coupled via completely attractive mean-field interactions. The stability of such a system depends on a number of parameters, in-particular the various scattering lengths, the number of atoms in each component and also the trapping geometry. To gain insight into the collapse dynamics of the binary system, consider the energy functional where the wave function of component j is Y º Y( ) r j , and the s-wave scattering length a jk is contained in the parameter where m is the atomic mass. Note that in this system there are only two scattering parameters depending on the various scattering lengths, g 11 and g 12 , while g 22 =0, and so the second component is linear and moves in the effective potential defined by the first component. The single-particle Hamiltonian H 0j appearing in equation (1) is defined as where ω ⊥ defines the transverse trapping frequency of the cloud, and = + y z r 2 2 2 defines the radial coordinate. Then, this problem contains three length scales, two associated with the two scattering lengths, as well as one from the harmonic trapping term appearing in equation (2). The collapse instability for the cylindrically symmetric single-component gas has been studied previously, including the effect of additional axial confinement [56]. To understand the nature of the collapse, we employ the cylindrically symmetric Gaussian variational ansatz where N j is the atom number in each component. Note that the ansatz of equation (3) is appropriate since both scattering lengths are attractive, so the system is miscible with both components spatially overlapping. It is also possible to consider the immiscible case, where one scattering length is repulsive and the other attractive, which has also been shown to support stable solitary wave structures [57] in a quasi two-dimensional scenario. We can then insert the ansatz equations (3) into (1), yielding  Figure 2 shows the numerically obtained solutions to equation (5). These solutions are obtained using an iterative procedure to procure the collapse point starting from a point in the parameter space with known analytical solution, in this case the point g 12 =0, from which the critical collapse point for a cylindrically symmetric trap is = -Na a 1 3 where N is the atom number and a s the s-wave scattering length. The mean-field collapse phase diagram is shown in figure 2(a), the volume enclosed by the a 11 , a 12 and N 2 /N 1 axis define the space of stable three dimensional solitons. Here the red lines show the boundary between stable and unstable regimes in each parameter plane. It can be seen that when N 2 /N 1 is small, corresponding to a small impurity population the collapse point is moved to larger values of a 12 . As the number of atoms in the impurity N 2 increases, the collapse point in the a 12 -a 11 plane moves to smaller values of a 12 . This result is intuitive, since one can interpret the additional attractive inter-species mean-field potential as providing an extra destabilizing contribution to the mean-field energy. shows a different cut through (a) for constant N 2 /N 1 . The green shaded region bounded by the blue dashed line indicates the stable region for N 2 /N 1 =0.1 in the a 11 -a 12 plane. The blue shaded region is stable for N 2 /N 1 = 0.01, but not N 2 /N 1 =0.1. Again, the white region is unstable to collapse. This rudimentary analysis shows that any experiment to realize a fully attractive two-component system would be favorable to a moderate mass imbalance, especially if one was interested in exploring the dynamics as a function of one of the scattering lengths of this system, as we will proceed to do in the following sections of this work. Since we consider a mean-field mass-imbalanced system it is worth considering when such a model is valid. It is known for example that on the repulsive side (a jk >0), this model undergoes composite fermionization [59]. We expect this attractive meanfield model to be suitable up to the collapse point, although it is conceivable that fluctuations could play an important role in this mass imbalanced system. However, this analysis lies beyond the scope of the current work.

Equations of motion
One of the characteristic attributes of solitary waves are their particle-like properties [60]. Consequentially, their inherent robustness leads to collision dynamics where they emerge unscathed, with the exception of a phase shift. For the single-component focussing nonlinear Schrödinger equation, the scattering of two bright solitons is always elastic, a consequence of the underlying integrability of the nonlinear Schrödinger equation. For the two component system the equations of motion for Y ( ) t r, j are found from the Lagrangian density and the associated Euler-Lagrange equations. We are interested in studying the soliton solutions which exist in the quasi one-dimensional limit. As such, we assume that there is tight radial confinement, such that any radial dynamics are effectively frozen out. Then the radial dynamics for both components can be factorized in the form defines the ground state of the radial trap. Proceeding, the dynamics in the quasi one-dimensional limit are captured by The equations of motion defined by equation (7) will form the work-horse for studying the binary attractive system. This mean-field model was originally studied by [61] who analyzed the localized solutions and their quantum fluctuations. We note that this model has also been studied recently in the context of repulsive meanfield interactions, where it was shown how the dark soliton solutions long lifetimes can be used to host qubits for quantum information applications [62]. Complementary to this the physics of polarons remains a topic of ongoing interest, with very recent theoretical work focussing on studying so-called Frölich polarons [63]; as well as the binding properties of trapped bosonic polarons [64].

Single polaron ground states
To understand the basic physics of the attractive binary condensate defined by equation (7), we begin by computing the ground state of this system as a function of the inter-component scattering parameter g 12 . This is shown in figure 3(a), which shows the density y | | . For large negative values of the inter-component scattering length, the impurity is well localized within the soliton. As  g 0 12 , the width of the impurity wave function starts to grow. This effect is investigated further in panel (b), where we compute the effective width (standard deviation) of the impurity, = á ñ ℓ x i 2 as a function of g 12 , for different mass ratios N 2 /N 1 . Each individual data set is scaled to the width of the impurity for the largest negative scattering length (ℓ i0 ) for ease of comparison. One can see that as

Binary soliton-impurity dynamics
To gain insight into the dynamics of the soliton impurity system, we simulate collisions between an 'empty' bright soliton, that is a soliton solution obtained for g 12 =0 from equation (7) with the soliton containing an impurity. As such, our initial condition takes the form where y j represents the numerical solution to equation (7) for a given finite choice of g 11 and g 12 . The first component ψ 1 contains the first bright soliton which hosts the impurity (ψ 2 ) localized at the origin x=0. The length scale (size) of each component depends critically on the scattering parameters g 11 and g 12 . The function y ( ) x S is the single soliton solution given by 0 and the length scale appearing in equation (10) is 11 1 , with N 1 giving the number of atoms in each soliton of the first component ψ 1 , while the scaled quasi-one-dimensional scattering parameter is = g 11  wâ 2 11 which describes the solitonic nonlinearity strength. The initial phase difference is given by δ. The parameter space associated with equation (7) contains two scattering lengths, two atom numbers, the initial velocity v 0 and position x 0 of the soliton and impurities as well as the initial phase difference, and as such is generally complicated to understand completely. To draw out the main features of the model, we simulate collisions for fixed g 11 and atom number, but vary the inter-component scattering length g 12 .
A useful measure for soliton collisions with non-integrable dynamics is the coefficient of restitution. This is a dimensionless quantity defined as the total kinetic energy of two particles after a collision to the total kinetic energy before the collision [65,66] Now, if η=1 the collision is elastic with conserved momenta before and after collisions, while h ¹ 1 indicates an inelastic collision between the solitons. The masses m i and velocities v i appearing in equation (11) are computed from Both quantities appearing in equation (12) are computed locally around the center of mass of each individual soliton.
In our simulations presented in figure 4 we have taken is the initial displacement of the empty soliton. The , while each simulated collision is run for t=500τ units of real time. The numerical simulations are handled using a spectral (split-operator) method, and we work in the so-called soliton units [30], where ℓ = ÿ/mv, τ = ℓ/v and E = mv 2 define the units of length, time and energy respectively. To understand how these units correspond to physical quantities, we can use the experimental parameters of [11], who produced a bright solitary wave with a 85 Rb condensate. Then one has = m u 85 , where u is the atomic mass unit, N 1 =2000 atoms in each soliton, and a transverse trapping frequency of ω ⊥ =27 Hz. Using these parameters one finds a natural length scale ℓ ; 11 μm, and small value of  -- ℓ N mg 10 1 2 12 2 4 , so in reality it would be necessary to use the powerful tool of Feshbach resonances in an experiment in order to bring the system into the regime described in this work.
In figure 4 we explore the binary dynamics of the soliton-impurity system. Here, a single empty soliton collides with the soliton-impurity system which is positioned initially at the origin. The coefficient of restitution, equation (11) is then computed as a function of the inter-component scattering length, g 12 for several initial phase differences, δ. It should be noted that although the phase difference is set initially in our simulations, this quantity evolves dynamically [67], so the initial value is not necessarily the phase difference at the point of collision. This phase evolution can be inferred from η as displayed in figure 4(a), where we plot η, as a function of g 12 , resulting from initial phase-differences δ=−π/2, 0, π/2, π. The dynamics of η can be roughly partitioned into two regimes, a 'weak' dimer phase (light-blue shading) that manifests for , and a second more conventional non-integrable regime with In the non-integrable regime the dependence of η on the solitons' relative phases is demonstrated by the phase-winding of the η oscillations associated with the different δs. Furthermore, η is seen to oscillate with an increasing amplitude and frequency as g 12 ʼs magnitude is increased. The frequency increases because increasing g 12 increases the chemical potential of the carrier soliton, and so its phase winds more quickly, effecting an additional phase-shift. Due to the nonintegrability of this system, energy is not conserved if h ¹ 1. The energy of the solitons is redistributed postcollision. The primary mechanism for this is the change of the solitons masses post-collision. As well as this, some of the kinetic energy initially carried by the moving soliton is redistributed into potential energy postcollision, affecting an additional perturbation to the systems two scattering lengths. In the 'weak' dimer phase we note that all four curves meet at g 12 =0 when η=1, where integrability is restored and the system is reduced to the single component focussing cubic  , the (dimensionless) initial velocity of the (empty) soliton, again for four equally spaced phase differences, δ=−π/2,0, π/2, π.
Sol. j 2 2 which are calculated by integrating the density of the impurity locally around each solitons center of mass. Then, for a system with n sol solitons the total impurity population is given by The panels on the right column figure 4(c) show the equivalent dynamics but for  - ℓ N mg N 0.24 1 12 2 1 . Clearly the dynamics of this system are highly non-integrable, showing a number of unusual dynamical effects. In particular, the second impurity component does not behave like a soliton, instead behaving like a quantum particle trapped by the potential generated by the first, solitonic component. Then, by exploring the scattering dynamics as a function of the inter-component scattering parameter g 12 , per figure 4(a) one can interpret the dynamics of the system. For larger negative values of g 12 , the impurity is localized deep within the soliton, and as such its length scale is typically less than that of the length scale (size) of the soliton in which it is initially localized. On the other hand, for smaller negative values of g 12 the potential felt by the impurity component is quite shallow, and the impurity in this regime is comparatively more weakly bound, having a length scale which can be significantly larger than that of its solitonic host.
Since the impurity component feels the solitonic component as an effective dynamical potential, its dynamics can show some unusual features. In particular, the impurity can transfer itself into the other, empty soliton. This is shown in the lower panel of figure 4(b), where the impurity population of each soliton is computed as a function of time. During dynamical evolution, ∼85% of the initial impurity population is smoothly transferred from the second to the first soliton. In the second example shown in figure 4(c) the initial population of the second soliton is almost completely transferred to the first, and then back again. We attribute the population transfer effect to quantum mechanical tunneling. It is also worth mentioning that there is an additional subtlety in the interpretation of these dynamics. In the integrable limit, each soliton can be identified by its amplitude and velocity [68], which means that the transfer of population between the two solitons in figure 4(b) could also be interpreted with the labels of the two solitons switched, post collision. We will instead keep the labeling of the solitons more in the style of two potentials that the second impurity potential feels. This choice makes a quantitative but not qualitative difference to the interpretation of our results. These dynamics could be useful for atomtronics applications [69], indeed these examples shown in figures 4(b) and (c) behave somewhat like an analog of a conventional transistor, where the first soliton that initially hosts the impurity can be interpreted as the 'source' while the second empty soliton can be labeled the 'drain', while the effective gate voltage is controlled by the inter-component scattering parameter, g 12 . , as a function of the initial velocity of the empty soliton,  ℓ v m 0 . For large initial velocities, the scattering is comparatively less sensitive to the scattering length g 12 . However, as the initial velocity is lowered, a prominent dip develops, whose depth and position on the g 12 axis depends on the initial phase difference δ. This effect is demonstrated in figure 5, which shows the data displayed from figure 4 reshaped. Each curve is taken for the fixed velocity Here, one can see that the position and depth of the 'dip' is quite sensitive to the initial phase difference, and seems to deepen for smaller scattering lengths as the phase difference is modulated. The lower row of panels in this The relative amount of kinetic to (attractive) potential energy in this system is crucial to the observed dynamics. Indeed, for collisions approaching zero inter-component scattering length, one has a highly delocalized impurity, which is weakly bound in its solitonic host. Coupled to this is the fact that the collision of the solitons in this system expels radiation in the form of small amounts of atomic density of the cloud. This can in turn interact with the solitons in this regime to further destabilize the observed dynamics. For scattering lengths slightly smaller in magnitude than figures 5(b) and (c), the post collision dynamics are found to be exceptionally sensitive to this radiation. This can partially be overcome by simulating collisions with increasingly larger numerical boxes, (in our work we typically use L box =200ℓ) however as  | | a 0 12 the size of the impurity will always be larger than one can realistically simulate, an unavoidable limitation inherent to this system.

Soliton molecules 4.1. Soliton trimers
Models of nonlinear systems can also play host to higher-order soliton states, in the form of soliton molecules ( i.e. several individual solitons forming bound objects) and also breathers, which are single solitonic entities that can be thought of as excited states of the focussing nonlinear Schrödinger equation, which have recently been engineered experimentally with matter waves for the first time [70] using an attractive gas of 85 Rb. Soliton molecules have been studied in various guises within the context of ultracold matter, for example the realization of degeneracy with atomic species possessing significant dipole-dipole interactions has led to the prediction of novel molecular states in these systems [71][72][73]. Related to this are the realisation of 'droplets' of both dipolar matter [74][75][76] and intriguingly, also light with a non-trivial angular momentum structure [77], as well as the prediction of soliton molecules in systems with nonlocal interactions [78].
Here, we consider a stationary spatially symmetric initial state to study the possibility of molecule-like states in the system described by equation (7). In the previous section it was found that the low velocity scattering of a pair of solitons leads to increasingly inelastic dynamics as the initial kinetic energy in the system approaches zero. In fact, if we try to form a simple soliton molecule with a pair of initially stationary solitons, one with an impurity, and one without the resulting molecular state rapidly destabilizes. This is due to the impurity that causes the phase of the soliton in the first component to wind, eventually breaking the molecule. Instead we focus on understanding molecules formed from three individual solitons. In order to create an initially symmetric state, we must place the impurity either in the center soliton with the outer two solitons initially empty, or visa-versa. In our simulations we have chosen the former, so that the initial state is S and x j are the centers of mass of the two outer solitons, ψ S (x) is defined per equation (10), and the x j are chosen symmetrically such that x 1 +x 2 =0. This initial configuration, built from the ground state and known exact solutions in the limit g 11 =0 is also shown in figure 3(c). Figure 6 shows example dynamics of the three soliton system. From left to right, panels  (c) respectively. The initial phase of the central soliton is δ=π/2. As the inter-component scattering length g 12 is increased, the dynamics of the system change quite drastically. This is reflected in that fact that in (a) the solitons move apart, with only one 'switch' of population occurring during the dynamics, as shown in the lower panel of (a). As g 12 is increased, the impurity is delocalized, promoting tunneling to the outer solitons, as shown in (b). Finally in (c) a molecular-like state is

Thermal fluctuations
Given the fragile nature of bright soliton states, it is important to understand when the predicted soliton trimer presented in figure 6 is stable to thermal fluctuations that are present in real systems. One way to understand the conditions under which the trimer is stable to thermal fluctuations is to compare the energy difference between the absolute ground state of the system and the trimer state with the thermal energy present in the system. We denote each of these quantities by E gnd and E tri respectively. Then the energy difference d = -E E E tri gnd we are interested in is given by which demonstrates that δE depends only on the outer solitons, and not the central soliton that carries the impurity. We can use the known analytical expression for the stationary bright soliton (equation (10)) profile to obtain an exact expression for δE. The total energy of each outer soliton is where N 1 and g 11 are the atom number and the quasi one-dimensional scattering parameter associated with the first component. To understand when thermal fluctuations play a role, we can form a dimensionless figure of merit as the ratio of the energy difference and the thermal energy present in the system at temperature T as   . This rudimentary argument suggests that producing stable dynamics requires both a reasonable atom number as well as a low temperature. Since equation (19) depends on the cube of the atom number, it should in principle not be too difficult to satisfy this condition. Alternatively, one could also calculate the thermal stability of the trimer state from the quantity E k T B tri alone. This could also give a deeper insight into the parameter regimes where this state is stable to thermal fluctuations and importantly how the inter-component scattering length g 12 affects this stability.

Coherent impurity dynamics
The dynamics of the impurity, presented in figure 6(c) are suggestive that the soliton molecule could host coherent population dynamics. To investigate this effect we perform a comparison of the dynamics of the impurity component with a simple three level system, modeled in terms of a 'vee' type atom. This model is chosen since the ground state energy (chemical potential) of the central soliton is slightly lower in energy; due to the presence of the impurity. Then, the equations of motion for the complex amplitudes c j (t) that determine the population of each soliton are We can connect the solutions c j (t) of equation (21) to the populations presented in figure 6 since = # ( ) | ( )| P t c t j Sol. j 2 . The dynamical system described by equation (21) introduce the 'Rabi' frequency Ω, which defines the frequency of population transfer between solitons, and the effective 'detuning' Δ. The total population is a conserved quantity given by å = | ( )| c t N j j 2 2 . Figure 7 shows comparisons of the Rabi model, equation (21) with Gross-Pitaevskii simulations. The analogous 'vee' atom level diagram is shown in figure 7(a), where the states ñ ñ | | 1 , 2 and ñ |3 represent the potential generated by the left, middle and right soliton felt by the impurity atoms. Then one can associate a state vector with equation The solutions given by equation (23) can be used to gain insight into the nature of the underlying tunneling effect responsible for the impurities transport inside the solitons. To do this, we calculate the tunneling current using the solutions for c j (t) from equation (23) giving and n is zero or a positive integer, which as can be seen from figures (7) (c) and (d) is exactly when the impurity in the inner soliton (soliton two) has a maximum in its impurities population. Likewise, the tunneling current  ( ) t t goes to zero when p W = t n d for even integer n; which corresponds to when the outer solitons (soliton one and two) have their maximum impurity population.
To numerically obtain the Rabi frequency Ω d , we take the Fourier transform  In both presented cases, figures 7(c) and (d) we find excellent agreement to the Rabi model. It is important to note that at much longer times, the outer solitons are attracted towards the central soliton, which causes the effective Rabi frequency Ω d to increase, but by sensibly choosing the system parameters such that the outer solitons are not initially too close to the central soliton, good agreement to equation(21) is obtained. The coherent oscillations presented in figure 7 could form the basis for future applications. In particular, the identification of these types of dynamics could find practical application in atomtronics [69,79] and quantum information processing [80], where the coherent dynamics of atomic systems are a required ingredient for many effects of interest in these fields.

Impurity localization transition
The dynamics of the impurity presented in figures 6 and 7 are suggestive of rich transport behavior. To understand the transport properties of the multi-soliton system further, we probe the dynamics of the three soliton system across the full parameter space. Obtaining a well-behaved, intuitive measure for the multiple soliton system is challenging. To understand the effect of the various system parameters on the dynamics of the impurity, we employ the inverse participation ratio (IPR) as a measure to quantify the dynamics of the system. The IPR provides a well-behaved measure of how localized a particular state is [81]. Related to this, recent work has also examined the effect of 'dynamical localization' in dynamical optical lattice potentials [82]. In particular we wish to calculate For non-interacting spatially localized states, the IPR takes a value of one such that  = -( ) t 1 1 , while delocalized states are found instead when  - ( ) t 1 1 . This definition is however strictly speaking only applicable to non-interacting systems, the introduction of mean-field interactions can yield value of the IPR that are greater than one. Nonetheless this quantity still provides a useful measure of the impurities spatial dynamics. Since we are dealing with a two component system where both components evolve dynamically, it is necessary to consider the time averaged version of equation (25) in order to make a meaningful analysis. The time average of equation (25) is defined as where T defines the length of the particular numerical simulation. To investigate the behavior of the IPR, equation (25), and in particular its time average given by equation (26), we perform numerical simulations using the initial state defined by equation (15), with the parameters N 1 =1000 and N 2 =1 giving the atom numbers for the soliton and impurity respectively, while  =ℓ N mg N 8  figure 8(a). Here, data is presented for several different initial phase differences: δ=0, ±π/2, ±π. The dynamics can be divided into three regions, a localized region (red gradient), a delocalized region, (blue gradient) and an intermediate region (white). The trend in (a) shows that for large negative g 12 the impurity is localized, since   1, here all data fall onto a common curve. The observed behavior of the IPR attaining values greater than one notably differs from its original definition where localized states are defined for  á ñ = ( ) t 1 only. We attribute this departure to the fact that we are considering an interacting, rather than non-interacting system. Then as the scattering length is increased, the impurity starts to delocalize across the three solitons and individual datum no longer follow a common trend, instead the particular value of  á ñ -1 one obtains is found to be sensitive to the initial phase δ. As the scattering parameter g 12 approaches zero, the data again fall onto a common curve, and the impurity is completely delocalized between the three solitons. In this region stable molecules are found that support this effect.
To understand the impurity dynamics in the intermediate region, (white region in figure 8) the fluctuations during dynamics of the IPR (equation (25)) are studied by calculating the standard deviation in figure 8(b). The standard deviation of the IPR is plotted with the average of the IPR (light blue shading and solid blue respectively). One can see that the fluctuations associated with equation (26) start to grow as  á ñ -1 falls below one. Indeed, it would seem within the mean-field model considered in this work one can attribute the point  á ñ = ( ) t 1 as the point in the parameter space where fluctuations of the IPR grow from zero and the impurity begins to delocalize between the outer solitons. The final row of figures shown in figure 8 shows example dynamics for each dynamical region. In particular figure 8(c) shows an example of a localized impurity. Then

Conclusions
In this work we investigated the scattering properties of a two-component Bose condensate with wholly attractive mean field interactions. By interpreting the second component as an impurity, this system was found to support unusual transport phenomena, including the appearance of a dimer like phase close to zero inter-component scattering length, where a pair of bright solitons in the first component can coherently transfer the impurity between each other many times. Such an effect could be useful for example in the emergent field of atomtronics, where atomic systems are used to build circuits analogous to their electronic counterparts. The ability to use solitary waves to coherently shuttle atomic density over macroscopic distances could form a novel tool in this endeavor.
It was also found that stable soliton molecules formed from three solitons can also be produced in parameter regimes where the equivalent single component system is unstable to the formation of molecular bound states. This stability was attributed to the nontrivial phase winding that occurs during dynamical evolution of the twocomponent system. Since the impurity that constitutes the second component can effectively delocalize itself across the whole system, the atom number of both components of the gas can change. Accompanying this change is a winding of the phase, which for a critical scattering length can be favorable to the formation of three soliton molecules. The population dynamics of the impurity was scrutinized using a simple three level atomic 'Rabi' model. For sensible choices of parameters excellent agreement was obtained with GPE simulations. Finally, the trimer-impurity system was analyzed using the tools of localization theory. It was found that the impurity undergoes a delocalization as a function of the inter-component scattering length.
It would be interesting to investigate the effect of trapping fermions in this physical setup, in a similar spirit to the experiment of [47]. The ability to build larger systems of solitons with this particular system opens a novel avenue in studying lattices formed from solitary waves, with the twist that one can have different numbers of impurities present, which could be used to study effects analogous to condensed matter, for example a soliton-Hubbard model could be potentially explored, as well as understanding the generalized Toda lattice that this system would constitute.