Interplay of Solitons and Radiation in One-Dimensional Bose Gases

Interplay of Solitons and Radiation in One-Dimensional Bose Gases

Yuan Miao    Enej Ilievski    Oleksandr Gamayun Institute for Theoretical Physics, Institute of Physics and Delta Institute for Theoretical Physics, Universiteit van Amsterdam, Science Park 904, 1098XH Amsterdam, the Netherlands
July 19, 2019

We study relaxation dynamics in one-dimensional Bose gases, formulated as an initial value problem for the classical non-linear Schrödinger equation. We propose an analytic technique which takes into account the exact spectrum of non-linear modes, that is both soliton excitations and dispersive continuum of radiation modes. Our method relies on the exact large-time asymptotics and uses the so-called dressing transformation to account for the solitons. The obtained results are quantitatively compared with the predictions of the linearized approach in the framework of the Bogoliubov–de Gennes theory. In the attractive regime, the interplay between solitons and radiation yields a damped oscillatory motion of the profile which resembles breathing. For the repulsive interaction, the solitons are confined in the sound cone region separated from the supersonic radiation.

I Introduction

Recent years have brought tremendous experimental and theoretical progresses in understanding various aspects of equilibrium and non-equilibrium physics in strongly correlated many-body systems, especially in the domain of quantum gases Bloch et al. (2008); Jaksch and Zoller (2005); Kinoshita et al. (2006); Hofferberth et al. (2007); Langen et al. (2015). There has been a particularly intense focus on studying relaxation phenomena and microscopic mechanisms responsible for thermalization in isolated many-body systems Deutsch (1991); Srednicki (1994); Rigol et al. (2008, 2007); Caux and Essler (2013); Eisert et al. (2015); Vidmar and Rigol (2016); Gogolin and Eisert (2016); DAlessio et al. (2016); Caux (2016); Cugliandolo et al. (2018); Deutsch (2018). Equilibration in generic chaotic (i.e. ergodic) systems is nowadays quite well understood, primarily using the arguments of Eigenstate Thermalization Hypothesis Deutsch (1991); Srednicki (1994); Rigol et al. (2008); Alba (2015); DAlessio et al. (2016); Deutsch (2018). On the other hand, it has been argued that non-ergodic systems, e.g. models which lie in the proximity of an integrable point, fail to thermalize in the conventional sense owing to an extensive amount of local conservation law which severely constraints the dynamics. The proposed generalized Gibbs ensembles Rigol et al. (2007); Vidmar and Rigol (2016); Essler and Fagotti (2016) have subsequently been scrutinized in a variety of non-interacting and interacting exactly solvable quantum many-body dynamics Caux and Konik (2012); Ilievski et al. (2015); Essler et al. (2015); Ilievski et al. (2016, 2017), and in the context of classical integrable systems by taking the classical limit of quantum fields De Luca and Mussardo (2016).

Theoretical studies of non-equilibrium phenomena in exactly solvable models are for the most part concerned with steady states and their properties, while much less is known about the relaxation dynamics at large (or intermediate) time-scales. Despite integrability, the latter represents a formidable task in both classical and quantum many-body systems with interacting degrees of freedom.

For instance, while the theoretical framework for solving integrable differential equations describing classical field theories is very well developed, a full-fledged analytic treatment of non-linear wave equations is unfortunately not tractable in full generality, explaining the scarcity of closed-form results in the literature. This work aims to partially fill this gap by presenting some non-trivial analytic results in a physically relevant setting. In particular, we consider interacting one-dimensional (1D) Bose gases which are, in the weakly-coupled regime, well described by the classical non-linear Schrödinger equation (NLSE) (also known as the Gross–Pitaevskii equation), one of the prime examples of exactly solvable non-linear wave equations Pitaevskii and Stringari (2016); Pethick and Smith (2001); Faddeev and Takhtajan (1987); Novikov et al. (1984a).

We study the initial value problem for the NLSE by implementing a ‘classical quench protocol’, i.e. initializing an inhomogeneous profile and letting it evolve under the non-linear evolution law. By exploiting integrability of the equation of motion, we devise an analytic technique by combining the (Darboux) dressing transformation and exact asymptotic formulae, which allows us to accurately approximate the evolution even on moderately short time scales. Our results are benchmarked against the standard linear approximation and numerical simulations.

We consider field configurations which decay towards a constant vacuum density at large distance. The setup we study therefore differs from the conventional setting in quantum quenches which address local equilibration in thermodynamic quantum gases at finite density. Nonetheless, it turns out that an extensive amount of conservation laws once again play a pivotal role in constraining the relaxation process.

Integrable classical field theories in general feature two types of solutions with a distinct character: (i) non-linear interacting particles known as solitons, representing non-dispersive localized field configurations, and (ii) a continuum of non-linear dispersive modes called radiation. It is common practice however to treat small fluctuations of a uniform background density within the linear theory of non-interacting (Bogoliubov) quasi-particles. Contrary to the solitons and radiation modes the latter are not the proper elementary modes of the NLSE, and it is thus natural to wonder whether they remain a meaningful concept in genuine far-from-equilibrium scenarios such as classical quenches considered in this paper. In an attempt to answer this question, we investigate the linearized dynamics at a qualitative and quantitative level, focusing on intermediate time-scales where the effects of non-linearity cannot be neglected.

Aside from the theoretical interest, it is worthwhile to briefly mention some experimental aspects of the quench protocol proposed here. For nearly two decades, the experimental realizations of Bose-Einstein condensation (BEC) Anderson et al. (1995); Davis et al. (1995); Bradley et al. (1995) offer new routes for investigating the many-body phenomena in a controllable and precise manner. In particular, 1D Bose gases with tunable s-wave scattering interaction (induced by Feshbach resonance) Burger et al. (1999); Muryshev et al. (2002); Khaykovich et al. (2002); Lepoutre et al. (2016); Schemmer et al. (2018) are of great interest, due to the simple and exactly solvable theoretical models which displaying rich physical behaviour. Theoretically speaking, 1D Bose gases with s-wave scattering can be modelled as the Lieb-Liniger model Lieb and Liniger (1963); Lieb (1963), which is, in the weak-coupling (Gross–Pitaevskii) limit, nothing but a 1D NLSE Pethick and Smith (2001); Pitaevskii and Stringari (2016); Carr et al. (2000a, b). Above all, solitonic excitations in both 1D attractive and repulsive NLSE have been observed experimentally Burger et al. (1999); Khaykovich et al. (2002); Becker et al. (2008); Weller et al. (2008); Stellmer et al. (2008); Lepoutre et al. (2016), making it promising to realize our quench proposal in ultracold atom experiments.

The paper is organized as follows. In Sec. II we consider quenches in the attractive 1D NLSE. In Sec. (II.1) we derive the long-time asymptotic solutions for a solitonless quench, while in Sec. (II.2) we describe quenches when both the soliton and radiation modes are present, and discuss the “soliton breathing” effect.

The quench in the repulsive NLSE with finite density is studied in Sec. III, revealing different physical phenomena such as e.g. the separation of the soliton sound cone from the supersonic radiation modes. We conclude in Sec. IV by summarizing the main results and listing some open questions. A concise introduction to the classical inverse scattering method and a detailed derivation of the complete spectrum of the linearized equations are presented in the appendices.

Ii Attractive Interaction

We consider an attractive (focusing) NLSE, which in dimensionless units takes the form


with non-linearity coefficient .

This equation can be solved by the inverse scattering method, which we briefly describe below using Refs. Faddeev and Takhtajan (1987); Novikov et al. (1984b). The essential idea is to reformulate the original non-linear equation as an auxiliary linear problem for the ‘wave-function’ ,




which is interpreted as a scattering problem: the field is viewed as the scattering potential, while the spectral parameter plays the role of energy. The transfer matrix for this problem can be presented as


with . One can compute the time dependence of these scattering data assuming that satisfies NLSE (1), obtaining a remarkably simple dependence Faddeev and Takhtajan (1987):


This way, to solve the initial value (Cauchy) problem of Eq. (1), one has to compute and diagonalize the transfer matrix for a given initial profile, evolve the transfer matrix according to Eq. (5) and, finally, retrieve the time-evolved potential from the scattering data. The last step is referred to as the “Inverse Scattering” and requires to solve the Gelfand–Levitan–Marchenko linear integral equation. In the very special cases of the reflectionless potentials, characterized by , this equation can be solved analytically and its solutions are called solitons. In particular, the one-soliton solution reads


describing a bell-shaped profile travelling with the velocity , usually named as a bright soliton. Without loss of generality, one can set the initial positions , , the inverse width as well as the velocity by an appropriate Galilean transformation. For generic initial conditions, i.e. , in addition to solitons the spectrum also contains radiation modes, which unlike solitons are dispersive.

In order to model various types of initial conditions while still being able to perform analytic computations, we consider a rescaled Satsuma-Yajima profile Satsuma and Yajima (1974); Miles (1981); Gamayun and Semenyakin (2016) with


The quench parameter can be removed from the profile by a suitable rescaling of Eq. (1) which will change the non-linearity coefficient (the coupling constant), . Hence, it is feasible to realize this profile in the ultracold atom experiments by preparing a one-soliton state and subsequently quenching parameters of the holding trap and external fields to induce the change in the coupling constant Gamayun et al. (2015); Gamayun and Semenyakin (2016); Caudrelier and Doyon (2016); Franchini et al. (2016). One can think of this as a classical analogue of the interaction (non-linearity coefficient) quench in the Lieb-Liniger model Kormos et al. (2013); De Nardis et al. (2014); Piroli et al. (2016).

The scattering data of potential (7) can be computed exactlyGamayun et al. (2015); Gamayun and Semenyakin (2016)


For integer quench parameters , corresponding to the reflectionless potential , the solution only involves solitons. The soliton parameters are identified with the zeroes of in the upper half plane, which has the form


corresponding to an -soliton solution. In this scenario, all solitons have zero velocity and the time-evolution exhibits a periodically oscillating behavior due to their mutual interaction which we refer to as “soliton breathing” Gamayun and Semenyakin (2016).

In the generic case, we put (the greatest integer smaller or equal to ), and the factor can be presented as




This representation indicates that the field profile involves solitons (described by parameters ), superimposed on a continuous background of radiation modes encoded by the non-zero reflection coefficient .

We focus our analysis on two representative regimes: (A) ; (B) . In the (A) case, solitons are absent and the evolution is governed solely by radiation modes, resulting in a ballistic widening of the initial profile similar to the wave-packet spreading in the linear Schrödinger equation. In the (B) case, the non-linearity becomes much more important since the profile contains a soliton. Therefore, in addition to the ballistic expansion there remains a “stable” part. The attractive interaction between the soliton and radiation will show up in the form of “breathing” of the whole profile, similarly to the two-soliton solution Gamayun and Semenyakin (2016). In the next section we examine these two regimes using an analytic expression for the time-asymptotic behavior and give a quantitative analysis of the soliton breathing phenomenon. Furthermore, we compare our analytic findings with numerical integration.

Figure 1: Time evolution of quenched profile (7) with and . The dashed lines denote the half width of the profile, indicating a ballistic (i.e. linear in time) expansion of the half width.

ii.1 Solitonless quench

A typical time evolution of a quenched profile (7) with is plotted in Fig. 1. One can observe a ballistic expansion, reminiscent of the wave-packet spreading governed by the linear Schrödinger equation.

Indeed, at the phenomenological level the situation is described as follows. At the initial moment of time a very narrow profile can be rescaled to the unit characteristic width by the rescaling and , resulting in the non-linearity coefficient proportional to and the dynamics dominated by the linear part of the equation. The evolution causes the spreading amplitude damping of the profile, thus suppressing the importance of the non-linear terms. The complete account of the non-linear terms needed for the analysis of the exact asymptotic expression requires sophisticated techniques based on the Riemann–Hilbert problem Zakharov and Manakov (1976); Deift et al. (1993). We shall instead proceed with a more physically transparent (but less rigorous) method, inspired by the above analogy with the linear case.

We begin by recalling the general solution of the linear Schrödinger equation and its large-time asymptotics,


Function is uniquely determined from the initial condition. This result motivates to look for the solution of Eq. (1) in the following form Segur and Ablowitz (1976); Segur (1976); Zakharov and Manakov (1976)


where , and are all functions of the scaling variable . Substituting this expression into Eq. (1), we find that all functions can be expressed in terms of the amplitude function and its derivative, and in particular


The aim now is to determine . In order to relate it to the initial profile, we borrow the logic of Ref. Segur and Ablowitz (1976), and compute the local conserved charges both on the initial profile and on the asymptotic expression (14).

According to Eq. (5), is a conserved quantity and the conventional integrals of motions are defined as coefficients in the asymptotic expansion, namely


Taking into account Eq. (11), the conserved charges can be found from the moment expansion


On the other hand, the local conserved charges can be presented as spatial integrals of local densities


which could be found from the recurrence relation Faddeev and Takhtajan (1987); Ablowitz and Segur (1981a)

Figure 2: Time evolved profile (7) for and at , obtained by numerical integration (solid line) and from the asymptotic solution given by Eq. (24) (dashed line). The inset plot shows the relative difference between the two integration methods .

Since asymptotically the non-linear part is suppressed in the inverse powers of , the conserved densities will be approximately the same as in the linear case. Evaluating them on the profile (14), one gets


After the change of variable, , the conserved charges read


Hence, taking into account Eq. (17), we readily identify


which in turn allows to completely determine the asymptotic solution from the extensive number of the conserved charges (i.e. the action variables for the 1D NLSE).

The determination of the phase of is more technical, and the result reads Zakharov and Manakov (1976)


The phase provides the angle variables for the 1D NLSE Its and Ustinov (1991). The full expression for the asymptotic solution reads


with the phase


Noticing that the density of the asymptotic profile is fully characterized by the modulus , i.e. local properties (charges) of individual initial configurations, our scenario closely resembles the problem of identifying a generalized Gibbs ensemble which corresponds to the reduced density matrix in the steady-state limit of a quantum quench. The information about the asymptotic phase (25) cannot be fully restored from the integrals of motion but requires extra information about the initial angle coordinates. Nevertheless, by discarding the information about the initial phase (e.g. by uniformly averaging over it), one can define a microcanonical ensemble of states whose long-time asymptotics is completely determined by the initial values of the conserved charges. In other words, the (microcanonical) ergodic average of local observables (e.g. observables proportional to field and its derivatives) only retains information about the initial action variables, and this can be extracted via the method presented here.

In Fig. 2 we compare cthe asymptotic expression with the results of numerical integration. It is worth stressing that the asymptotic solution (24) to Eq. (1) only becomes exact when all higher order terms in Eq. (14) are taken into account. Remarkably, we observe no dramatic effect of non-linearity coefficient on the half-width of the quenched profile, which means that the linearized solution offers reasonably good results for all times if the initial profile is small enough (of course, the initial profile is proportional to , but the half-width grows strictly ballistically as in the linear equation). Still, the non-linearity has to be accounted for in order to give the correct values of the (infinitely many) conserved charges (cf. dependence in Eq. (22)). The presented method for the analytic determination of the asymptotic solution can be easily generalized to the other classical quench protocols, as long as the class of initial conditions permits to extract the scattering data, which enables computation of the conserved charges. We have demonstrated that the asymptotic solution (24) accurately describes the situation of the solitonless initial profiles.

In the presence of solitons, the outlined method cannot be straightforwardly generalized despite the fact that the phase-space expressions for the local conserved charges remain exactly the same as in solitonless case (see Eq. (18)). However, from the asymptotic expansion of Eq. (16) using general presentation (11) one can clearly see that the conserved charges acquire soliton corrections. Therefore, instead of (14), a new ansatz for the asymptotic field profile is needed. The latter can be, in principle, deduced from the large-time asymptotic solution of the NLSE linearized on the soliton background. But even in this case, it is not clear how to effectively evaluate the densities and compute the corresponding charges. Similar approaches developed in Ref. Segur (1976) give unsatisfactory results. Thus, to describe the (B) regime, , we pursue a different strategy presented in the next section. Firstly, we compute the full time dependence of the NLSE linearized on the soliton background, which yields reasonably good results on short and intermediate times. Secondly, we employ the Darboux transformation, which permits us to satisfy condition (43) exactly, but requires explicit knowledge of the asymptotics for the radiative part.

ii.2 Interplay of solitons and radiation

The time evolution of the quenched profile undergoes a qualitative change in the regime where the non-linear effects become more pronounced. For , according to Eq. (11), the spectral data contains a single static soliton whose form can be explicitly computed from the scattering data Gamayun et al. (2015); Gamayun and Semenyakin (2016)


The results of numerical integration are shown in Fig. 3. The effect of the attractive interaction is that ballistic spreading of the radiation modes is now accompanied by a “breathing motion”, i.e. the peak of the soliton oscillates, which, in particular, can be seen in the time-dependence of the profile’s half-width, see Fig. 4. The period of oscillations can be approximated by the period of the soliton,


given that near the origin () the evolved profile is a sum of the soliton and slowly varying component consisting of the radiation modes.

Figure 3: Time evolution of the intensity of the quenched profile (7) with and , showing persistent oscillations referred to as the “soliton breathing”. The red lines denote the half width of the profile.

Bogoliubov–de Gennes theory.

In order to present a heuristic picture of this phenomenon, we consider the linearized theory with respect to the soliton solution (26). Namely, we split


and account only for the linear terms in the equation of motion for , yielding


This equation, presented in the matrix form, is nothing but the Bogoliubov–de Gennes (BdG) equation Bogoljubov (1958); Pitaevskii and Stringari (2016); Walczak and Anglin (2011)


where and . In order to get rid of the time dependence which is due to the “potential” in Eq. (30), we introduce


so the corresponding BdG equation in the rotating frame reads


where now and are independent of time. General solutions of BdG equation (32) can be presented in the form


where and stand for the continuous and discrete spectrum of linear modes, respectively. It is natural to further decompose them into the real and imaginary parts, namely for the real part we have


and similarly for the imaginary part


The explicit expressions for the complete spectrum of modes are given in Appendix C. Time dependence of the expansion coefficients of the continuous part takes a simple form


with the dispersion


The dependence of the discrete part is at most linear in time


The initial values , and are determined from the initial profile.

Note that the presence of a soliton (non-zero ) opens a gap in the continuous spectrum of Bogoliubov quasiparticles (39). Moreover, the dispersion (39) coincides with magnon excitations (linear spin waves) in the presence of a magnetic field (). This is can be explained by the gauge similarity of the focusing NLS to the isotopic Landau–Lifshitz ferromagnet Faddeev and Takhtajan (1987), and the fact that the field intensity, , corresponds to the gradient of spin field.

We wish to stress that the discrete modes are essential to satisfy the completeness of the solutions of BdG equation (30). They can be obtained, in particular, by the observation that if represents a solution of the 1D NLSE (1) depending on some parameter , then the parametric derivative


satisfies BdG equation (29) (see also Refs. Walczak and Anglin (2011); Dziarmaga (2004)). In our case, can be any parameter in the one-soliton solution (6), e.g. , , , etc. Therefore, the discrete modes for are proportional to which corresponds to the discrete zero-energy solutions (also known as zero-modes). Notice however that the discrete modes generated this way are not all linearly independent. In Appendix C we carefully check the completeness and orthogonality relations.

The initial profile for BdG equation (29) is obtained by subtracting the soliton (26) from the quenched profile (7),


The initial values of the expansion coefficients of both continuous and discrete modes are calculated numerically from this profile. For the discrete modes we find that , while , which means that the imaginary part of the evolved profile grows linearly with time. This, in particular, results in the divergence of the conserved charges of the original non-linear problem as , signalling a breakdown of the linearization approach and its inadequacy for describing non-equilibrium scenarios such as the quenches presented here.

In Fig. (5) we present a comparison with the exact numerics, demonstrating that in relatively short amount of time the discrepancy blows up. It is worth mentioning at this stage that a general analysis of non-linear equations containing soliton-like solutions suggests that the zero-mode contributions are absent in the asymptotics Buslaev and Perelman (1995). Hence, the results of the BdG theory can be trusted only on sufficiently short time-scales. Curiously enough, if we voluntarily discard the contributions of discrete modes, i.e. putting and retaining only the finite-frequency continuum, the approximation improves quite noticeably, see Fig. (4). In order to describe the large-time asymptotics, one has to be able to properly take into account the infinite set of the conserved charges and find their values on the initial profile. We achieve this in the next section using the dressing method.

Figure 4: Time dependence of the profile’s half-width obtained from the linearized (BdG) equation with parameter and , with (diamonds) and without (dashed line) the zero-mode components, compared to the exact numerical integration.

Darboux transformation.

Figure 5: Solution to BdG equation (30) (diamonds) versus the solution constructed from the asymptotics profile dressed by a soliton (24) (dashed), shown for () at times , , , and , respectively. The results are benchmarked against the numerical integration (solid line).

The conserved charges can be found using the asymptotic expansion (16) and the usual parametrization (11). Namely, they can be split into soliton and radiation contributions


where for one can easily deduce that


whereas the radiation part remains the same as in the solitonless case (17) (with the corresponding not depending on the soliton parameters). Hence, if we were only to include the constraints of , we would have arrived at the same solution as previously, cf. Eq. (24). To overcome this issue with only minor modifications, we propose to employ the Darboux (dressing) transformation.

The very idea of the dressing transformation Matveev and Salle (1991) is to start from some reference solution and construct new solutions to the equation of motion by constricting a suitable non-linear transformation. This will ensure that (43) is satisfied exactly Gu et al. (2005); Its et al. (1988). In distinction to the linearization procedure outlined previously, the present construct heavily relies on the underlying integrability, in particular on the existence of the auxiliary linear problem (2), which in the short hand notations can be presented as


with connection


and with and defined in Eq. (3). Applying a gauge transformation with the matrix ,


we find that satisfies the linear problem (45) with the dressed connection


We should in addition demand that is such that has the same structure as Eq. (46), but with the new field corresponding to another solution of the NLSE. The gauge can be sought as a polynomial in the spectral parameter , but for our purposes it is enough to consider a linear function


The equation for then reads


while the dressed potential (connection) takes the form


To solve (50), let denote a general solution of the linear problem (2) with and the potential of being . A formal expression is given by Eq. (114). The solution for the complex-conjugate value is given by . These two solutions can be combined in the matrix , satisfying


with . This way, the solution to Eq. (50) is given by


The dressed field (cf. Eqs. (51), and (3)) takes the form


The scattering data for the dressed potential can be computed from the asymptotics of , which is determined by the asymptotics of (the undressed scattering data), and by the asymptotics of , which for reads


Finally, the dressed scattering data is given by


The obtained expression is consistent with the general form (11), and moreover manifestly respects the form of Eq. (43) for the conserved charges. In other words, we have constructed a new solution to the equation of motion, which contains the original field (described by in the spectral space) dressed by a soliton (with the spectral parameter ).

The most renowned application of the Darboux transformation is the construction of a one-soliton profile by dressing the trivial vacuum solution . In this case , and one gets a one-soliton solution with the corresponding given by Eqs. (6),(8), respectively. For the purposes of our application, we should instead apply the dressing to the asymptotic solution (24). The condition (43) is guaranteed to be satisfied exactly. To get a symmetric profile we use a specific solution for (see Eq. (117) and appendix A).

The comparison with numerical integration is shown in Fig. (5). We see that the dressing method stays very close to the exact result for all times, similarly as in the solitonless case. The origin of a small discrepancy is due to the fact that we have not used the exact solitonless solution but rather the asymptotic one given by Eq. (117).

Iii Repulsive interaction

Now we consider the repulsive interaction, also allowing for soliton modes if the system is initialized at finite density. We consider the finite density repulsive (defocusing) NLSE,


with non-linearity coefficient , and denoting the asymptotic background density, as . With no loss of generality we subsequently put .

The corresponding auxiliary linear problem is now of the form Faddeev and Takhtajan (1987)




As the initial profile we consider the “rescaled dark soliton” with quench parameter


The corresponding scattering data can be computed directly from the Eq. (58), reading Gamayun et al. (2015); Gamayun and Semenyakin (2016)


where and . Parameter is the speed of sound, corresponding to the upper limit for the velocity of dark/gray solitons Pitaevskii and Stringari (2016); Faddeev and Takhtajan (1987). The existence of such an upper bound is contrastingly different from the case of bright solitons in the 1D attractive NLSE (1) where the the range of soliton velocities is not restricted. Another important difference with respect to the attractive case is that for repulsive interactions the soliton part is present for all . Recall that solitons correspond to zeroes of in the upper-half spectral plane, so there are solitons coexisting with the radiative part (which exactly vanishes for the integer ). In particular, for we can present


or, alternatively, using the uniformization spectral variable, defined by


in an equivalent canonical form


The soliton part now describes a dark soliton with spectral parameter . The corresponding one-soliton profile is of the form


Similarly to the attractive case (cf. Eq. (11)), the remaining analytic part of can be seen as the radiation modes, characterized by non-zero reflection coefficient .

A physical consequence of repulsive interactions is that the radiation effectively decouples from the solitons. This means that to obtain the asymptotic profile it will be sufficient to consider their contributions independently. We proceed by computing the asymptotics of the generic solitonless profile, following the same lines as in the attractive case (cf. Secs. II.1). Subsequently we shall analyze the Bogoliubov-de Gennes theory describing linear fluctuations on the soliton background. We show that the decoupling between the radiation and soliton components is quite profound and in fact happens exponentially fast. The main benefit of this observation is that in order to describe asymptotic solutions in the presence of both the solitons and radiation there is no longer any need of resorting to the Darboux transformation.

iii.1 Solitonless asymptotic

Similarly to the attractive case, we consider linear perturbations of the trivial vacuum, which we choose to be . We put to obtain from Eq. (57) a linear equation


which has a general solution of the form



Figure 6: The asymptotic solution of Eq. (76) at (dotted line), compared with the numerical solution of the exact equation of motion (57) (solid line), for parameter and .

with being arbitrary functions with the following reflection property . For a smooth initial profile, the asymptotic behavior of , can be calculated by the steepest descent method. This requires to find the critical points of the phase factor ,


One can immediately notice that there are no solutions for , meaning that is exponentially small in that region. This provides an explanation for the supersonic nature of the radiation modes, as discussed later in Sec. III.2. For , the critical points are contained only in the “ part” of Eq. (68). To describe them explicitly we introduce and , . Then the -dependence of the critical points can be found from the relation


Notice also that , and have the same form as for a suitable parametrization for the scattering data (64)


however, in this case is not a spectral parameter but a function of and . (cf. with and for the attractive case (II.1)).

The critical points correspond to


which means that


This way, the asymptotics of for reads




Similar expressions can be found for . Utilizing the same approach as in Eq. (14), we can write down the following ansatz for the asymptotic expansion




In analogy to the attractive case, we fix by evaluating the local conserved charges on the asymptotic profile (68) and matching it to the initial values obtained from the spectral presentation. The calculations are slightly more involved and thus are relegated to Appendix (B). The end result is however analogous to (22), reading


The logarithmic phase correction is obtained after substituting the asymptotic ansatz (80) into the 1D NLSE (57),


This method allows us to fix the “action part” of the asymptotic solution and write


In the same manner as in the attactive case (cf. (23)) the phase of , encoding the angle variables Its and Ustinov (1991), can be fixed from the corresponding Riemann–Hilbert problem Deift et al. (1993); Its and Ustinov (1991); Vartanian (2000).

The comparison of asymptotic results against the numerics is displayed in Fig. (6). There we consider the quench which involves a dark soliton, but primarily focus on the region with no solitons. Next we shall discuss the separation between radiation modes and solitons in more detail.

iii.2 Solitons versus radiation at finite density

To gain some insight into how the presence of a soliton influences the analysis of the previous section, we follow the logic of the attractive case and consider a linearized theory with respect to the soliton background. Putting , we obtain a linearized equation


describing the BdG equation of the following matrix form


Here , and the soliton solution is given by (66). Unlike in the vacuum case (68), the general solution of Eq. (87) now comprises both the continuous and discrete modes.

Figure 7: Time evolution of of the quenched profile (60) with coupling and quench parameter . The radiation modes propagate faster than any dark/gray solitons, with lower velocity bound . All the radiation modes are confined outside of the sound cone, marked with dashed lines.

Formally, it can be split as


where the continuous part is further decomposed into normal modes, which for convenience we present separately for real and imaginary parts,


The exact form of and , as well as their orthogonality relations, are presented in Appendix (D). The time dependence of the expansion coefficients is given by


with the dispersion relation


One can observe that this dispersion is exactly the same as for excitations on the solitonless background. For small momenta it describes linear sound excitations, indeed as . Notice that the group velocity of propagation always exceeds the sound velocity


meaning that the radiation modes are “supersonic”, while the solitons on the contrary have “subsonic” velocities . Therefore, irrespective of the initial profile, the overlap between radiation and solitons becomes exponentially small in a short amount of time. This can be clearly observed in the numerical simulations of the exact non-linear dynamics of 1D NLSE (57) shown in Fig. 7. We would like to emphasize that all these effects can be attributed to the finite density boundary condition which (irrespective of the sign of the coupling constant) changes the dispersion relation of the Bogoliubov quasi-particles from the magnonic (i.e. quadratic) relation (39) to (93), with linear dependence at long wave-lengths as for acoustic phonons.

In addition to the continuum part, we have the discrete modes , which still can be regarded as zero-modes even though there is no gap in the spectrum. Their general form reads


with real constants and . These modes are essential to prove the completeness relations discussed in Appendix D (cf. Eqs. (185), (186) and (191)).

Figure 8: Solution of BdG equation (87) (diamonds) compared to the exact numerical integration (solid line) for the 1D NLSE with finite asymptotic density (57), with parameters at . The imaginary part of the solution ( and modes) is responsible for the observed excess in the asymptotic density (as shown by the dotted line).

The initial values of and , as well as, and are determined from the initial condition, which for the quench (60) can be written as


This profile is antisymmetric which implies no contribution from the zero-modes, , meaning that the linear evolution is governed entirely by the quasi-particle continuum.

The linearized evolution is compared with the exact numerics in Fig. (8), where the agreement at sufficiently large distances is quite good. It is particularly interesting to notice that the imaginary part (i.e. the contributions of and ) results in a higher asymptotic density near the origin (see red curve in Fig. 8). This happens because the full density () is not a conserved quantity of the linearized BdG equation (87). However, if we again permit ourselves to discard this contribution, we get a much better approximation (see dashed black curve in Fig. 8).

For a generic profile with , the BdG approach is a reasonable approximation only for short times. As commented earlier, unlike in the attractive case now there is not need of performing a Darboux transformation to merge the solitons and radiation. Indeed, due the to the exponentially small overlap the conserved charges separate automatically