Nonlinear management of topological solitons in a spin-orbit-coupled system

We consider possibilities to control dynamics of solitons of two types, maintained by the combination of cubic attraction and spin-orbit coupling (SOC) in a two-component system, namely, semi-dipoles (SDs) and mixed modes (MMs), by making the relative strength of the cross-attraction, gamma, a function of time periodically oscillating around the critical value, gamma = 1, which is an SD/MM stability boundary in the static system. The structure of SDs is represented by the combination of a fundamental soliton in one component and localized dipole mode in the other, while MMs combine fundamental and dipole terms in each component. Systematic numerical analysis reveals a finite bistability region for the SDs and MMs around gamma = 1, which does not exist in the absence of the periodic temporal modulation ("management"), as well as emergence of specific instability troughs and stability tongues for the solitons of both types, which may be explained as manifestations of resonances between the time-periodic modulation and intrinsic modes of the solitons. The system can be implemented in Bose-Einstein condensates, and emulated in nonlinear optical waveguides.


I. INTRODUCTION
An active direction in the current work with Bose-Einstein condensates (BECs) in atomic gases is using them as a testbed for emulation of various effects originating in condensed-matter physics, which may be reproduced in a clean and easy-to-control form in ultracold bosonic gases [1][2][3]. In particular, a binary gas, with a pseudo-spinor twocomponent wave function, may emulate spin-orbit coupling (SOC) in semiconductors, i.e., the interaction between the electron's spin and its motion across the underlying ionic lattice [4,5], as first demonstrated in Ref. [6], see also reviews [7]- [11]. While most experimental works on the BEC simulation of SOC dealt with effectively one-dimensional (1D) settings, implementation of SOC in the quasi-2D geometry was reported too [12], making it relevant to consider 2D (and 3D) systems coupled by the spin-orbit interaction. In this way, SOC opens a straightforward way to the creation of topological modes characterized by vorticity, because linear operators accounting for the coupling of two components in the corresponding system of Gross-Pitaevskii equations (GPEs), see Eq. (1) below, generate vorticity in one component if the other one is taken in the zero-vorticity form.
The SOC effect, being linear by itself, may be naturally combined with the intrinsic nonlinearity of bosonic gases, represented by cubic terms in the respective GPEs [13] (and/or by nonlocal cubic terms accounting for long-range interactions in BEC built of dipole atoms [14]). The interplay of SOC and nonlinearity makes it possible to predict a great variety of stable modes, including 1D and 2D solitons [15][16][17][18] and various nonlinear topological states in 2D, such as vortices and vortex lattices [19]- [28] and skyrmions [29]. In fact, the 2D and 3D SOC systems is one of the most prolific sources of nonlinear states with intrinsic topological structures.
A majority of works addressing nonlinear dynamics of SOC systems investigated the case of self-repulsion, which is relevant to the current experiments with 87 Rb [6]. Nevertheless, the interplay of SOC with intrinsic attraction is possible too. Theoretical considerations predict that the latter setting gives rise to 2D states with intrinsic topological structures and very unusual dynamical properties: until recently, it was commonly assumed that any 2D model with cubic self-attraction may only generate unstable self-trapped states, such as Townes solitons [30] and their vortical counterparts [31,32]. The fundamental (zero-vorticity) Townes solitons are destabilized by the critical collapse (or by the supercritical collapse in 3D), while vortex solitons are subject to a still stronger splitting instability [35,36]. A new paradigm was revealed by the analysis of the 2D SOC system with cubic self-and cross-attractive interactions [37][38][39]: the linear SOC terms lift the specific conformal invariance of the cubic GPEs in 2D, which is responsible for the instability, as it makes norms of the entire soliton family degenerate, allowing them to take a single valueexactly the critical one which launches the 2D collapse. As a result of lifting the degeneracy, the norm of the 2D solitons falls below the critical value, thus protecting them against the onset of the collapse [37]. The specific form of SOC creates two distinct species of topological solitons in this case, namely, semi-vortices (SVs), which combine vorticities S = 0 and S = ±1 in the two components, and mixed modes (MMs), which juxtapose zero-vorticity and vortex terms in both components [37]. Their stability is determined by relative strength γ of the cross-attraction between the components and self-attraction (or the XPM/SPM (cross/self-phase modulation) ratio, in terms of optics [40]): the SVs and MMs are stable (actually, realizing the system's ground state) at γ ≤ 1 and γ ≥ 1, respectively. The stability boundary shifts to γ > 1 under the action of the Zeeman splitting [38]. On the other hand, SVs and MMs are stable at all values of γ in a model where the self-trapping is provided not by attractive interactions, but by repulsion, with the local strength growing fast enough from the center to periphery [41,42]; moreover, even excited states of SVs and MMs, which are completely unstable in the case of the self-attraction [37], are partly stable in the latter case.
SOC implemented in 2D BEC as the emulation of the solid-state phenomenology may, in turn, be emulated in optical media, in terms of the spatiotemporal propagation in dual-core planar waveguides, which simulates the pseudo-spinor (two-component) structure of the wave field, while the SOC proper is simulated by temporal [43] or spatial [44] shift of the linear coupling between the parallel waveguiding cores. The former possibility is provided by the known effect of temporal dispersion of the coupling [45], while the latter scheme may be supported by a skewed structure of the medium separating the cores [44]. Combining these settings with the natural Kerr self-focusing in the dielectric material opens the way to predict stable 2D optical solitons with an intrinsic topological structure [43,44], which may be construed as counterparts of the above-mentioned SVs and MMs.
A known method which makes it possible to additionally stabilize soliton modes which are, otherwise, vulnerable to instabilities, is the nonlinearity management, i.e., periodic modulation of the self-focusing strength in time, in the case of BEC [46][47][48][49][50][51], or along the propagation distance in optics [52,53]. In the 2D SOC system, the application of the management technique is an especially interesting possibility, as, by means of the Feshbach resonance controlled by a low-frequency ac magnetic field [54], one can apply time-periodic modulation to the above-mentioned XPM/SPM ratio γ. As a result, the system will periodically pass from the SV-stability region, γ < 1, to the MM-stability one, γ > 1. Such a possibility poses the problem of the existence and stability of the two species of the 2D solitons in such dynamical states. The situation is somewhat similar to the earlier studied situation with the Feshbach-resonance management periodically alternating the sign of the nonlinearity in the single-component case, which drives periodic transformations between regular and gap-type solitons [55].
The present work addresses this dynamical problem by means of simulations of GPEs including the SOC terms and the time-modulated coefficient γ. However, performing such systematic simulations in the 2D model with many control parameters is a challenging numerical problem. On the other hand, it may be efficiently emulated by the similar 1D system, where a counterpart of the SV is a semi-dipole (SD), with a fundamental (spatially even) structure in one component, and a dipole structure (a spatially odd localized state with zero at the central point) in the other (see, e.g., Ref. [17]). MM states are possible in 1D as well, with the fundamental and dipole terms mixed in both components. A crucially important property of the 1D system, which suggests to use it for the emulation of the 2D prototype which is critically sensitive to the sign of γ − 1, is the fact that, exactly like in 2D, the one-dimensional SDs and MMs are stable, respectively, at γ < 1 and γ > 1, and the nonlinearity-controlling techniques, such as the Feshbach resonance, apply even easier in the 1D settings.
The rest of the paper is organized as follows. The model is introduced in Section II, along with some analytical results which help to illustrate the structure of static SD and MM solitons. Results of the numerical analysis of the management model are collected in Section III. In particular, a region of the SD-MM bistability is found, and qualitative explanations are presented for specific dynamical features (instability troughs and stability tongues) induced by the management. The paper is concluded by Section IV.

II. THE MODEL AND ANALYTICAL RESULTS
In scaled units, the system of GPEs for the 1D binary BEC under the action of the Rashba-type SOC and attractive nonlinearity, takes the form of [38,39] where φ + and φ − are two components of the pseudo-spinor wave function, λ is the SOC strength, self-attraction coefficients are also scaled to be 1, and γ is the above-mentioned relative strength of the nonlinear cross-interaction. By means of scaling transformation, which does not affect γ, we further fix λ = 1 in Eq. (1) (λ is kept as a free parameter in analytical expressions given below by Eqs. (7)- (12), to make the structure of those expressions clearer). Of course, transformation (2) cannot be applied in the absence of SOC, λ = 0. Stationary solutions to Eq. (1) of the SD type, with chemical potential µ < 0, have the form of where f even,odd (x) are even and odd functions, respectively (see, e.g., Eq. (12) below). For the MMs, the appropriate ansatz is Stationary modes are characterized by their norm, As mentioned above, for constant γ = γ 0 , SDs and MMs are stable, severally, at γ 0 < 1 and γ 0 > 1 (in Ref. [37], the latter result was originally obtained for SOC of the Rashba type, but it was then established that the same is true for a general Rashba-Dresselhaus combination [38]). At γ 0 = 1, both SDs and MMs are marginally stable.
We introduce the nonlinearity management by making γ a periodic function of time, the corresponding frequency, ω, being a free parameter of the management: Thus, in our model, with λ = 1 fixed by the scaling, there are four free parameters: N , γ 0 , γ 1 and ω. We report systematic numerical results for N = 3, checking that this value adequately represents the generic case. Then, there remain three free parameters to vary: γ 0 , γ 1 and ω in Eq. (6). While the governing equations are cast here in the scaled form, the rescaling can be undone to estimate characteristic values of control parameters in physical units. In particular, in terms of the BEC realization for light atoms, such as 7 Li, and the assuming the size of the soliton ∼ 3 µm, the time unit in the scaled equations corresponds to t 0 ∼ 10 ms. Then, characteristic scaled values of the management frequencies, ω ∼ 1 (see below) correspond, roughly, to ω ∼ 2π × 20 Hz.
A more specific situation corresponds to the limit of N ≪ 1, i.e., broad small-amplitude solitons filled by the striped phase [56], which is represented below by factors cos (λx) and sin (λx) in Eqs. (7) and (8). In the absence of the management, they take the following approximate forms, for the SD and MM species, respectively (here, for the clarity's sake, we keep the SOC strength, λ, as a free parameter, rather than fixing λ = 1): the chemical potential being µ ≈ −(1/2) λ 2 + (3 + γ) 2 (N/8) 2 for both species (µ ≤ −λ 2 /2 is the semi-infinite spectral gap which may be populated by soliton states).
Another specific case is one corresponding to large N , i.e., narrow solitons. In particular, the respective SD state has a large fundamental component, while a relatively small dipole one, with y ≡ N x/2 and real function u(y) determined by a linearized equation, It is worthy to note that Eq. (11) admits an exact solution precisely in the case of γ = 1, when both SD and MM are stable: The application of the nonlinearity management to these specific cases should be considered elsewhere. Simulations of Eq. (1) were run starting with the initial state built as a stable stationary soliton of the SD or MM type (the ground state), produced by means of the imaginary-time-integration method applied to Eq. (1), with γ 1 = 0 in Eq. (6). Figures 1(a) and (b) display typical profiles of the corresponding static MM and SD states, obtained at γ 0 = 1.1 and γ 0 = 0.9, respectively. In the simulations of the full management model, with γ 1 = 0 in Eq. (6) for various values of ω, stable solitons were identified as those which keep their integrity and initial structure (SD or MM) in the course of long real-time simulations until t = 1000. Note that, in terms of the above-mentioned estimates for physical parameters of the BEC setting, this corresponds to times 10 s, which definitely covers the range of times that may be realized in any experiment.

III. RESULTS: STABILITY REGIONS FOR SOLITONS UNDER THE ACTION OF THE MANAGEMENT
Systematic simulations of Eq. (1), with γ taken as per Eq. (6), are summarized in Fig. 2, which displays stability regions in the parameter plane of (γ 0 , ω) for the MM-and SD-type solitons, in panels (a) and (b), respectively, with a fixed value of the management amplitude, γ 1 = 0.05 in Eq. (6). One conclusion is that the application of the management somewhat expands the MM and SD stability areas, from the above-mentioned half-planes in the absence of the management (γ 1 = 0), i.e., γ 0 ≥ 1 and γ 0 ≤ 1, respectively, to values which may be smaller than γ 0 = 1 for the MMs, and larger than 1 for the SDs. The expansion gives rise to an MM-SD bistability region, which is presented in detail at the end of this section.
Conspicuous features revealed by Fig. 2 are instability troughs created by the management in the originally stable half-planes. These features may be understood as manifestations of resonant interaction of the ac drive, supplied by the management, and an eigenmode of intrinsic oscillations of stable solitons, existing in the absence of the management. To consider this possibility, the MM's intrinsic eigenmode and respective eigenfrequency were identified by direct real-time simulations of Eq. (1), with γ 1 = 0 in Eq. (6), adding a small initial perturbation to the numerically exact MM soliton, as follows: where {φ 0+ , φ 0− } are the components of the MM state displayed in Fig. 1(a) (i.e., with γ 0 = 1.1). For the perturbation strength δ = 0.01 in Eq. (13), the ensuing evolution of an essential characteristic of the perturbed soliton, which we define as the share of the total norm staying in one component, is presented in Fig. 3(a). The evolution of R(t) clearly exhibits eigenfrequency ω 0 ≈ 0.34 of the MM's intrinsic mode. The most conspicuous feature in Fig. 2(a) is a relatively wide diagonal instability trough. In particular, at γ 0 = 1.1 its size is 0.585 < ω < 0.785. (15) To consider a possible explanation of this feature in terms of the resonance, in Fig. 3(b) we display a typical example of the instability for γ 1 = 0.05 at the point taken in the middle of interval (15), The instability is shown by means of the corresponding time dependence for R(t) and the MM's center-of-mass position, together with the modulation format, γ(t). The initial condition is taken as Eq. (13) with δ = 0.001. The instability manifests itself, in Fig. 3(b), by oscillations of R(t) and x with a growing amplitude (the spontaneously emerging oscillatory motion of unstable solitons is illustrated by an example displayed below in Fig. 5 (d)). Figure 3(b) shows the initial stage of the perturbation growth. At larger t, R(t) exhibits irregular oscillations with an amplitude ≃ 0.22 around R = 0.5. The situation observed in Fig. 3(b) is a typical picture of mechanical instability caused by the parametric resonance with the frequency ratio 1 : 2 [57], in agreement with relation (16). Further, Fig. 3(c) shows the same dynamical characteristics, R(t), x and γ(t) for γ 0 = 1.1, γ 1 = 0.05, and ω = 0.335, which is a point belonging to the second (much more narrow) instability trough in Fig. 2(a). In this case, the instability is again manifested by the growth of R(t) and x , although the instability is weaker than in Fig. 3(b). Comparing the current value of ω with ω 0 (see Eq. (16)), we conclude that this instability may be interpreted as caused by a direct resonance, with frequency ratio 1 : 1. A plausible explanation of additional small "notches" in Fig.  2(a) is the presence of very weak higher-order (subharmonic) resonances. Similarly, the instability-trough pattern observed in Fig. 2(a) can be construed as manifestations of the parametric, direct, and subharmonic resonances with an intrinsic mode of the SD soliton. Another relevant summary of the stability results is presented in Fig. 4 for the solitons of the MM (a) and SD (b) types in parameter plane (γ 1 , ω), fixing the value of the constant term in the nonlinearity coefficient (6), viz., γ 0 = 1.05 > 1 in (a), and γ 0 = 0.95 < 1 in (b). The solitons are stable at sufficiently large ω, where the high-frequency modulation is effectively averaged out, hence it does not produce a conspicuous effect. It is also natural that stability areas tend to shrink as γ 1 increases. Nevertheless, narrow stability tongues are found too. For example, the SD soliton is stable at 0.265 < ω < 0.295 for γ 1 = 0.12. Examples of the SD dynamical states, taken at the same fixed value of the management amplitude, γ 1 = 0.12, are shown in Fig. 5 (c) ω = 0.25 (unstable). The stable SD in panel (b) belongs to the narrow "tongue" in Fig. 4(b). It was found that the instability region above the "tongue" is similar to the main instability trough shown in Fig. 2(a). While accurate explanation of the nature of the stability tongues requires a more detailed analysis, it is plausible that they also originate from a nonlinear resonance, which, as it is known, gives rise to both unstable and stable solution branches [57].
Lastly, as mentioned above, a noteworthy feature observed in Fig. 2 is bistable coexistence of the SD and MM states in a small but finite region near γ 0 = 1 (in the absence of the management, γ 1 = 0, the bistability is only possible strictly at γ 0 = 1 [37]). Figure 6 shows stability boundaries for the SD solitons (solid lines) and MMs (dashed lines) in the parameter space of (γ 0 , ω) (the same plane which is displayed in Fig. 2), at three fixed values of the management amplitude: (a) γ 1 = 0.05; (b) γ 1 = 0.1; (c) γ 1 = 0.15. The SD and MM solitons are stable, severally, in regions bounded by solid lines and dashed ones. Accordingly, the bistability takes place in finite areas between the dashed and solid lines, which are designated by vertical shading in Fig. 6. The bistability region is small, but its size increases with the growth of γ 1 .

IV. CONCLUSION
In this work, we address the possibility to control stability and dynamics of two types of two-component spin-orbitcoupled solitons, SDs (semi-dipoles) and MMs (mixed-modes), which represent distinct species of 1D topological modes, by means of the nonlinearity management, i.e., making the relative strength of the cross-attraction, γ, a function of time which periodically oscillates around γ = 1. This value is the boundary between stability regions of the static SDs and MMs. By means of systematic simulations, we have found a finite bistability region around γ = 1, which expands with the increase of the management amplitude. In the usual SDs' and MMs' stability domains (γ < 1 and γ > 1, respectively), the analysis reveals new features generated by the management, in the form of long instability troughs and, on the other hand, stability tongues penetrating into instability domains. These features may be explained as manifestations of resonances between the time-periodic management and excitation modes of stable 2D solitons, that exist in the absence of the management.
As said above, the 1D soliton species of the SD and MM types emulate their 2D spin-orbit-coupled counterparts, in the form of SVs (semi-vortices) and two-dimensional MMs, respectively. It may be quite interesting, although somewhat challenging, in terms of collecting systematical numerical data, to identify stability charts for the 2D topological modes of both types under the action of the nonlinearity management. Another relevant extension may be realization of the nonlinearity management for the two-component system combining SOC and long-range dipoledipole interactions [58]. It may be realized by periodically varying the direction of the external magnetic field which determines the orientation of the atomic magnetic moments.