1. Introduction
Spin–orbit coupling (SOC) plays an essential role in understanding several physical phenomena which include quantum Hall effects [1, 2], topological insulators [3, 4] and topological superconductors [5, 6]. It promises applications in spintronics and topological quantum computing. In solid-state materials, the presence of disorder and contaminant, however, makes it difficult to observe such phenomena. The remarkable experimental breakthroughs in the realization of SOC in Bose–Einstein condensates (BECs) of neutral atoms provides a platform for exploring SOC in a defect-free medium with great flexibility [7, 8]. The possibility to generate relatively stronger SOC in BECs leads to exotic phases due to the non-conservation of the spin part of the particles during the motion. In addition, the BEC with angular SOC supports half-skyrmions in the ground state [9]. The breathing mode oscillation frequency of this skyrmion state is not universal but it depends on the strengths of SOC and nonlinear interaction [10].
The inter-atomic interaction in BECs allows us to explore nonlinear phenomena such as rouge waves, dark solitons, gap solitons and bright solitons. Under the action of Raman-induced SOC, the dynamics of solitons change appreciably. The matter–wave soliton in a SOC-BEC suffers from the lack of Galilean invariance and induces shape change with the increase of velocity [11]. A SOC-BEC supports the so-called stripe solitons characterized by density modulations in the form of stripes. The spatial motion of the soliton is coupled to the spin degree of freedom since the Raman transition can cause both momentum transfer and spin flipping. It has recently been demonstrated experimentally that the center of mass motion of a BEC in a harmonic trap is significantly affected by the spin degree of freedom [12, 13].
In this work, we comprehensively analyze the spin dynamics and center-of-mass motion of bright–bright solitons in spin–orbit coupled BECs in presence of optically induced periodic potential. We pay special attention to the influence of spin degrees of freedom on the center-of-mass motion. It is seen that the spin part of the solution in the presence of an optical lattice (OL) can be represented on the Bloch sphere and the effective spin dynamics can be described by the so-called nonlinear Bloch equation [14, 15]. Understandably, an effective magnetic field is produced in the spinor reference frame. This field causes the spinors to precess. The precessional frequency is found to depend on the lattice parameters and, hence, coupling between the spatial motion of soliton and its spin degrees of freedom is enriched due to OLs. Specifically, the amplitude and frequency of population oscillation due to spin flipping between the spin–orbit coupled states are modified in presence of OLs.
SOC creates anharmonic characteristics in the collective spin dynamics in BECs with OLs. We employ a fixed point analysis both in the presence and absence of OLs and find that the effective spin dynamics of the system are periodic with a single period in the absence of OLs and periodic multiplicity occurs in the presence of OLs. This ultimately leads to chaos in the effective spin dynamics. Chaotic dynamics of BECs in a one-dimensional quasi-periodic OL and random disordered potentials were studied in [16]. Recently, it has been demonstrated that SOC may create anharmonic characteristics in collective dipole oscillations beyond the effective mass approximation [17, 18].
The paper is organized as follows. In section 2 , we introduce a set of coupled Gross–Pitaevskii equations (GPEs) to describe the spin dynamics of two solitons in the presence of the spin–orbit interaction and OLs. Based on variational calculations we derive equations for different parameters of solitons and hence find equations for the spin dynamics in the presence of OLs using a similarity transformation. We study the dynamics of effective spin components and the center-of-mass from the variational calculations and also justify the variational predictions by the direct numerical simulation of the GPEs in section 3 . We envisage nonlinear spin dynamics of the soliton pair in section 4 through the analysis of phase-space trajectories and Lypunov exponents in the parameter space [19–21]. Finally, we summarize the results in section 5 .
2. Theoretical formulation
The spin–orbit interaction generates coupling between the two pseudo-spin states created by using proper detuned Raman lasers in BECs. The two-component BECs thus created can be described by the coupled mean-field GPE. The GPEs for spin–orbit coupled BECs in quasi-one-dimension is given by2 ), β and α are the strengths of the spin–orbit interaction and Raman coupling, strengths of inter and intra-component interactions are denoted by cj and cij, respectively. We replace t → t/ω⊥, x → x a⊥ and V0 → V0ℏω⊥ in equation (2 ) such that it becomes dimensionless.
$\begin{eqnarray}i\left(\begin{array}{c}{{\rm{\Psi }}}_{t}\\ {{\rm{\Phi }}}_{t}\end{array}\right)=\left(\begin{array}{cc}{H}_{11} & {H}_{12}\\ {H}_{21} & {H}_{22}\end{array}\right)\left(\begin{array}{c}{\rm{\Psi }}\\ {\rm{\Phi }}\end{array}\right).\end{eqnarray}$
Here $\begin{eqnarray*}\begin{array}{rcl}{H}_{11} & = & -\displaystyle \frac{1}{2}{\partial }_{x}^{2}-{\rm{i}}\beta {\partial }_{x}+V(x)+{c}_{1}| {\rm{\Psi }}{| }^{2}+{c}_{12}| {\rm{\Phi }}{| }^{2},\\ {H}_{22} & = & -\displaystyle \frac{1}{2}{\partial }_{x}^{2}+{\rm{i}}\beta {\partial }_{x}+V(x)+{c}_{12}| {\rm{\Psi }}{| }^{2}+{c}_{2}| {\rm{\Phi }}{| }^{2}\end{array}\end{eqnarray*}$
and $\begin{eqnarray}{H}_{12}={H}_{21}=\alpha \end{eqnarray}$
with the effective potential $\begin{eqnarray}V(x)=\displaystyle \frac{1}{2}{\lambda }_{\perp }^{2}{x}^{2}+{V}_{0}\cos (2{k}_{\mathrm{lat}}x+{\phi }_{L}).\end{eqnarray}$
Here the first term stands for harmonic trap with λ⊥ = ωx/ω⊥, the ratio of longitudinal and transverse frequencies, and the second term gives an OL potential with amplitude V0, wave number klat and phase φL. In equation (The SOC-BEC with an attractive atomic interaction supports the bright soliton solution in the absence of OLs and this soliton may contain nodes. In view of this, we consider the following trial solution of equation (1 ) [14]Appendix ).7a )–(7d ) describe the dynamics of different parameters of soliton solutions of the system. Particularly, equation (7d ) shows how the population imbalance between the soliton components varies with phase difference. In this context, we note that this equation is also useful to study Josephson-type oscillation in SOC-BECs [22, 23].
$\begin{eqnarray}\left(\begin{array}{c}{\rm{\Phi }}\\ {\rm{\Psi }}\end{array}\right)=\sqrt{a/2}\left(\begin{array}{c}\sin \theta \,{\rm{sech}} ({ax}+{x}_{c}){{\rm{e}}}^{{\rm{i}}[{p}_{1x}x+{\phi }_{1x}]}\\ \cos \theta \,{\rm{sech}} ({ax}+{x}_{c}){{\rm{e}}}^{{\rm{i}}[{p}_{2x}x+{\phi }_{2x}]}\end{array}\right),\end{eqnarray}$
where θ, a, xc, pjx and φjx are variational parameters. Here θ and a−1 determine, respectively, the population imbalance between two components and their width, xc is the center-of-mass, pjx and φjx represent the wave number and phase of the jth component, respectively. The trial solution is normalized to N = a/4. We obtain an averaged Lagrangian density involving the variational parameters (a, xc, p1x, p2x, φ1x and φ2x) using $\left\langle { \mathcal L }\right\rangle ={\int }_{-\infty }^{+\infty }{ \mathcal L }\,{\rm{d}}{x}$. It is given by $\begin{eqnarray}\left\langle { \mathcal L }\right\rangle =\left\langle {{ \mathcal L }}_{1}\right\rangle +\left\langle {{ \mathcal L }}_{2}\right\rangle +\left\langle {{ \mathcal L }}_{3}\right\rangle +\left\langle {{ \mathcal L }}_{4},,\right\rangle \end{eqnarray}$
where $\begin{eqnarray}\left\langle {{ \mathcal L }}_{1}\right\rangle =-\pi {V}_{0}\displaystyle \frac{\cos (2{k}_{{\rm{lat}}}{x}_{c}/a+{\phi }_{L})}{\sinh (\pi /a)},\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}\left\langle {{ \mathcal L }}_{2}\right\rangle & = & -\displaystyle \frac{\pi \alpha }{2}\cos \left[\displaystyle \frac{{x}_{c}}{a}{p}_{x}+\phi \right]\displaystyle \frac{{p}_{x}\sin 2\theta }{\sinh (\pi {p}_{x}/2a)}\\ & & +\,{\beta }_{1}\left({p}_{1x}{\sin }^{2}\theta -{p}_{2x}{\cos }^{2}\theta \right),\end{array}\end{eqnarray}$
$\begin{eqnarray}\begin{array}{rcl}\left\langle {{ \mathcal L }}_{3}\right\rangle & = & \left(\displaystyle \frac{1}{2}\,{p}_{2x}^{2}-\displaystyle \frac{{x}_{c}}{a}{\dot{p}}_{2x}+a{\dot{\phi }}_{2x}\right){\cos }^{2}\theta \\ & & -\,\left(\displaystyle \frac{1}{2}{p}_{1x}^{2}-\displaystyle \frac{{x}_{c}}{a}{\dot{p}}_{1x}+a{\dot{\phi }}_{1x}\right){\sin }^{2}\theta ,\end{array}\end{eqnarray}$
$\begin{eqnarray}\left\langle {{ \mathcal L }}_{4}\right\rangle =-\displaystyle \frac{{a}^{2}}{6}(a+{c}_{1}{\cos }^{4}\theta +{c}_{2}{\sin }^{4}\theta +\displaystyle \frac{{c}_{12}}{2}{\sin }^{2}2\theta ),\end{eqnarray}$
where kn = (p1x − p2x)/2, kp = (p1x + p2x)/2 and φn = φ1x − φ2x. We make use of the Ritz optimization procedure and obtain the following equations for the variational parameters (see the $\begin{eqnarray}\displaystyle \frac{{\rm{d}}\left\langle z\right\rangle }{{\rm{d}}{t}}={k}_{p},\end{eqnarray}$
$\begin{eqnarray}{\dot{k}}_{p}=\displaystyle \frac{2\pi \alpha {\lambda }^{2}\sin \phi }{a\sinh (\pi \lambda /a)}\sin 2\theta +\displaystyle \frac{2\pi {V}_{0}{\lambda }^{2}\sin \phi }{a\sinh (\pi \lambda /a)},\end{eqnarray}$
$\begin{eqnarray}\dot{\phi }=2\lambda {k}_{p}-\displaystyle \frac{2\pi \alpha \lambda \cos \phi }{a\sinh (\pi \lambda /a)}\cot 2\theta ,\end{eqnarray}$
$\begin{eqnarray}\dot{\theta }=-\displaystyle \frac{\pi \alpha \lambda }{a\sinh (\pi \lambda /a)}\sin \phi .\end{eqnarray}$
Here φ = 2kn⟨z⟩ + 2φn, ⟨z⟩ = − xc/a, β ≈ kn, kn = λ, klat ≈ kn, φL ≈ 2φn, c1 = c2 = c12 and ${\dot{k}}_{n}=0$. Equations (2.1. Effective equations for the spin components
Let us consider normalized complex-valued spinors: ${{\rm{\Psi }}}_{j}=\sqrt{\rho (x,t)}{\chi }_{j}$ with χ = (χ↑, χ↓) and $\Psi$j = ($\Psi$, Φ) such that ρ = ∣$\Psi$∣2 + ∣Φ∣2 and ∣χ↑∣2 + ∣χ↓∣2 = 1, where ${\chi }_{\uparrow }=\sin \theta \,{{\rm{e}}}^{{\rm{i}}{\phi }_{1x}}$ and ${\chi }_{\downarrow }=\cos \theta \,{{\rm{e}}}^{{\rm{i}}{\phi }_{2x}}$. Spin expectation value can be defined through the transformation9 ) as
$\begin{eqnarray}{\bf{S}}={\chi }^{{\rm{T}}}\sigma \chi .\end{eqnarray}$
Here σ ≡ {σx, σy, σz} is the set of Pauli spin operators. Expectation values of different spin components are given by $\begin{eqnarray}{S}_{x}={\chi }^{{\rm{T}}}{\sigma }_{x}\chi ={\chi }_{\uparrow }^{* }{\chi }_{\downarrow }+{\chi }_{\downarrow }^{* }{\chi }_{\uparrow }=\sin 2\theta \cos \phi ,\end{eqnarray}$
$\begin{eqnarray}{S}_{y}={\chi }^{{\rm{T}}}{\sigma }_{y}\chi ={\chi }_{\downarrow }^{* }{\chi }_{\uparrow }-{\chi }_{\uparrow }^{* }{\chi }_{\downarrow }=-\sin 2\theta \sin \phi ,\end{eqnarray}$
$\begin{eqnarray}{S}_{z}={\chi }^{{\rm{T}}}{\sigma }_{x}\chi =| {\chi }_{\uparrow }{| }^{2}-| {\chi }_{\downarrow }{| }^{2}=-\cos 2\theta .\end{eqnarray}$
Equations for the dynamics of spin components are obtained from equation ( $\begin{eqnarray}{\dot{S}}_{z}=2\tilde{{\rm{\Omega }}}\,{S}_{y},\end{eqnarray}$
$\begin{eqnarray}{\dot{S}}_{y}=-2\tilde{{\rm{\Omega }}}{S}_{z}-2\lambda {k}_{p}{S}_{x},\end{eqnarray}$
$\begin{eqnarray}{\dot{S}}_{x}=2\lambda {k}_{p}{S}_{y},\end{eqnarray}$
where $\tilde{{\rm{\Omega }}}=\tfrac{\pi \alpha \lambda }{a\sinh (\pi \lambda /a)}$ and $\begin{eqnarray}{k}_{p}=-\lambda {S}_{z}-\displaystyle \frac{2\lambda {V}_{0}}{\alpha }{\sin }^{-1}({\tilde{S}}_{z})+{c}_{0}\end{eqnarray}$
with ${\tilde{S}}_{z}=\sqrt{(1+{S}_{z})/2}$ and ${c}_{0}=\lambda {S}_{z0}+\tfrac{2\lambda {V}_{0}}{\alpha }{\sin }^{-1}({\tilde{S}}_{z0})$. Here we take $\left\langle \dot{z}(0)\right\rangle ={k}_{p}(0)=0$. Note that a = (Sx, Sy, Sz) represents a vector on the Bloch sphere and it satisfies, ρ = (1 + a · σ).Equations (10a )–(10c ) for the spin dynamics can be expressed in the form
$\begin{eqnarray}\dot{{\boldsymbol{S}}}={\boldsymbol{S}}\times {\boldsymbol{B}}\end{eqnarray}$
with ${\boldsymbol{B}}=\left(-2\tilde{{\rm{\Omega }}},0,2\lambda {k}_{p}\right)$. This is the so-called nonlinear Bloch equation. Thus, we see that spinors χ↑ and χ↓ face an effective magnetic field B [24]. The field causes spin precession and the frequency of which depends on the lattice and SOC parameters. Understandably, the frequency of precession is ${\omega }_{p}=\gamma | {\boldsymbol{B}}| ={(4{\lambda }^{2}{k}_{p}^{2}+4{\tilde{{\rm{\Omega }}}}^{2})}^{1/2}$, where the effective gyromagnetic ratio γ = 1. Therefore, the spin parts of the solutions can be expressed as ${\chi }_{\uparrow }=\sin \theta \,{{\rm{e}}}^{{\rm{i}}({\phi }_{1x}+{\omega }_{p}t)}$ and ${\chi }_{\downarrow }=\cos \theta \,{{\rm{e}}}^{{\rm{i}}({\phi }_{2x}-{\omega }_{p}t)}$.In order to study the dynamics of the spin component, we combine (10a )–(10c ) and (11 ), and write14 ), ${v}_{1}=2{c}_{1}{V}_{0}{\lambda }^{2}\tilde{{\rm{\Omega }}}/\alpha $, v2 = − 8V0λ3(c0 + λ)/α and ${v}_{3}=-8{V}_{0}^{2}{\lambda }^{4}/{\alpha }^{2}$. We see that for V0 = 0, the frequency of oscillation of Sz under harmonic approximation is $\sqrt{{m}_{1}}$. In the presence of the OLs, the oscillation frequency further changes to $\sqrt{{m}_{1}+{v}_{1}}$ in the linear limit. In the nonlinear limit, one can, however, expect richer spin dynamics due to the anharmonic response of the system. In the following, we consider both regular and nonlinear spin-dynamics in the presence of OLs.
$\begin{eqnarray}{\ddot{S}}_{z}+{m}_{1}{S}_{z}+{m}_{2}{S}_{z}^{2}+{m}_{3}{S}_{z}^{3}+{m}_{4}={V}_{L}({S}_{z})\end{eqnarray}$
with $\begin{eqnarray}{V}_{L}({S}_{z})={\sin }^{-1}({\tilde{S}}_{z})\left[{v}_{1}+{v}_{2}{S}_{z}+{v}_{3}{S}_{z}{\sin }^{-1}({\tilde{S}}_{z})\right].\end{eqnarray}$
Here ${m}_{1}=-2{c}_{1}\tilde{{\rm{\Omega }}}{\lambda }^{2}+4{\tilde{{\rm{\Omega }}}}^{2}+2{\lambda }^{2}{c}_{0}^{2}$, m2 = 4λ3c0, m3 = 2λ4 and ${m}_{4}=2\lambda \tilde{{\rm{\Omega }}}{c}_{0}{c}_{1}$ with ${c}_{1}={S}_{x0}-{S}_{z0}(\lambda {k}_{p}/\tilde{{\rm{\Omega }}})$. In equation (3. Regular spin dynamics of BB-type solitons in OLs
Note that the center-of-mass dynamics of the solitons pair is governed by equation (7a ) where the parameter kp is a function of spin expectation value. From equation (11 ) and equation (7a ), we find an effective equation15 ) describes a dependence of center-of-mass motion on the expectation value of spin components. It is seen that the center-of-mass is accelerated due to the force $-2\lambda \tilde{{\rm{\Omega }}}{S}_{y}$ produced by the SOC. The solitons also encounter a linear restoring force due to the OL potential. To illustrate the dynamics of the spin expectation value and the center-of-mass in detail, we numerically solve equations (10a )–(10c ) and (7a ) for θ(t = 0) = π/4 and φ(t = 0) = π/4. The results obtained for the spin dynamics are shown in the left panel of figure 1 while in the right panel, we show the time evolution of the center-of-mass. Here, we fix the initial values of all the spin components of the system on the Bloch sphere such that ${\sum }_{j=x,y,z}{S}_{j}^{2}=1$. The figure clearly shows that expectation values of the spin components oscillate periodically with time due to SOC. However, the periodicity and amplitudes of oscillations of different spin components are not the same. We see that the center-of-mass of the system also oscillates periodically.
$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}}^{2}\left\langle z\right\rangle }{{{\rm{d}}{t}}^{2}}=-2\lambda \tilde{{\rm{\Omega }}}{S}_{y}+\displaystyle \frac{2\lambda \tilde{{\rm{\Omega }}}{V}_{0}}{\alpha }\sin (\phi ).\end{eqnarray}$
Clearly, equation (
Figure 1. Top panel: Spin components of the solitons for the initial conditions θ(0) = π/4 and φ(t = 0) = π/4. In the left figure, the red line represents Sx, the blue line represents Sy and the green line represents Sz. The right panel gives the evolution of center-of-mass coordinate of the solitons. Here the initial state is chosen for α = 0.5, $\beta =0.5\sqrt{\alpha }$, c = −10, V0 = 0. Bottom panel: Same as those shown in the top panel but in the presence of OLs with V0 = −5. In both the panels, solid lines are obtained from variational calculations and circles are generated by direct numerical simulations of GPE. |
It is seen that the center-of-mass of solitons oscillates due to SO coupling and the expectation values of spin components move periodically in a closed orbit on the Bloch sphere. Since the spin dynamics can influence the center-of-mass motion and the latter is expected to be affected by the OLs, the oscillations of spin components are also expected to be affected by the external periodic potential. In order to find the effect of OL on the spin dynamics and hence also on the center-of-mass coordinate of the solitons, we take V0 = −5. In the bottom panel of figure 1, we show the variations of the expectation value of spin components (left panel) and of center-of-mass (right panel) for the interaction strength c = −10. It is seen that the amplitude and frequency of center-of-mass motion increase with the introduction of OL. The OL also causes a phase shift in the center-of-mass dynamics. The periodic dynamics of spin components are affected by lattice potential. Particularly, the amplitude of oscillation of only the x-component of spin changes significantly due to the OL.
In order to authenticate the results obtained from the variational approximation, we numerically solve the GPE in equation (1 ). The time-dependent GPE is discretized in space and time using the split-step Crank–Nicolson scheme. Real-time propagation is employed to solve the discretized equation to study the BEC dynamics by using sufficiently small space and time steps [25, 26]. Taking the ansatz of BB-type solitons in equation (4 ) as an initial condition (t = 0) for p1x = − p2x = λ and xc = 0 we calculate the time evolution of Sx, Sy, Sz and ⟨z⟩ both in absence and presence of OLs.
A comparison of results obtained from variational and direct numerical simulations is presented in figure 1 where the circles are generated by direct GPE numerical simulations. We see that the direct numerical simulations of the GPE agree very well with those obtained from variational calculations. Both the studies infer that the dynamics of spin expectation values and center-of-mass of the soliton pair are periodic in time. The numerical simulation of the GPE with SOC in presence of OLs also shows good agreement with the variational results. More specifically, modifications of amplitudes, phase and frequency of spin components and center-of-mass in OLs are clearly reflected in both cases.
We calculate the spatial variation of density distributions with time for both components from direct numerical simulation of the GPE. It is seen that the density varies periodically with time. The numerical simulation confirms that the density variation occurs due to the periodic interchange of particles between the components (figure 2). As a result, the value of spatial densities of two components alternately disappear and revive keeping the total density constant. Interestingly, the variation of total density with time is also periodic. We also note that, for a comparatively weaker atomic interaction, the solitons decay under the action of SOC.
Figure 2. Top panel: Evolution of densities of both components produced by the GPE simulations with the initial conditions θ(t = 0) = π/4, φ(t = 0) = π/4, α = 0.5, $\beta =0.5\sqrt{\alpha }$, c = −10 and V0 = 0. Here the panel marked by (a) gives the density (∣ψ∣2) of the first component while the panel (b) represents the density (∣φ∣2) of the second component. The total density (∣ψ∣2 + ∣φ∣2) variation is shown in the panel (c). Bottom panel: The same as that shown in the top panel but V0 = − 5. |
In presence of OLs, the density of the solitonary wave oscillates in space with a relatively larger frequency. During the motion densities of both the components vary with time and the variations are not independent of each other. Particularly, if the density is large in the first half period in one component then it is large in the second half period in the second component. In this way, the density of the two components alternately decays and revives with time. However, the total spatial density remains constant. This implies that spin flip-flop takes place periodically and causes variation in the population density of each component. The density varies largely with a larger frequency as the strength of the lattice increases. A striking feature is that the variation of spatial density shows an additional collective periodic variation. From equation (14 ) it is evident that the center-of-mass directly depends on the OL parameter. Therefore, an appropriate strength of lattice can also result in such an additional collective variation in the density of each component.
4. Nonlinear spin dynamics of BB-type solitons in OLs
The precession of soliton spin is described by the nonlinear Bloch equations, in which the nonlinear terms mostly result from SOC. We can see that the spin precession couples to the center-of-mass motion under the action of SOC. A stronger lattice potential results in an increase of effective magnetic field in the z-direction and thus augments the precession frequency. OL potential acts as a small perturbation to this nonlinearity which may have a significant impact on the nonlinear spin dynamics.
In order to study the influence of the lattice potential on the nonlinear spin dynamics of BB-type solitons, we solve the coupled nonlinear Bloch equation (12 ) numerically taking the initial condition at the fixed point Sc = (Sxc, Syc, Szc) of the system. The fixed point is determined by solving the Bloch equation at a steady state [27]. More specifically, we make a fixed point analysis for the spin components in the absence (V0 = 0) and presence (V0 ≠ 0) of OLs and display the outcome in figures 3–5. The phase plot of x, y and z components of the effective spin S for zero lattice strength is periodic with period 1. The period of phase-space trajectory increases with the increase of lattice strength. It roughly resembles the famous Rössler system where periodic doubling occurs with the increase of a control parameter [28]. In the present case, the lattice strength may be considered a control parameter. We see that the OL introduces periodic multiplicity and it increases with the increase of lattice strength. This ultimately results in chaos in the dynamics of spin-components. We also notice that prominent chaotic dynamics can be observed for a weaker lattice strength if the inter-atomic interaction is taken relatively stronger. In contrast, the trajectory of the chaotic pattern diverges dramatically with the rise of lattice potential while maintaining the same interaction strength.
Figure 3. Structure of phase space trajectory of spin component Sx for (i) c = −2, V0 = 0 (top left panel), (ii) c = −2, V0 = −1 (top right panel) (iii) c = −2, V0 = −2 (bottom left panel) and (iv) c = −4, V0 = −1 (bottom right panel). |
Figure 4. Structure of phase space of spin component Sy for (i) c = −2, V0 = 0 (top left panel), (ii) c = −2, V0 = −1 (top right panel), (iii) c = −2, V0 = −2 (bottom left panel) and (iv) c = −4, V0 = −1 (bottom right panel). |
Figure 5. Structure of phase-space of spin component Sz for (i) c = −2, V0 = 0 (top left panel), (ii) c = −2, V0 = −1 (top right panel), (iii) c = −2, V0 = −2 (bottom left panel) and (iv) c = −4, V0 = −1 (bottom right panel). |
The variation of spin expectation values, Sx, Sy and Sz in phase-space for different values of lattice and interaction parameters is displayed in figure 6. We see that the relative change of spin values in phase-space is double periodic in the absence of OLs (V0 = 0) and the periodic multiplicity increases with the increase of lattice strength (V0 ≠ 0). For a relatively stronger value of nonlinear interaction, the time evolution of spin projections in phase-space also becomes multi-periodic. This gives an indication to the chaotic dynamics of spin components.
Figure 6. Phase space trajectory of three spin components with the initial conditions Sx = −0.12, Sy = −1, Sz = −1. for (i) c = −2, V0 = 0 (top left panel), (ii) c = −2, V0 = −1 (top right panel), (iii) c = −2, V0 = −2 (bottom left panel) and (iv) c = −4, V0 = −1 (bottom right panel). |
In order to check our prediction from the phase-space pattern of spin dynamics, we calculate the Lyapunov exponent (LE) at different times both in the presence and absence of lattice potential for the given values of interaction strengths. The results are shown in figure 7. It shows that the LEs for all the spin components remain close to zero for V0 = 0 and thus phase-space trajectory is non-chaotic for given values of inter-atomic interaction. But in the presence of a weak lattice potential, the LE associated with the Sx and Sy components becomes positive. However, for the Sz component, the LE is always negative. Thus, we say that the effective spin dynamics become chaotic if the LE of any one or more spin components is positive.
Figure 7. Lyapunov exponent (Lλ) of spin components. The blue curve represents Sx while the red and black curves give Sy and Sz, respectively. The top left panel gives the result for (i) c = −2, V0 = −1 (solid curves) (ii) c = −2, V0 = 0 (dashed curves). The top right panel gives the results for c = −2, V0 = −2. The bottom left panel shows Lλ for c = −4, V0 = −1. The bottom right panel gives (Lλ) for c = −4, V0 = −2. |
It may be interesting to note that the temporal variation of the positive value of the LE corresponding to Sx and Sy grows rapidly after a certain evolution for relatively higher values of lattice depth and interaction strength. Thus the system becomes more chaotic.
5. Conclusion
We consider spin–orbit coupled BECs in OLs and study the effective spin dynamics within the framework of mean-field equations. More specifically, we take bright–bright solitons solution consisting of spin up and down particles and study both regular and nonlinear effective spin dynamics of the coupled system using variational and numerical approaches. We extract the spinor part of the solution and write a nonlinear Bloch equation for the effective spin and thus find an effective magnetic field in presence of OLs.
The spin angular momentum of the soliton changes by an effective magnetic field in accordance with the Bloch equation and generates a force that modifies the motion of the soliton's center-of-mass coordinate. We have demonstrated how the interplay among nonlinearity, OL, and SOC influences the soliton's spin expectation value and center-of-mass motion of solitons. It is seen that the OL significantly affects the dynamics of spin components and the center-of-mass motion. More specifically, the OL is found to generate an extra collective oscillation with a greater frequency and amplitude in addition to the movements of solitons generated by the SOC.
We have checked the effects of periodic variation of spin expectation value on the density profile through direct numerical simulation. It is seen that spin flipping and/or exchange occur periodically. As a result one of the states remains highly populated during half of the period while the other state remains less populated and vice versa in the next half of the period. The amplitude and frequency of population imbalance are found to be sensitively affected by the OL.
We have also examined the effective nonlinear spin dynamics in the presence of the OL potential using fixed-point analysis. It is seen that the lattice strength can be considered as a control parameter for changing the dynamics from regular to chaos. The spin dynamics exhibit regular periodic oscillation in the absence of the lattice potential. However, with the introduction of lattice potential, the dynamics show chaotic behavior. This is confirmed by the structure of the phase trajectory of spin components and the values of LEs.
Appendix
Equations for the dynamics of different parameters are obtained from the vanishing conditions of variational derivatives $\left\langle { \mathcal L }\right\rangle $ with respect to different parameters. From $\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta {\phi }_{1}}-\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta {\phi }_{2}}=0$ we get
$\begin{eqnarray}\dot{\theta }=-\displaystyle \frac{\alpha \pi \lambda }{a\sinh (\pi \lambda /a)}\sin \phi .\end{eqnarray}$
Again from $\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta {p}_{1x}}+\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta {p}_{2x}}=0$ we get $\begin{eqnarray}{k}_{p}=({k}_{n}-\beta )\cos 2\theta +\displaystyle \frac{{\rm{d}}}{{\rm{d}}{t}}(-{x}_{c}/a).\end{eqnarray}$
and the use of $\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta {p}_{1x}}-\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta {p}_{2x}}=0$ gives $\begin{eqnarray}\displaystyle \frac{\pi \lambda }{a}\coth \left(\displaystyle \frac{\pi \lambda }{a}\right)-1=0.\end{eqnarray}$
Similarly, from $\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta {x}_{c}}=0,$
$\begin{eqnarray}{\dot{k}}_{p}=\displaystyle \frac{2\pi \alpha {\lambda }^{2}}{a\sinh (\pi \lambda /a)}\sin 2\theta \sin \phi +\displaystyle \frac{2\pi {k}_{{\rm{lat}}}^{2}{V}_{0}}{a\sinh (\pi {k}_{{\rm{lat}}}/a)}\sin \phi \end{eqnarray}$
and from $\tfrac{\delta \left\langle { \mathcal L }\right\rangle }{\delta \theta }=0,$ $\begin{eqnarray}\dot{\phi }=2{k}_{n}{k}_{p}-\displaystyle \frac{2\pi \alpha {k}_{n}\cos \phi }{a\sinh (\tfrac{\pi {k}_{n}}{a})}\cot 2\theta .\end{eqnarray}$