Generic transverse stability of kink structures in atomic and optical nonlinear media with competing attractive and repulsive interactions

Kink solitons (alias domain walls) are nonlinear excitations encountered in disparate disciplines ranging from optical Kartashov et al. (2006), and magnetic media Pérez-Junquera et al. (2008); Buijnsters et al. (2014), living cellular structures Ivancevic and Ivancevic (2013); Yomosa (1983); Caspi and Ben-Jacob (1999), folding protein chains Hu et al. (2011a, b) to atomic gases Tylutki et al. (2020); Katsimiga et al. (2023a) and cosmology Vilenkin and Shellard (2000). Their topological character has been recently unveiled in two-dimensional (2D) van der Waals materials Lee et al. (2023), opening up the possibility of robust computations Wang et al. (2022). Stabilization of multidimensional spatially localized states is a fundamental challenge of scientific interest in physics and beyond Kartashov et al. (2019). It is well-known that 1D topological (e.g. dark solitons) and non-topological (bright solitons) defects suffer from the so-called snake-instability Tikhonenko et al. (1996); Mamaev et al. (1996); Kivshar and Pelinovsky (2000); Anderson et al. (2001) and infrared catastrophe Bergé (1998); Fibich (2015), respectively, once embedded in 2D and three-dimensional (3D) geometries. Mechanisms of suppression of the ensuing instabilities of defects have also been discussed. These mainly consider nonuniform media, i.e., incorporating external trapping geometries Proukakis et al. (2004); Wen et al. (2013); Kevrekidis et al. (2004), or are accompanied by nonlocal interactions Skupin et al. (2006); Lin et al. (2008); Nath et al. (2008), or accounting for the active competition of higher-order nonlinearities Zeng and Zeng (2020).

Cold atoms are ideal quantum many-body simulators, i.e., platforms featuring remarkable tunability of system parameters, such as e.g. interactions and external traps Bloch et al. (2008); Lewenstein et al. (2012). In this context, ultradilute and incompressible self-bound states of matter, referred to as quantum droplets Luo et al. (2021); Mistakidis et al. (2023) were recently experimentally detected in Bose mixtures Cheiney et al. (2018); Semeghini et al. (2018); Cabrera et al. (2018); D’Errico et al. (2019); Burchianti et al. (2020) and dipolar gases Chomaz et al. (2023); Böttcher et al. (2020). Stabilization of these states is achieved through the counterbalance of attractive and repulsive interactions modeled by mean-field nonlinear couplings and quantum fluctuations. The latter are incorporated perturbatively through the famous Lee-Huang-Yang (LHY) Lee et al. (1957); Larsen (1963) correction whose sign and form depend on dimensionality Ilg et al. (2018) leading to an extended Gross-Pitaevskii equation (eGPE) Petrov (2015); Petrov and Astrakharchik (2016). Such environments sustain also different kinds of nonlinear excitations such as dark solitons Edmonds (2023); Katsimiga et al. (2023b); Saqlain et al. (2023), bubbles Katsimiga et al. (2023a); Edmonds (2023), kinks Kartashov et al. (2022); Tylutki et al. (2020); Katsimiga et al. (2023a) and vortices Li et al. (2018); Yoğurt et al. (2023).

From a complementary nonlinear optics-seeded perspective, the Cubic-Quintic (CQ) model is an extension of the famous nonlinear Schrödinger model Killip et al. (2017) constituing a universal mathematical toolbox involving higher-order nonlinearities. In the optical setting, CQ is utilized to study the propagation of electromagnetic waves Falcão Filho et al. (2013); Reyna et al. (2020) in photorefractive materials. In this context, the quintic term accounts for higher-order corrections to the nonlinear refractive index Akhmediev and Afanasjev (1995), i.e., incorporating susceptibilities up to fifth-order. Experimentally, it is often also relevant to include the nontrivial (dissipative) absorption effects. The broad applicability of CQ manifests by the fact that it has been deployed to describe various phenomena, for instance, liquid waveguides Falcão Filho et al. (2013), special types of glasses Boudebs et al. (2003), and colloids containing metallic nanoparticles Falcão-Filho et al. (2007). In atomic gases, the cubic (quintic) coupling is related to the two- (three-) body interactions Luckins and Van Gorder (2018); Cardoso et al. (2011); Bulgac (2002) (see also, e.g., Kolomeisky et al. (2000)).

In this Letter we explore the dynamical transverse stability of kinks in a wide range of models from nonlinear optics and atomic physics, such as the 2D and 3D eGPE settings, as well as the 2D CQ model. A key feature of these setups is the competition between attractive nonlinearities (dominant at low density) and repulsive ones (prevailing at high density) whose presence is proved crucial for the stability of kinks. Specifically, the stability and robustness of these entities is demonstrated in a twofold manner. Firstly, we analytically extract and numerically evaluate the relevant for each model Bogoliubov-de-Gennes (BdG) spectrum. Secondly, we subject the relevant one-dimensional (1D) kink solutions to external, spatially distributed fluctuations (customarily present in experimental settings) and subsequently monitor their remarkable persistence for evolution times of the order of seconds. This demonstrates the relevance of these 1D topological defects in current state-of-the-art cold atom and nonlinear optics experiments and their structural robustness as potential information carriers.

To highlight the universal features of the existence and stability of kink solutions in higher spatial dimensions, three different nonlinear models characterized by competing attractive (at low density) and repulsive (at high density) interactions are investigated. These correspond to the eGPEs describing, for instance, self-bound quantum droplets in 2D Petrov and Astrakharchik (2016); Luo et al. (2021) and 3D Petrov (2015); Luo et al. (2021), and the CQ equation. Experimentally, the droplet bearing systems can be emulated e.g. by considering two hyperfine states of ³⁹K Semeghini et al. (2018). In optics the CQ model is typically used to monitor, for instance, optical beam profiles of topological defects in liquid $\rm{CS}_{2}$ Reyna et al. (2020). In dimensionless form (see also supplemental material (SM) sup ) these equations read

Here, Eqs. (1a) and (1c) refer to the 2D and 3D eGPEs, while Eq. (1b) is the 2D CQ model with $g_{CQ}>0$. In the 2D eGPE, the logarithmic nonlinearity encompasses both mean-field and first-order quantum fluctuation effects Petrov and Astrakharchik (2016); Luo et al. (2021), while its 3D counterpart features a cubic (resp. quartic) attractive with $g_{1}<0$ (resp. repulsive) mean-field (LHY) contribution Petrov (2015), see also SM sup . Furthermore, in the CQ model there is an attractive (repulsive) cubic (quintic) nonlinear coupling commonly accounting for two- (three-) body interactions in the realm of cold gases Luckins and Van Gorder (2018); Bulgac (2002) or higher order corrections to the Kerr effect in nonlinear optics Akhmediev and Afanasjev (1995); Reyna et al. (2020).

To infer the existence of planar kinks (i.e., one dimensional ones embedded in higher dimensions), the following ansatz is employed, $\Psi(x,r_{\perp},t)=e^{-i\mu t}u(x)$. Here, $\mu$ denotes the chemical potential and $r_{\perp}=y$ ($r_{\perp}=(y,z)$) are the transverse coordinates of the wave function in 2D (3D). A uniform background is assumed along the transverse directions, and the $u(x)$ waveform is taken to be real without loss of generality. Introducing this ansatz into Eqs. (1a)-(1c) results in their reduced effective 1D time-independent standing wave analogues. These can be cast into Newtonian equations of motion, $d^{2}u(x)/dx^{2}=-dV(u)/du$, subjected to the relevant for each model effective potential $V(u)$. These turn out to be

Integration of the Newtonian equations of motion yields the effective energy, $E$, of the 1D reduced system Katsimiga et al. (2023b). While the resulting potentials feature multiple maxima, the frequency parameter $\mu$ can be generically tuned to render these maxima equi-energetic, thus enabling a genuine heteroclinic orbit, i.e., a kink, between $u=0$ and $u=u_{*}$ with $E=0$.

To systematically determine the finite kink backgrounds and chemical potentials, two conditions need to be satisfied, $V(u_{*})=0$ and $dV/du|_{u_{*}}=0$, namely, conditions for the maximum at $u=u_{*}$ being equi-energetic with the one at $u=0$. For the 2D eGPE (CQ) these values lead analytically to $\mu=-g/(2\sqrt{e})$, $u_{*}=e^{-1/4}$ ($\mu=-3/(16g_{CQ})$ and $u_{*}=\sqrt{3/(4g_{CQ})}$) [Fig. 1(a)]. The $u_{*}$ boundary for the 2D eGPE separates the parameter region of 2D droplets from that of 2D bubbles Luo et al. (2021); Petrov and Astrakharchik (2016), analogously to what is the case also in 1D Katsimiga et al. (2023a). For the 3D eGPE, the chemical potential and extremum $u_{*}$ are $\mu=25g_{1}^{3}/216$ and $u_{*}=-5g_{1}/6$ respectively [Fig. 1(b)]. Notice that while such kink-type states have been discussed for CQ settings previously (see, e.g., Jin-Liang et al. (2006)), they are unprecedented, to our knowledge, for the emergent atomic theme of 2D and 3D quantum droplet settings. Interestingly, the asymmetry of the effective potential $V_{3D}(u)$ [Fig. 1(b)] suggests the concurrent existence of a heretofore unexplored quantum droplet (for values of $u<0$) whose analysis is deferred for future studies.

Within the effective potential framework, unique kink solutions exist for each coupling strength when competing interactions exist, i.e., when $g_{CQ}\neq 0\neq g_{1}$. In what follows, without loss of generality $g$ and $g_{CQ}$ are set to unity in order to inspect the stability properties of these topological defects. On the other hand, the mean-field interaction $g_{1}$ in the 3D eGPE is routinely tunable in experiments via Feshbach resonances, see Semeghini et al. (2018) and SM sup .

According to the aforementioned effective potential description the 1D planar kink solution occurs for $\mu=-1/(2\sqrt{e})$, $\mu=-3/16$, and e.g. $\mu=-0.1$ where $g_{1}=-5/6$, within the 2D eGPE, 2D CQ and 3D eGPE model, respectively. The relevant effective 1D but also the higher dimensional stationary states are then obtained upon utilizing heteroclinic (tanh-shaped) initial guesses, bearing the correct asymptotics described above. In the CQ case, an analytical form exists as $u=\sqrt{\frac{3}{8}(1+\tanh(\frac{\sqrt{3}x}{2\sqrt{2}}))}$ Jin-Liang et al. (2006); Birnbaum and Malomed (2008). It is then numerically confirmed that such 1D planar waveforms satisfy the corresponding 2D and 3D time-independent version of Eqs. (1a)-(1c).

Then, in order to extract the stability properties of these configurations, small perturbations of the numerically exact (up to a prescribed tolerance) kink solutions are considered having the form $\Psi=\left(\Psi_{0}+\epsilon ae^{-i\omega t}+\epsilon b^{\star}e^{i\omega^{\star}t}\right)e^{-i\mu t}$. Here, $``\star"$ denotes complex conjugation, $a,b$ and $\omega$ refer to the corresponding eigenvectors, and eigenfrequencies, whilst $\epsilon\ll 1$ is a small perturbation parameter. Utilizing the above wave function ansatz and keeping terms of order $\mathcal{O}(\epsilon)$ one arrives at an eigenvalue problem that, irrespectively of the system under consideration [Eqs. (1a)-(1c)], can be written as

The model-dependent diagonal matrix elements are found to be $\hat{L}_{11}=-\frac{1}{2}\nabla^{2}-\mu+\Psi^{2}_{0}+2\Psi_{0}^{2}\rm{ln}(\Psi^{2}_{0})$, $\hat{L}_{11}=-\frac{1}{2}\nabla^{2}-\mu-2\Psi^{2}_{0}+3\Psi_{0}^{4}$, and $\hat{L}_{11}=-\frac{1}{2}\nabla^{2}-\mu+2g_{1}\Psi_{0}^{2}+\frac{5}{2}\Psi_{0}^{3}$, for the 2D eGPE, CQ, and 3D eGPE respectively. Their relevant off-diagonal counterparts read $\hat{L}_{12}=\Psi^{2}_{0}+\Psi^{2}_{0}\rm{ln}(\Psi^{2}_{0})$, $\hat{L}_{12}=-\Psi^{2}_{0}+2\Psi_{0}^{4}$, and $\hat{L}_{12}=g_{1}\Psi^{2}_{0}+\frac{3}{2}\Psi_{0}^{3}$, respectively. In all cases, $\hat{L}_{11}=-\hat{L}_{22}$, $\hat{L}_{12}=-\hat{L}_{21}$ and $\nabla^{2}$ denotes the $\mathcal{N}$-dimensional Laplace operator. Additionally, $\Psi_{0}$ designates the stationary, real kink state for each distinct model.

The stability analysis outcome for all three models investigated herein, is depicted in Figs. 2(a)-(c). In all scenarios, spectral stability of the kink state can be inferred from the absence of a finite imaginary part, $\Im(\omega)\sim 10^{-12}$, in the pertinent BdG spectrum. For all settings, the first hundred eigenvalues of the continuous quasi-1D BdG spectra are illustrated, while we note that we verified that their relevant higher dimensional analogues produce the same stability result. By quasi-1D BdG spectra, what is meant here is that the transverse direction is decomposed into Fourier modes (with the “quantization” thereof from a finite transverse computational domain accounted for) and the stability of (a large number of) effective 1D spectra for different transverse wavenumbers is confirmed.

To corroborate our BdG findings we next consider the dynamical evolution of the stationary kink solution for each of the models under study. Specifically, we aim not only to mimic unavoidable experimental imperfections but also to exemplify the robust nature of the kink structure in the presence of transverse perturbations. To that effect, the kink initial state is significantly perturbed via a random normal distribution generator, $\delta(x,r_{\perp})$. The ensuing wave function acquires the form $\psi(x,r_{\perp},t=0)=\Psi_{0}\left[1+\varepsilon\delta(x,r_{\perp})\right]$. In this context, $\delta(x,r_{\perp})$ is characterized by zero mean and unit variance, while $\varepsilon=\varepsilon_{R}+i\varepsilon_{I}$ such that also phase disturbances are encountered. In particular, the considered amplitudes of the perturbation range lie in the interval $[10\%,50\%]$ of the solution amplitude. Typical evolution times (in the dimensionless units adopted herein) are of the order of $\sim 500$ and $\sim 10^{3}$ for the 3D and 2D settings, respectively, confirming the longevity of these planar 1D defects upon their exposure to transverse excitations.

Evidently, despite the significantly excited nature of the initial state the kink remains intact preserving its shape throughout the evolution, while “shedding” some of the relevant “radiation” in the form of dispersive wavepackets, as it can be seen in Figs. 3(a)-(c). Particularly here, distinct plates correspond to density, $|\Psi(x,y,t)|^{2}$, contours in the $x-y$ plane taken at different times during the propagation of the perturbed kink configuration. For the 3D dynamics the $z$-direction is integrated out, while we note that the robustness of the kink configurations has also been verified for different values of the chemical potential and hence distinct interaction strengths.

A binary bosonic mixture in the ultracold regime, where interspecies attraction is balanced by intraspecies repulsion, can host self-bound states Petrov (2015); Cheiney et al. (2018); Semeghini et al. (2018). The two hyperfine levels $\ket{F=1,m_{F}=0}$, $\ket{F=1,m_{F}=-1}$ of ³⁹K in the vicinity of the $59\,\rm{G}$ intraspecies Feshbach resonance Chin et al. (2010); D’Errico et al. (2007) were proven to be suitable candidates satisfying the above criteria Cabrera et al. (2018); Semeghini et al. (2018). These two states are experimentally populated using radiofrequency spectroscopy Cabrera et al. (2018) and their dynamics can be adequately described by the eGPE models [Eqs. (1a), (1c)] in 2D and 3D. For such a setup with averaged mean-field scattering length $\delta a=-0.24\,a_{0}$, the propagation times of the kink defect in 3D ($\sim 500$ in dimensionless units) refer approximately to $800\,ms$ sup . A corresponding pancake geometry confining such a 2D mixture can be realized by utilizing a tight harmonic trap along one spatial dimension Hadzibabic and Dalibard (2011). For an experimentally relevant trapping frequency $\omega_{z}=2\pi\times 2\,\rm{kHz}$ Cabrera et al. (2018); Saint-Jalm et al. (2019) and $\delta a=-1.5\,a_{0}$ [see SM sup ], the relevant timescale ($\sim 1000$ in dimensionless units) for the kink evolution is $2\,$ secs. These long evolution times guarantee the resilience of the kink structure for its realization in current state-of-the-art experiments that suffer in addition by three-body recombination processes Semeghini et al. (2018).

On the other hand, the 2D CQ equation is used to describe, for instance, the propagation of electromagnetic waves in optical media incorporating higher-order nonlinear refractive indices Sutherland (2003). Characteristic materials featuring such properties are chalcogenide glasses Boudebs et al. (2003), such as liquid $\rm{CS}_{2}$ Falcão Filho et al. (2013); Reyna et al. (2020). In such setups (where optical vortices and their destabilization has been monitored Desyatnikov et al. (2005)), the time-evolution can be probed by measuring the transmission of the kink along the axial direction of the nonlinear medium. Considering a light beam of $900\,nm$ Falcão Filho et al. (2013), a propagation distance of $10^{3}$ corresponds to $0.14\,mm$ axial depth for the kink transmission. In the above-discussed setups, namely either cold bosonic mixtures or nonlinear optical media, kink defects can be imprinted in the bulk by means of the well-established technique of phase masking Denschlag et al. (2000); Pedaci et al. (2006) or utilizing digital micromirror device patterned optical traps Saint-Jalm et al. (2019); Navon et al. (2021); Gauthier et al. (2016).

The existence and corresponding spectral and dynamical stability of 1D kink solutions embedded in higher-dimensional gases or nonlinear optical media featuring competing attractive lower order and repulsive higher order interactions is unveiled. This is exemplarily demonstrated for three distinct, representative models referring to the 2D and 3D eGPE droplet bearing systems and the 2D CQ setting relevant in nonlinear optics, in order to establish the breadth and generality of results. We explicated the existence of these 1D topological defects upon analytically extracting, in each case, a quasi-1D effective potential picture, subsequently supplemented by their spectral stability inferred through linearization analysis. Additionally, monitoring their dynamical evolution in the presence of external fluctuations, with amplitudes up to 50% of the actual solution, revealed their profound robustness and resilience, for times up to few seconds.

These results suggest the use of kink defects as promising candidates for topological quantum computing. Indeed, contrary to the generically (in such nonlinear Schrödinger models) transversely unstable dark soliton states, which have been studied and reviewed extensively Kivshar and Pelinovsky (2000), kinks appear to defy the proneness to transverse modulations. Importantly, also, these results pave the way for significant further both physical and mathematical studies. On the latter side, a rigorous proof of such stability is a key open direction. On the other hand, from a physical perspective, our effective potential analysis suggests the existence of further lower-dimensional nonlinear excitations (including droplets, bubbles, etc.), whose transverse dynamics and potential dynamical features may be of further theoretical and experimental interest in its own right. Such studies are currently in progress and will be reported in future publications.

We acknowledge fruitful discussions with B. Malomed. This material is based upon work supported by the U.S. National Science Foundation under the awards PHY-2110030 and DMS-2204702 (PGK).