Spontaneous generation of dark-bright and dark-antidark solitons upon quenching a particle-imbalanced bosonic mixture

TL;DR

This paper investigates how quenching a particle-imbalanced Bose-Bose mixture leads to the spontaneous formation of dark-bright and dark-antidark solitons, revealing the mechanisms, control parameters, and robustness of these localized structures.

Contribution

It uncovers the dynamical generation mechanisms of dark-bright and dark-antidark solitons during interaction quenches in imbalanced Bose mixtures, including effects of particle imbalance and trap geometry.

Findings

More solitons form with larger particle imbalance.

Stronger postquench interactions increase soliton number.

Tighter traps suppress soliton formation.

Abstract

We unveil the dynamical formation of multiple localized structures in the form of dark-bright and dark-antidark solitary waves that emerge upon quenching a one-dimensional particle-imbalanced Bose-Bose mixture. Interspecies interaction quenches drive the system out-of-equilibrium while the so-called miscible/immiscible threshold is crossed in a two directional manner. Dark-bright entities are spontaneously generated for quenches towards the phase separated regime and dark-antidark states are formed in the reverse process. The distinct mechanisms of creation of the aforementioned states are discussed in detail and their controlled generation is showcased. In both processes, it is found that the number of solitary waves generated is larger for larger particle imbalances, a result that is enhanced for stronger postquench interspecies interactions. Additionally the confining geometry highly affects the production of both types of states with a decaying solitary wave formation occurring for tighter traps. Furthermore, in both of the aforementioned transitions, the breathing frequencies measured for the species differ significantly for highly imbalanced mixtures. Finally, the robustness of the dynamical formation of dark-bright and dark-antidark solitons is also demonstrated in quasi one-dimensional setups.

Figures

Tables

Table 1. Table 1: Fitted DB soliton characteristics referring to the different collisional events that take place during the spatio-temporal evolution of both species shown in Figs.  1 (c) and (d). Recall that A ​ cos ⁡ ϕ 𝐴 italic-ϕ A\cos\phi [B] refers to the amplitude of the dark [bright] solitons generated in the A 𝐴 A [ B 𝐵 B ] species. The selected time instants, i.e. t = 35 𝑡 35 t=35 , t = 106 𝑡 106 t=106 , t = 175 𝑡 175 t=175 , and t = 250 𝑡 250 t=250 , for which we have performed the fits in the density profiles are indicated in Figs.  1 (c) and (d) by vertical dashed lines, and the collision events during propagation are marked with coloured circles. The number of particles of each bright solitary wave, N b subscript 𝑁 𝑏 N_{b} , prior and after each collision point is depicted in boldface. Note also that in all cases the accuracy of the fitted values is of about 0.98 0.98 0.98 .

DB characteristics between collisions
		$t=35$	$t=106$	$t=175$	$t=250$
DB soliton 1	$A\cos \varphi$	$1.1683$	$1.1434$	$1.2908$	$1.2301$
	$B$	$0.9087$	$1.0299$	$0.7435$	$1.1371$
	$D$	$0.5684$	$0.5372$	$0.4484$	$0.4195$
	$x_{0_{B_{\mathrm{𝟏}}}}$	$27.3425$	$27.3687$	$30.8764$	$22.6316$
	$x_{0_{D_{\mathrm{𝟏}}}}$	$27.3715$	$27.3613$	$31.1916$	$22.4180$
	$N_b$	$2.9053$	3.9491	2.4656	6.1651
DB soliton 2	$A\cos \varphi$	$1.1023$	$1.2351$	$1.2136$	$-$
	$B$	$0.9788$	$1.0675$	$1.1102$	$-$
	$D$	$0.5604$	$0.7378$	$0.5444$	$-$
	$x_{0_{B_{\mathrm{𝟐}}}}$	$17.8409$	$20.4372$	$19.7581$	$-$
	$x_{0_{D_{\mathrm{𝟐}}}}$	$17.8281$	$20.5576$	$19.7446$	$-$
	$N_b$	$3.4191$	3.0893	4.5278	$-$

Table 2. Table 2: Controlled DB solitary wave formation for a miscible-to-immiscible transition with g A ​ B i = 0.2 subscript superscript 𝑔 𝑖 𝐴 𝐵 0.2 g^{i}_{AB}=0.2 , upon significantly enlarging the system size. In all cases N A = 10 4 subscript 𝑁 𝐴 superscript 10 4 N_{A}=10^{4} is fixed while the bold-faced quantities are the ones that are varied in each distinct simulation. Notice the significantly larger DB soliton count as the system size increases.

Miscible to Immiscible transition
	$g_{AB}^f$	$\omega$	$N_A$	$N_B$	$N_{DB}$
	$1.2$	$0.05$	${10}^4$	${\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$47$
	$1.2$	$0.05$	${10}^4$	$\mathrm{𝟓}\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$114$
	$1.2$	$0.05$	${10}^4$	$\mathrm{𝟖}\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$136$
	$1.5$	$0.05$	${10}^4$	$5\times {10}^3$	$160$
	$2.0$	$0.05$	${10}^4$	$5\times {10}^3$	$237$
	$1.2$	$0.10$	${10}^4$	$5\times {10}^3$	$96$
	$1.2$	$0.15$	${10}^4$	$5\times {10}^3$	$85$

Table 3. Table 3: Controlled DAD formation for an immiscible-to-miscible transition with g A ​ B i = 1.5 subscript superscript 𝑔 𝑖 𝐴 𝐵 1.5 g^{i}_{AB}=1.5 , upon significantly enlarging the system size. In all cases N A = 10 4 subscript 𝑁 𝐴 superscript 10 4 N_{A}=10^{4} is fixed while the bold-faced quantities are the ones that are varied in each distinct simulation. Notice the significantly larger DAD soliton generation as the system size increases.

Immiscible to Miscible transition
	$g_{AB}^f$	$\omega$	$N_A$	$N_B$	$N_{DAD}$
	$0.75$	$0.05$	${10}^4$	${\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$40$
	$0.75$	$0.05$	${10}^4$	$\mathrm{𝟓}\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$160$
	$0.75$	$0.05$	${10}^4$	$\mathrm{𝟖}\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}$	$234$
	$0.80$	$0.05$	${10}^4$	$5\times {10}^3$	$142$
	$0.40$	$0.05$	${10}^4$	$5\times {10}^3$	$334$
	$0.75$	$0.10$	${10}^4$	$5\times {10}^3$	$130$
	$0.75$	$0.15$	${10}^4$	$5\times {10}^3$	$120$

Equations

$$i\hslash {\partial }_{\tilde{t}}{\tilde{\Psi }}_{\sigma }=\left(-\frac{{\hslash }^2}{2m}{\partial }_{\tilde{x}}^2+\tilde{V}(\tilde{x})-{\mu }_{\sigma }+\sum _{k=A,B}{\tilde{g}}_{\sigma k}|{\tilde{\Psi }}_k{|}^2\right){\tilde{\Psi }}_{\sigma }.$$

$$i\frac{\partial {\Psi }_A}{\partial t}$$

$$i\frac{\partial {\Psi }_B}{\partial t}$$

$$R_{\sigma }^{TF}$$

$$R_{\bar{\sigma }}^{TF}$$

$$\nu =\frac{n_{max}-n_{min}}{n_{max}+n_{min}}.$$

$$\begin{array}{rl}& i\hslash \frac{\partial {\Psi }_{\sigma }(\text{r},t)}{\partial t}\\ & =\left[-\frac{{\hslash }^2{\nabla }^2}{2m_{\sigma }}+V_{\sigma }(\text{r})+\sum _{{\sigma }^{\prime }=A,B}{g}_{\sigma {\sigma }^{\prime }}|{\Psi }_{{\sigma }^{\prime }}(\text{r},t){|}^2\right]{\Psi }_{\sigma }(\text{r},t).\end{array}$$

$${\Psi }_{\sigma }(x,y,z,t)={\psi }_{\sigma }(x,t){\varphi }_{\sigma }(y){\varphi }_{\sigma }(z).$$

Full text

I Introduction

Over the past two decades, nonlinear wave phenomena Kevrekidis et al. (2007) in Bose-Einstein condensates (BECs) Pitaevskii and Stringari (2016); Pethick and Smith (2002) have attracted considerable interest. Early experiments on single component BECs Anderson et al. (1995); Davis et al. (1995); Bradley et al. (1997) triggered a vast amount of studies devoted in investigating the properties of nonlinear excitations that arise in them, such as dark Burger et al. (1999); Frantzeskakis (2010); Becker et al. (2008); Weller et al. (2008); Anderson et al. (2001); Shomroni (2009) and bright Strecker et al. (2002); Khaykovich et al. (2002); Cornish et al. (2006); Dabrowska-Wüster et al. (2007) solitons. These structures exist as stable configurations in highly elongated quasi one-dimensional (1D) BECs featuring repulsive and attractive interatomic interactions respectively. Besides single component BECs, the experimental realization of two-component mixtures Myatt et al. (1997); Hall et al. (1998); Stamper-Kurn et al. (1998); Modugno et al. (2001); Bloch et al. (2001) enabled the investigation of more complex compounds, the so-called vector solitons. These include multiple dark and/or bright states such as dark-dark Hoefer et al. (2011); Yan et al. (2012), and dark-bright (DB) Rajendran et al. (2009); Yin et al. (2011); Yan et al. (2011); Katsimiga et al. (2017a); Busch and Anglin (2001); Nistazakis et al. (2008a); Vijayajayanthi et al. (2008); Brazhnyi and Pérez-García (2011); Álvarez et al. (2011); Achilleos et al. (2011); Álvarez et al. (2013); Achilleos et al. (2012); Middelkamp et al. (2011); Katsimiga et al. (2017b) (see also the review Kevrekidis and Frantzeskakis (2016)), as well as more peculiar entities like dark-antidark (DAD) solitons Danaila et al. (2016); Kevrekidis et al. (2004); Mistakidis et al. (2018); Katsimiga et al. (2017b). The latter antidark state consists of a density hump on top of the matter-wave background. Initially, such solitons were discovered and broadly studied in the field of nonlinear optics Chen et al. (1997); Ostrovskaya et al. (1999); Christodoulides (1988); Afanasyev et al. (1989); Kivshar et al. (1996); Radhakrishnan and Lakshmanan (1995); Buryak et al. (1996); Sheppard and Kivshar (1997); Park and Shin (2000), but they were later on shown to appear also in BECs Becker et al. (2008); Hamner et al. (2011). Focusing on DB solitons multiples of these types of entities have been experimentally realized Hamner et al. (2011). Moreover, considerable attention has been paid in understanding their properties and interactions not only in the so-called integrable Manakov limit Manakov (1974); Prinari et al. (2015) but also upon breaking integrability. The latter breaking can be achieved e.g. by the presence of external trapping Yan et al. (2015); Karamatskos et al. (2015), or by considering different interactions between the same or different atomic species Yan et al. (2015); Katsimiga et al. (2017a); Karamatskos et al. (2015); Katsimiga et al. (2017c, 2018a).

The generation of DB states is a topic of ongoing research activity for a multitude of reasons. Remarkably these states appear in binary BECs with repulsive interparticle interactions, an environment into which bright solitons, being localized density peaks, would disperse if dark solitons were absent. This symbiotic nature of DB solitons leads to their remarkably prolonged lifetimes compared to their single component analogues. Thus DB states constitute ideal candidates for applications as e.g. coherent storage as well as processing of optical fields, multi-channel signals and their switching Rajendran et al. (2009); Folman et al. (2000); Petrov et al. (2009); Rajendran et al. (2011).

DB solitons exist both in weakly miscible Hamner et al. (2011) and immiscible Hall et al. (1998) BEC mixtures. The miscibility condition within mean-field theory is determined by $a_{AB}^{2}<a_{AA}a_{BB}$ Ho and Shenoy (1996); Timmermans (1998); Ao and Chui (1998); Esry and Greene (1998) where $a_{jj}$ is the intra- and $a_{ij}$ denotes the interspecies scattering lengths. In turn these scattering lengths can be experimentally tuned via Feshbach resonances Inouye et al. (1998); Vogels et al. (1997); Thalhammer et al. (2008) allowing the realization of a transition between the miscible and immiscible regime of interactions Papp et al. (2008); Wang et al. (2016); Tojo et al. (2010). The tunability of the scattering lengths opened a new path for studying, in a controllable manner, the properties of binary BECs and the solitonic structures that form in them. In this way, vector solitons have been identified upon considering e.g. temporal Kanna et al. (2013); Manikandan et al. (2016); Yu (2016) and spatial Cheng (2009); Belmonte-Beitia et al. (2011); Mareeswaran et al. (2014) variations of the scattering lengths. In the same context control over the system’s parameters has been examined by performing temporal modulations of the confinement frequencies Kanna et al. (2013); Manikandan et al. (2016); Yu (2016) and also by adding appropriate time-dependent gain and loss terms Rajendran et al. (2009, 2011). Moreover, control over the system’s dynamical evolution and pattern formation has been reported in Eto et al. (2016), upon quenching the binary BEC towards miscibility. Also very recently the spontaneous generation of vector solitons upon quenching the interspecies interactions crossing the miscibility/immiscibility threshold both within mean-field but also considering the many-body treatment of this problem has been reported Mistakidis et al. (2018). Here, robustly propagating DAD states were shown to dynamically emerge. The latter observation rises a natural question on how different variations of the systems’ parameters can lead to the controllable creation of these states at least within mean-field theory.

Motivated by the aforementioned studies, in the present work we take a closer look on the multi-soliton formation upon considering the interspecies interaction-quenched dynamics of a particle imbalanced BEC mixture while crossing the miscibility/immiscibility threshold in both directions. In this way we aim at deepening our understanding regarding vector soliton generation in an integrability broken scenario. In this context dynamical instabilities lead to peculiar evolution of the resulting states including mass redistribution between the solitons Katsimiga et al. (2017c, 2018a), and thus shape changing collisions Rajendran et al. (2011); Manikandan et al. (2014). In particular, we investigate how the number of vector solitons formed can be controlled under different variations of the binary systems’ parameters as e.g. the trapping frequency and the associated with each species particle number.

In the miscible-to-immiscible quench scenario, the unstable dynamics leads to the filamentation of the density of both species in line with our earlier findings Mistakidis et al. (2018). This filamentation process entails the formation of DB solitary waves. The formation of the latter is verified upon fitting the waveforms corresponding to a dark and a bright soliton stemming from the exact (in the Manakov limit) single DB state Yan et al. (2011); Katsimiga et al. (2017a, 2018a). Focusing on the initial stages of the nonequilibrium dynamics, an increase of the number of the DB states generated is found upon varying the particle number of the minority species until the particle balanced limit is reached. A result that is more enhanced for transitions that enter deeper in the immiscible phase Mistakidis et al. (2018). Lastly, a decaying formation of DB states is observed as the trapping frequency is increased. In all cases both clouds undergo a collective breathing motion. Moreover, the generated DB states are seen to oscillate within the parabolic trap, featuring multiple collisions with one another. However, due to the non-integrable nature of the system, these collisions are also accompanied by a transfer of mass between the soliton constituents Katsimiga et al. (2017c, 2018a). The latter effect reduces, during evolution, the number of states generated with significantly fewer DB solitary waves remaining and robustly propagating at large evolution times. Finally, inspecting the corresponding breathing frequency of each bosonic cloud we demonstrate that significant differences occur for highly particle imbalanced mixtures Sartori and Recati (2013).

Following the reverse process, namely an immiscible-to-miscible transition, robust DAD solitons are formed Mistakidis et al. (2018). The binary system, is the 1D analogue of the so-called “ball and shell” configuration Pu and Bigelow (1998); Mertes et al. (2007), with one species occupying the center of the trap and the other species residing on the two sides of the harmonic potential. After the quench the outside species is allowed to “fall” towards the central region inducing a counterflow. The latter leads in turn to the appearance of interference fringes Andrews et al. (1997); Naraschewski et al. (1996); Hoston and You (1996) and dark solitons emerge in the regions of destructive interference Hamner et al. (2011). It is the effective potential created by the emergent dark states that causes the other species to split and fill the newly formed density dips, leading in turn to the formation of DAD structures. An almost linear increase of the number of DAD states formed occurs either as the particle number of the inner species is increased or entering deeper in the miscible phase. Also, tighter trapping results to fewer DAD solitons occurring during the nonequilibrium dynamics for larger postquench interspecies interactions. Finally, different breathing frequencies are measured for the two species independently of the postquench interaction strength and the particle imbalance. A remark is appropriate here. Note that the quench-induced dynamics studied in the present work is restricted to a 1D geometry. As such, the dynamical formation of the DB and the DAD solitons identified herein, can be tested in current quasi-1D experiments as it will be demonstrated later on. However, we must stress at this point that alterations of the observed interaction quench dynamics can take place in higher dimensional settings, e.g. due to the presence of transversal instabilities of dark solitons.

Our presentation is structured as follows. In section II we provide the theoretical framework of our mean-field approach in the form of a set of two coupled Gross-Pitaevskii equations (GPEs). We also briefly comment on the numerical methods utilized herein. In section III we present our results for miscible-to-immiscible quenches while section IV contains our results for the reverse process. The robustness of our results in a quasi 1D harmonic trap is demonstrated in section V. Finally section VI summarizes our findings and also provides future perspectives.

II Theoretical Framework

II.1 Mean-Field ansatz for a binary Bose-mixture

We consider a binary mixture of repulsively interacting BECs composed of two different hyperfine states, namely $\ket{F=1,m_{F}=-1}$ and $\ket{F=2,m_{F}=1}$, of 87Rb Egorov et al. (2013) being confined in a one-dimensional (1D) harmonic oscillator potential. Such a cigar-shaped geometry can be realized experimentally Becker et al. (2008); Hoefer et al. (2011); Middelkamp et al. (2011) upon considering a highly anisotropic trap with the longitudinal and transverse trapping frequencies obeying $\omega_{x}\ll\omega_{\perp}$. Within mean-field theory the dynamics of this binary mixture can be described by the following system of coupled GPEs Kevrekidis et al. (2007); Pitaevskii and Stringari (2016):

[TABLE]

In the above expression $\tilde{\Psi}_{\sigma}(\tilde{x},\tilde{t})$ ($\sigma=A,B$) denotes the wavefunction for the $A$- and $B$-species respectively. Each $\tilde{\Psi}_{\sigma}(\tilde{x},\tilde{t})$ is normalized to the corresponding number of atoms i.e. $N_{\sigma}=\int_{-\infty}^{+\infty}|\tilde{\Psi}_{\sigma}|^{2}d\tilde{x}$. Also $\tilde{V}(\tilde{x})$ represents the external trapping potential, while $m_{A}=m_{B}=m$ and $\mu_{\sigma}$ refer to the atomic mass and chemical potential for each of the species respectively. The effective 1D coupling constants are given by $\tilde{g}_{\sigma k}=2\hbar\omega_{\perp}a_{\sigma k}$, where $a_{\sigma k}$ are the three $s$-wave scattering lengths (with $a_{AB}=a_{BA}$) accounting for collisions between atoms that belong to the same ($a_{\sigma\sigma}$) or different ($a_{\sigma k},\sigma\neq k$) species. We note that both the intra and interspecies scattering lengths can be manipulated experimentally by means of Feshbach Köhler et al. (2006); Chin et al. (2010) or confinement induced resonances Olshanii (1998); Kim et al. (2006).

For the numerical findings to be presented below, we choose to express the system of Eqs. (1) in the following dimensionless form:

[TABLE]

In the above Eqs. (2)-(3), the interaction coefficients are normalized to the intraspecies scattering length, $a_{AA}$, of the $A$-species, namely $g_{AB}\equiv\tilde{g}_{AB}/\tilde{g}_{AA}$, and $g_{BB}\equiv\tilde{g}_{BB}/\tilde{g}_{AA}$. Moreover, densities $|\Psi_{\sigma}|^{2}$, length, energy and time are measured in units of $2a_{AA}$, $\sqrt{\frac{\hbar}{m\omega_{\perp}}}$, $\hbar\omega_{\perp}$ and $\omega_{\perp}^{-1}$ respectively. Finally, $V(x)=\frac{1}{2}\omega^{2}x^{2}$ is the dimensionless trapping potential with $\omega_{A}=\omega_{B}\equiv\omega=\omega_{x}/\omega_{\perp}$.

In the present work we aim at revealing the out-of-equilibrium dynamics of the aforementioned binary system [see Eqs. (2)-(3)] upon performing an interspecies interaction quench from the miscible to the immiscible regime of interactions and vice versa. To this end we fix the intraspecies coefficients to the experimentally relevant values for a 87Rb mixture, i.e. $g_{AA}=1$, $g_{BB}=0.95$ Egorov et al. (2013). To induce the dynamics an interspecies interaction quench is performed with the system being initially relaxed to its ground state configuration lying within the miscible regime of interactions. The latter is characterized by an initial (prequench) interaction $g^{i}_{AB}$ subject to the miscibility condition $g^{i}_{AB}<\sqrt{g_{AA}g_{BB}}$ Ao and Chui (1998). Consecutively the system is quenched towards the immiscible (alias phase separated) regime having the desired postquench amplitude $g^{f}_{AB}$ [see Sec.III]. The same procedure is followed for the inverse process, namely considering a relaxed stated for immiscible components ($g^{i}_{AB}>\sqrt{g_{AA}g_{BB}}$) and quenching the system towards the miscible domain [see Sec.IV].

In both quench scenarios under investigation in order to obtain the ground state configuration a fixed-point numerical iteration scheme is employed Kelley (2003). To simulate the nonequilibrium dynamics of the binary mixture governed by Eqs. (2)-(3) a fourth-order Runge-Kutta integrator is employed and a second-order finite differences method is used for the spatial derivatives. The spatial and time discretization are $dx=0.1$ and $dt=0.005$ respectively. Moreover our numerical computations are restricted to a finite region by employing hard-wall boundary conditions. The latter are chosen wide enough to avoid boundary effects. In particular, in the dimensionless units adopted herein and also for the majority of our simulations, the hard-walls are located at $x_{\pm}=\pm 80$ and we do not observe any appreciable density for $|x|>60$.

III Quench Dynamics from the Miscible to the Immiscible Regime

In the following we discuss the vector soliton formation considering quenches from the miscible to the immiscible regime, namely for $g_{AB}^{f}>\sqrt{g_{AA}g_{BB}}$. Before proceeding, it is worth mentioning at this point that the latter relation is valid for the homogeneous case (V=0). Nonetheless, corrections due to the presence of the trap are only significant for strong trapping frequencies Navarro et al. (2009) that are not considered in the present work. In what follows the trapping frequency is fixed to $\omega=0.05$ unless it is stated otherwise. Additionally, since the miscibility depends on the interaction strength between the species, for a 87Rb mixture considered herein the transition between the two regimes occurs at $g^{th}_{AB}=\sqrt{0.95}\approx 0.974$. Thus for quenches with $g^{f}_{AB}<g^{th}_{AB}$ it is expected that the two species prefer overlapping thereby residing around the trap center. However, for postquench interspecies interactions having $g^{f}_{AB}>g^{th}_{AB}$ component separation takes place.

III.1 Density evolution within the miscible domain and breathing mode

Figures 1 (a)-(f) illustrate the quenched density evolution upon increasing $g_{AB}^{f}$. In all cases the system is initially prepared in the miscible regime with $g_{AB}^{i}=0.2$, while $N_{A}=65$ and $N_{B}=25$. As per our discussion above when the interaction is quenched e.g to $g_{AB}^{f}=0.8$ depicted in Figs. 1 (a)-(b) for the $A$-species and the $B$-species respectively, component separation is absent. Instead the two species are overlapping while residing around the trap center for all evolution times. Notice also that since $g_{BB}<g_{AA}$, the $B$-species is localized closer to the trap center when compared to the $A$-species. Moreover, a collective breathing motion consisting of a periodic expansion and contraction of each bosonic cloud is observed during evolution.

Let us now investigate the effect that different variations of the binary systems’ parameters have on the breathing frequency. In all cases illustrated in Figs. 2 (a)-(c) in order to obtain the breathing frequency, $\omega^{\sigma}_{br}$, we start from the non-interacting limit ($g^{i}_{AB}=0$) and we measure the Fourier transform of $\braket{x^{2}}=\int^{+\infty}_{-\infty}x^{2}|\Psi_{\sigma}|^{2}dx$ for each species $\sigma=A,B$. In this way we deduce $\omega^{\sigma}_{br}$ upon varying either $g^{f}_{AB}$ [see Figs. 2 (a) and (b)] or the particle number, $N_{B}$, of the $B$-species [see Fig. 2 (c)]. It is found that for particle number imbalances $N_{B}/N_{A}=0.5$, $\omega_{br}^{\sigma}$ is independent of the $g^{f}_{AB}$ variation acquiring the constant value of $\omega_{br}\approx\sqrt{3}\omega$ Sartori and Recati (2013), see Fig. 2 (a). However, for larger particle imbalances as e.g. $N_{B}/N_{A}=0.1$ depicted in Fig. 2 (b) a completely different behavior is observed, with the breathing frequency of the minority $B$-species being drastically affected by the postquench variation. In particular as $g^{f}_{AB}$ increases towards the miscibility/immiscibility threshold, a monotonic decrease of $\omega^{B}_{br}$ is observed Sartori and Recati (2013), e.g. with $\omega^{B}_{br}=0.5$ for $g^{f}_{AB}=0.9$. This decrease of $\omega^{B}_{br}$ suggests that instead of the harmonic trapping a modified frequency comes into play. Since the $A$-species is the majority component, we assume that it creates an effective potential Mistakidis et al. (2019); Ferrier-Barbut et al. (2014) of the form $V_{eff}(x)=V(x)+g_{AB}^{f}|\Psi_{A}(x,0)|^{2}$ [see Eqs. (2)-(3)] into which the $B$-species is trapped and that $|\Psi_{A}(x,0)|^{2}$ has roughly the form of a Thomas-Fermi (TF) profile. As a consequence we obtain the following effective trapping frequency $\omega_{eff}=\omega\sqrt{1-\frac{g_{AB}^{f}}{g_{AA}}}$. Thus the corresponding effective breathing frequency is $\omega^{B,eff}_{br}=\sqrt{3}\omega_{eff}$. A remarkable agreement between $\omega^{B,eff}_{br}$ and $\omega^{B}_{br}$ is also shown in Fig. 2 (b) (see dashed black line). Finally, in line with the above discussion, deep in the miscible regime $\omega^{B}_{br}$ shows an increasing tendency for larger $N_{B}$ and acquires the constant value of $\omega_{br}\approx\sqrt{3}\omega$ for imbalances $N_{B}/N_{A}\gtrsim 0.3$ [see Fig. 2 (c)].

III.2 Density evolution towards the immiscible regime and identification of DB states

In contrast to the above-observed dynamics within the miscible domain, for values of $g^{f}_{AB}$ that are above a certain threshold the miscible-to-immiscible phase transition takes place. As it is evident in Figs. 1 (c)-(f) the dynamics in this regime is unstable and the initially localized, for each species, configuration breaks into several filaments Mistakidis et al. (2018). The filaments refer to the individual density branches appearing in $\left|\Psi_{\sigma}\right|^{2}$ during evolution. For further details regarding this dynamical instability we refer the reader e.g. to Ref. Tommasini et al. (2003) and references therein where the homogeneous case is investigated or the recent work of Ref. Mistakidis et al. (2018) that includes the presence of the trap. By comparing the filamentation of each density for two different postquench interactions, we observe that for values of $g_{AB}^{f}$ that enter deeper into the immiscible phase more filaments are formed [compare Figs. 1 (c) and (d) to Figs. 1 (e) and (f) respectively]. Yet another important observation is that in all cases the number of filaments formed at the initial stages of the out-of-equilibrium dynamics is significantly larger compared to the remaining number of filaments robustly propagating at larger evolution times. In particular e.g. in Figs. 1 (c) and (d) around $t=35$ eight filaments are found to be spontaneously generated in both species [see also Figs. 3 (a), (b)] while only four and three of them are present in the $A$ and $B$-species respectively for $t\geqslant 200$.

Focusing on the initial stages of the dynamics, let us closer inspect the filaments formed. Figures 3 (a) and (b) present profile snapshots of the densities $|\Psi_{A}(x)|^{2}$ and $|\Psi_{B}(x)|^{2}$ of the majority and the minority species respectively at $t=35$ and for $g_{AB}^{f}=1.5$. Indeed, eight filaments are formed in each of the species with the density dips of the majority one being filled by density peaks of the minority species. Since such a filling mechanism resembles the formation of DB solitons in defocusing media below we attempt an identification of the structures formed to these symbiotic states.

As a first step towards confirmation we invoke the exact, at the integrable Manakov limit, single DB soliton solution Yan et al. (2011); Katsimiga et al. (2017a, 2018a). The latter reads $\Psi_{A}(x,t)=A\left(\cos\phi\tanh\left[D\left(x-x_{0}(t)\right)\right]+i\sin\phi\right)$ and $\Psi_{B}(x,t)=B\operatorname{sech}\left[D\left(x-x_{0}(t)\right)\right]e^{ikx+i\theta(t)+i(\mu_{B}-\mu_{A})t}$. Here, $\phi$ denotes the soliton’s phase angle, while $A\cos\phi$, and $B$ refer to the amplitude of the dark and the bright soliton respectively. Furthermore, $D$ and $x_{0}(t)$ correspond to the common inverse width and the soliton’s center respectively. Finally, $k=D\tan\phi$ is the constant wave-number of the bright soliton and $\theta(t)$ its phase. The above expressions are used for the density profile fits shown with dashed lines in Figs. 3 (a) and (b). Evidently a remarkably good agreement between our numerical findings and the fitted profiles is seen [see also Table 1]. The latter verifies that indeed the filamentation process leads to structures that possess a DB solitary wave character. We note that these emergent dark states are characterized by a phase jump being a multiple of $\pi$ (results not shown here for brevity). The same overall phenomenology is observed for even higher values of $g^{f}_{AB}$ such as the one illustrated in the profile snapshots of Figs. 3 (c) and (d) for $g^{f}_{AB}=2.5$. Here, the number of DB states is significantly increased compared to the $g^{f}_{AB}=1.5$ scenario but again an excellent agreement between the numerical data and the fitted DB waveforms can be inferred [see also Table 1]. Notice also that since in this case the number of solitary waves generated is large, already from these early stages of the dynamics collisions between the DB states are clearly captured by this snapshot. Indeed, the collisions leave their fingerprints in the deformed densities, see the dashed circles which mark the two most inner pairs symmetrically placed around $x=0$ depicted in Figs. 3 (c) and (d).

After the initial stages of formation the DB solitary waves are left to dynamically evolve within the parabolic trap [see Figs. 1 (c)-(d) for $t>50$]. A closer inspection of the evolution of the density of each species reveals a collective breathing motion with four contraction or focus points where the DB states are closest. It is this breathing motion that leads in turn to the collision events that occur during evolution. Several collision events take place both around the center of the trap as well as at the edges of each bosonic cloud. These events are much more pronounced for smaller $g^{f}_{AB}$ [compare Figs. 1 (c)-(d) with Figs. 1 (e)-(f)] with the DB states that are less mobile, i.e. those closer to the trap center, featuring more dramatic collisions. Recalling now that we operate in the nonintegrable limit these collisions can be rather asymmetric Katsimiga et al. (2018a) leading to a transfer of mass between the solitary wave constituents and thus to shape changing collisions Rajendran et al. (2011); Manikandan et al. (2014); Katsimiga et al. (2018a). It is this mass redistribution that results in a decreasing number of DB solitons but renders them wider and robustly propagating for large evolution times.

In order to monitor such events and shed light on the observed dynamics we consider several collision points at different time instants during evolution [see the circles in Figs. 1 (c), (d)]. For these selected events we perform fittings prior and after its collision vertex. In doing so we immediately have access to the relevant amplitudes and widths of the states and we can further compare the number of particles, $N_{b}$, contained in each bright solitary wave prior and after a collision. In order to perform the latter calculation we make use of the analytical expression for the single DB state Yan et al. (2011); Katsimiga et al. (2017a, 2018a), i.e. $N_{b}=\int~{}dx~{}|\Psi_{B}|^{2}=2B^{2}/D$, that is exact at the integrable limit but approximately holds also in our setup 111For the numerical evaluation of $N_{b}^{num}$ for a selected DB pair we used $N_{b}^{num}=(1/N_{B})\int_{c_{1}}^{c_{2}}dx\left|\Psi_{B}\right|^{2}$. Here, $N_{B}$ is the total number of particles in the $B$-species and the integration interval $[c_{1}$, $c_{2}]$ refers to the spatial area around the center of the selected bright soliton. The observed deviation between $N_{b}^{num}$ and $N_{b}$ is smaller than $1\%$.. The corresponding outcomes of such a process are summarized in Table 1 for e.g. postquench amplitude $g^{f}_{AB}=1.5$. In particular, in this Table we present the fitted solitary wave parameters when monitoring the two most outer DB states labeled as “1” and “2” in Figs. 1 (c)-(d). These DB waves participate in two distinct collision events that take place during evolution one at around $t\approx 140$ and a second one around $t\approx 210$ [see the light blue and green circles in Figs. 1 (c)-(d)]. By measuring the number of particles, $N_{b}$, contained in each bright solitary wave prior and after the first collision namely at $t=106$ and $t=175$ respectively (see Table 1) it is found that indeed a mass redistribution between the solitary waves occurs (see the corresponding boldfaced $N_{b}$ values in Table 1). Notice that after the first collision the number of particles contained in “2” is almost two times that of state “1”. Remarkably, at later propagation times, i.e. after the second collision takes place, an almost complete transfer of mass results in a single but wider bright solitary wave (see Table 1 at $t=250$) that remains and robustly propagates in the BEC medium. Notice that its corresponding dark counterpart is also wider [see Fig. 1 (c) at $t=250$]. We remark that the symmetric, with respect to $x=0$, DB states have exactly the same dynamical evolution.

Finally, we perform the relevant fittings for the huge merger that appears around the trap center after the collision of at least three DB pairs that takes place at around $t\approx 140$ [see the blue circle in Fig. 1 (d)]. From this fitting we can verify via the particle number conservation prior and after the collision point, that indeed a total transfer of mass of the two outer DB pairs to the central one occurs. In particular, $N_{b_{I}}=4$, $N_{b_{II}}=5$ and $N_{b_{III}}=4$ prior the collision leading to the single DB solitary wave with $N_{b}=13$.

III.3 Controlling the number of DB states

Having identified the DB solitary waves generated via the filamentation process, in what follows we aim at controlling their formation. Thus we investigate how the number of DB states, $N_{DB}$, changes upon considering different variations of the binary systems’ parameters.

In particular, Fig. 4 $(a)$ illustrates $N_{DB}$ upon increasing the particle number, $N_{A}$, of the $A$-species for $N_{B}=50$. As it can be seen, solitary wave formation occurs already from small particle numbers ($N_{A}>10$) in the $A$-species. As $N_{A}$ increases $N_{DB}$ increases in a square root fashion with more DB states being generated for larger postquench interactions $g^{f}_{AB}$. This behaviour can be intuitively understood by considering the size of the two BECs. To estimate the system’s size, we start from Eqs. (2)-(3) and by making use of the TF approximation we derive, within the miscible regime, the following expressions for the corresponding radii of the two clouds Pethick and Smith (2002)

[TABLE]

In the above expressions, $\sigma=A,B$ ($\bar{\sigma}=B,A$) denotes the majority (minority) component. It is now important to note that DB solitary waves will only form in regions where both species are present. The latter sets the scale for solitary wave formation as the region defined by the smaller of the two aforementioned radii. In the inset of Fig. 4 $(a)$ $R^{TF}_{A}$, and $R^{TF}_{B}$ are shown as a function of $N_{A}$ for $g^{f}_{AB}=1.5$. For $N_{A}<N_{B}$ the TF radius of the $A$-species is the smaller one and as such is provided by Eq. (5). Indeed, $R^{TF}_{A}$ increases until the particle balanced limit is reached ($N_{A}=N_{B}=50$) where $R^{TF}_{A}=R^{TF}_{B}$ 222Note that these equalities are exact at the integrable limit, but approximately hold also here with a deviation of about $2\%$. It is in this region, i.e. $N_{A}<N_{B}$, that the growth rate of $N_{DB}$ is significant. On the other hand when $N_{A}>N_{B}$, $R^{TF}_{B}$ becomes the smaller radius. In this case, unlike $R_{A}^{TF}$, $R^{TF}_{B}$ remains almost constant for increasing $N_{A}$. Thus, since $R^{TF}_{B}$ remains constant as $N_{A}$ is increased, a plateau of almost constant DB soliton production appears.

The same qualitative picture is expected to hold also upon varying the number of particles, $N_{B}$, of the second species, since fixing $N_{A}$ and varying $N_{B}$ simply interchanges the role of the two components. Indeed as illustrated in Fig. 4 $(b)$, solely focusing on the case where $N_{B}<N_{A}$, again a square-root-like increase of $N_{DB}$ is observed leading to more solitary waves being generated as we enter deeper into the immiscible regime of interactions. Notice that for $N_{B}\geq N_{A}$ again an almost constant production of DB states is observed. However since in this case $N_{A}$ is fixed to a larger value compared to the previous variation, $N_{DB}$ acquires larger values before the particle balanced limit is reached.

The influence of the postquench interspecies interaction $g_{AB}^{f}$ on the DB formation is shown in Fig. 4 $(c)$. Here and as per our discussion above, for $g^{f}_{AB}<g^{th}_{AB}\equiv 0.974$ (i.e. within the miscible domain) $N_{DB}=0$. Solitary wave formation takes place only for values of $g^{f}_{AB}>g^{th}_{AB}$, with $N_{DB}$ increasing in an almost square root manner for increasing postquench interactions. Additionally here, $N_{DB}$ is larger even for slightly larger particle number imbalances. Recall that the observed dynamics shown in Figs. 1 (c)-(f) dictates also an increase of $N_{DB}$ being generated faster for larger $g_{AB}^{f}$. The latter is consistent with earlier predictions Mistakidis et al. (2018) regarding particle balanced mixtures.

Furthermore, in order to examine the effect of the trapping geometry on the solitary wave formation, we next calculate $N_{DB}$ upon varying the trapping frequency $\omega$. The outcome of this parametric variation is depicted in Fig. 4 $(d)$. $N_{DB}$ decreases rapidly as $\omega$ increases but acquires larger values for larger $g^{f}_{AB}$. Note also that as we approach the homogeneous limit, $\omega\rightarrow 0$, $N_{DB}$ tends to infinity a result that is explained by the increasing size of the condensates. Indeed the radii of the condensates [see Eqs. (4)-(5) depend on $\omega^{-2/3}$ which is in accordance to what we observe.

Finally, let us stress at this point that the observed dynamical evolution of the binary mixture can be experimentally monitored by performing in-situ single-shot absorption measurements Mistakidis et al. (2018); Bloch and Zwerger (2008); Katsimiga et al. (2017d). Indeed, the DB solitary waves building upon the densities of the different species can be captured by such images. For this purpose we further consider mixtures with larger particle numbers in both species showcasing the generalization of our findings. Evidently, as can be seen in Table 2 the same overall controlled generation of DB states persists upon varying distinct parameters of the system.

IV Quench Dynamics from the Immiscible to the Miscible Regime

Up to now we only considered quenches from the miscible to the immiscible regime of interactions. Next, immiscible-to-miscible transitions will be investigated. For this reverse process, the binary system is initially prepared in its ground state having a fixed $g_{AB}^{i}=1.5$. Then, the $A$-species resides at the edges of the harmonic trap since it possesses a larger intraspecies coefficient ($g_{AA}$) while the $B$-species occupies the central region. Such an initial state configuration represents the 1D analogue of the so-called “ball and shell” structure found in higher dimensional mixtures Pu and Bigelow (1998); Mertes et al. (2007). After this initial state preparation the system is abruptly quenched to a lower $g_{AB}^{f}$ value lying within the miscible regime of interactions.

IV.1 Effective counterflow dynamics and emergence of DAD states

Case examples of the spatio-temporal evolution of the densities of both species for two different particle number imbalances, namely for $N_{B}/N_{A}=1/3$ and $N_{B}/N_{A}=3$, but for the same postquench interspecies repulsion are shown in Figs. 5 (a), (b) and (c), (d) respectively. In particular, Figs. 5 (a) and (c) depict the evolution of the density of the $A$-species and Figs. 5 (b) and (d) the corresponding evolution of the $B$-species. In both cases pattern formation takes place from the early stages of the nonequilibrium dynamics but with a dramatically larger number of states formed when the inner $B$-species is the majority component of the bosonic mixture. The latter outcome can be inferred by inspecting the corresponding profile snapshots illustrated in Figs. 6 (a)-(d). Here, Figs. 6 (a), (b) and (c), (d) show the density profiles of both species at the initial stages of formation, i.e. at $t=45$ and $t=60$ respectively, for $N_{B}/N_{A}=1/3$ and $N_{B}/N_{A}=3$. Evidently, at least twelve localized entities appear in the former snapshots, namely when the outer $A$-species is the majority component, while at least thirty occur when the inner $B$-species is the majority component. In this latter case the localized structures closely resemble interference fringes. Notice that in both cases, since component mixing is favored in this regime, the $B$-species tend to fill the quench induced minima of the $A$-species. However, when the $B$-species is the majority component an overfill of the corresponding $A$-species minima occurs [see the relevant profiles in Figs. 6 (b), (d)] and thus no robust localized configuration can be realized. Before delving in the associated dynamics, we should emphasize that in contrast to the previous quench scenario, the structures formed in this case are not only different in nature when compared to the DB solitary waves previously generated, but appear also to be more robust during evolution [compare Figs. 5 (a)-(d) to Figs. 1 (c)-(f)].

Thus the breaking of integrability within the miscible domain doesn’t seem to affect the evolution of the states generated, a result that is in line with earlier predictions Mistakidis et al. (2018).

In particular, a closer inspection of the spatio-temporal evolution of the densities [see Figs. 5 (a)-(d) and Figs. 6 (a)-(d)] reveals that dark soliton entities emerge during evolution in the $A$-species, while antidark solitons develop in the corresponding $B$-species. Note that the antidark solitons are states consisting of a density peak on top of the BEC background Danaila et al. (2016); Kevrekidis et al. (2004); Mistakidis et al. (2018); Katsimiga et al. (2017b). Thus in contrast to the DB solitary waves formed in the previous quench scenario, here multiple DAD solitons are spontaneously generated. The formation of these DAD states has been reported also in the particle balanced case Mistakidis et al. (2018). Below, and in contrast to earlier findings we will explore in detail the mechanism of the formation of these DAD states and we will also attempt to control their generation upon considering different variations of the binary systems’ parameters.

To unravel the underlying mechanism of formation of the aforementioned DAD states we must first recall that the initial state configuration consists of two spatially separated components. As the system is dynamically quenched towards the species miscible regime, the interspecies interaction energy, $\Delta E_{int}=2\left|\Psi_{A}\right|^{2}\left|\Psi_{B}\right|^{2}(g_{AB}^{f}-g_{AB}^{i})$, decreases the more the smaller the $g_{AB}^{f}$ becomes. Then, the $B$-species is no longer deemed to reside at the central region of the parabolic potential but instead it can expand. More importantly the $A$-species residing previously symmetrically around the edges of the trap is now allowed to “fall” towards the center. It is this counterflow of the $A$-species that causes the formation of the dark solitons right at the trap center where the two parts of the gas destructively interfere Andrews et al. (1997). Note also that in our initial state preparation we start with a constant zero phase in the $A$-species and thus we expect and indeed observe an even number of dark solitons to emerge Zakharov and Shabat (1973) [see also here the profiles shown in Figs. 6 (a), (c)]. This dark soliton formation leads in turn, due to the effective potential that these states create, to the breaking of the $B$-species into several density peaks each of which filling the newly formed density dips. Essentially it is energetically beneficial for the $B$-species to avoid regions where the $A$-species has a high density, since the interaction between the species is repulsive ($g^{f}_{AB}>0$). Thus the $B$-species tends to fill the $A$-species minima. This effect is more severe for weaker repulsions that enter deeper in the miscible regime, with the number of the emergent DAD states found to increase for smaller $g_{AB}^{f}$ (see our discussion below). However since the binary system is in the miscible regime where a mixed state is favorable the generated density peaks are formed on top of their BEC background. The above-discussed mechanism of formation of the observed DAD states is fundamentally different in nature, and is to be contrasted with the unstable dynamics that leads to the formation of the DB structures identified in the miscible-to-immiscible transition.

To testify the formation of dark solitons in our setup, we next rely on predictions stemming from counterflow experiments Weller et al. (2008). It is known that for single component BECs a critical distance between the two initially separated BECs exists below which soliton formation takes place Scott et al. (1998). This critical distance in the dimensionless units adopted herein reads $D_{c}=\pi\left(\frac{6\pi N_{A}g_{AA}}{\omega_{eff}^{A}}\right)^{1/3}$. Here, $\omega_{eff}^{A}=\omega\sqrt{1-\frac{g_{AB}^{f}}{g_{BB}}}$ is the modified trapping frequency that incorporates the postquench cross interaction term $g_{AB}^{f}$. This effective potential picture Mistakidis et al. (2019); Ferrier-Barbut et al. (2014), $V_{eff}=V(x)+g_{AB}^{f}|\Psi_{B}(x,0)|^{2}$, is analogous to the one introduced in section III.

We remark here that in contrast to the previous quench scenario [section III] the $A$-species performs the counterflow dynamics with the $B$-species residing at the trap center. Thus, we consider the dynamics of the $A$-species into the effective potential formed by the $B$-species and the external trap. For the case example of $N_{B}/N_{A}=1/3$ presented in Fig. 5 (a) $D_{c}=197$ while $D=20.8$, verifying that the observed structures building upon the majority $A$-species are indeed dark solitons.

IV.2 Breathing mode

Besides the formation of the DAD structures the two species also undergo a collective breathing motion as it can be readily seen in Fig. 5. Results for the corresponding normalized breathing frequency, $\omega_{br}^{\sigma}/\omega$, are shown in Figs. 7 (a), (b), upon considering the effect of different particle number imbalances and postquench interspecies interaction strengths respectively. Referring to a fixed postquench interaction strength $\omega_{br}^{A}$ and $\omega_{br}^{B}$ are significantly different especially in the region where $N_{A}>N_{B}$, see Fig. 7 (a). This deviation is much more pronounced for smaller $g_{AB}^{f}$. Within the region $N_{A}<N_{B}$, the frequencies $\omega_{br}^{A}$ and $\omega_{br}^{B}$ tend to approach each other but they never coincide. For fixed particle number imbalances, i.e. $N_{B}/N_{A}=1/3$, an increase of $\omega_{br}^{\sigma}/\omega$ is observed as the system is quenched to lower values of $g^{f}_{AB}$ that enter deeper in the miscible regime of interactions [see Fig. 7 (b)]. More specifically, both species appear to have almost the same breathing frequency for quenches near the miscibility/immiscibility threshold but deviate significantly for values of $g^{f}_{AB}$ within the interval $[0.6,0.1]$. In this later case, $\omega_{br}^{B}/\omega$ appears to be higher for all values of the postquench interaction $g^{f}_{AB}$ until the noninteracting limit is reached, with both species acquiring the value $\omega_{br}^{\sigma}/\omega\approx\sqrt{3}$. Moreover, the effective breathing frequency, $\omega_{br}^{A,eff}$ (see also section III), seems to underestimate $\omega_{br}^{\sigma}/\omega$ the more, the closest we are to the miscibility/immiscibility threshold. The observed deviations are mainly attributed to the fact that the effective model assumes a TF profile for the minority $B$-species which differs considerably from the actual density of this species due to the dynamical generation of antidark solitons.

IV.3 Controlled creation of DADs

Having discussed the mechanism of formation of DAD states, we next investigate whether their number can be tuned via parametric variations. Figures 8 (a)-(d) summarize our findings. We find that the number of DAD states, $N_{DAD}$, strongly depends on variations of the particle number either of the $A$- or the $B$-species [see Figs. 8 (a), (b)]. In particular, in both cases illustrated in Fig. 8 (a) a linear increase of $N_{DAD}$ for increasing $N_{A}$ is observed. This increase is followed by a plateau of constant soliton count until the particle balanced limit is reached, with $N_{DAD}$ being larger for larger $N_{B}$. Further increase of $N_{A}$, namely for $N_{B}/N_{A}<1$ leads to a gradual decrease of $N_{DAD}$ with less entities formed for smaller ratios. This behavior can be intuitively explained by recalling that increasing $N_{A}$ also increases the size of the $A$-species. The latter initially reside at the edges of the trap squeezing more and more the central $B$-species and thus reducing the initial separation of the outer BECs. Smaller initial spatial separation leads in turn to less kinetic energy for the counterflowing parts of the $A$-species and thus to less solitons being spontaneously generated. To inspect more carefully the case of $N_{B}/N_{A}<1$, Fig. 8 (b) shows $N_{DAD}$ for varying $N_{B}$ and for different values of $g_{AB}^{f}$ but fixed $N_{A}$. We observe that increasing $N_{B}$ leads to a linear increase of $N_{DAD}$ with the latter being larger for smaller $g_{AB}^{f}$. Also $N_{DAD}$ is larger for smaller ratios of $N_{B}/N_{A}$ and fixed $g_{AB}^{f}$.

The impact of $g_{AB}^{f}$ on $N_{DAD}$ is also examined in Fig. 8 (c) for constant particle imbalances. As it can be seen, the number of solitons increases in an almost linear manner for smaller $g_{AB}^{f}$, namely deeper in the miscible domain. The above-mentioned tendency can be understood as follows. For smaller values of $g_{AB}^{f}$ the two species featuring less repulsion tend to be closer to one another favoring a phase mixed state. As such the $A$-species acquires more kinetic energy during the dynamics caused by a quench to weaker $g_{AB}^{f}$ and thus resulting in a larger $N_{DAD}$. However, it is important to note that as the noninteracting limit is approached the DAD states become gradually less robust.

Focusing on the antidarks formed in the $B$-species it is found that as $g^{f}_{AB}$ decreases the amplitude of these states also decreases. To appreciate this behaviour we introduce the visibility

[TABLE]

Here, $n_{max}$ is the density at the peak of an antidark and $n_{min}$ is the corresponding background density in the minimum. In the case that the minimum densities are notably different on the left/right side of the peak, we define $n_{min}$ as the average of the two minima. As an example the maximum visibility of the central antidark soliton, is $\nu=80.5\%$ in Fig. 6 (a) and grows to $\nu=82.8\%$ in Fig. 6 (c), while is found to be $\nu=20.4\%$ and $\nu=25.7\%$ in Figs. 6 (b), (d) respectively. Moreover, the maximum visibility measured for the corresponding central antidark state for $g_{AB}^{f}=0.1$ is $\nu=10\%$, while it approaches zero for $g_{AB}^{f}\rightarrow 0$. Thus, in the following, we define a visibility threshold for DAD formation. In particular, we consider as DAD states the ones whose central antidark structure exhibits a visibility $\nu\geq 10\%$. For instance in the case of $N_{B}/N_{A}=1/3$, according to our simulations this condition is fulfilled only for quenches with $g_{AB}^{f}>0.1$ [see also Fig. 8 (c)].

The role of the trapping geometry is examined in Fig. 8 (d). Here, a hyperbolic-like increase of $N_{DAD}$ can be seen for a decreasing trapping frequency $\omega$, leading to more DAD structures being spontaneously generated for larger $N_{B}/N_{A}$ ratios. The observed increase in $N_{DAD}$ can be easily explained by the overall size of the system which for lower $\omega$’s increases providing in this way more space for solitons to form.

Finally, Table 3 illustrates that the same overall phenomenology regarding the controlled formation of multiple DAD solitons holds for even larger bosonic ensembles.

V Dynamics in a Quasi One-Dimensional Harmonic Trap

To further showcase that our above-described results regarding the spontaneous generation of solitonic structures can be detected in contemporary ultracold atom experiments we simulate the corresponding interaction dynamics in a quasi 1D geometry. In particular, we consider a binary BEC consisting of 87Rb atoms prepared in its hyperfine states $\ket{F=1,m_{F}=-1}$ (species A) and $\ket{F=2,m_{F}=1}$ (species B) and being confined in a quasi 1D harmonic trapping potential Egorov et al. (2013); Yan et al. (2012). The underlying dynamics of the system is governed by the following set of three-dimensional coupled GPE equations

[TABLE]

In this expression, the spatial coordinates are denoted by $\textbf{r}=(x,y,z)$, the species index $\sigma=A,B$ while $m_{\sigma}$ and $\Psi_{\sigma}(\textbf{r},t)$ denote the mass and the wavefunction of the $\sigma$-species respectively. Also, the intraspecies interaction strenghts $g_{\sigma\sigma}=2\pi\hbar a_{\sigma\sigma}/m_{\sigma}$ involve the scattering lengths $a_{AA}$ and $a_{BB}$, whereas the interspecies coupling is $g_{AB}=2\pi\hbar a_{AB}/m_{AB}$ where $m_{AB}=m_{A}m_{B}/(m_{A}+m_{B})$ is the reduced mass. Furthermore, we assume a three-dimensional harmonic oscillator potential, namely $V_{\sigma}(\textbf{r})=V_{\sigma}(x,y,z)=m_{\sigma}\omega^{2}_{\sigma}(x^{2}+\alpha_{\sigma}^{2}y^{2}+\lambda_{\sigma}^{2}z^{2})/2$ with $\omega_{\sigma}$ being the trapping frequency along the $x$ direction, while $\alpha_{\sigma}$ and $\lambda_{\sigma}$ are two anisotropy parameters. For our purposes, the atoms of both components possess the same mass (two different hyperfine states of the same atomic species) and are confined in the harmonic potential $V_{\sigma}(\textbf{r})$. As a consequence $m_{A}=m_{B}=m$, $\omega_{A}=\omega_{B}=\omega_{x}$, $\alpha_{A}=\alpha_{B}=\omega_{y}/\omega_{x}$, and $\lambda_{A}=\lambda_{B}=\omega_{z}/\omega_{x}$. The length and the energy scales of the system are expressed in terms of the harmonic oscillator length $a_{\rm osc}=\sqrt{\hbar/m\omega_{x}}$ and the energy quanta $\hbar\omega_{x}$ respectively. For the convenience in our numerical simulations the above set of coupled GPE equations [see Eq. (7)] is expressed into a dimensionless form by rescaling the spatial and time coordinates as $x^{\prime}=x/a_{\rm osc}$, $y^{\prime}=y/a_{\rm osc}$, $z^{\prime}=z/a_{\rm osc}$ and $t^{\prime}=\omega_{x}t$ repsectively. Additionally, the species wavefunction is $\Psi^{\prime}_{\sigma}(x^{\prime},y^{\prime},z^{\prime})=\sqrt{a_{\rm osc}^{3}/N_{\sigma}}\Psi_{\sigma}(x,y,z,t)$ and we also set $\hbar=1$. For simplicity below we omit the primes.

Moreover, due to the quasi 1D geometry of the external trap $\omega_{z},\omega_{y}\gg\omega_{x}$, the wavefunction of each species is decomposed as follows

[TABLE]

Here, $\phi_{\sigma}(y)$, $\phi_{\sigma}(z)$ refer to the normalized ground state wavefunctions in the $y$ and $z$ spatial direction respectively. Note also that for all the simulations to be presented below we use $N_{A}=5000$ and $N_{B}=3000$ atoms in the corresponding species and the values $a_{AA}=100.4a_{0}$, $a_{BB}=95.44a_{0}$ for the intraspecies scattering lengths Egorov et al. (2013). The latter correspond to the experimentally relevant scattering lengths of 87Rb, where $a_{0}$ is the Bohr radius. Additionally, for the harmonic oscillator potential we use the trapping frequencies $(\omega_{x},\omega_{y},\omega_{z})=2\pi\times(1.5,140,178)Hz$ as per the experiment of Yan et al. (2011). Notice that in this section we chose to present our findings in dimensional units in order to provide a direct connection with current state-of-the-art experimental settings. Our system is initially prepared in its ground state configuration being characterized by an initial interspecies scattering length $a_{AB}^{i}$. To trigger the dynamics an interaction quench is performed (as in Sections III, IV) towards $a_{AB}^{f}$ and we monitor the time-evolution of the binary mixture up to $600$ ms utilizing Eq. (7), see Fig. 9. Evidently, following an interspecies interaction quench within the miscible phase e.g. from $a_{AB}^{i}=20a_{0}$ towards $a_{AB}^{f}=80a_{0}$ both species undergo a breathing motion manifested by the expansion and contraction of their one-body densities [Figs. 9 (a), (b)]. However, for a quench from the miscible towards the immiscible regime of interactions the spontaneous generation of DB soliton structures is expected and indeed observed [Figs. 9 (c), (d)]. Notice that dark solitons build upon the majority component and the bright ones develop in the minority species following a quench from $a_{AB}^{i}=20a_{0}$ to $a_{AB}^{f}=150a_{0}$. Finally, upon considering an interspecies interaction quench from the immiscible to the miscible phase DAD solitonic structures appear in the corresponding one-body density evolution as a result of the quench-induced counterflow dynamics as shown in Figs. 9 (e), (f) for a quench from $a_{AB}^{i}=150a_{0}$ to $a_{AB}^{f}=75a_{0}$. Indeed, we observe that dark states build upon the majority species and antidark structures emerge in the minority species. From the above discussion we can conclude that the dynamical formation of the DB and the DAD solitons discussed in the preceding Sections III, IV appears also in quasi 1D setups and as such can be potentially observed in current state-of-the-art ultracold atom experiments. Indeed, the waveforms identified herein are known to be robust in quasi 1D geometries Yan et al. (2011). Exposure of these structures to higher-dimensional setups may lead to transversal instabilities and thus the formation of nonlinear excitations of a different kind such as vortex-bright ones in two spatial dimensions Law et al. (2010); Pola et al. (2012); Mukherjee et al. (2019) and vortex-line-bright solitons or vortex-ring-bright solitons in three dimensions respectively Charalampidis et al. (2016).

VI Conclusions

In the present work the interaction quench dynamics of a 1D harmonically confined particle imbalanced BEC mixture has been investigated within a mean-field theoretical approach. In particular upon considering quenches crossing the miscibility/immiscibility threshold in both directions the dynamical formation of multiple solitary waves in this non-integrable model has been unraveled. For miscible-to-immiscible transitions the unstable dynamics leads to the filamentation of the density of both species in line with our recent findings regarding particle balanced mixtures Mistakidis et al. (2018). However for the particle imbalanced scenario studied herein we were able to identify that via the filamentation process DB solitary waves emerge. The number of the latter is significantly reduced during the dynamics. Indeed, due to the non-integrable nature of the system, the emergent DB structures feature highly asymmetric collisions with mass redistribution Katsimiga et al. (2018a) between the solitary wave constituents. A controlled generation of DB solitons is achieved by considering distinct variations of the experimentally tractable system parameters such as the particle numbers, the interspecies repulsion coefficient, and the trapping frequency. Finally, and since both clouds undergo a collective breathing motion, we further evaluate the relevant breathing frequency of each species. We find that the two breathing frequencies are substantially different for strong particle imbalances. A deviation that is more dramatic for larger postquench repulsions Sartori and Recati (2013).

On the other hand, for immiscible-to-miscible transitions the spontaneous creation of multiple DAD solitons is observed linking to previous studies Mistakidis et al. (2018). Here, we extend our investigations shedding light on the underlying mechanism that leads to the formation of these states. It is showcased that these states emerge due to a counterflow process that takes place as the system is quenched towards miscibility. Additionally, the generated DAD solitons are found to be robust throughout the evolution. This result is to be contrasted to the observed decrease of the DB states formed in the previous transition, indicating that the breaking of integrability in this regime does not seem to affect the dynamics. Measuring the breathing frequency of each bosonic cloud for such a transition, we find that it differs from one component to the other independently of the postquench interaction strength or the particle number imbalance. The above behavior of the breathing frequency for this immiscible-to-miscible transition, is to the best of our knowledge, reported herein for the first time. Also here, the controlled creation of multiple DAD states is revealed by manipulating the particle numbers, the interspecies interaction, and the trapping frequency.

Finally, we should emphasize that in both quench processes the observed dynamics persists for large bosonic ensembles with particle numbers being of the order of $10^{4}$ suggesting a possible observation of our findings in current state-of-the-art experiments.

There are several possible extensions of the present work that can be considered in future endeavours. Of particular interest would be to examine the quench-induced dynamics of 1D spinor BECs Bersano et al. (2018). In such settings dark-dark-bright and dark-bright-bright soliton complexes are known to form Bersano et al. (2018); Nistazakis et al. (2008b). Therefore, attempting to control the generation of multiple such states upon quenching would be a natural generalization of our current effort. Also for this spinorial BEC environment one could further go beyond the mean-field approximation Katsimiga et al. (2017b, d, 2018b), addressing the fate of multiple solitonic entities in the presence of quantum fluctuations. Yet another fruitful direction would be to investigate the controlled dynamical formation of other nonlinear excitations such as the beating dark-dark solitons. It is known that these states can be experimentally created in miscible binary mixtures upon considering the fast counterflow of both components Hoefer et al. (2011). Finally, a generalization of the quench-induced dynamics of binary BEC mixtures to higher dimensions represents an intriguing perspective. It is well-known that e.g. in two-dimensions there exist vortex-bright soliton states Law et al. (2010); Pola et al. (2012); Mukherjee et al. (2019), while in three-dimensions structures consisting of vortex-line-bright and vortex-ring-bright solitons have been identified Charalampidis et al. (2016). It would be interesting to examine if such structures can be dynamically produced under the influence of an interspecies interaction quench.

Acknowledgements

S. I. M and G. C. K would like to thank P. G. Kevrekidis for fruitful discussions. S. I. M and G. C. K gratefully acknowledge K. Mukherjee for insightful discussions and assistance regarding the quasi one-dimensional simulations.