Dark-bright solitons in Bose-Einstein condensates at finite temperatures

We study the dynamics of dark-bright solitons in binary mixtures of Bose gases at finite temperature using a system of two coupled dissipative Gross-Pitaevskii equations. We develop a perturbation theory for the two-component system to derive an equation of motion for the soliton centers and identify different temperature-dependent damping regimes. We show that the effect of the bright ("filling") soliton component is to partially stabilize"bare"dark solitons against temperature-induced dissipation, thus providing longer lifetimes. We also study analytically thermal effects on dark-bright soliton"molecules"(i.e., two in- and out-of-phase dark-bright solitons), showing that they undergo expanding oscillations while interacting. Our analytical findings are in good agreement with results obtained via a Bogoliubov-de Gennes analysis and direct numerical simulations.


Introduction
Macroscopic nonlinear excitations of atomic Bose-Einstein condensates (BECs) [1,2] have been a subject of intense theoretical and experimental research over the last few years [3]. More specifically, matter-wave dark and bright solitons, that can be formed in single-component BECs with repulsive or attractive interatomic interactions respectively, have been observed in a series of experiments while their statics and dynamics have been extensively studied theoretically in various settings (see, e.g., [4][5][6] for recent reviews). Of particular interest are coupled dark-bright (DB) solitons that may exist in binary mixtures of BECs with repulsive interatomic interactions (such as ones composed by different hyperfine states of 87 Rb atoms [7,8]): these solitons are frequently called symbiotic ones, as the bright soliton component (which does not exist in the system with repulsive interactions [4]) can be supported due to the nonlinear coupling with the dark soliton component. Such structures have recently been observed experimentally in a 87 Rb BEC mixture using a phase-imprinting method [9] or in two counter-flowing 87 Rb BECs [10,11], while they have also been studied in various theoretical works in continuum [10][11][12][13] and discrete [14] settings.
The above theoretical studies on atomic DB solitons have been performed in the ideal case of zero temperature: in fact, finite-temperature induced dissipation of matterwave solitons have basically been studied, so far, in the simpler case of dark solitons in single-component BECs [15][16][17][18][19][20]. In particular, this problem was first addressed in Ref. [15] (see also Ref. [16]), where a kinetic-equation approach, together with a study of the Bogoliubov-de Gennes (BdG) equations, was used. In that work, it was found that the dark soliton center obeys an equation of motion of a harmonic oscillator, which incorporates an anti-damping term accounting for the finite temperature effect. The presence of this term alters the soliton trajectories so that the experimentally observed dark soliton dynamics can be qualitatively understood: solitons either decay fast at the rims of the BEC (for high temperatures) [21][22][23] or perform oscillations of growing amplitude (for low temperature) [9,[24][25][26] and eventually decay. A similar equation of motion for the dark soliton center was also derived in Ref. [17] by applying the Hamiltonian approach of the perturbation theory for dark matter-wave solitons [6] to the so-called dissipative Gross-Pitaevskii equation (DGPE). This model incorporates a damping term (accounting for finite temperature), first introduced phenomenologically by Pitaevskii [27], and later shown to be relevant from a microscopic perspective (see, e.g., the review [28]). It is important to note that, as shown in Ref. [17], the analytical results obtained in the framework of the DGPE were found to be in very good agreement with numerical results obtained in the framework of the stochastic Gross-Pitaevskii equation (SGPE); see, e.g., Ref. [29] for a review on the SGPE model. It should also be mentioned that while the above works chiefly considered finite temperature effects for the case of a single dark soliton, the DGPE model and the anti-damping-incorporating ordinary differential equations (ODEs) for the soliton center were also examined in the case of multiple dark solitons. In particular, the cases of two and three oscillating and interacting, anti-damped dark solitons were considered in Ref. [30].
In the present work, we study finite-temperature dynamics of DB solitons in harmonically confined Bose gases. In particular, we adopt an effective meanfield description and analyze theoretically and numerically a system of two coupled DGPEs, describing the evolution of a binary quasi-one-dimensional (1D) BEC at finite temperature. We extend the considerations of Ref. [17] and develop a Hamiltonian perturbation theory for the two-component system at hand. This way, we obtain an equation of motion for the DB soliton center, similar to the one derived in Refs. [15,17]. This equation, which includes an anti-damping term accounting for finite temperature, provides a characteristic eigenvalue pair (i.e., a pair of solutions of the characteristic equation associated with the linear equation of motion), which is connected to the eigenvalue associated with the anomalous mode of the DB soliton. Performing a Bogoliubov-de Gennes (BdG) analysis, we show that the anomalous mode eigenvalue becomes complex as the dissipation (temperature-dependent) parameter is introduced, leading to an instability of the DB soliton pair. The temperature-dependence of the eigenvalues (determined analytically) is found to be in good agreement with the one of the anomalous mode eigenvalue (determined numerically).
Furthermore, these considerations are generalized in the case of a DB soliton "molecule", composed by two-DB-solitons. In the latter setting, both configurations featuring in-phase and out-of-phase bright components can be obtained in the trap [31]. We illustrate their dynamical instabilities as a function of temperature and capture them analytically by means of coupled nonlinear ODEs accounting for the three ingredients (trap restoring force, interaction between DB solitons and thermally induced antidamping). We show that, due to finite temperature, the nature of their interaction (and collisions) changes: for short times individual solitons behave as repelling particles, while for longer times they gain kinetic energy and completely overlap at the collision point. Our analytical considerations and numerical results reveal a fundamental effect: the partial stabilization that the bright ("filling") soliton component offers to the corresponding "bare" dark soliton against temperature-induced anti-damping. This way, a significantly longer lifetime of the symbiotic (dark-bright) structure can be achieved, in comparison to its bare dark soliton counterpart.
The paper is structured as follows. In section II we present the model and study some of its basic properties such as the evolution to the equilibrium state. In section III we develop the perturbation theory to derive and solve the equation of motion for the single DB soliton; we also compare our analytical findings to numerical results. In section IV we generalize relevant considerations to the case of multiple DB solitons, and in section V we present our conclusions. 3

The system of dissipative Gross-Pitaevskii equations
We consider a two-component elongated (along the x-direction) repulsive Bose gas, composed of two different hyperfine states of the same alkali isotope, and confined in a highly anisotropic trap (such that the longitudinal and transverse trapping frequencies are ω x ≪ ω ⊥ ). In such a case, the system can be considered as quasi-1D and, hence, the coupling constants take their effectively 1D form, namely g jk = 2 ω ⊥ a jk , where a jk denote the three s-wave scattering lengths (note that a 12 = a 21 ) which account for collisions between atoms belonging to the same (a jj ) or different (a jk , j = k) species. Let us now focus on the experimentally relevant case of a two-component BEC consisting of two different hyperfine states of 87 Rb, such as the states |1, −1 and |2, 1 used in the experiment of Ref. [8], or the states |1, −1 and |2, −2 used in the experiments of Refs. [10,11]. In the first case, the scattering lengths take the values a 11 = 100.4a 0 , a 12 = 97.66a 0 and a 22 = 95.00a 0 , while in the second case the respective values are a 11 = 100.4a 0 , a 12 = 98.98a 0 and a 22 = 98.98a 0 (where a 0 is the Bohr radius). In either case, it is clear that the scattering lengths and, accordingly, the effectively 1D coupling constants take approximately the same values, say a ij ≈ a and g ij ≈ g = 2 ω ⊥ a, respectively, which is what we will assume henceforth.
We now consider the case where the two-component Bose gas under consideration is at finite temperature. In particular, we assume that the thermal modes of energies > ω ⊥ are at equilibrium, accounting for a heat bath in contact with the axial part of the gas, while the modes in the x-direction are highly occupied so that the classical field approximation is valid [32,33]. Then, extending considerations pertinent to singlecomponent Bose gases [28,29,32,33] to the two-component case, we may use the following set of two coupled 1D SGPEs to describe the axial modes of the system: Here, ψ j (x, t) (j = 1, 2) are complex order parameters characterizing each component of the binary Bose gas, m is the atomic mass, µ j are the chemical potentials, while V (x) = (1/2)mω 2 x x 2 is the external trapping potential. Furthermore, η j (z, t) are complex Gaussian noise terms with correlations of the form η , where brackets denote averaging over different realizations of the noise. The strength of the latter can be calculated ab initio by the Keldysh selfenergy [32]; for thermal clouds close to equilibrium, the relevant integrals determining the dissipation γ j (x, t) can be expressed as follows: n are the energies of the n-th excited states, N n and k n = 2mE n / 2 denote, respectively, the kinetic energies and momenta of single particles in the n-th excited state. Physically speaking, Eq. (2) describes the exchange of atoms between the thermal clouds and the condensates due to elastic collisions; notice that in the above description we have taken into regard exchanges up to the third excited state while, to leading order approximation, we have omitted exchanges between the different hyperfine states (in other words, we have considered the simplest situation where each condensate component interacts with its own thermal cloud).
Under the above assumptions, the dissipation terms γ j (x) may in principle be calculated numerically, for several temperatures, as was done in the case of a singlecomponent Bose gas in Refs. [17]. In this work, it was shown that, sufficiently close to the trap center (i.e., in the interval [−R/2, R/2], where R is the Thomas-Fermi radius), the dissipation takes approximately constant values for a relatively wide range of temperatures. Furthermore, as shown in Refs. [15,16,19] (see also the discussion in Ref. [17,19] and, more recently, in Ref. [20]), the value of γ-which determines the dark soliton's life time-scales with temperature as γ ∝ T α , with 1 < α < 4; note that the case γ ∝ T 4 corresponds to the regime k B T ≪ µ, while the case γ ∝ T corresponds to the regime k B T ≫ µ (where µ is the chemical potential of the background Bose liquid).
Taking into regard the above findings, below we will consider the situation where both dissipative terms γ j are constant: such an assumption is consistent with our scope, i.e., to analyze the dynamics of the DB-soliton near the center of the trap. Furthermore, based on the fact that simulations investigating soliton dynamics in the framework of the SGPE model were found to be in fairly good agreement with analytical and numerical results relying on the respective DGPE model, below we will omit the noise terms η j (x, t); this way, we will use the following system of two coupled DGPEs to describe the DB soliton dynamics in the two-component Bose gas at finite temperatures: Note that the above model was recently used in Ref. [34], where the quantum Kelvin-Helmholtz instability of a two-component BEC was studied.
The system of Eqs. (3) can be expressed in dimensionless form as follows. Measuring the densities |ψ j | 2 , length, time and energy in units of 2a, a ⊥ = /ω ⊥ , ω −1 ⊥ and ω ⊥ , respectively, Eqs. (3) become: where we have used the notation ψ 1 = u d and ψ 2 = u b , indicating that the component 1 (2) is supposed to support a dark (bright) soliton, and the respective chemical potentials are now µ 1 = µ d = µ and µ 2 = µ b = µ + ∆; in our considerations below we assume that where Ω = ω x /ω ⊥ ≪ 1 is the normalized trap strength; the latter, along with the thermally induced damping parameters γ d,b , are considered to be small parameters of the system (these will be treated as formal perturbation parameters in our analytical approximation -see below). We should add a comment here about the relevant range of values of the parameter γ. A number of recent experiments, including the ones in Hamburg [9,25], Heidelberg [24,26] and Pullman [10,11], have focused on regimes of very low temperature where the effect of the term associated with γ is imperceptible (over the experimentally relevant time scales). The focus of these experiments was on the soliton dynamics and an effort was made (by operating at T /T c ≤ 0.1) to correspondingly minimize the thermal effects. It is easier to appreciate the latter features in the context of the earlier experiments of the Hannover group [21,23], which were conducted in the regime of T /T c ≈ 0.5. In that realm, the relevant values of γ can be estimated to be up to 10 −2 [35]. In what follows, we will treat γ generally as a free parameter, in order to illustrate the available wealth of bifurcation and dynamical phenomena of this system. Nevertheless, the reader more keen on the physical applications of the model to the physics of finite-temperature BECs should keep in mind the above values as a guideline towards the parameter regimes pertinent therein. We finally note that our analysis may also be used as a theoretical basis for understanding results of future experiments on dark and dark-bright solitons exploring finite-temperature effects (see, e.g., discussion in the Supplemental Material of Ref. [10]).

Relaxation to the ground state of the system
Since our purpose is to study the dissipative dynamics of DB solitons in this setting, it is natural to consider at first the dynamics of the pertinent background wave functions, namely a Thomas-Fermi (TF) wave function for the u d component and a zero wave function for the u b component. In particular, we will show that the coupled DGPEs Eqs. (4)-(5), similarly to their one-component counterpart (see, e.g., discussion in Ref. [36]), describe a relaxation process. Namely, as a result of the finite temperature, the two components, starting (at t = 0) from suitable initial conditions, will evolve so that, at sufficiently large times, u d will converge towards a TF cloud with the prescribed value of the chemical potential µ, while u b will vanish.
To show that this is the case indeed, we examine the peak amplitudes U d,b (t) of the wave functions u d,b (x = 0, t), corresponding to their (absolute) values at the center of the trap (i.e., at x = 0, where V (x) = 0 as well), and assume respective phases θ d,b (t). The evolution equations for U d,b (t) and θ d,b (t), which can directly be obtained by introducing the ansatz into Eqs. (4)- (5), are of the form: where overdots denote time derivatives. Next, utilizing Eqs. (6), we obtain from Eqs. (7)-(8) the following system: . It is clear that that the system of Eqs. (9)-(10) has a fixed point (U d0 , U b0 ) = ( √ µ, 0) [a similar analysis can be done for the fixed point The evolution of small perturbations U d1,b1 around this fixed point can then readily be found introducing the ansatz and U b0 = U b1 (t) into Eqs. (9)-(10) and linearizing with respect to U d1,b1 ; this way, we can easily solve the equations for U d1,b1 and finally obtain the following approximate expressions for the peak amplitudes of the wave functions: where U d,b (0) are initial conditions. Thus, at sufficiently large times, the peak amplitude of u d will decay to the value √ µ, while the one of u b will become zero. Accordingly, during the relaxation to equilibrium process, one may expect the following type of evolution towards relaxation. If the u d component is initially a Thomas-Fermi (TF) cloud of amplitude U d (0), its density will evolve as, On the other hand, if the u b component has initially the form of an arbitrary localized function, e.g., a Gaussian, of amplitude U b (0), it will asymptotically approach the trivial stationary state. The above predictions can be directly compared to numerical simulations. In particular, in Fig. 1 we show the evolution of a state characterized by the initial densities

Analytical results
Having studied the relaxation process described by Eqs. (4)-(5), we will now proceed to investigate, in the same framework, the dissipative dynamics of DB solitons. We  will assume that the dark soliton is on top of an already formed TF cloud with the equilibrium density |u d,T F | 2 = µ − V (x); this way, the density |u d | 2 in Eqs.
, and cast Eqs. (4)-(5) into the following form: (14)- (15) can be viewed as a system of two coupled perturbed nonlinear Schrodinger (NLS) equations, with perturbations given by Eqs. (16)- (17). In the absence of the perturbations, i.e., at zero temperature (γ b = γ d = 0) and for the homogeneous system (V (x) = 0) subject to the boundary conditions |u d | 2 → 1 and |u b | 2 → 0 as |x| → ∞, the NLS Eqs. (14)-(15) possess an exact analytical one-DB-soliton solution of the following form: where φ is the dark soliton's phase angle, cos φ and η represent the amplitudes of the dark and bright solitons, and D and x 0 (t) are associated with the inverse width and the center position of the DB soliton. Furthermore, k = Dtanφ = const and θ(t) are the wavenumber and phase of the bright soliton, respectively. The above parameters of the DB-soliton are connected through the following equations: withẋ 0 denoting the DB soliton velocity. Notice that the amplitude η of the bright soliton, the chemical potential µ of the dark soliton, as well as the (inverse) width parameter D of the DB soliton are connected to the number of atoms of the bright soliton by means of the following equation: Let us now employ the Hamiltonian approach of the perturbation theory for the matter-wave solitons to study the dissipative dynamics of DB solitons. We start by considering the Hamiltonian (total energy) of the system of Eqs. (14)- (15), in the absence of the perturbations (i.e., for R b = R d = 0), namely, The energy of the system, when calculated for the DB soliton solution of Eqs. (18)- (19), takes the following form: We now consider an adiabatic evolution of the DB soliton and, particularly, we assume that, in the presence of the perturbations of Eqs.
where we have used Eq. (23). The evolution of the parameters φ(t), D(t) and x 0 (t) can be found by means of the evolution of the DB soliton energy. In particular, employing Eq. (25), it is readily found that On the other hand, using Eqs. (14)- (15) and their complex conjugates, it can be found that the evolution of the DB soliton energy, due to the presence of the perturbations, is given by: where the asterisk denotes complex conjugation. Substituting R d and R b into Eq. (29) and evaluating the integrals, we finally obtain from Eqs. (28)-(29) the following result: Equation (30), together with Eqs. (26)-(27), constitute a system of equations for the unknown soliton parameters φ(t), D(t) and x 0 (t). In the case of a DB soliton near the center of the trap with an almost "black" dark-soliton-component (i.e., x 0 ≈ 0 and cos φ ≈ 1), the above system has a fixed point x 0,eq = 0 and φ eq = 0, and Considering now small perturbations around the fixed points, i.e., x 0 → 0+x 0 , φ → 0+φ and D → D eq + D 1 , we linearize Eqs. (26)- (27) and Eq. (30) with respect to x 0 , φ and D 1 , and obtain the following results: x 0 = D eq φ.
Differentiating Eq. (34) with respect to time once, and using Eqs. (33)-(34), we obtain after some straightforward algebraic manipulations the following equation of the motion for the DB soliton center x 0 : where the oscillation frequency ω osc and the anti-damping parameter a are respectively given by: We note that in the absence of dissipation [a = 0 in Eq. It is clear that the nature of the soliton trajectories x 0 (t) as predicted by Eq. (35) depend on whether the roots of the auxiliary equation s 2 − as + ω 2 osc = 0 are real or complex. The roots are given by with the discriminant D ≡ a 2 −a 2 cr determining the type of the motion. In particular, we identify different temperature-dependent damping regimes: the subcritical weak antidamping regime (D < 0, a < a cr ), the critical regime (D = 0, a = a cr ), and the super-critical strong antidamping regime (D > 0, a > a cr ). In the first regime the soliton performs oscillations of growing amplitude, with x 0 (t) ∝ exp(at) cos(ω osc t), while in the latter two regimes the soliton follows an exponentially growing trajectory, i.e., x 0 (t) ∝ exp(s 1,2 t) (with s 1,2 ∈ R), and decays at the rims of the condensate cloud (see also below).

Numerical results
We now turn to a numerical examination of the above findings. First, we will show that our analytical predictions are supported by a linear stability analysis around the stationary DB soliton, say u 0 ≡ (u, v) T [see Eqs. (18)- (19) for φ = 0 and x 0 = 0]. For such a state, the right-hand side of Eqs. (4)-(5) still vanishes and, thus, stationary DB solitons are exact solutions of the problem with γ d,b = 0. We obtain this solution by means of a fixed point algorithm and then find the linearization spectrum around the stationary DB soliton state as follows. We introduce the ansatz into the DGPEs Eqs. (4)-(5) (here {λ, (a, b)} define an eigenvalue-eigenvector pair, and ǫ is a formal small parameter), and then solve the ensuing BdG eigenvalue problem. In Fig. 2, we observe a prototypical realization of a stationary DB soliton in a trap of strength Ω = 0.1 (for simplicity, we consider the case with γ d = γ b = γ). Notice that upon the variations of γ (and hence of temperature) considered in the figure, the solution profile does not change, as mentioned above; however, the linearization problem and its eigenvalues significantly depend on the value of γ, as is shown in the four bottom panels of Fig. 2. In the zero-temperature (Hamiltonian) case of γ = 0, all eigenvalues are imaginary. Furthermore, the oscillatory motion of a single DB soliton in the trap [12] (see also recent work in Refs. [10,11,31]) is spectrally associated with the existence of a single anomalous (alias negative Krein sign, "translational") mode in the linearization around the stationary soliton. In analogy to the case of dark solitons (see e.g. Refs. [15,26,30]), this anomalous mode possesses a frequency identical to the frequency of the DB soliton oscillation, i.e., ω AM ≡ Im(λ AM ) = ω osc .
Our analytical approximation for ω osc is tested against the numerical results for Im(λ AM ), both in the case examples of Fig. 2, as well as in the parametric dependence results of Fig. 3. It is clear from the spectral plots (middle and bottom panels of Fig. 2) that, as soon as γ = 0, the relevant anomalous eigenmode (indicated by red stars in Fig. 3) becomes complex, leading to soliton oscillations of growing amplitude; this behavior, which corresponds to the "subcritical" regime mentioned above, is similar to the case of dark solitons [17,30] and in accordance with rigorous results pertaining to dissipative NLS systems [37]. If γ is increased beyond a critical point, namely γ cr ≈ 0.141, the relevant eigenvalue pair collides with the real axis, leading to the emergence of a pair of real eigenvalues (cf. bottom right panel of Fig. 2 and Fig. 3). This corresponds to the "super-critical" regime mentioned above (see also Refs. [17,30]), where the divergence of the soliton from its center equilibrium is purely exponential. Notice that the analytical predictions for the relevant unstable eigenvalue (and the oscillatory or purely exponential divergence from the equilibrium position) in Figs. 2 and 3 are generally fairly accurate, although their accuracy is decreasing as γ gets larger; this can be understood by the fact that our analytical approximation relies on the smallness of γ which was treated as a small parameter of the problem within our perturbation theory approach.
Last but not least, the role of the bright-soliton component in the dynamics should be highlighted in connection to the case of a dark soliton in a single-component condensate (where the bright soliton is absent). It can be directly seen from Eq. (37) that the anti-damping effect is always weaker for the DB soliton in comparison to the dark one (at least in the case γ d = γ b = γ we consider herein). Hence, the lifetime of the DB soliton is always longer than that of the dark soliton and, in fact, it becomes larger,  as the bright-soliton component "filling" of the dark one becomes stronger. This partial stabilization of the dark soliton evolution by means of its symbiotic second component is clearly illustrated in Fig. 4, where the bifurcation diagrams for the DB soliton are directly compared to the ones corresponding to the "bare" dark soliton. It is clear that the whole bifurcation diagram for the DB soliton is "drifted" towards larger values of γ (e.g., γ cr = 0.212 for the DB soliton and γ cr = 0.155 for the dark soliton), acquiring also smaller values of the instability growth rate λ r as compared to the ones found for the dark soliton for the same values of the temperature parameter γ. This is a clear indication that the "filled" dark soliton in a two-component BEC is more robust in the presence of finite temperature than a "bare" dark soliton in a single-component BEC. This is one of the principal findings of the present work.
Our analytical predictions were also tested against direct numerical simulations illustrating the evolution of the single DB soliton, both for the sub-critical case of oscillatory growth (see Fig. 5), and for the supercritical case of purely exponential growth (see Fig. 6). In both cases it can be seen that the dashed line corresponding to the analytical solution of the ODE (35) accurately tracks the evolution of the center of the DB soliton, which progressively loses its contrast and eventually disappears in the condensate background, with the system converging to its ground state (see section 2).
It is worth noting that in the results of Figs. 5 and 6 we have used, as initial condition, a TF cloud with a density at the trap center equal to the chemical potential µ appearing in the DGPEs (4)-(5). Nevertheless, we have also briefly studied a case where the density at x = 0 of the TF cloud was different from µ. The evolution in such a far from equilibrium scenario, is shown in Fig. 7, where the parameters are as in the subcritical case of Fig. 5, but with U d (0) = 0.8. It is readily observed that apart from the transient period towards equilibrium (i.e., when the density is rearranged so that it properly corresponds to the relevant value of the chemical potential µ = 1.5), the agreement between analytical and numerical results is fairly good. I.e., the fast scale of the background relaxation does not substantially affect the evolution of the DB wave on the slower time scale of the oscillatory decay of the latter.

Two Dark-Bright Soliton States
We now focus on the study of DB soliton "molecules" composed by two-DB-soliton states. Let us first consider the homogeneous case (Ω = 0), and use the following ansatz to describe a two-DB-soliton state composed by a pair of two equal-amplitude, oppositely located (at x = ±x 0 ) single DB solitons: where X ± = D (x ± x 0 (t)), 2x 0 is the relative distance between the two solitons, and ∆θ is the relative phase between the two bright solitons (assumed to be constant); below we will consider both the out-of-phase case, with ∆θ = π, as well as the in-phase case, corresponding to ∆θ = 0. Note that, similarly to the case of a single-DB soliton, the number of atoms N b of the bright-soliton component in the above two-DB-soliton state can be used to connect the DB-soliton parameters; in particular, if the two DB solitons are well-separated then N b is approximately twice as large compared to the result of Eq. (23), namely, N b ≈ 4η 2 √ µ/D. As was recently shown in Ref. [31], at zero temperature (i.e., γ d = γ b = 0), the evolution equation for the DB soliton center (for Ω = 0) reads: In the above equations, F int is the interaction force between the two DB solitons, which consists of three different components: the interaction forces F DD and F BB between the two dark and two bright components, respectively, as well as the interaction force F DB of the dark soliton of the one soliton pair with the bright soliton of the other pair (and vice-versa). These forces depend on the soliton coordinate x 0 , as well as on the DB soliton parameters, as follows [31]: where we have assumed thatḊ(t) ≈ 0 and, thus, D(t) ≈ D eq .
In order to complete the consideration of the case at hand, we will finally study the finite-temperature effect on a two DB-soliton state in the trap. To do so, we will combine the thermal effect on each DB soliton in the trap, represented by Eq. (35), and interaction effects included in Eq. (42). This way, we may use the following approximation to describe the motion of the centers of the two DB solitons: The equilibrium points x eq , can easily be found as solutions of the transcendental equation resulting from Eq. (48) lettingẋ 0 =ẍ 0 = 0 in both the in-and out-of-phase cases. To study the stability of these equilibrium points in the framework of Eq. (48), we use the ansatz x 0 (t) = x eq + δ(t), and obtain a linear equation for the small-amplitude perturbation δ(t), namely:δ − aδ + ω 2 1 δ = 0, where the frequency ω 1 is given by, where ω 2 osc and a are respectively given by Eq. (36) and Eq. (37). We now test the relevant predictions against BdG simulations, first for the inphase case in Figs. 8-9 and then for the out-of-phase case in Figs. 10-11. As expected, in the case of the in-phase configuration the BdG analysis reveals the existence of two anomalous modes: the one with the smaller (larger) eigenvalue-in the zero-temperature case-corresponds to an in-(out-of-) phase motion of the two DB solitons, similarly to the case of a two-dark-soliton state in a single-component BEC [26,30]. These two anomalous mode pairs lead to complex eigenfrequencies for γ = 0, and the two-DBsoliton state performs oscillations of growing amplitude. Similarly to the case of the single DB soliton, as γ is increased these pairs collide pairwise on the real axis in two critical points, namely γ 1 ≈ 0.153 and γ 2 ≈ 0.38. Beyond the second critical point γ 2 , the growth of the trajectory of the DB soliton center becomes purely exponential. The theoretical approximation of the relevant complex (and subsequently real) eigenvalues depicted by dashed line in Fig. 9 is again fairly accurate, becoming progressively worse as γ increases.
A similar phenomenology arises in the case of out-of-phase two-DB-soliton states, as shown in Figs. 10-11. However, there exists a rather nontrivial twist in comparison to the previous case. In particular, a third pair of complex eigenvalues emerges due to the fact that a third anomalous mode exists for γ = 0. This mode is no longer a translational one associated with the in-or out-of-phase motion of the two soliton centers (as before and as shown in the bottom left and bottom right eigenmodes of Fig. 10). It is instead a mode associated with the π relative phase of the peaks: if we add the eigenvector of this unstable (for γ = 0) mode to the two-DB-soliton out-of-phase solution, we observe that while the center location of the state remains intact, the relative heights of the two solitons are affected, leading to a symmetry breaking of the configuration. We will not consider this unstable mode further since its induced instability is weaker than those of the (in-phase and out-of-phase) translations. Nevertheless, we note that all three pairs of modes eventually collide on the real axis, eventually leading to pairs of purely real eigenvalues. Finally, we turn to direct numerical simulations for both the in-phase two-DBsoliton state in Fig. 12, and for the out-of-phase two-DB state in Fig. 13. In both cases, we show only the low-γ, oscillatory growth (subcritical) regime. Despite the complexity of the resulting system and of the DB soliton interactions, it can still be clearly observed that the ODE (48) can be used to capture fairly accurately the relevant dynamics even for the long time evolutions considered in these figures. Here, it should be mentioned that the temperature-induced dissipation results in an interesting effect. Particularly, as observed in Figs. 12 and 13, for short times, the individual DB solitons clearly behave like repelling particles, which can always be characterized by two individual density minima -even at the collision point. Nevertheless, for longer times, the nature of their interaction changes: due to dissipation, they gain kinetic energy and completely overlap at the collision point.

Conclusions
In the present work, we presented a systematic analysis of a prototypical model (the so-called dissipative Gross-Pitaevskii equation) incorporating the effects of temperature on the dynamics of dark-bright (DB) solitons. This was done both in the case of a single DB soliton, as well as in the case of DB soliton "molecules", composed by multiple (inor out-of-phase) DB solitons. We have developed a perturbation theory for the two-component system to analytically show the following: similarly to dark solitons, dark-bright ones execute anti-damped oscillations of growing amplitude for sufficiently low temperatures, while if the relevant parameter becomes sufficiently large, then the decay of the contrast of the solitons (and their disappearance in the background) becomes exponential.
A fundamental effect revealed by our analysis is that the presence of the bright ("filling") component hinders the temperature-induced dissipation associated with the dark soliton, and offers a significant partial stabilization (i.e., a significantly longer lifetime) to the corresponding symbiotic DB soliton structure, in comparison to its "bare" dark soliton counterpart. The above effect relies on the fact that the critical value of the relevant parameter (labeling the different damping regimes) is increased, while the instability growth rate is decreased, for the DB solitons. Similar conclusions were reached in the case of two dark-bright entities, with the added twist that their relative phase may introduce (in the out-of-phase case) additional anomalous modes and instability sources in the system. The latter are not associated with in-or out-of-phase translational motion of the solitons but rather with a symmetry-breaking in their relative amplitudes.
As concerns the relevance of our findings with pertinent experimental efforts we note the following. First, all relevant recent experiments for dark and darkbright solitons were conducted at extremely low temperatures, aiming to minimize corresponding anti-damping effects. In our setting this corresponds to subcritical dynamics of small γ. Nevertheless, we believe that our findings may be relevant to future experiments exploring in more detail finite-temperature effects (see Supplemental Material of Ref. [10]).
A natural direction to extend the present studies is to consider the higher dimensional setting of vortices [38] and of their two-component generalizations, namely the vortex-bright solitons [39]. Understanding the thermally induced dynamics and the modifications of the corresponding precessional motion, especially in the presence of multiple coherent structures would constitute an interesting topic for future study. On the other hand, it would certainly be relevant to extend the present studies to more complex models that provide coupled dynamical equations for the condensate and the thermal cloud [28] (rather than use a single equation directly incorporating the effects of the thermal cloud on the condensate without the possibility of "feedback"). Such studies are in progress and pertinent results will be reported in future publications.