Nonlinearity managed vector solitons

The evolution of vector solitons under nonlinearity management is studied. The averaged over strong and rapid modulations in time of the inter-species interactions vector Gross-Pitaevskii equation (GPE) is derived. The averaging gives the appearance of the effective nonlinear quantum pressure depending on the population of the other component. Using this system of equations, the existence and stability of the vector solitons under the action of the strong nonlinearity management (NM) is investigated. Using a variational approach the parameters of NM vector solitons are found. The numerical simulations of the full time-dependent coupled GPE confirms the theoretical predictions.


Introduction
The properties of an atomic Bose-Einstein condensate (BEC) under the modulation in time of a scattering length have attracted large attention the last years [1,2].The existence of two-dimensional bright solitons, in one and two component attractive condensates [3,4,5,6], in media with competing nonlinearities [7,8,9], long lived Bloch oscillations of gap solitons [10], Faraday waves [11], and other interesting phenomena, has been predicted and observed.Recently it has been shown, that the strong and rapid periodic modulation in time of the scattering length (so called Feshbach resonance management) in a BEC can lead to the existence of stable compactons [12].
We are here interested in investigating if such structures also can exists in two-component BECs with equal inter-species interactions.In distinction from the scalar case, effects of Feshbach resonance management can then be on inter-and intra-species scattering lengths, leading to richer varieties of Faraday waves and localized states in such systems.
Mathematically the problem is reduced to the investigation of modulational instability and localized states in two coupled nonlinear Schrödinger equations (NLSE) with nonlinearity management of the self-and cross-phase modulation terms.
This paper is devoted to the investigation of these problems.First we consider the modulational instability of nonlinear plane waves under periodic modulations in time of the nonlinearity.The analysis will be performed for the case of rapid and strong modulations.We derive the averaged vector-GPE, which is called vector-NLSE in optics, with effectively nonlinear dispersion terms.Based on this averaged equation, we analyse the conditions of the existence of bright solitonic states and study the properties of these solutions.By the numerical simulations of the original vector-GPE with nonlinearity management, we verify the analytical predictions, including the averaged equation approach, and analyze the results further than the analytical predictions.

Model
Two-component BECs are here described by the system of coupled Gross-Pitaevskii equations: Here gij = 2 ω ⊥ a ij , where ω ⊥ is the transverse frequency of the trap, and a ij are the inter-and intra-species atomic scattering lengths.ψ 1,2 represent the wavefunctions of the individual components of a BEC.Considering the case of equal masses m = m 1 = m 2 , for the atoms of both components, and introducing the dimensionless variables according to: we obtain the system of a coupled GPE, without an external potential, that describe the vector solitons: With the coefficients g 12 = γ, we have the corresponding Hamiltonian: In the field of optics, when studying spatial solitons, t represents the propagation coordinate, while the u and v variables denote two mutually incoherent beams.We consider the case of rapidly changing time dependent coefficients, that is To obtain averaged equations, we use the following transformations [12,13,14,15], which removes the rapidly and strongly varying terms from the GPE, and allows to obtain the averaged over rapid modulations system of equations: where Γ t (t) = γ 1 (t), i.e., Γ(t) = γ 1 ω sin(ωt).By inserting the new field transformations (4) into Eqs.(2), we get the following equations Here and hereafter we omit the bar signs for simplicity.It follows from Eqs. (5) that After inserting Eqs.(6) into the transformed equations (5), and averaging over the period of rapid oscillations Λ = 2π/ω, we derive the following set of coupled averaged equations: where σ 2 = Γ(t) 2 = γ 2 1 /(2ω 2 ).Hence, the standard Hamiltonian (3) can now be written in the following averaged form: The averaged system (7) shows that the effective nonlinear quantum pressure appears in addition to the original linear quantum pressure.The magnitude of this correction depends on the population of the other component.The appearance of the nonlinear dispersive terms will lead to new effects in the dynamics of the matter waves.In particular to the existence of nonlinearity managed vector solitons.Standard vector solitons, i.e. in one component, exists due to the balance between the second-order linear dispersion (the quantum pressure) and the mean-field cubic nonlinearity.Now the balance between linear dispersion, nonlinear dispersion and the mean-field nonlinearity can give rise to the existence of NM vector solitons.
In the next sections we will consider the modulational instability in the NM system, which is important for defining the region of parameters where solitons can be generated.The existence and stability of NM vector solitons will also be considered.
To assess the accuracy of the averaged model, we conducted numerical simulations of both the original system (2) and the averaged equations (7).The simulation utilized the same initial condition that was obtained by numerically solving for the stationary states of Eqs.(7).An example of the comparative analysis of the stationary solutions' evolution is depicted in Fig. 1.There we used the parameters ω = 10π, γ 0 = −2, and γ 1 = 15, for which we illustrate the comparison.We have also successfully verified the consistency between the results obtained from the original and averaged models for other sets of parameters with moderate values of γ 1 .In all the following graphs, the numerical widths are calculated using the relation width = fwhm/ √ 8 ln 2, which represents the relationship between the actual width of a Gaussian beam and the full width at half maximum, denoted as fwhm.

Modulational Instability
In this section, we consider modulational instability of matter waves in two-component BEC under nonlinearity management.Analytical considerations will be performed using the averaged system and then compared to numerical simulations of the full time dependent model.To do so, we first consider constant amplitude solutions of Eqs. ( 7) of the form where A and B are constants, and q u , q v , ω u , and ω v satisfy the dispersion relations We now study the stability of the constant amplitude solutions against small perturbations.Hence, we look for solutions in the form where δA(x, t) and δB(x, t) are complex functions that represents small perturbations.Substituting Eqs.(11) into the averaged equations (7), taking into account the dispersion relations (10) and linearizing the resulting equations, we obtain the following two coupled equations: Decomposing the perturbations into real and imaginary parts and looking for solutions of these unknown functions in the form of plane waves exp(iΩt + ikx), we obtain the following dispersion relation: where and The modulational instability (MI) occurs when Ω 2 < 0 for any wavenumbers k.For small k we then get the following condition from Eq. ( 13): The growth rate G, also called the gain of the modulational instability, is related to the imaginary part of Ω, and is given by When considering small values of k, the growth rate ( 15) exhibits an absolute maximum, which corresponds to the maximum gain.By applying Fermat's theorem dG(k)/dk = 0, we can determine the maximum gain (G max ) and its corresponding wavenumber (k 0 ): From Eqs. ( 16) and (17), it becomes evident that as the parameter σ 2 increase, which is related to the modulation, both the maximum gain G max and the corresponding wavenumber k 0 decrease in accordance with the factor 1/ √ 1 + 4σ 2 A 2 B 2 .The curves plotted in Fig. 2 illustrates the relationship between the wave number of the perturbation k, and the gain, as defined by Eq. (15).We used the values A = 1.5 and B = 1, with three different modulational parameter values (σ 2 ).Studying Fig. 2, we can verify that the relationship between modulation strength and maximum gain adheres to the ratio 1.0196 : 0.9450 : 0.8450 = 1 + 9σ 2  1 : 1 + 9σ 2 2 : 1 + 9σ 2 3 , i.e., as the last factor in ( 16).We have performed numerical simulations of the original model ( 2) using the same parameters as in Fig. 2. We compare the analytically found growth rates, which are marked at different points with different styles (see the legend of Fig. 2), to growth rates obtained through numerical solutions of the original model (2).Specifically, we use perturbed plane waves with the frequencies (k) corresponding to the marked points in Fig. 2 as initial conditions u 0 and v 0 .In Fig. 3, it is shown that the analytical predictions corresponding to the (blue) square in Fig. 2 is confirmed by full simulations of the original model (2).For the rest of the marked points in Fig. 2, numerically calculated growth rates, corresponding to spatial frequencies k = 1.18 (red triangle), and 1.2 (red dot) are G = 0.85, and 0.98, respectively.Where k = 1.6 corresponds to the (blue) square depicted in Fig. 2. The inset shows the exponential growth rates of the amplitudes (≈ 0.8 for both A u and A v ) found from solving Eq. ( 2) numerically, which is very close to the analytical prediction (≈ 0.84, see (blue) square on the γ 1 = 6 curve in Fig. 2).Parameters used are: g 11 = 1.2, g 22 = 1.4,γ 0 = −2, γ 1 = 6, and ω = 10π.

Variational approach
Let us consider the existence and stability of solitons in the nonlinearity managed vector-GPE.The existence follows from the results of the existence of the modulational instability, considered in the Section 3. It is useful to employ the variational approach to Eqs. (7) for the description of the vector solitons [16].Eqs. ( 7) can be obtained by the following Lagrangian density We use trial Gaussian functions as ansatz functions with real time-dependent amplitudes A i (t), widths w i (t), chirps b i (t) and phases φ i (t) (i = 1, 2) A set of coupled ordinary differential equations for the eight ansatz parameters (A i , w i , b i , and φ i ) can be obtained from the following variational principle: where L denotes the averaged Lagrangian density, which is obtained by integrating the result of inserting the Gaussian functions given by ( 19) into the Lagrangian density (18).By performing the integration, we obtain where The equations for the Gaussian ansatz parameters The Euler-Lagrange equations provide conservation of norms (N u = const, and N v = const), for the phases (φ i ), and they result in the following system of equations for the widths (w i ) and chirps (b i ): respectively.The following coupled differential equations for the widths can be derived from Eqs. ( 22)-( 25): 2 ) 3 2 2 ) By identifying the fixed points of this system, one can determine the stationary waveforms.The symmetric configuration with w = w 1 = w 2 , N = N u = N v , and g = g 11 = g 22 , is more suitable for analytical analysis, and we now carry out the analysis for this simplified case.In this scenario, the system described by equations ( 26)-( 27) can be simplified into the single equation One can determine the width of the stationary localized state by solving the algebraic equation The real and positive root of this equation is which correspond to the possible value of w that yields a physically valid width for the state, and an amplitude A = N/( √ πw).In Fig. 4, comparisons of the stationary solution found by Eq. (30), and a numerically calculated solution is shown in the case of N = 3, g = 1.5, γ 0 = −2, γ 1 = 10, and ω = 10π.For experiments with 39K spin mixtures BEC, see for example [17,18], we have for the parameters γ 1 = 10, ω = 10π, the scattering lengths a 11 = a 22 ∼ 20a 0 , a 12 ∼ 200a 0 and ω ∼ 10 3 Hz when the trapping frequency is ω ⊥ ∼ 10 2 Hz.Dynamical evolution of the stationary solutions corresponding to Fig. 4 are presented in Fig. 5.We find that the solution becomes unstable over time, such type of instability for the scalar NM model has been observed in [19,20].For other values of the norms (N), comparisons of widths and amplitudes found by the variational and the numerical methods are shown in Fig. 6.The results shows good agreement for the variational and the numerical analysis.2).Frames c) and d) corresponds to the amplitudes and widths of the solitons, respectively.In both cases the stationary solutions depicted in Fig. 4 are used as initial conditions.Especially noteworthy is the observation in the comparison that the averaged equation ( 7) starts to deviate from the dynamics observed over long time-scales here.Parameters used are: g = 1.5, γ 0 = −2, γ 1 = 10, and ω = 10π.The norms were both taken as N = 3.
One can also approximate the bound states of the original model (2) using asymptotic approximations based on the solution of the averaged model (7), as discussed in [13].For the sake of analytical simplicity, we now consider the symmetric case, which transforms the averaged model into the following conventional NLSE with modified coefficients: To approximate the bound states, we employ a standard ansatz for the stationary solutions: u(x, t) = φ(x)e −iµt .This leads to an ordinary differential equation (ODE) for the stationary solution φ(x): The first integral of this equation is given by: Using the same procedure as described in [13], it is possible to determine the existence of stationary bright solitons φ(x) in the open quadrant where µ < 0 and g + γ 0 < 0. These solitons can be obtained by solving the equation: with an amplitude of In Fig. 6 we compare the above analytic amplitude with the variational and numerical results.

Extended numerical simulations
We have also examined how the dynamics of soliton solutions, obtained for specific parameters, respond to small changes in those parameters.As an illustration, Fig. 7 demonstrates the solution dynamics for the case of N = 2 when the selected parameters are slightly modified.with variations in the parameters.The red curve, derived from the variational analysis, see Eq. (30), reflects the soliton's amplitude and width under the parameters g = 1.5, γ 0 = −2, γ 1 = 10, and ω = 10π.The other curves demonstrates how altering one parameter, while keeping the rest unchanged, affects the width and amplitude: blue for a slight change in the intra-species interaction (g = 1.485), green for a minor change in the inter-species interaction (γ 0 = −2.02),and black for a shift in the modulation frequency (ω = 33).
In Fig. 7, the red curve represents the evolution of the amplitude and width of the soliton, as obtained from the variational analysis, see Eq. (30), with the following parameter values: g = 1.5, γ 0 = −2, γ 1 = 10, and ω = 10π.The other curves depict how the width and amplitude evolve when only one of these parameters is slightly modified, while keeping the remaining parameters constant.
Furthermore, the comparative Fig. 5 revealed a specific limit to the soliton's lifetime, beyond which the averaged equation diverges from the observed long time-scale dynamics.Extended numerical simulations, carried out over longer durations, revealed a strong dependence of the soliton's lifetime on both its norm (N), and on the modulation strength (σ 2 ). Figure 8 displays the relationship between the soliton's lifetime and the modulation parameter strength (σ 2 ) under constant norms (top two subfigures), and the relationship between the soliton's lifetime and the norm (N) with a fixed modulation strength (bottom subfigure).
We have numerically found specific threshold values for both the norm and the modulation strength.Below these threshold values, the lifetime remains independent of both the norm and the modulation strength.To illustrate, in the case of N = 1 (as shown in the top left subfigure in Fig. 8), this threshold value is σ 2 ≈ 0.146 (γ 1 = 17).In the case of N = 2 (as depicted in the top right subfigure in Fig. 8), this value is σ 2 ≈ 0.033 (γ 1 = 8).At last, in the case when the modulation strength is σ 2 = 0.18 (as indicated in the lower subfigure in Fig. 8), the threshold value for the norm is N ≈ 0.9.As an illustration, in Fig. 9, the behavior of the soliton is depicted when the parameters reach their threshold values.It is evident from this figure that the soliton remains intact over a long time-scale without breaking apart.

Conclusions
In conclusion, our investigation has provided several important results: The first result is the derivation of a system of equations obtained by averaging over rapid and strong modulations of the GPE.The second result reveals that the averaged system exhibits the presence of effective nonlinear quantum pressure, which is more pronounced compared to an unmodulated system.Furthermore, the process of modulation instability (MI) has been analytically investigated using the averaged system, and the results are confirmed by numerical simulations of the full time-dependent coupled GPE.Moreover, it is shown that the combined action of the linear and nonlinear dispersion from one side, and the mean-field nonlinearities from the other side, leads to the existence of nonlinearity managed (NM) vector solitons.Theoretical analysis based on the variational approach was performed, and the predictions of the theory regarding soliton parameters were confirmed through numerical simulations of the full system, including the sensitivity to variations in the parameter values, and the lifetime of the solitons.

Figure 1 :
Figure 1: Comparison of the averaged Eqs.(7) and the original model (2).The top frames display the temporal evolution of soliton components (u and v), with the averaged model shown as blue solid curves and the original model as red dashed curves.The middle frames show the amplitude and width of the u component, while the bottom frames depict the amplitude and width of the v component.In all frames, the dashed red curves correspond to the averaged model, while the solid blue curves represent the original model.Norms of the corresponding components are: N u = 1 and N v = 2. Parameters used are: g 11 = 1.4,g 22 = 1.1, γ 0 = −2, γ 1 = 15, and ω = 10π.

Figure 5 :
Figure 5: Dynamical evolution of the stationary solutions found by the variational approach.Frame a) represents the numerical solution of the averaged model of Eq. (7).Frame b) corresponds to the numerical solution of the original model of Eq. (2).Frames c) and d) corresponds to the amplitudes and widths of the solitons, respectively.In both cases the stationary solutions depicted in Fig.4are used as initial conditions.Especially noteworthy is the observation in the comparison that the averaged equation (7) starts to deviate from the dynamics observed over long time-scales here.Parameters used are: g = 1.5, γ 0 = −2, γ 1 = 10, and ω = 10π.The norms were both taken as N = 3.

Figure 7 :
Figure7: Evolution of amplitude (A) and width (w) of stationary solutions for N = 2 with variations in the parameters.The red curve, derived from the variational analysis, see Eq. (30), reflects the soliton's amplitude and width under the parameters g = 1.5, γ 0 = −2, γ 1 = 10, and ω = 10π.The other curves demonstrates how altering one parameter, while keeping the rest unchanged, affects the width and amplitude: blue for a slight change in the intra-species interaction (g = 1.485), green for a minor change in the inter-species interaction (γ 0 = −2.02),and black for a shift in the modulation frequency (ω = 33).

Figure 8 :
Figure 8: Relations of the soliton lifetime on the modulation strength (σ 2 ) under the norms N = 1, and N = 2 (upper two subfigures).The relation to the variation of the norm (N ) at σ 2 = 0.18 (lower subfigure).