Solitons and Josephson-type oscillations in Bose-Einstein condensates with spin-orbit coupling and time-varying Raman frequency

# Solitons and Josephson-type oscillations in Bose-Einstein condensates with spin-orbit coupling and time-varying Raman frequency

F. Kh. Abdullaev, M. Brtka, A. Gammal, Lauro Tomio Physical-Technical Institute, Uzbekistan Academy of Sciences, Tashkent, Bodomzor yuli, 2-b, Uzbekistan CCNH e CMCC, Universidade Federal do ABC, 09210-170, Santo André, Brazil Instituto de Física, Universidade de São Paulo, 05508-090, São Paulo, Brazil Instituto Tecnológico de Aeronáutica, CTA, 12228-900, São José dos Campos, Brazil Instituto de Física Teórica, Universidade Estadual Paulista, 01140-070, São Paulo, Brazil
July 15, 2019
###### Abstract

The existence and dynamics of solitons in quasi-one-dimensional Bose-Einstein condensates (BEC) with spin-orbit coupling (SOC) and attractive two-body interactions are described for two coupled atomic pseudo-spin components with slowly and rapidly varying time-dependent Raman frequency. By varying the Raman frequency linearly in time, it was shown that ordinary nonlinear Schrödinger-type bright solitons can be converted to striped bright solitons and vice versa. The internal Josephson oscillations between atom-number of the coupled soliton components, and the corresponding center-of-mass motion, are studied for different parameter configurations. In this case, a mechanism to control the soliton parameters is proposed by considering parametric resonances, which can emerge when using time-varying Raman frequencies. Full numerical simulations confirm variational analysis predictions when applied to the region where regular solitons are expected. In the limit of high frequencies, the system is described by a time-averaged Gross-Pitaevskii formalism with renormalized nonlinear and SOC parameters and modified phase-dependent nonlinearities. By comparing full-numerical simulations with averaged results, we have also studied the lower limits for the frequency of Raman oscillations in order to obtain stable soliton solutions.

###### pacs:
42.65.-k, 42.65.Sf, 42.81.Dp

## I Introduction

A progressively growing interest in the physics of Bose-Einstein condensate (BEC) with spin-orbit coupling (SOC) has been observed in recent years 2015Zhai (); 2011LJS (); LC2009 (); 2011Ho (); 2017Li (). For the coupling, different forms of Rashba BR84 () and Dresselhaus Dresselhaus55 (), as well as mixture of them, have been realized 2013Galitskii (). Important forms for the nonlinear excitations have been verified with structure of stable solitons, for condensed systems having attractive and repulsive interactions. In this regard, we can mention that the existence of solitons in BEC with SOC was investigated in Ref. Merkl10 (), for the case that we have repulsive interactions between atoms; and, in Refs. ZPZ13 (); AFKP (); DKM (); SM13 (), when the interactions are attractive. Gap solitons are predicted in Ref. KKA-PRL13 () for BEC with SOC in a spatially periodic Zeeman field, corresponding to a linear optical lattice. Experimentally, it is not an easy task to control the SOC parameter, having recent suggestions to tune it by applying rapid time variations of the Raman frequency 2013Zhang (); Salerno2016 (). The experimental observations reported in Ref. 2015Spielman () are shown that the spin-orbit coupling can be tuned in this way. In principle, the periodic variation in time of the condensate parameters can lead to new phenomena as the generation of new quantum phases Struck (), artificial gauge fields 2016FE (), compacton matter waves AKS (), etc. Therefore, actually it is quite relevant and of interest to investigate how the periodic time variations of the Raman frequency can affect the nonlinear modes of the condensate, such as solitons and vortices. In the limit of high frequencies, as shown in Ref. 2013Zhang (), the averaged Hamiltonian contains the nontrivial renormalization of the spin-orbit coupling, as well as the new effective nonlinear phase, which is sensitive to the interaction strengths SKMA ().

Our main task in the present work is to investigate the dynamics of solitons and Josephson-type oscillations between solitonic components, considering BEC with tunable spin-orbit coupling, with attractive interactions between the atoms, under slow and rapid time modulations of the Raman frequency. In this regard, related to Josephson oscillations in BEC, we should mention some previous studies in Refs. albiez2005 (); levy2007 (); abbarchi2013 (); 2014zou (). In particular, when considering the Raman frequency modulated in time, together with changes in other parameters of the system, one should expect to observe parametric resonance phenomena to occur in the internal Josephson oscillations, which have been introduced between the atom-number fraction existent in each of the components of the condensate with SOC. The parametric resonances in this case are introduced by the time-dependence of the Raman frequency and its corresponding relation with spin-orbit coupling of the two components of the condensate. Such study can be useful for possible experimental investigations, which can help the control of BEC parameters through resonance phenomena observations. We should point out that previous studies on parametric resonances in BEC are mainly concerned with time variations in trap configurations, optical lattices, as well as nonlinear parameters, looking for direct interference effects manifested in the densities 1999ripoll (); 2000kevrekidis (); 2002salasnich (); 2005tozzo (); 2014cairncross (); 2016abdullaev (). In our present study, by introducing a time modulation in the Raman frequency, the main focus is the oscillatory behavior between the internal atom-number population of the two components, during the time evolution of the condensate.

We start our study by first considering the case that we have defined spin-orbit parameter and constant Raman frequency, in order to verify the characteristics of the existent soliton solutions, which can be regular or striped solitons. Next, we introduce an adiabatic time modulation in the Raman frequency (growing and decreasing linearly), such that we can study how to switch between different kind of solitons. We follow our study by considering the Josephson oscillations between the components for both the cases that we have regular and striped soliton solutions. Parametric resonance effects in the oscillations are verified by considering periodically time variation of the Raman frequency at some given SOC parameters. The case of rapid time modulations of the Raman frequency can be treated by using a time-averaging approach, which is implemented over the time-dependent coupled system, implying in a renormalization of the SOC and nonlinear parameters. In this way, the interactions are effectively time independent. This case is being discussed in the final part (section IV) of the present work.

We are considering exact numerical simulations in all the cases, complemented by theoretical analysis, using variational approaches, whenever simplified solutions can be performed. As shown, the predictions derived by using the variational approach (VA) are verified to be fully consistent with the given numerical results in the region where regular bright solitons are obtained. In other cases, where the solutions are striped ones, demanding more parameters in the ansatz, the VA is quite helpful to indicate regions of stability, as well as initial starting profiles for the full numerical computation. By using the multi-scale expansion method for the averaged system, solitonic solutions are also found and confirmed by our numerical simulations of the full coupled system with time-dependent Raman frequency.

In the next, the basic formalism of the model is being presented in section II. For reference to other sections, we add two subsections: The first (A) where we provide some details on the linear energy dispersion, which is defining two regions for the kind of soliton solutions that we can obtain (regular or striped). In the subsection B, we add already some results to exemplify the two possible solutions and how to transform solitons between the two regions by considering adiabatic linear time variation of the Raman frequency. In section III, by considering Raman frequency modulated in time and Josephson oscillations, we analyze the possibility to obtain resonant responses. The case with Raman having high frequency modulations is presented in section IV. Finally, in section IV, we present our main conclusion.

## Ii Model formalism

In our approach, we consider a spin-orbit coupled BEC with equal Rashba and Dresselhaus contributions for the spin-orbit coupling terms, as in Ref. AFKP (), which can be described by a one-dimensional (1D) coupled equation for the two pseudo-spin components. For that, let us consider an harmonic trap where the frequency along one direction, , is much smaller than the frequency in the perpendicular direction, . In this case, given the units of energy, length and time, respectively, by , and (where is the mass of both atomic components), we can write in dimensionless units the corresponding SOC formalism for the two pseudo-spin components, and , of the total wave function

 ψ≡ψ(x,t)≡(uv). (1)

The corresponding matrix-formatted non-linear Schrödinger type coupled equation can be written as

 i∂∂t(uv) = [−12∂2∂x2−ikLσz∂∂x+Vtr+Ω(t)σx](uv) (2) − (|u|2+β|v|200β|u|2+γ|v|2)(uv),

where is the trap potential, assumed to be zero () in the present study. are the usual Pauli matrices, with being the strength of the spin-orbit coupling and the time-dependent Raman frequency (also given in units of the trap frequency ). In the non-linear terms we have the dimensionless parameters and , which are given by the ratio between the two-body scattering lengths, , of the two atomic components, such that . In the present work, we are going to consider attractive two-body interactions, such that we have an overall minus signal for the non-linear interaction. From the symmetry of the coupled equations (2) for , we can extract a simple relation between the two components and : Let us consider that, for a given parameter we are identifying the solution for by . In this case, a particular solution decoupling the equations can be verified with .

### ii.1 Linear energy spectrum

For a constant Raman frequency , the linear energy spectrum can be derived by considering a plane wave-function with momentum , , which will give us the following dispersion relation:

 w±(k)=12k2±√k2Lk2+Ω20. (3)

This relation, also shown in Ref. AFKP (), is plotted as a function of in Fig. 1, where two different regions (I and II) can be distinguished, according to the choice of parameters we have for the spin-orbit coupling and Raman frequency .

In region I, which happens when , we can only obtain two single minima in the dispersion relation; when . For , we are in region II, where has just one minimum (at , as in the case of region I). However, in this case, the dispersion relation for presents a local maximum at with two minima at , where ; both with .

As already verified in Ref. AFKP (), bright soliton solutions of the NLSE are obtained in region I. However, for (region II), the solutions are striped-type bright solitons. In Fig. 1, the two regions are represented, exemplified for particular values of , (region I) and (region II). The bright solitonic solutions can be found by the multiple scale analysis, which will be used in section IV to investigate the soliton dynamics under rapid modulations of parameters. Here, we should remark that observation of striped phases have been recently reported in Ref. 2017Li () for BEC with SOC. In the following subsection, we are exemplifying the kind of soliton solutions we obtain in both regions, using constant and linearly time-varying Raman frequencies.

### ii.2 Regular and striped soliton interchanged by adiabatic time variation of Ω(t)

With adiabatic time variation of the Raman frequency, , we can transform solitons from one region to another region, for a given fixed . By considering a regular soliton obtained in region I for , it can be transformed to a striped soliton by decreasing ; or, the other way, if started in region II. For that, let us consider a variation of the form

 Ω(t)=Ω0(1±Δt). (4)

Recently such variation has been used in numerical simulations for dark soliton generations in BEC with SOC josa17 (). The transition of a soliton solution obtained in region II (striped soliton) to the region I (regular soliton), and back from region I to region II, is illustrated in Fig. 2, by considering several panels, obtained for different values of a the time-dependent Raman frequency, such that we have an adiabatic transition. In terms of the step-function (=0 or 1, respectively for or ), we can write the time varying Raman frequency as . In these cases, the results are obtained with the same values of the nonlinear parameters, which are related to the scattering length ratios, , implying that the inter- and intra-species scattering lengths remain the same.

In all the cases considered in the present work, we are using exact full-numerical solutions of the coupled Gross-Pitaevskii (GP) formalism (2), by applying an imaginary-timer propagation method brtka (), using the Crank-Nicolson algorithm, followed by real-time evolution of the soliton profiles. The time steps, as well as the total space interval and corresponding discretization, have been adapted to obtain convergent and accurate results. In most of the cases, considering our dimensionless units, we found the time step to be sufficient in the imaginary-time relaxation procedure, and in the real-time evolution. In this regard, we can mention that particular care has to be considered in the time-evolution of striped solitons, where stable solutions demand a large enough number of grid points within a large interval (to avoid border effects). In view of that, for some results we decrease the time step to . Along the text, we are providing some theoretical analysis by using variational ansatzes, which are verified to be more efficient in region I, where the solutions are regular solitons. For the case that rapid modulations are used for the Raman frequency, in region II, it was also shown to be quite useful to employ a multi-scale expansion for the averaged coupled system.

## Iii Time-modulated Raman frequency and internal Josephson oscillations

An interesting case that can be explored is the influence of time-varying Raman frequency on the oscillations in atomic populations, which can occur between the components of the soliton solutions (the internal Josephson effect) in the regions I and II. The time-periodic modulation of the Raman frequency may lead to resonant responses of solitons in BEC with SOC. This phenomenon is possible to occur for the imbalanced populations between soliton components, which are produced initially with different phases. To study the dynamics of this kind of process, in both the cases of region I () and region II (), we first implement the Josephson oscillations for constant Raman frequencies , by introducing a phase between the two soliton components, before starting the time evolution of the coupled system.

For the time modulated Raman frequency, we consider the following expression:

 Ω(t)=Ω0+Ω1cos(ωt), (5)

where is the amplitude, with the frequency of the oscillations. Different regimes are possible in the dependence of the values of the modulating frequency , such that we can have slow (), resonant () or rapid () modulations. In this section, we are mainly concerned with the intermediate regime, where we can have resonant responses, such that we are going to assume small amplitude for the oscillations, . The case of very slow frequency can be reported to the previous study presented in section II-B. The other regime for the time-perturbed Raman, with high frequencies and large amplitudes , we are considering in section IV.

The studies in this section are done by a full numerical simulation, as well as by some analytical considerations through a variational procedure. In order to study the interference effect on the Josephson oscillation due to a time-modulated Raman frequency, we employ a variational approach in region I (where , and regular soliton solutions are obtained), considering the ansatz

 (6)

where , ,, and are time-dependent parameters, where we have the assumption that the solitons have the same width and center-of-mass (i.e., both components overlap), which are confirmed by numerical simulations. By considering the Lagrangian density for Eq. (2) and the above variational ansatz, we have

 L(x,t) = [i2(u∗dudt+v∗dvdt)+ikL2(u∗dudx−v∗dvdx)+cc] (7) − 12∣∣∣dudx∣∣∣2−12∣∣∣dvdx∣∣∣2−Ωu∗v−Ωv∗u + 12|u|4+γ2|v|4+β|u|2|v|2,

with the Lagrangian given by :

 L = −2∑i=1Ni[dϕidt+dkidtx0+14a2+k2i2+(−)ikikL] (8) − 2Ω(t)√N1N2e−a2k2−cos(2k−x0+ϕ) + 12√2πa(N21+γN22+2βN1N2).

Here, and in the following, we are using the definitions and . The number of atoms for each component is given by , with the total number being conserved. By assuming weak SOC parameter, , following arguments given in Ref. Malomed2017 (), the corresponding Euler-Lagrange equations for the parameters are given by

 dx0dt=k+,dk+dt=2e−k2La2Ω(t)kL√1−Z2sin(φ),dZdt=−2e−k2La2Ω(t)√1−Z2sin(φ),dφdt=2kLk++ΛZ+2e−k2La2Ω(t)Z√1−Z2cos(φ),⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (9)

where , , and, for simplificity, we fix . In case of a hyperbolic-type ansatz, the exponential factor in the above equations, derived from a Gaussian ansatz, is changed as (where is the soliton amplitude). In the weak SOC limit, when or , this factor reduces to one.

This system is analogous to the one considered in Malomed2017 () for a constant , where it was shown that in the Euler-Lagrange system for the parameters, the equation for the center-of-mass, , can be approximately solved as Thus, the center-of-mass oscillations are small, with its amplitude being of the order of By considering our numerical simulations, as discussed in more detail in the next subsection, with results given in Figs. 3 and 4, for the case that and , we have confirmed that the center-of-mass oscillations agree with the estimated value, being , as shown in Fig. 4.

Then, for small values of (as compared with ), we can consider the following coupled system:

 dϕdt = ΛZ+2e−k2La2Ω(t)Z√1−Z2cos(ϕ), (10) dZdt = −2e−k2La2Ω(t)√1−Z2sin(ϕ), (11)

which describes the internal Josephson oscillations of atomic populations between two pseudo-spin components. For a constant Raman frequency, these oscillations have been studied in Ref. Zhang2012 (). It appeared when investigating the macroscopic quantum tunneling obtained in a double-well potential having a barrier between the wells with constant height Shenoy (). This was also studied in Refs. AK (); Milburn (); Boukobza () for the case that the barrier between wells have their heights oscillating in time.

At some frequencies of the modulations, it is possible to verify parametric resonance in the Josephson oscillations. In the following we will discuss the full numerical results, by considering the two possible regions defined in Fig. 1 by the relations between the Raman and SOC parameters. The results obtained in region I, where we have bright type solitons, are shown to be fully compatible with a variational approach. However, for the region II, where we have striped soliton solutions, the coupled system is not so amenable to simplified variational analysis, such that we can provide only rough estimates for some limiting situations. Therefore, in the case of time evolution for the striped solitons, with constant and time-dependent Raman frequencies, our study relies mostly in full-numerical simulations, which are shown to provide convergent results with high numerical precision.

### iii.1 Results for Josephson oscillations in region I

The results in this case are for , considering the Raman frequency constant, as well as time-modulated. In the case of time-modulated Raman frequency, we have also introduced a variational analysis to provide an estimate for the localization of the modulation frequency leading to the resonant behavior.

#### iii.1.1 Constant Raman frequency: Ω(t)=Ω0

When is constant, the frequency of the free oscillations is given by

 ωJ=√2Ω0(2Ω0+Λ)(=2Ω0,forβ=1). (12)

From numerical simulations for free Josephson oscillations of the full system of the GP equation, the results obtained in the region I () are represented in Figs. 3 and 4, by considering the spin-orbit coupling parameter , with two constant Raman frequencies, given by (upper panel) and (lower panel). In Fig. 4 we show the density plots, for , and , corresponding to the case that , for the time interval .

The purpose, in this case, is to verify the dependence of the oscillating behavior on different values of the initial phase introduced between components when starting the evolution. By considering three values for the initial phase , it is shown that the maximum of the periodic atom transfer occurs at , reaching almost 100% of atoms in the case that . It is also shown that the phase affects only the amplitude of the oscillations, but not the frequency. The constant Raman parameter will determine the frequency of the oscillations, which is given by , confirming the theoretical prediction (12).

#### iii.1.2 Time-modulated Raman frequency: Ω(t)=Ω0+Ω1cos(ωt)

Now, let us analyze the case when is modulated in time, such that in order to verify the localization of possible resonant behaviors. For that, we can consider two limiting conditions for the Eq. (10): one applied for , when ; the other for the regime of macroscopic quantum localization, as follows.

Let us consider the linear regime case, when and . Within these conditions in Eq. (10), using the second derivative for , we obtain a modified Mathieu differential equation, with the main term in oscillating with a frequency , where For that, when considering small values of and , the standard analysis as given in Ref. LL () can be applied, which leads to a resonance at (in case , so ). In this case, parametric resonances are also expected to occur for , where is the frequency of free Josephson oscillations given in Eq. (12).

Results of our investigations on resonant interferences which can occur in region I are being presented in Figs. 5 to 8. They are obtained from numerical simulations of the full coupled system (2), with the Raman and SOC parameters and , respectively. For the time modulation of the Raman parameter, given by Eq. (5), we assume the amplitude given by . The choice of these parameters in region I are to distinguish more clearly the manifestation of resonant interferences in the Josephson oscillations between the atom numbers of the two components, , . We found also illustrative to provide some density plots, for the components and total profiles, corresponding to the results presented in Fig. 5. Therefore, in Fig. 6 we are showing the case where we have in Fig. 5. Here, we should remark that, the perturbed results obtained when we are not close to the resonant interference regions are shown to be almost identical to the unperturbed results (as observed by comparing the first column of panels of Fig. 6 with the third column). Indeed, quite small center-of-mass oscillations are verified near the initial localization. Another point that can be observed from these results is that the center of mass of the soliton is strongly affected by the resonant behavior, such that, it can be verified in the central panels of Fig. 6 that the central position is moving from at to at . There is no change observed in the center of mass for the other cases, outside the resonant region. The numerical simulations of the variational system (8) confirms this behavior of the motion of the center of mass for the solitonic components and the oscillations of the atomic imbalance at the resonance.

We should comment that, for the values of the frequency , the resonant position (“window”) is quite sharp, verified in our numerical simulation, such that the resonant perturbations are being confirmed only for very close to 2 and 4 , which makes the simulations quite time demanding. As shown in Fig. 5 and in the left panel of Fig. 7, for a large time interval going up to , one of the resonant position are detected for . By a slight larger deviation of this frequency the results for the oscillations are about the same as given for the non-perturbed case () shown in the left panel of Fig. 7. The other resonant position, as shown in the right panel of Fig. 7, for , is found in an even smaller range of , given by , with fluctuation start appearing when we have exact . In Fig. 5 we are also showing how the resonant behaviors are affected by changes in the nonlinear parameter . As shown it is enhanced in case that , with the variation having peaks with maxima close to 0.9. The case of , for , is also presented in the left frame of Fig. 7, for comparison with the unperturbed results of the Josephson oscillations.

With Fig. 8, we conclude the analysis of the results shown in Figs. 5 Figs. 6 and 7, by presenting the behaviors of density profiles (total and for each component) in two-dimensional plots, for different time positions of the evolution.

### iii.2 Results for Josephson oscillations in region II

Now, let us consider the Josephson oscillations between components of the striped soliton solutions, corresponding to the region II in the dispersion relation, which are given by . We perform this study by considering full numerical simulations of the corresponding GP coupled equations. First, we provide some results obtained for Josephson oscillations in the case that we have constant Raman frequency. Next, we consider the more general case, where the Raman frequency is modulated in time, and we can have resonant results at some particular values of the modulating frequency .

We start the study of this section by considering a variational analysis, where we need to introduce the momentum , which provides the momentum position of the minima shown in the lower panel of Fig. 1. By observing that the coupled equations for the imbalanced populations and relative phase are not easy to be derived in explicit form, in a more general case, let us assume that the tunneling between components occurs for the same sign of . The time modulations for are not inducing transitions between oppositely propagating modes with . To have such transitions, so called momentum Josephson oscillations, we need parameters with periodic modulation in space stripeJoseph (); 2016-Luo (). Therefore, by assuming for the components the same ansatz as given in Eq. (6), but with and considering the center fixed at , we arrive to the same coupled expressions (10) and (11), except that the equation for contains an additional term . By linearizing the system relative to , we obtain

 dϕdt = 2kLk0+[Λ+2Ω(t)]Zcos(ϕ), (13) dZdt = −2Ω(t)sin(ϕ). (14)

For a constant and with , we have

 ϕ≈ϕ0+(2kLk0)t, (15)

implying that the population imbalance is oscillating with the frequency . Therefore, we should expect a different behavior of the results, in comparison with the region I, in the initial stage (defined by ). For larger time of the propagation, the frequency for the oscillations should approach the same ones as verified for the region I.

#### iii.2.1 Constant Raman frequency: Ω(t)=Ω0

For the numerical simulation of the Josephson oscillations obtained in region II, we first select some results obtained for constant values of , which are given in Figs. 9 and 10, by considering an initial phase difference between components given by . In these cases, by considering (Fig. 9) and (Fig. 10), with several values of , we can verify clearly that we have an initial stage of the oscillations, where the frequencies are not depending on , but only on the values of , being for and for . The values of affects only the amplitude of the initial oscillations. However, for larger times, the behavior of the frequencies are similar as in region I. Another observation from these results is that, for long-time interval the Josephson oscillations are being damped as we increase the difference . In Fig. 10, we present two inset panels from where we can verify the initial and intermediate oscillation patterns. In the inset with , just after starting the evolution, the frequency is about the same for all the three cases, , not depending on . In the other inset, for , after a transient time interval, the frequency change to (as in case of region I).

#### iii.2.2 Time-modulated Raman frequency: Ω(t)=Ω0+Ω1cos(ωt)

When studying the phase-dependence of the atom-number oscillations for striped soliton solutions, we first observe that the amplitude of the oscillations depends on the initial phase difference, as already verified in the case that we have constant Raman frequency parameter, given by . Therefore, before considering the case where we have the Raman frequency perturbed in time, we have studied the phase-dependence of the atom-number oscillations for striped soliton solutions during time evolution. In this numerical study, we have verified that for arbitrary initial fixed phase (from 0.01 to ) introduced between components, only the amplitude of the oscillations are being affected, which are being verified by the transient time just after starting the evolution of the solutions. The frequency of the oscillations does not depend on the strength of the Raman frequency , at least during the transient time till the oscillations become stable. In a longer time interval, after the transient time, the frequency of the oscillations will correspond to the Raman frequency, being given by , as discussed for the case of regular soliton solutions. As the initial phase between the components can be arbitrary and will not affect the natural frequency of the oscillations between and , when studying time-perturbed Raman frequency, in general we choose this phase to be , such that the amplitude of the natural oscillations is not too large, as well as not too small.

In the case of periodic modulations of in the region II () the results of our numerical simulations to identify parametric resonances is first being exemplified by the Figs. 11 and 12. In Fig. 11, for , we present two panels considering 40 (upper) and 60 (lower). In both we are plotting the perturbed case considering the amplitude of the oscillations given by and , with the frequency the linear frequency verified in the transient time interval (about 126127). The non-perturbed case () is also shown for comparison, in both the cases. We should emphasize that in general the Josephson oscillating behavior is about the same as it happens for the unperturbed case , except close to the specific values for and where resonant interference behaviors are detected. With Fig. 12, for a time interval of half-period of the Josephson oscillations, we are representing the profiles of the two component densities ( and ), for the case shown in the upper panel of Fig. 11 with and (the other parameters are the same). The oscillation dynamics is represented in two panels, given for (upper panel) and (lower panel), considering that a complete period is close to . The panels are indicating (through the corresponding densities) the atom-number oscillation between the components

For a long-time interval, resonant behaviors are expected to occur when considering cases where the natural frequency of the oscillations are still surviving in the unperturbed case. For that, in Fig. 13, we are showing results for a simulation with and , with and initial phase . In this case, by taking , we can observe a resonant interference that occurs for . For this case, the striped soliton profiles of both components are also being represented in the two lower panels of the figure. In the left panel we have them at , and in the right panel for .

To conclude our study related to striped solitons and resonant interference effects, we present results obtained in longer time intervals, for the case that the unperturbed Raman is given by , with , as in Fig. 13, but with a much smaller amplitude of the modulations, such that . The investigation of the interval of where interferences can be found is shown in Fig. 14, considering a small initial phase of oscillations between components given by . As shown by the set of five panels (for ) with given values of , resonant interference effects due to the perturbation are verified only in the interval , with maxima interferences occurring for (the middle panel, where we have also included with dashed line the unperturbed case, for comparison). The results for and are almost identical with the non-perturbed case, . Therefore, we select the case with to show in more detail the oscillating behavior, which is presented in Figs. 15 and 16. In the lower panel of Fig. 15, we consider a larger time interval with (lower panel). The middle panel () serves to show the change in the frequency of the oscillations, such that for each two cycles another cycle is emerging, which can be verified for . In all the three given panels, for comparison we are including in dashed line the unperturbed case.

## Iv High frequency modulations. Averaged GP equations

In the case that we have rapidly and strongly varying Raman oscillations , it is useful to derive the corresponding averaged GP equation, such that one can reduce the time-dependent modulated Raman frequency to the constant one , by renormalizing the spin-orbit coupling and the non-linear parameters, as we are going to show in this section. By matching the averaged results with the full-numerical ones, obtained with real-time evolution, we are also verifying numerically how fast and strong should be the time oscillations in order to validate the averaging approach. In order to derive the averaged over rapid modulations system of equations, we first apply the following unitary transformation SKMA (); 2013Zhang () in Eq. (1):

 Φ≡(UV) = (cos(q)isin(q)isin(q)cos(q))(uv), (16)

where is given by the requirement that the time-dependent part of the Raman frequency does not appear explicitly in the coupled equation for . When performing the time averaging of the interactions, together with the SOC parameter , the parameters of the non-linear interaction have also to be renormalized. They are replaced by parameters that contains zero-order Bessel function, considering that

 12π∫2π0d(ωt)exp(inΩ1ωsin(ωt))=J0(nΩ1ω). (17)

Then, by defining , the coupled equation, averaged over the period of rapid modulations, with , can be written as2013Zhang ():

 i∂Φ∂t =[−12∂2∂x2−ikLJ0(χ)σz∂∂x+Ω0σx](UV) − (α+|U|2+α1|V|2α0U∗Vα0V∗Uα−|V|2+α1|U|2)(UV),

where

 α0≡(β−1+γ2)1−J0(2χ)4,α±≡α0+1+γ2±1−γ2J0(χ),α1≡β−2α0.⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (19)

In the case of gauge symmetry, with , we have and ; i.e., the nonlinear part of the above coupling equation for () is exactly the same is the ones obtained for , such that the time averaging is only renormalizing the SOC parameter to

 κeff=kLJ0(2Ω1ω), (20)

as one can verify by comparing the coupled Eqs. (IV) with (2). This approach for tuning of the SOC parameter has been confirmed recently in an experiment reported in Ref. 2015Spielman (). Therefore, when considering rapid variations of the Raman frequency, the spin-orbit coupling can be tuned in order to control the solitons in a BEC with SOC. In particular, it can be quite useful to transform striped solitons to regular solitons, and vice versa. With the appropriate ratio between amplitude and frequency of the Raman oscillations, a given value of for region II, where , can be changed to , where we obtain regular soliton solutions, such that all the theory developed before (in sections II and III) for constant Raman frequency, can be applied.

The above can be exemplified by the results shown in Fig. 3, which are for regular soliton solutions, with = 4 and 80 and 20, respectively. These results are for region I, but can also be applied to the case that we have originally larger than , if the time modulation of the Raman frequency , given by Eq. (5), is such that the ratio between and will give us 4. We could take initially, , for example, as it is larger than in both the cases shown in Fig. 3 , with the parameters of the time-modulating Raman such that .

When considering other values for the nonlinear parameters, as a general remark we noticed that stable soliton solutions are obtained for attractive two-body interactions. Another remark, when considering the averaged approach, is that for , we can also have conditions with zero in the off-diagonal terms of the non-linear two-body matrix, which brings the Eq. (IV) to the same format as Eq. (2). This happens for , with , , and . The more general cases, as for , or , new terms appear, corresponding to the effective four-wave mixing ( ). These terms can lead to new possibilities, such as a way to control the atom number oscillations between two components (internal Josephson effectShenoy (); Zhang2012 ()).

#### iv.0.1 The solitonic solutions

The solitonic solutions for the averaged GP equations can be found by applying the multi-scale method AFKP () to the two regions defined by the linear spectrum, which are given by Eq(3). By using this multi-scale method, for values of the chemical potential near the bottom of the dispersive curve, with (), is the free parameter, in region I (see Fig. 1), we obtain

 u(I)s = ε√2w0α++α1+α0sech(ε√2w0Δeffx), (21) v(I)s = −u(I)s,Δeff=1−κ2effΩ0. (22)

In the region II, where and two minima exist in the momentum space, we can take the chemical potential as (see Fig. 1), with

 wmin=12κ20−κ2eff,κ0≡±√κ2eff−Ω20/κ2eff (23)

and look for solutions of the form . For the result, we obtain a bright soliton solution with the form given by

 (u(II)sv(II)s)=(Ω0−κeff(κeff±κ0))εf(x)e±iκ0x−iμt)√|κeff±κ0| (24)

where

 f(x) = √2w0κeff√α+(κ4eff+κ2effκ20)+(α1+α0)Ω20×, (25) × sech⎛⎜⎝ε ⎷2w0κ2effκ20x⎞⎟⎠.

Analogically, the striped soliton solution can be found as linear superpositions of solutions represented by Eq.(24). As already known, these solutions are used to describe the longitudinal and transversal spin polarizations of the solitons AFKP (), with

 ⟨σz⟩ = 1N∫∞−∞dx(|u|2−|v|2), ⟨σx⟩ = 1N∫∞−∞dx(u∗v+uv∗), (26) N ≡ ∫∞−∞dx(|u|2+|v|2).

In the region I, where , the solitons are fully polarized along the axis. The same approach is valid for stripe solitons in the region II. However, the polarization along is not zero for solitons with momentum . From Eqs. (21)-(25), we obtain

 ⟨σz⟩(II)=−√1−Ω20/κ4eff,⟨σx⟩(II)=−Ω0κ2eff. (27)

Thus, by varying the ratio , and so , we can observe quantum phase transition in the pseudo-spin polarization of the soliton. These results for the soliton polarization are analogous to the ones obtained for the repulsive BEC in the framework of the Dicke model in 2013Zhang ().

With the understanding that the results obtained in this section are valid in a more general context for constant values of the Raman frequency, with Fig. 17 we are showing the dependence of the energy and chemical potential on the number of atoms for the case that , , , when considering (with both and very large), which give us . Note that, in this simple case we have for the dispersion relation (3) , with the signal of the give SOC moving from the original to a negative one, . Threfore, after considering the time averaging, in this particular case with , we obtain , such that both and have the same shape as shown in Fig. 1, but with minima given at (for ; and (for ). As we are in region II, even after the averaging, the soliton solutions are not regular ones, being expected to have shapes with some oscillations. In Fig. 18, we are illustrating the kind of solutions we obtain, by presenting the real and imaginary parts of the wave-function components when considering the particular case with , and .

#### iv.0.2 Full numerical versus averaged results

To conclude this section, we are comparing the time evolution results obtained with the effective time-averaging approach (where the SOC parameter is ), with the ones obtained in real time, with the SOC parameter and explicit Raman frequency modulated by .

With Figs. 19, 20 and 21, we are exemplifying our results for the comparison of averaged results with real-time dependent numerical simulations. All the results for the time-averaged formalism that are shown in these examples are verified to be numerically very stable in the time evolution.

The results given in Fig. 19 are for the region I, with 120 (), by considering the SOC parameter . For the time-dependent Raman frequency we assume and such that , implying that . Therefore, in this case, the averaged SOC parameter is . As shown in the four panels, the averaged results present good agreement with the real-time simulations when is about 10 times or more larger than .

For the region II, we are illustrating with Figs. 20 and 21, for two quite different combinations of Raman frequency and SOC parameters.

In Fig. 20, we present our results obtained for and 12, with and such that . As and , we are in region II (). The results are shown in three panels, for different values of (lower panel), (middle panel) and (upper panel). The parameters used in this case correspond to one of the examples presented in Fig. 9 for Josephson oscillations, where we have constant , with . We should note that the striped solitons shown in Fig. 20 have the main maximum at the center, with only one pair of maxima visible in each side, due to the choice of parameters which are close to the border between the regions for striped and regular solitons.

In Fig. 21, we are considering a case where the effective SOC becomes negative, and we are more deeply inside the region II. By departing from a large value of