The experimental realization of the artificial Abelian or non-Abelian gauge potential in neutral atoms^[1,2] gives rise to the spin–orbit coupling (SOC) effect in the ultracold atoms. The one-dimensional (1D) and two-dimensional (2D) SOC effects have been experimentally realized in Bose–Einstein condensates (BECs) or ultracold Fermi gases.^[3–8] It is well-known that the plane wave phase, the stripe phase, and the meron with dipole–dipole repulsion^[9–13] can be generated in spin–orbit coupled spin-1/2 BECs depending on the relative magnitude of intraspecies and interspecies interactions. The realization of SOC in experiment also brings a completely new avenue for the study of the many-body dynamics of the spin–orbit coupled BECs. The dark and bright solitons, quantized vortices dynamics, and Josephson dynamics are also studied in spin–orbit coupled spin-1/2 BECs.^[14–20]

In the spinor BECs, the complex atomic interactions give rise to many exotic ground states. For the spin-1 BECs, the mean-field ground state is a polar state with repulsive spin–spin interaction (such as ²³Na) and ferromagnetic state with attractive spin–spin interaction (such as ⁸⁷Rb).^[21] For the spin-2 BECs, due to the density–density, spin–spin, and spin–singlet pair terms, the mean-field ground state exhibits a rich phase diagram mainly comprised of the ferromagnetic, uniaxial nematic (UN), biaxial nematic (BN), and cyclic states.^[22–24] The spin–orbit coupled BECs with spin larger than 1/2 also have been studied in theory and experiment.^[8,25] Recently, the spin–orbit coupled spinor BECs have been attained through Raman coupling among the hyperfine states or the gradient magnetic field in experiment,^[26,27] which provides an interesting avenue to explore the high-spin systems with the SOC effect. The exotic ground states, phase transitions, spin squeezing, and solitons are revealed in the spin–orbit coupled spin-1 BECs.^[28–38] The square lattice phase, triangular lattice phase, and vector solitons are found in the spin–orbit coupled spin-2 BECs.^[39–44] In the spin–orbit coupled spin-1 BECs, most researches restricted on the SU(2) type SOC, i.e., the internal states are coupled to their momenta via the SU(2) Pauli matrices. However, the SU(2) spin matrices cannot couple the internal states |1〉 and |–1〉 in the spin-1 BECs directly. The SU(3) SOC with the spin operator spanned by the Gell–Mann matrices can completely describe the internal couplings of the spin–orbit coupled spin-1 Bose gases. The spiral spin textures, the double-quantum spin vortices, solitons, and threefold-degenerate plane wave with nontrivial spin textures were predicted in the SU(3) spin–orbit coupled Bose gases.^[45–49]

In this paper, we investigate the ground-state phases of isotropic and anisotropic SU(3) spin–orbit coupled spin-1 BECs in a 2D harmonic trap respectively. The competition between the SU(3) SOC and the spin–spin interaction results in a rich variety of lattice configurations. The ground-state phase diagram spanned by the isotropic SU(3) SOC and the spin–spin interaction is presented. Five ground-state phases can be identified on the phase diagram, including the plane wave (PW) phase, the stripe (ST) phase, the kagome lattice (KL) phase, the stripe-honeycomb lattice (SHL) phase, and the honeycomb hexagonal lattice (HHL) phase. The system undergoes a sequence of phase transitions from the rectangular lattice (RL) phase to the HHL phase, and to the triangular lattice (TL) phase in the spin-1 BECs with anisotrpic SU(3) SOC.

The paper is organized as follows. In Section 2, we introduce the model of 2D SU(3) spin–orbit coupled spin-1 BECs in a harmonic trap. In Section 3, we display the ground-state phases of 2D spin-1 BECs with isotropic and anisotropic SU(3) SOCs, respectively. The PW phase, the ST phase, the KL phase, the SHL phase, and the HHL phase are found with isotropic SU(3) SOC in Subsection 3.1. The RL phase, the HHL phase, and the TL phase are found with anisotropic SU(3) SOC in Subsection 3.2. A summary is included in Section 4.

We study the ground-state phases of the SU(3) spin–orbit coupled spin-1 BECs in a 2D harmonic trap. The expectation value of the Hamiltonian is given as

$$ \begin{eqnarray}\begin{array}{ll}E[\varPsi ] & \equiv \langle \hat{H}\rangle =\displaystyle \int {\rm{d}}{\boldsymbol{r}}{\varPsi }_{j}^{\ast }\left\{\left(-\displaystyle \frac{{\hslash }^{2}}{2m}{\nabla }^{2}+V({\boldsymbol{r}})+{\upsilon }_{{\rm{soc}}}\right){\varPsi }_{j}\right.\\ & \quad \left.+\displaystyle \frac{1}{2}{c}_{0}{n}^{2}+\displaystyle \frac{1}{2}{c}_{2}|{\boldsymbol{F}}{|}^{2}\right\},\end{array}\end{eqnarray}$$ (1)

The time evolution of the mean field is governed by

$$ \begin{eqnarray}{\rm{i}}\hslash \displaystyle \frac{\partial {\varPsi }_{j}}{\partial t}=\displaystyle \frac{\delta E}{\delta {\varPsi }_{j}^{* }}.\end{eqnarray}$$ (3)

$$ \begin{eqnarray}\begin{array}{lll}{\rm{i}}\hslash \displaystyle \frac{\partial {\varPsi }_{-1}}{\partial t} & = & \left(-\displaystyle \frac{{\hslash }^{2}}{2m}{\nabla }^{2}+V({\boldsymbol{r}})\right){\varPsi }_{-1}+({c}_{0}n-{c}_{2}{F}_{z}){\varPsi }_{-1}\\ & & +\displaystyle \frac{{c}_{2}}{\sqrt{2}}{F}_{+}{\varPsi }_{0}-\hslash \left({\gamma }_{y}\displaystyle \frac{\partial }{\partial y}+{\rm{i}}{\hslash }_{x}\displaystyle \frac{\partial }{\partial x}\right){\varPsi }_{0}\\ & & +\hslash \left({\gamma }_{y}\displaystyle \frac{\partial }{\partial y}-{\rm{i}}{\hslash }_{x}\displaystyle \frac{\partial }{\partial x}\right){\varPsi }_{1};\\ {\rm{i}}\hslash \displaystyle \frac{\partial {\varPsi }_{0}}{\partial t} & = & \left(-\displaystyle \frac{{\hslash }^{2}}{2m}{\nabla }^{2}+V({\boldsymbol{r}})\right){\varPsi }_{0}+({c}_{0}n){\varPsi }_{0}+\displaystyle \frac{{c}_{2}}{\sqrt{2}}{F}_{-}{\varPsi }_{-1}\\ & & +\displaystyle \frac{{c}_{2}}{\sqrt{2}}{F}_{+}{\varPsi }_{1}+\hslash \left({\gamma }_{y}\displaystyle \frac{\partial }{\partial y}-{\rm{i}}{\hslash }_{x}\displaystyle \frac{\partial }{\partial x}\right){\varPsi }_{-1}\\ & & -\hslash \left({\gamma }_{y}\displaystyle \frac{\partial }{\partial y}+{\rm{i}}{\hslash }_{x}\displaystyle \frac{\partial }{\partial x}\right){\varPsi }_{1};\\ {\rm{i}}\hslash \displaystyle \frac{\partial {\varPsi }_{1}}{\partial t} & = & \left(-\displaystyle \frac{{\hslash }^{2}}{2m}{\nabla }^{2}+V({\boldsymbol{r}})\right){\varPsi }_{1}+({c}_{0}n+{c}_{2}{F}_{z}){\varPsi }_{1}\\ & & +\displaystyle \frac{{c}_{2}}{\sqrt{2}}{F}_{-}{\varPsi }_{0}-\hslash \left({\gamma }_{y}\displaystyle \frac{\partial }{\partial y}+{\rm{i}}{\hslash }_{x}\displaystyle \frac{\partial }{\partial x}\right){\varPsi }_{-1}\\ & & +\hslash \left({\gamma }_{y}\displaystyle \frac{\partial }{\partial y}-{\rm{i}}{\hslash }_{x}\displaystyle \frac{\partial }{\partial x}\right){\varPsi }_{0},\end{array}\end{eqnarray}$$ (4)

The momentum distribution of the ground-state phase can be derived from the Fourier transformation

$$ \begin{eqnarray}|{\varPsi }_{j}({\boldsymbol{k}}){|}^{2}={\left|\displaystyle \int {\rm{d}}{\boldsymbol{r}}\exp (-{\rm{i}}{\boldsymbol{k}}\cdot {\boldsymbol{r}}){\varPsi }_{j}({\boldsymbol{r}})\right|}^{2}.\end{eqnarray}$$ (5)

In our paper, we introduce a 2D harmonic trap V(r) = mω²/2(x² + y²) to simulate the realistic devices of ultracold atomic experiments. We set γ_x = γ and γ_y = ζγ_x = ζγ, where the coefficient ζ is the SOC strength ratio between γ_y and γ_x. The length (x, y) is normalized by l=ℏ/(mω). In this unit, the units of the SOC (γ) and interaction (c₀, c₂) are lω and ℏω l³/N, respectively. The density–density interaction is fixed at c₀ = 50.

We study the ground-state phases of spin-1 BECs in a 2D harmonic trap with isotropic and anisoropic SU(3) SOCs, respectively. The ground-state phases are obtained by using the time-splitting Fourier pseudospectral method with the imaginary time propagation (t → –i t).^[51–55] The wave function in Eq. (4) can be written as Ψj(r,t)=∑k=1nΨjk(r)exp(−itEk) (j = –1,0,1), where Ψ_jk(r) is the eigenfunction of the eigenvalue E_k. Each eigenfunction is attenuated by exp(– i tE_k). For the different eigenfunctions, the eigenvalues E_k are different. The ground-state has the lowest energy, the decay rate is the slowest. Therefore, after each iteration, the proportion of the ground-state becomes larger and larger. The returned wave function evolves to the ground-state wave function, and its limit is also the ground-state wave function. However, in the numerical calculation, the imaginary time propagation neither ensures the normalization nor conserves the spin-density in spin–orbit coupled spinor BECs. In order to solve this question, a backward-forward Euler Fourier-pseudospectral discretization is implemented to ensure the normalization of the total density and conservation of the spin-density. For the time discretization, we use the backward or forward Euler scheme for linear or nonlinear terms in the time derivatives. For the spatial discretization, we take fast Fourier transform in spatial derivatives. We adopt three conserved quantities to calculate the normalization constants in this paper. These quantities ensure the normalization of the total density and conservation of the spin-density after each iteration in imaginary time, the continuous normalized gradient flow has been discussed in Refs. [56,57]. In each iteration, the wave function in Eq. (4) is transformed as Ψ_j (r,t + Δt) = d_jΨ_j (r,t), where d_j are the normalization constants. So the constraint on the total number of atoms can be written as N=Σjdj2Nj, where Nj=Σ−∞∞|Ψj|2dr. We introduce the constraint on the spin-density as M=Σ−∞∞(d12|Ψ1|2−d−12|Ψ−1|2)dr. The method is discussed and used in the study of the spin–orbit coupled spin-1 and spin-2 BECs.^{[43,44,58–60]} In our numerical simulation, the initial wave function is the normalized Gaussian wave packet, i.e., Ψj(r)=(πωr2)1/4exp(−r2/2ωr2), where ω_r is the width along the r direction. The spatial and time steps employed are Δr = 0.05 and Δt = 0.0001, respectively.

In the SU(3) spin–orbit coupled spin-1 BECs, the competition between the spin–spin interaction and the isotropic SU(3) SOC plays an important role in determining the ground-state phases, a rich variety of lattice configurations are present, i.e., the PW phase, the ST phase, the KL phase, the SHL phase, and the HHL phase. The system undergoes a sequence of phase transitions from the RL phase to the HHL phase, and to the TL phase with anisotrpic SU(3) SOC.

In this section, we study the ground-state phases of the ferromagnetic and antiferromagnetic spin-1 BECs with isotropic SU(3) SOC, respectively. We first consider the ferromagnetic spin-1 BECs, i.e., c₂ < 0. The ground-state densities of spin-1 BECs with the different SOC strengths γ = 0, 1, 3, and 5 are shown in Figs. 1(a)–1(d), respectively. The spin–spin interaction c₂ = –50. The particles favor the m_F = 0 component without SOC γ = 0 in Fig. 1(a). With the SOC strength increasing, the particles begin to prefer the m_F = ± 1 components. As the SOC strength increases further, the densities of the three components are uniformly distributed, i.e., |Ψ₋₁|² = |Ψ₀|² = |Ψ₁|² = N/3, as shown in Fig. 1(d) with γ = 5. For the SU(3) spin–orbit coupled ferromagnetic spin-1 BECs, the ground state is a PW phase, and one of the three minima of the single-particle energy spectrum is occupied.

Figures 2–4 exhibit the ground-state phases of antiferromagnetic spin-1 BECs with isotropic SU(3) SOC, the phenomena become interesting. The ground-state phases with the weak spin–spin interaction (c₂ = 100) are shown in Fig. 2. The system prefers the state that particles locate in the m_F = ±1 components evenly, i.e., |Ψ₋₁|² = |Ψ₁|² = N/2, one can see in Fig. 2(a) without the SOC γ = 0. When considering the SOC, the translational symmetries of each component along the x-direction and y-direction are broken by the SOC, the densities of the three components are immiscible. The density of each component forms a kagome lattice structure^[47,48] in Figs. 2(b)–2(d), we call this ground-state phase as the KL phase. Due to the presence of the strong SU(3) SOC, the momentum distribution of the KL phase has three discrete minima. The three discrete minima are occupied with the equal weights, as shown in Fig. 5(a). With the SOC strength increasing, the period of the density modulation decreases. We next discuss the ground-state phases of the SU(3) spin–orbit coupled spin-1 BECs with the large spin–spin interaction c₂ = 2000, a rich variety of lattice configurations are displayed. The translational symmetry of each component is broken only along one direction with the weak SOC γ = 0.75, 1.10, 1.65, as shown in Figs. 3(b)–3(d). The ground-state phases consist of alternating spin domains between m_F = 0 component and m_F = ±1 components, which realizes the ST phase. The ST phase is consistent with the ground-state phase of the SU(2) spin–orbit coupled spin-1 BECs. If the SOC strength is beyond a critical value (γ ≃ 1.70), the translational symmetries along both the x-direction and the y-direction are broken in each component. The double vortices structure in the total density |Ψ|² is shown in Fig. 3(e4) with the SOC strength γ = 1.75. Figures 4(a)–4(e) show the ground-state phases with the strong SOC strengths γ = 2, 2.5, 3.5, 4, and 4.4, respectively. We can find that the densities of the three components are miscible. The density of the m_F = 0 component has the honeycomb lattice structure and the m_F = ±1 components have the rectangular lattice structures, as shown in Fig. 4(a). The total density realizes the honeycomb lattice with a vortex core structure, which is different from the phases shown in Fig. 3. As the SOC strength increases further, the vortex core structure is replaced by the honeycomb hexagonal lattice structures in the edge regime of total densities in Figs. 4(b)–4(d). The honeycomb hexagonal lattice structure can exist stably when the SOC strength γ = 4.4 in Fig. 4(e), we take it as the HHL phase. The ground-state phases in Figs. 3(e) and 4(a)–4(d) are the transition phases between the ST phase and the HHL phase, we call them as the SHL phases.

The momentum distributions of each component and the spin texture are shown in Fig. 5. The parameters of Figs. 5(a)–5(d) are the same as those of Figs. 2(d), 4(a), 4(b), and 4(e), respectively. The momentum distributions of each component of the KL phase (see Fig. 2(d)) show three discrete minima occupied with the equal weights in Fig. 5(a). When the SOC strength γ = 2, four discrete minima in the momentum distributions of each component are shown in the SHL phase of Fig. 5(b). The three minima of the single-particle energy spectrum are broken by the interaction and SOC, and the SHL phase with four discrete minima is a metastable state. As the SOC strength increases, one of the minima gradually weakens and eventually disappears. We find that the SHL phase with four discrete minima is confined to a very small SOC regime 1.95 ≤ γ ≤ 2.25 in our numerical calculation. The SHL phase with three discrete minima is shown in Fig. 5(c) with the SOC strength γ = 2.5. For the HHL phase (see Fig. 4(e)), three discrete minima of the equilateral triangle have unequal weights, as shown in Fig. 5(d). The spin textures of the KL phase, the SHL phase, and the HHL phase are shown in Figs. 5(a4), 5(b4), 5(c4), and 5(d4). The spin textures show a spontaneous magnetic ordering in the form of crystals of the meron pairs and antimeron pairs with the strong SU(3) SOC strength. Previous studies indicated that stable meron-pair lattice can be obtained in two-component BECs in a periodic potential^[47] or with rotation.^[61] Our results show that meron-pair lattice can also be stabilized in alternating spin domains by the SU(3) SOC of spin-1 BECs.

The ground-state phase diagram spanned by the SU(3) SOC strength γ and the spin–spin interaction strength c₂ is shown in Fig. 6. Five ground-state phases can be identified on this phase diagram, including the PW phase, the KL phase, the SHL phase, and the HHL phase. The ferromagnetic SU(3) spin–orbit coupled spin-1 BECs (c₂ < 0) only have the PW phase. The PW phase is also shown in the antiferromagnetic SU(3) spin–orbit coupled spin-1 BECs (c₂ > 0) with the weak SOC strength γ ≤ 0.7. When the spin–spin interaction strength is in the regime of 0 < c₂ ≤ 600, the system undergoes the phase transition from the ST phase to the KL phase as the SOC strength increases. With the spin–spin interaction strength c₂ increasing, the KL phase is replaced by the STL and HHL phases. The system undergoes a sequence of phase transitions from the ST phase to the STL phase, and to the HHL phase. From the phase diagram, we can find that both the spin–spin interaction and the SU(3) SOC play an important role on the ground-state phases of the SU(3) spin–orbit coupled spin-1 BECs.

In order to better investigate the effect of the SU(3) SOC on the ground-state phases of spin-1 BECs, we study the ground-state phases of antiferromagnetic spin-1 BECs with anisoropic SU(3) SOC in a harmonic trap. Figure 6 exhibits the ground-state phases of anisotrpic SU(3) spin–orbit coupled spin-1 BECs with the spin–spin interaction c₂ = 2000 and the SOC strength along the x direction γ_x = 5. The SOC ratios in Figs. 6(a)–6(c) are ζ = 0.9, 0.99, and 1.1, respectively. Three types of lattice phases are found, i.e., the RL phase, the HHL phase, and the TL phase. The translational symmetries of the total density along both the x-direction and y-direction are broken. The vortices in the neighboring chains are parallel, which form a RL phase, as shown in Fig. 6(a) with ζ = 0.90. With the SOC ratio increasing, the HHL phase is shown. The HHL phase only can be found in the regime of |ζ – 1| ≤ 0.02. which shows that the HHL structure is the unique solution of ζ = 1 in the numerical calculation. As the SOC ratio increases further, the translational symmetries of the total density along the x-direction and y-direction are also broken, the vortices of the neighboring chains are stagger. The TL phase is found in Fig. 6(c) with ζ = 1.1. The translational symmetry is broken by the SOC ratio, and the system undergoes a sequence of phase transitions from the RL phase to the HHL phase, and to the TL phase in spin-1 BECs with anisotrpic SU(3) SOC.

We have investigated the ground-state phases of 2D isotropic and anisotropic SU(3) spin–orbit coupled spin-1 BECs in a harmonic trap respectively. The competition between the SU(3) SOC and spin–spin interaction results in a rich variety of lattice configurations. Five ground-state phases, i.e., the PW phase, the ST phase, the KL phase, the SHL phase, and the HHL phase, are identified on the phase diagram with isotropic SU(3) SOC. The system undergoes a sequence of phase transitions from the RL phase to the HHL phase, and to the TL phase in spin-1 BECs with anisotrpic SU(3) SOC.