Soliton lattices in the Gross-Pitaevskii equation with nonlocal and repulsive coupling
Hidetsugu Sakaguchi



Department of Applied Science for Electronics and Materials,
Interdisciplinary Graduate School of Engineering Sciences, Kyushu
University, Kasuga, Fukuoka 816-8580, Japan

Spatially-periodic patterns are studied in nonlocally coupled Gross-Pitaevskii equation.  We show first that spatially periodic patterns appear in a model with the dipole-dipole interaction. Next, we study a model with a finite-range coupling, and show that a repulsively coupled system is closely related with an attractively coupled system and its soliton solution becomes a building block of the spatially-periodic structure.  That is, the spatially-periodic structure can be interpreted as a soliton lattice. An approximate form of the soliton is given by a variational method. Furthermore, the effects of the rotating harmonic potential and spin-orbit coupling are numerically studied. 

Soliton lattices in the Gross-Pitaevskii equation with nonlocal and repulsive coupling
    Hidetsugu Sakaguchi
    December 30, 2023
======================================================================================



§ INTRODUCTION

Supersolid is an interesting material which is characterized with both the superfluidity and crystal structure <cit.>. The superfluidity in solid ^4He has been studied as an example of supersolid. Recently, the supersolid phases  are studied in Bose-Einstein condensates (BECs). The Gross-Pitaevskii equation is  a mathematical model of matter waves of BECs <cit.>. The GP equation is derived by a mean-field approximation of weakly interacting bosons at zero temperature. The effect of multi-body correlation and the fluctuation by finite temperature are neglected in the description by the GP equation. Various types of ground states  were successfully described by the GP equations. 
For example, a vortex lattice state was observed in rotating BECs with repulsive interaction <cit.> and was described by the GP equation <cit.>. Localized matter waves called solitons were found in experiments of BECs <cit.> and the solitons were expressed as a ground state in the GP equation for Bose-Einstein condensate with attractive interaction <cit.>. Pomeau and Rica studied a nonlocally coupled Gross-Pitaevskii (GP) equation as a model of supersolid <cit.>.
They showed the excitation spectrum exhibits the so-called roton minimum when a parameter exceeds a critical value which leads to the spatially periodic structure. Henkel, Nath, and Pohl studied the supersolid in Rydberg-dressed BECs <cit.>. Since then, the supersolid in the Rydberg-dressed BECs has been studied by many authors <cit.>. The GP equation with the finite-range interaction was also studied by several authors <cit.>. Prestipino et al. studied the spatially-periodic solution of the GP equation with a variational method <cit.>. Various other effects can be incorporated in the nonlocal GP equation. The effect of rotation to the Rydberg-dressed BECs was studied by Henkel et. al. <cit.>, and the effect of spin-orbit coupling was studied by Han et al. <cit.>. Cinti showed the supersolid phase using a Monte-Carlo simulation of N soft-core bosons <cit.>.

For a certain parameter sets, a strongly localized structure appears in the spatially periodic solution. In this paper, we would like to propose a viewpoint that the localized structure can be interpreted as a soliton to a GP equation with an attractive interaction, and the spatially periodic solution can be interpreted as a soliton lattice. We will show that the localized solution can be approximated as a soliton to the nonlinear Schrödinger equation at a certain parameter range in a one-dimensional system. 
In Sec.2, we show spatially periodic patterns in a model with dipole-dipole interaction which might be more easily realized in experiments, although similar results were already shown in the GP equation for the Rydberg-dresses BECS. In Secs.3 and 4, we study spatially-periodic patterns in a simpler finite-range interaction model. The model equation with the repulsive interaction is closely related with an attractively coupled system. A soliton appears naturally in the attractively coupled system. The soliton is a building block of the spatially-periodic structure in the repulsive system.  Furthermore, an approximate solution of the soliton is given by a variational method for the attractive GP equation. In Sec.5, spin-orbit coupled systems are studied, in which semi-vortex lattices are created. We show that the vortex structure is easily broken in a rotating harmonic potential. This type of spin-orbit coupled system under rotation has not been studied before. 


§ SPATIALLY PERIODIC PATTERN IN THE ONE-DIMENSIONAL GROSS-PITAEVSKII EQUATIONS WITH DIPOLE-DIPOLE INTERACTION

The one-dimensional nonlocally coupled Gross-Pitaevskii equation is expressed as

    i∂ϕ/∂ t =-β/2∇^2ϕ+{∫ G(x-x^')|ϕ(x^')|^2dx^'}ϕ,

where the Planck constant is set to be 1, m_0=1/β denotes a parameter proportional to particle mass. The second term on the right-hand side represents the repulsive and nonlocal interaction, where G(x-x^') is an integral kernel characterizing the nonlocal coupling. 
For the Rydberg-dressed BECs, the integral kernel takes a form G(x-x^')=g/[ϵ+(x-x^')^6/R^6]. For the dipole-dipole interaction, G(x-x^')=g/[ϵ+|x-x^'|^3/R^3].  
There is a uniform solution to Eq. (1): ϕ(x)=Ae^-iμ t where μ=A^2∫ G(x)dx. The modulational instability of the uniform solution can be investigated by assuming the wave function as ϕ=A(1+δ(x))e^iμ t+iθ(x)+λ t. The eigenvalue λ(k) for the perturbation with wavevector k is given by 

    λ^2(k)=-β k^2/2 (β k^2/2+G(k)A^2 ),

where  G(k)=∫ G(x)e^ikxdx.
If λ^2(k)<0, λ is pure imaginary and the uniform state is stable against the perturbation, however, if λ^2(k)>0, λ takes a real positive value and the uniform state becomes unstable. This corresponds to the roton softening in the excitation spectrum <cit.>.

As an example, we show the instability in the case of the dipole-dipole interaction, since the nonlocal GP equation has been mainly studied for the Rydberg-dressed BECs and the finite-range interaction. Figure 1(a) shows λ^2(k) for β=0.1 and the dipole-dipole interaction with g=1, ϵ=1, and R=5.  The eigenvalue becomes positive for k∼ 1. Figure 1(b) shows a ground state of Eq. (2) in one dimension obtained by the imaginary-time simulation method. In the imaginary-time evolution,  the numerical simulation of the equation obtained by replacing t with -iτ from Eq. (1) is performed under the condition of the constant norm N=∫_0^L|ϕ|^2dx. Generally, the ground state is obtained by the imaginary-time evolution method. The total norm  N is fixed to be 30. The system size L is 30 and the periodic boundary conditions are imposed at the boundaries x=-L/2 and L/2. A periodic pattern of wavenumber L/5=6 appears as a result of the instability. Figure 1(c) shows the real time evolution of |ϕ| by Eq. (1) after t=1000 starting from the initial condition ϕ(x)=1+r(x) where r(x) is a random number between -0.1 and 0.1. A spatially periodic pattern with five peaks appears from the randomly perturbed uniform initial condition even in the real time evolution of Eq. (1).  



§ SOLITON LATTICES IN THE GROSS-PITAEVSKII EQUATION WITH FINITE-RANGE COUPLING

Hereafter, we consider a simple finite-range coupling expressed as G(x-x^')=g for |x-x^'|≤ L_g and 0 for |x-x^'|> L_g.  This function can be expressed as a limit of g/(1+|x-x^'|^n/L_g^n) including both the dipole-dipole interaction and the interaction in the Rydberg dressed BECs, i.e., G(x-x^')=lim_n→∞g/[1+|x-x^'|^n/L_g^n]. Furthermore, the parameter β is set to 1. The model equation is rewritten as 
    i∂ϕ/∂ t =-1/2∂^2ϕ/∂ x^2+g{∫_x-L_g^x+L_g|ϕ(x^')|^2dx^'}ϕ.

Figure 2(a) shows a spatially-periodic structure obtained by the imaginary-time evolution method at g=5, L=20, L_g=5.5, and N=1. Rhombi in Fig. 2(b) show the peak amplitudes of periodic solutions as a function of L_g for g=5,L=20, and N=1. Spatially-periodic solutions of wavelength L/3 appears at L_g=3.8 when L_g increases.  The periodic solution of mode m=3 changes to a solution of mode m=2 at L_g=5.6, and the one-hump solution of m=1 appears at L_g=10.1. 
  
Equation (3) can be rewritten as

    i∂ϕ/∂ t   =   -1/2∂^2ϕ/∂ x^2+g{∫_x-L^x+L|ϕ(x^')|^2dx^'-∫_x-L^x-L_g|ϕ(x^')|^2dx^'-∫_x+L_g^x+L|ϕ(x^')|^2dx^'}ϕ
       =   -1/2∂^2ϕ/∂ x^2+2gNϕ-g{∫_x-(L-L_g)^x+(L-L_g)|ϕ(x^')|^2dx^'}ϕ,

where the total norm is N=∫_0^L|ϕ|^2dx and the spatial periodicity ϕ(x+L)=ϕ(x-L)=ϕ(x) is used.
Note that the nonlocal repulsive interaction is changed into an attractive interaction. The BECs tend to be localized owing to the attractive interaction. 
This is a reason why a soliton-like localized structure appears.  The soliton-like structure for mode m can be expressed as a solution to  

    i∂ϕ/∂ t=-1/2∂^2ϕ/∂ x^2-g{∫_x-(L^'-L_g)^x+(L^'-L_g)|ϕ(x^')|^2dx^'}ϕ,

in a system of size L^'=L/m, and the norm is N^'=N/m.  Here, the term 2gNϕ in Eq. (4) is neglected because the term can be removed by the transformation ϕ e^-i 2gNt→ϕ, and the term is irrelevant to the spatial structure of the solution. The solid line in Fig. 2(b) shows the peak amplitude of the one-hump solution to the nonlocally-coupled attractive model equation Eq. (5) for m=1,2,, and 3. Good agreement is seen between the numerical results (rhombi) by the direct numerical simulations of Eq. (3) and those for Eq. (5) (solid lines).      

The spatially periodic solution to Eq. (3) can be expressed as m trains of  one-hump solutions to Eq. (5). 
When L_g is changed as a control parameter, the peak amplitude of the soliton takes the maximum at a value slightly smaller than L/m and decreases rapidly to zero near L/m for each m. A strongly localized structure appears near the maximum point. We call the localized structure a soliton in this paper. The soliton generally appears in the GP equation with attractive interaction <cit.>. The soliton appears because the originally repulsive GP equation (3) can be transformed to an attractive GP equation (5). 

If (L-L_g) is sufficiently small, Eq. (5) for m=1 is approximated at 

    i∂ϕ/∂ t=-1/2∂^2ϕ/∂ x^2-2g(L-L_g)|ϕ|^2ϕ.

This is equivalent to the nonlinear Schrödinger equation. The soliton solution to Eq. (6) is expressed as

    ϕ=A/cosh Bx,

where 2A^2/B=N and 2g(L-L_g)A^2=B^2.
The solid line in Fig. 2(c) shows one-hump solution to Eq. (5) for m=1 at L_g=19.6. The dashed line in Fig. 2(c) denotes the soliton solution 1/cosh(2x) to Eq. (6). The pulse solution is well approximated by the soliton solution. 
The solid line in  Fig. 2(d) shows the peak amplitude for the one-hump solution of Eq. (5) for m=1 as a function of L_g. The dotted line denotes A=√(2.5(L-L_g)) for the soliton solution. Good agreement is observed for L_g>19.4. However, the peak amplitude increases with L_g for L_g<19.4. 

The total energy E of the solution to Eq. (5) for m=1 is expressed as

    E=1/2∫ |∂ϕ/∂ x |^2dx-g/2∫_0^L∫_x-(L-L_g)^x+(L-L_g)|ϕ(x)|^2|ϕ(x^')|^2dxdx^'.

If the Gaussian function is assumed for ϕ, i.e., ϕ=Ae^-α x^2, E is expressed as 

    E(α)=N/2α-2/√(π)gN^2∫_0^√(2α)L_ge^-z^2dz,

where N=√(π/(2α))A^2. The variational approximation gives α from the minimization of E(α) and A is obtained from A=(N√(2α/π))^1/2. The dashed line in Fig. 2(d) shows the relationship between A and L_g obtained by the variational method, which is a good approximation for L_g<18.2.  




§ TWO-DIMENSIONAL GROSS-PITAEVSKII EQUATION
 
Next, we study the two dimensional nonlocally coupled GP equation. The two-dimensional integral kernel is expressed as  G(x-x^',y-y^')=g only for |x-x^'|≤ L_g and |y-y^'|≤ L_g. The model equation is written as

    i∂ϕ/∂ t=-1/2∇^2ϕ+g{∫∫_S|ϕ(x^',y^')|^2dx^'dy^'}ϕ,
 
where the integral region S satisfies  |x-x^'|≤ L_g and |y-y^'|≤ L_g. Figure 3(a) shows a spatially periodic stationary solution of wavelength L/3 at g=20, L=20, L_g=5, and N=16. The region satisfying |ϕ|>0.3 is plotted in this figure. The rhombi in Fig. 3(b) shows the peak amplitude as a function of L_g for  g=20, L=20, and N=16. Spatially-periodic states with wavelength L/m (m=1,2,3, and 4) appears by changing L_g. 
The solid line shows the peak amplitude of the corresponding attractive model equations in a system of (L/m)× (L/m) for m=1,2,3, and 4:

    i∂ϕ/∂ t=-1/2∇^2ϕ-g{∫∫_S^'|ϕ(x^',y^')|^2dx^'dy^'}ϕ,
 
where S^' is a region satisfying |x-x^'|>L_g/m or |y-y^'|>L_g/m, and the norm is N^'=N/m^2.  
Good agreement is seen in numerical results for Eqs. (10) and (11), although discontinuous mode transitions occur in the direct numerical simulation of Eq. (10). Intermediate states sometimes appear near the transition of modes. An example of the intermediate state between m=3 and 4 at L_g=4.25 is shown in Fig. 3(c).  

The total energy E of Eq. (11) for m=1 is expressed as

    E=1/2∫ |∇ϕ |^2dx-1/2∫_0^L∫_0^L ∫∫_S^' |ϕ(x,y)|^2|ϕ(x^',y^')|^2dx^'dy^'dxdy.

If the Gaussian function is assumed for ϕ, i.e., ϕ=Ae^-α (x^2+y^2), E is expressed as 

    E(α)=Nα-2/πgN^2{I_S^2+2(0.5-I_S)I_S},

where I_S=∫_0^√(α)L_ge^-z^2dz, and N=π A^2/(2α). The variational approximation gives α from the minimization of E(α) and A is obtained from A=(2Nα/π)^1/2. The solid line in Fig. 3(d) is direct numerical results of Eq. (11) for m=1. The dashed line in Fig. 3(d) shows the relationship between L_g and A obtained by the variational method, which is a good approximation for L_g<18. The peak amplitude takes the maximum value at a certain L_g slightly smaller than L/m, and the peak amplitude decreases rapidly and becomes 0 near L/m. That is, a very strongly localized structure or a two-dimensional soliton of large amplitude appears near L_g=L/m. The behavior is similar to the one-dimensional system, however, the change of the peak amplitude near L/m is even larger in the two-dimensional system. 
However, there is no collapse phenomenon in the nonlocally coupled equation Eq. (11) with attractive interaction in contrast to the locally coupled two-dimensional nonlinear Schrödinger equation.   


If the BECs are trapped in a rotating harmonic potential, the GP equation is expressed with

    i∂ϕ/∂ t =-1/2∇^2ϕ+g{∫∫_S|ϕ(x^',y^')|^2dx^'dy^'}ϕ+1/2K(x^2+y^2)ϕ-Ω L_zϕ,

where L_z=-i(x∂/∂ y-y∂/∂ x) and Ω is the frequency of rotation. When K=Ω^2, Eq. (14) becomes

    i∂ϕ/∂ t =- [1/2 (i∂/∂ x-Ω y )^2+1/2 (i∂/∂ y+Ω x )^2 ]ϕ+g{∫∫_S|ϕ(x^',y^')|^2dx^'dy^'}ϕ.

If g=0, this equation is equivalent to the 2D Schrödinger equation of a charged particle in uniform magnetic field Ω directed perpendicular to the (x,y) plane. For the locally coupled Gross-Pitaevskii equation with g>0, a vortex lattice state appears <cit.>. Figure 4(a) shows a spatially-periodic state of wavelength L/4  at Ω=1, g=20, L=20, L_g=3.5, and N=16. In Fig. 4(a), |ϕ|>0.3 is satisfied in the green region and Reϕ>0.3 is satisfied in the blue region. The wave function ϕ changes in the azimuthal direction like ϕ(x,y)∼ e^ik(r)θ where θ is the angle around the center (0,0).  Figure 4(b) shows the numerically obtained wavenumber k(r) at several positions near the peaks of |ϕ|  as a function of the distance from the center at Ω=1. The dashed lines denote k(r)=Ω r. The wavenumber k(r) increases with r and the proportionality coefficient is equal to Ω. The term of -Ω L_z in Eq. (15) induces the clockwise rotation of frequency Ω, which is compensated by the anti-clockwise rotation induced by the phase factor e^ik(r)θ, and a stationary state appears.   
We can also interpret that this state is a vortex lattice state where vortices are expelled from the regions of the soliton lattice to the region of |ϕ|∼ 0 and only the phase factor survives. 




§ SPIN-ORBIT COUPLED BOSE-EINSTEIN CONDENSATES

Spin-orbit coupled BECs have been intensively since Lin et al. realized the spin-orbit coupling in BECs experimentally <cit.>. In the spin-orbit coupled BECs, the spin-up and spin-down components interact with each other through the spin-orbit coupling k_yσ_x-k_xσ_y where k_x and k_y are the momentum operator in the x and y directions and σ_x,y are the Pauli matrix. The spin-orbit coupled GP equations have been studied by several authors <cit.>. Many remarkable nonlinear phenomena such as stripe patterns and solitons have been found in the spin-orbit coupled GP equations. Sakaguchi et al. found a vortex soliton called a semi-vortex in two-dimensional spin-orbit coupled GP equation <cit.>. In the semi-vortex, the spin-up component takes a form of a two-dimensional soliton and the spin-down component takes a form of localized vortex.   
 A nonlocally coupled system with spin-orbit interaction in two dimensions is expressed as

    i∂ϕ_+/∂ t   =   -1/2∇^2ϕ_++g{∫∫_S(|ϕ_+|^2+|ϕ_-|^2)dx^'dy^'}ϕ_++λ (∂ϕ_-/∂ x-i∂ϕ_-/∂ y ),
    
    i∂ϕ_-/∂ t   =   -1/2∇^2ϕ_-+g{∫∫_S(|ϕ_-|^2+|ϕ _+|^2)dx^'dy^'}ϕ_--λ (∂ϕ_+/∂ x+i∂ϕ_+/∂ y ),

where λ is the strength of the spin-orbit coupling and S is the integral region satisfying |x-x^'|<L_g and |y-y^'|<L_g. 


Figures 5(a) and (b) show the regions satisfying |ϕ_+|>0.2 and |ϕ_-|>0.2 for g=20, L=20, L_g=5, N=16, and λ=1. A lattice of semi-vortices with wavelength L/3 appears <cit.>. A vortex exists at each center of annular structure of |ϕ_-|. Similarly to the previous section, each semi-vortex can be expressed by the attractively interacting model equation:

    i∂ϕ_+/∂ t   =   -1/2∇^2ϕ_+-g{∫∫_S^' (|ϕ_+|^2+|ϕ_-|^2)dx^'dy^'}ϕ_++λ (∂ϕ_-/∂ x-i∂ϕ_-/∂ y ),
    
    i∂ϕ_-/∂ t   =   -1/2∇^2ϕ_--g{∫∫_S^'(|ϕ_-|^2+|ϕ _+|^2)dx^'dy^'}ϕ_--λ (∂ϕ_+/∂ x+i∂ϕ_+/∂ y ),

where S^' is a region satisfying |x-x^'|>L_g/m or |y-y^'|>L_g/m. 
Figure 5(c) compares the peak amplitudes A of |ϕ_+| of the semi-vortex lattice in Eq. (16) with the peak amplitudes  of |ϕ_+|  obtained by Eq. (17). Good agreement is seen. It is noted that discontinuous transition from m=3 to m=2 modes occurs near L_g=6.2 in Eq. (16). These results are similar to the ones shown in Fig. 4. 


If the spin-orbit coupled BECs are trapped in the rotating harmonic potential satisfying K=Ω^2, the coupled GP equations are expressed as <cit.>

    i∂ϕ_+/∂ t   =   - [1/2 (i∂/∂ x-Ω y )^2+1/2 (i∂/∂ y+Ω x )^2 ]ϕ_+
          +g{∫∫_S (|ϕ_+|^2+|ϕ_-|^2)dx^'dy^'}ϕ_++λ (∂ϕ_-/∂ x-i∂ϕ_-/∂ y ),
    
    i∂ϕ_-/∂ t   =   - [1/2 (i∂/∂ x-Ω y )^2+1/2 (i∂/∂ y+Ω x )^2 ]ϕ_-
          +g{∫∫_S(|ϕ_-|^2+|ϕ _+|^2)dx^'dy^'}ϕ_--λ (∂ϕ_+/∂ x+i∂ϕ_+/∂ y ).

The semi-vortex lattice disappears when Ω increases.
Figures 6(a)-(c) show regions satisfying |ϕ_-|>0.1 at (a) Ω=0.001, (b) 0.002, and (c) 0.004 for L=20, g=20, L_g=3.5, λ=1, and N=16. 
Spatially-periodic semi-vortex lattices of wavelength L/4 appear at L_g=3.5. When Ω=0, the vortices of ϕ_- locate at the center of each localized structure. The vortex positions move toward the origin (0,0) of the whole system. The movement is larger for the semi-vortices which is more distant from the origin. At Ω=0.001, the vortices near the four corners go out of the shaded regions representing the localized structure. At Ω=0.002, only the four vortices near the center survive.
 At Ω=0.004, all vortices go out from the shaded regions. That is, the semi-vortex structure is broken at relatively small Ω. The peak amplitude of |ϕ_+| and |ϕ_-| take a similar value at Ω=0.02.    



§ SUMMARY
 
First, we have shown that a spatially-periodic structure can be created in a system with the dipole-dipole type nonlocal coupling. BECs with the dipole-dipole interaction might be more easily realized than the Rydberg-dressed Bose-Einstein condensates studied previously. Next, we have studied a model system with a finite-range nonlocal interaction more in detail. A spatially-periodic pattern like a square lattice appears in the repulsively coupled system in two dimensions owing to the anisotropy. Dynamical transitions between spatially periodic structures with different wavenumber occur when the interaction range L_g is changed. 
We have found that the peak amplitude takes a sharp maximum slightly below L_g=L/m (m integer) when L_g is changed as a control parameter. The spatially periodic structure can be interpreted as a soliton lattice, in that the soliton structure is expressed by a solution to an attractively coupled system and works as a building block of the spatially-periodic pattern. Approximate solutions were derived by the variational method. In a rotating harmonic potential, a phase wave is induced in the azimuthal direction as ϕ∼ e^ik(r)θ.  The wavevector along the azimuthal direction increases in proportion to the distance from the center of the harmonic potential, which implies the rigid rotation. In the spin-orbit coupled nonlocal model, semi-vortex lattices are created, however, the vortex structure is easily broken in the rotating harmonic potential.      99
Andreev A. F. Andreev and I. M. Lifshitz, Sov. Phys. JETP 29 (1969) 1107.
Chester G. V. Chester, Phys. Rev. A 2 (1970) 256.
Pethick C. J. Pethick and H. Smith, Bose-Einstein Condensation in Dilute Gases (Cambridge University Press, Cambridge, 2002)
Stringari L. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Oxford University Press, Oxford, 2003),
Williams J. E. Williams and M. J. Holland, Nature 401 (1999) 568.
Matthews M. R. Matthews, B. P. Anderson, P. C. Haljan, D. S, Hall, C. E. Wieman, and E. A. Cornell, Phys. Rev. Lett. 83, 2498 (1999).
Tsubota M. Tsubota, K. Kasamatsu, and M. Ueda, Phys. Rev. A 65 (2002) 023603.
Sakaguchi1 H. Sakaguchi and H. Takeshita, J. Phys. Soc. Jpn. 77  (2008) 054003.
Khaykovich L. Kaykovich, F. Schreck, G. Ferrari, T. Boudel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, Science 296 (2002) 1290.
Malomed B. Malomed, L. Toner, F. Wise, and D. Mihalache, J. Phys. B 49 (2016) 170502.
Pomeau Y. Pomeau and S. Rica, Phys. Rev. Letts. 72 (1994) 2426.
Henkel N. Henkel, R. Nath, and T. Pohl, Phys. Rev. Letts. 104 (2010) 195302.
Maucher F. Maucher, N. Henkel, M. Saffman, W. Królikowski, S. Skupin, and T. Pohl, Phy. Rev. Letts. 106 (2011) 17401.
Hsueh C. H. Hsueh, T.-C. Lin, T.-L. Horng, and W. C. Wu, Phys. Rev. A 86 (2012) 013619.
Sepulveda N. Sepúlveda, C. Joosserand, and S. Rica, Eur. Phys. J. B 78 (2010) 439.
Kunimi M. Kunimi and Y. Kato, Phys. Rev. B 86 (2012) 060510R.
Ancilotto F. Ancilotto, M. Rossi, abd F. Toigo, Phys. Rev. A 88  (2013) 033618.
Macri T. Macri, F. Maucher, F. Cinti, and T. Pohl, Phys, Rev. A 87 (2013) 061602R.
Prestipino S. Prestipino, A. Sergi, and W. Bruno, Phys. Rev. B 98 (2018) 104104.
Henkel2 N. Henkel, F. Cinti, P. Jain, G. Pupillo, and T. Pohl, Phys. Rev. Letts. 105 (2012) 265301. 
Han W. Han, X F .Zhang, D.-S. Wang, H .F. Jian, W. Zhang, and S -G. Zhang, Phys. Rev. Letts 121 (2018) 0330404.
Cinti1 F. Cinti, P. Jain, M. Boninsegni, A. Micheli, P. Zoller, and G. Pupillo, Phys. Rev. Lett. 105 (2010) 135301.
Cinti2 F. Cinti, T. Marci, W. Lechner, G. Pupillo, and T. Tohl, Nature Commum. 5 (2014) 3235.
Lin Y. J. Lin, Jimenez-Garcia, and I. B. Spielman, Nature 471 (2011) 83. 
Zhai H. Zhai, Int. J. Mod. Phys. B 26 (2012) 1230001. 
Malomed2 B. A. Malomed, EPL 122 (2018) 36001. 
Sakaguchi2 H. Sakaguchi, B. Li, and B. A. Malomed, Phys. Rev. E 89 (2014) 032920.
Sakaguchi3 H. Sakaguchi and K. Umeda, J. Phys. Soc. Jpn. 85 (2016) 064402.

