Visible stripe phases in spin–orbital-angular-momentum coupled Bose–Einstein condensates

Recently, stripe phases in spin–orbit coupled Bose–Einstein condensates (BECs) have attracted much attention since they are identified as supersolid phases. In this paper, we exploit experimentally reachable parameters and show that annular stripe phases with large stripe spacing and high stripe contrast can be achieved in spin–orbital-angular-momentum coupled (SOAMC) BECs. In addition to using Gross–Pitaevskii numerical simulations, we develop a variational ansatz that captures the essential interaction effects to first order, which are not present in the ansatz employed in previous literature. Our work should open the possibility towards directly observing stripe phases in SOAMC BECs in experiments.


Introduction
The realization of synthetic gauge fields and spin-orbit coupling for ultracold atoms has opened new opportunities for creating and investigating topological matters in a clean and easy-to-manipulate environment [1][2][3][4]. In the spin-orbit coupled Bose-Einstein condensates (BEC) realized in early works [5][6][7], the internal spin states are coupled to the center-of-mass linear momentum of the atoms via Raman laser dressing. There, the Raman beams transfer photon momentum to the atoms as the spin state changes. By using a similar method with Laguerre-Gaussian (LG) Raman beams which transfer orbital-angular-momentum (OAM) between atomic spin states, physicists recently demonstrated a coupling between internal spin states and the center-of-mass OAM [8][9][10]. In the following, we refer to the former as spin-linear-momentum coupling (SLMC) and the latter as spin-orbital-angular-momentum coupling (SOAMC) 4 .
The interplay between interactions and SLMC leads to interesting quantum phases [11][12][13][14][15]. For a pseudospin 1/2 system, by tuning the Raman coupling strength from large to small values, the energy-versus-momentum dispersion transforms from a single minimum to double minima (see figure 1(b)). For the latter case, whether the atoms occupy one of the minima or both minima is determined by a competition between inter-and intra-species interactions, where the two species refer to atoms which occupy the two respective minima. When the atoms occupy both minima associated with different quasimomentum, the interference results in density modulations in the position space, which is known as a stripe phase. This is a spatially mixed phase, i.e., miscible phase, of the two quasimomentum states. When only one of the minima is occupied, it is the separated phase, i.e., the plane-wave phase with single quasimomentum, and is immiscible for the two quasimomentum states. When the Raman coupling is below (above) the critical value, the ground state is the stripe phase (separated phase); see figure 1(b). The miscibility is indicated by the critical Raman coupling strength and the magnitude of stripe contrast, where both are determined by the spin-dependent interaction. The effective two-level scheme can be realized by using a large quadratic Zeeman energy ω q such that |m F = 1 is off-resonance. (b) Schematic SLMC phase diagram (c) SOAMC phase diagram of atoms in a harmonic trap. Δ = 20, r M = 5 μm, R TF = 12.5 μm, and μ/h = 926 Hz. For δ = 0, the ground state is the stripe (separated) phase with Ω M < (>)Ω c . For small nonzero detuning, the critical coupling is smaller than the Ω c (δ = 0) and decreases with increasing |δ|. Near Ω M = 0, the stripe phase exists for |δ|/2π 1 Hz. At δ = 0, the magnetization σ z = N ↑ − N ↓ is zero; the stripe contrast at δ = 0 is larger than that at any other nonzero δ. The grey bar at δ = 0 and Ω M > Ω c indicates the separated phase can have either σ z = 1 or −1.
The stripe phase in SLMC BECs is intriguing since it spontaneously breaks the translational symmetry (being a solid) and the U(1) gauge symmetry (being a superfluid) simultaneously, leading to a so-called supersolid [16]. Analogously, the ground state of SOAMC BECs also has an annular stripe phase and separated phases, which are theoretically studied in references [17][18][19][20]. The annular stripe phase of SOAMC BECs corresponds to occupying both energy minima with different quasiangular-momentum. The stripe spatial period is then ≈2πR/Δ , where R is a typical length scale smaller than the BEC size R BEC , and Δ is the transferred OAM between spin states in units of . Since R is the order of micrometers, the spatial period can be made larger than that in SLMC, which is λ/2 with λ being the optical wavelength of the Raman laser. The submicron stripe period of SLMC BECs is difficult to resolve even with the state-of-the-art quantum gas microscope [21].
Due to both the small spatial period and small contrast resulting from small miscibility, direct observations of stripe phases in position space remain elusive to date. Recent experimental works have demonstrated signatures of stripe phases in spin-linear-momentum-coupled BECs [22][23][24], where Bragg spectroscopy is employed to detect the stripe density modulation in references [22,23]. In the experiment with Raman-coupled internal spin states [22], the spatial phase coherence of both the stripe and separated phase is demonstrated interferometrically. In reference [23], atoms localized within each side of a double well serve as two pseudospin states. This circumvents the problem of detuning noises owing to the magnetic field noises for internal spin states, and enhances the miscibility. The observed stripe contrast is ∼8% limited by heatings from the Raman driving fields which create SLMC.
In this paper, we exploit the advantages of SOAMC systems and demonstrate the feasibility to directly observe annular stripe phases in situ with practical experimental parameters. We observe that interactions reduce the stripe density contrast. Here the interaction strength is ε int /E L < 1, where ε int is the mean field interaction energy and E L = 2 Δ 2 /2mR 2 is the characteristic energy scale of SOAMC systems. The effects of interactions discussed in previous papers [14,18,20,25] are based on the wave function ansatz that is not fully self-consistent in the presence of interaction. Even within first order in interaction strength, we find that the results of references [14,18,20,25] are subject to significant corrections. We use an improved ansatz and obtain results that are correct to first order in interaction. While in SLMC systems, the analogous interaction strength is ε int /4E r and is typically small, where the photon recoil energy E r is larger than E L . We investigate how the stripe density contrast depends on experimentally accessible parameters: the transferred angular momentum Δ , the size of the OAM-carrying LG Raman beam, the BEC cloud size, and the mean field energy. By optimizing these parameters, we achieve a stripe period of ∼2 μm and a 30% contrast of density modulations. This is detectable using high-resolution imaging with about 1 μm resolutions [21]. Further, the contrast can be made larger than 30% by increasing the BEC cloud size. Finally, we point out that by using synthetic clock states [26], the stripe phase of the thermodynamic ground state can be stable against external magnetic field noises despite the narrow detuning window within which the stripe phase exits.

Formalism
We consider pseudospin 1/2 atoms tightly confined along z in a quasi-two-dimensional (quasi-2D) geometry, where ω z > μ with ω z being the trap frequency along z and μ the chemical potential. Two Raman beams couple the two spin states with a transfer of OAM Δ in unit of , and the frequency difference between the two beams is Δω L . In the rotating frame at frequency Δω L with rotating wave approximation, the single-particle Hamiltonian iŝ where L z = −i ∂ φ is the angular momentum operator, V(r) is the spin-independent trapping potential, δ = Δω L − ω 0 is the Raman detuning, and ω 0 = E ↓ − E ↑ is the energy splitting between |↓ and |↑ . The Raman beams are two Laguerre-Gaussian beams of order Δ /2 and −Δ /2, and the coupling strength is where the peak coupling Ω M is at r = r M , and the waist of each beam is w = 2r M / √ Δ . In addition toĤ 0 , we have the mean field energy where |ψ ↑ | 2 , |ψ ↓ | 2 are the 2D density of |↑ , |↓ . The wave functions are normalized as drr dφn(r, φ) = N where n = |ψ ↑ | 2 + |ψ ↓ | 2 and N is the atom number. The 2D interaction strengths are g = g ↑↑ = g ↓↓ and g ↑↓ . We define g 1 = (g + g ↑↓ )/2, g 2 = (g − g ↑↓ )/2, g 2 being the spin-dependent interaction strength. We use real experimental parameters by taking the pseudospin states as |↑ = |F = 1, m F = 0 and |↓ = |1, −1 of 87 Rb atoms, for which g = (g 00 + g −1,−1 )/2, g ↑↓ = g 0,−1 with The scattering lengths are a 00 = 100.86a B and a −1,−1 = 100.40a B , where a B is the Bohr radius [27]. This gives g > g ↑↓ and positive g 2 /g 1 = 0.00114. As compared to the realistic case with g ↑↑ = g ↓↓ , here our simplification of using g = g ↑↑ = g ↓↓ is based on the results of uniform SLMC systems in the absence of trapping potentials in reference [14], which is just a shift in detuning for the ground state. We show that this is a good approximation for the trapped atoms with inhomogeneous n(r) under SOAMC.
For δ = 0 in the non-interacting limit, the ground state may be expressed as where > 0, |C + | 2 + |C − | 2 = 1, andn(r) is the density after azimuthal average with dr2πrn = N. Since the Raman coupling has a phase winding number Δ , i.e., an OAM of light, the Raman beams couple |↑, ↑ to |↓, ↓ where the OAM difference between the spin states is ↑ − ↓ = Δ . By introducing the quasiangular momentum , ↑ and ↓ are rewritten as ↑ = + Δ /2 and ↓ = − Δ /2. Then, equation (4) is referred as 'two-quasiangular-momentum ansatz', which have two running wave components along φ with quasiangular momentum ± . For sufficiently small Raman coupling, the ground state has = Δ /2 [17,18], and we focus on this regime throughout this paper. With = Δ /2, there is a critical coupling Ω c below which the ground state has |C + C − | > 0. For Ω M < Ω c , equation (4) shows that the spin component |↑ |↓ has an OAM superposition of ↑ = Δ and 0 ( ↓ = 0 and −Δ ), leading to a density modulation along φ, which is then called stripe phase (see figure 2(a)). Here, where ϕ is the relative phase between C + and C − . With Ω M > Ω c , the ground state is the separated phase with |C + C − | = 0, i.e., |C + | = 1, |C − | = 0 or |C − | = 1, |C + | = 0 (see figure 2(b)), which are equivalent for δ = 0. For the stripe phase with |C + C − | > 0, the ground state has |C + | 2 = |C − | 2 = 1/2 and |C + | = |C − | where |C + | 2 |C − | 2 is maximized. Note that at |C + | = |C − | the wave function equation (4) is an eigenstate of the time-reversal operator T =σ x K with K being the complex-conjugate operator, which is possible because the Hamiltonian commutes with T. At radial position r, the contrast of the azimuthal density modulation is where n max , n min and n avg are the maximum, minimum and average of the density along φ, respectively. Then, the contrast of both |↑ and |↓ from equation (5) is and the spatial period of the density stripe is 2πr/Δ . Now we consider the general form of the spinor wave function for all interaction strengths, which is where p = 0, ±1, . . . . are integers. The multiple OAM components differing by Δ are due to the nonlinear interaction term, and equation (8) is similar to that in reference [28] for spin-linear-momentum coupling.
In the two-quasiangular-momentum ansatz equation (4) with = Δ /2, it has only a 0 , a Δ and b 0 , b −Δ . For small θ, and small interactions in the first order perturbation regime, equation (8) can be simplified as (see appendix A) having only a 0 , a ±Δ and b 0 , b ±Δ , This is referred as 'four-quasiangular-momentum ansatz', which has quasiangular momentum of ±Δ /2, ±3Δ /2. Since the ground state energy is independent of the relative phase between C + and C − , we can take C + < 0, C − > 0 as real values without loss of generality. Similar to the case with two-quasiangular-momentum ansatz, we consider states that are eigenstates of the time-reversal operator T, which will be confirmed by numerical simulations. The variational parameters then satisfy where for |↑ , and −C + cos θ(−C − sin θ + A + )e iΔ φ + c.c. for |↓ . Thus, for the ground state with minimized density modulation, the relative phase between C + sin θ and A * − and between −C − sin θ and A + is π. This gives real and positive A + , A − . We then have Similarly, if we choose C + , C − > 0, the condition becomes as real numbers, the densities of |↑ and |↓ are The normalization is C 2 and for small A + , Thus, the density contrast of the cos(Δ φ) term in equation (12) is by using equation (6) with (n max − n min )/2 equal to the amplitude of the cos(Δ φ) term and n avg ≈ 1/2 for each spin component. With equation (11), the contrast is Since the density modulation of the cos(2Δ φ) term is much smaller than η ↑,↓ , we use equation (15) as the contrast in our simulations.

Simulations methods
We perform both the Gross-Pitaevskii (GP) simulations and the variational calculations to find the ground state in the SOAMC system. The GP simulation gives the ground state with the full Hamiltonian including bothĤ 0 and the interaction energy. Additionally, we perform the variational calculations with a simplified picture: we neglect the radial kinetic energy associated with ∂ r inĤ 0 , thus the rest of all the energy terms are functions of radial position r. We find the results of GP and variational methods have good agreements.

Gross-Pitaevskii ground state
We use the GP simulations to find the ground state by numerically solving the Gross-Pitaevskii equation (GPE). We perform imaginary time propagations, where the initial state of the stripe phase for the imaginary time propagation is and for the separated phase it is The initial state of the stripe phase has a superposition of OAM differing by Δ in either spin up and down, such that all the OAM components differing by Δ [see equation (8)] can be reached in the final ground state. The initial state of the separated phase corresponds to C + = 1, C − = 0. After the numerical computation, we compare the energy differences between the two phases and determine the ground state from the lower energy state.

Variational method
We adopt a variational method to minimize the energy E var for δ = 0 and obtain the variational ground state. E var includes the single particle Hamiltonian in equation (1) but excluding the radial kinetic energy from ∂ r , and the mean field interaction equation (3). This gives E var = drr2πn(r)ε(r) where ε(r) is the energy per atom after the azimuthal average. We discuss calculations based on the two-quasiangular-momentum ansatz, equation (4), and four-quasiangular-momentum ansatz, equation (9), respectively. The variational ground state from the simple equation (4) agrees with our GP simulation in the non-interacting limit. This variational form, equation (4), is used in earlier papers [14,18,20,25]; [14] for SLMC and [18,20,25] for SOAMC BECs. In our simulations, we find the variational ground state from equation (4) is inconsistent with the GP result in the non-negligible interaction regime, where the additional OAM components must be taken into account as the ansatz equation (9).

Two-quasiangular-momentum ansatz
With the ansatz of equation (4), the variational energy per atom ε var0 is given by 0 is the single-particle energy arising from the Raman coupling and the centrifugal potential L 2 z /2mr 2 , where the latter is characterized by E L = 2 (Δ ) 2 /2mr 2 at position r. Here we exclude the trap energy V(r) inĤ 0 since V(r) does not depend on any variational parameters and is simply an offset. ε var0 int is the mean field interaction energy with β = |C + | 2 |C − | 2 satisfying 0 β 1/4. Given the local energy ε var0 (r) at a radial position r with the averaged densityn(r), we take θ(r) and β as variational parameters. β = 0 for the separated phase with |C + | = 1, |C − | = 0 or |C + | = 0, |C − | = 1, and β = 1/4 for the stripe phase with Within 0 β 1/4, the energy difference between the stripe and the separated phase is lowest at either β = 0 or β = 1/4 (see appendix A). At a given β, minimizing ε var0 (r) with respect to θ determines θ as a solution of the following equation: In the non-interacting case, the solution of equation (19) is tan 2θ 0 = Ω/E L , or equivalently which can be approximated for small Ω/E L as With interactions, the solution θ var0 of the stripe phase with β = 1/4 is smaller than θ 0 , given by from which the contrast is given by as derived in equation (7). We usen(r) obtained from the GP simulation, which is the same for the stripe phase and separated phase and is well approximated by the Thomas-Fermi (TF) profile except for small r M . By expanding to first order inng 1 /E L and Ω/E L , For the separated phase with β = 0, θ sep is well approximated with θ 0 owing tong 2 E L .

Four-quasiangular-momentum ansatz
With the ansatz of equation (9), the single-particle part of the variational energy is given by where C + , C − , A + , A − are real. The interaction energy ε var int is also a function of θ, C + , C − , A + , A − ; by using equation (11) for the stripe phase, we plug in (24) Then we minimize ε var = ε var 0 + ε var int with respect to (θ, A ± ), respectively, giving the numerical solutions for the stripe phase, θ var and A var ± > 0, where the sign of A var ± agrees with equation (11). The contrast of the for both |↑ , |↓ following equation (15). By expanding A var ± to first order inng 1 /E L and Ω/E L , we obtain A var ± ≈ng 1 θ/ √ 2E L , which also agrees with the result using perturbation (see appendix A). After plugging it into equation (25) and expand θ var , we have Comparing to η var0 in equation (22), we find the coefficient ofng 1 That is, including the additional OAM ±3Δ /2 in equation (9) is necessary for correct results to first order inng 1 /E L . For the separated phase with C + = 1, C − = 0 and A + = A − = 0, it is identical to that using the two-quasiangular-momentum ansatz. We comment on earlier theoretical papers on SOAMC systems [18,20,25]. We examine the dimensionless interaction strengthng 1 /E L in these papers. In reference [20],ng 1 /E L is 0.01, and single-particle eigenstates are taken as the basis of the variational method, i.e., θ = θ 0 . References [18,25] use variational methods with the wave function ansatz equation (4), where reference [18] has ring traps withng 1 /E L > 100, and the interactionng 1 is not specified in reference [25].

Results and discussions
We consider practical experimental parameters to maximize the density contrast of the stripe phase. We first discuss BECs in harmonic traps in the Thomas-Fermi regime along the radial direction with the Thomas-Fermi radius R TF . We study how the GP stripe phase contrast depends on (Δ , r M , R TF , μ); μ is the chemical potential and the peak mean field energy in the harmonic trap.
For comparisons, we also consider atoms in a ring trap with radius r 0 . Here r M = r 0 is the only length scale, unlike the harmonically trapped systems where there are two relevant length scales, (r M , R TF ).

Harmonic traps
We first obtain the GP ground state phase diagram as shown in figure 1(c). We then focus on the GP stripe phase at δ = 0, setting Ω M = Ω c . We run simulations for Δ between 2 and 30, all with r M = 17 μm, R TF = 46 μm, and μ = h × 21 Hz. Δ = 30 corresponds to the LG beam with phase winding number of ±15, which can be achieved experimentally (in reference [29], LG beams with phase winding number of 45 are realized). From the GP wave function ψ(r, φ), we evaluate the density contrast η GP (r) from the normalized Fourier componentsψ of ψ(r, φ) following equation (15). For |↓ ,ψ ↓ ,↓ are given by [see equation (9) from which the contrast η GP (15) can be described as 4ψ 0,↓ψ−20,↓ (the first term) and 4ψ 0,↓ψ20,↓ (the second term). We then compare η GP (r) to the variational solutions of the contrast, which are sin 2θ 0 (r) for the non-interacting case, η var0 (r) for using the two-quasiangular-momentum ansatz [equation (21) from the ansatz equation (4)] and η var (r) for using the four-quasiangular-momentum ansatz [equation (25) from the ansatz equation (9)]. In figure 3(a), we plot sin 2θ 0 (r), η var0 (r), η var (r) and η GP (r) for the example value Δ = 20; their maxima are at r r M . In figure 3(b), We plot the peak values of sin 2θ 0 (r), η var0 (r), η var (r) and η GP (r) versus Δ , which are denoted as sin 2θ 0 , η var0 , η var and η GP , respectively. We observe that the single-particle contrast sin 2θ 0 is significantly larger than η GP for small Δ , while sin 2θ 0 and η GP are close for Δ 20. As the dimensionless interaction ng 1 /E L increases with decreasing Δ , the contrast η GP decreases. When the interaction is taken into . From variational solutions: sin 2θ 0 (r) for non-interacting atoms (grey curve), η var0 (r) for using two-quasiangular-momentum ansatz (green), and η var (r) for using four-quasiangular-momentum ansatz (red). Blue curve indicates η GP (r) for the GP. (b) Peak values of the contrast vs Δ . Grey, green, red, and blue symbols denote sin 2θ 0 , η var0 , η var , η GP , respectively. (c) Annular Fourier power spectrum after integration along r for |↓ (dark grey) and |↑ (light grey). consideration using the ansatz equation (4), the resulting η var0 overestimates η GP , indicating that equation (4) is insufficient. We can understand this from the annular Fourier transform of the GP wave function ψ ↓ for Δ = 20. Figure 3(c) shows the power spectrum of the normalized Fourier components ψ ↓ ,↓ , where there are ↓ = 0, ±20 components, and the ↓ = −40 is negligible since its power is about 10 −3 of that of ↓ = 20. The appearance of the ↓ = 20 component signifies that the more general equation (9) should be used in the variation method with the A + term accounting for ↓ = 20, while A + , A − are absent in the simple equation (4). The spectrum for |↑ is also displayed in figure 3(c). The GP results have C ± ≈ ∓1/ √ 2 and A + = A − = A GP ± , which confirms the time-reversal symmetry condition, equation (11). [From the GP results, the signs of C ± ≈ ∓1/ √ 2 are applied in equation (11) and in the variational method using equation (9).] We also show the peak values, η var of η var (r), in figure 3(b), where η var0 > η var η GP and η var fits well with η GP .
From the above studies, we find the maximum of the stripe contrast is at r = r peak ≈ r M , where the spatial period is ≈2πr M /Δ . (r peak is only slightly larger than r M for all the contrasts, e.g., r peak = 17.6 μm for η GP in figure 3(a).) The peak value of the contrast increases with Δ , and thus a larger contrast corresponds to a small stripe period ∝Δ −1 . To observe the stripe phase in experiments, we note that with state-of-the-art imaging techniques in ultracold atoms, e.g. those using quantum gas microscopes, one can resolve as small as 0.5 μm with λ = 0.78 μm for 87 Rb [21]. This sets the lower bound on 2πr M /Δ in our simulations. Since the peak contrast is the signal we optimize, r peak ≈ r M and thus E L (r peak ) ≈ E L (r M ) are the relevant length and energy scale for SOAMC, respectively.
Next, we fix Δ = 20 and μ = h × 93 Hz, and vary (r M , R TF ). We study r M < R TF where there is sufficient atomic density at r = r M for various combinations of (r M , R TF ). Here we set the smallest r M = 5 μm for Δ = 20, where the spatial period of the stripe is ≈1.6 μm and is larger than the diffraction limit of the imaging, 0.5 μm. With Ω M = Ω c , the GP results have the peak contrast η GP increasing with decreasing r M and increasing R TF , as shown in figure 4(b). We then compare the GP to the variational calculations with four-quasiangular-momentum ansatz. In figure 4(a), we plot the density contrast contributed from θ, A ± and the sum, respectively. These are −sin 2θ GP , 2 √ 2 cos θ GP A GP ± , and η GP versus r for GP, and −sin 2θ var , 2 √ 2 cos θ var A var ± , and η var versus r for the variational calculations. We find the contrast of GP obtained from equation (6) versus r agrees well with η GP (r), showing that η GP from the cos(Δ φ) term dominates the contrast in equation (12), where the second harmonics is negligible. We display the peak values η var versus r M for all R TF in figure 4(b) and compare them to the peak values η GP , where η var overestimates η GP by 3-8%. For all (r M , R TF ), the Fourier spectrum of GP has ↑,↓ = 0, ±20 components and ↑ = 40, ↓ = −40 are negligible, being consistent with equation (9). We then compare the peak values (θ GP , A GP ± ) to (θ var , A var ± ) versus r M for R TF = 50 μm in figure 4(c). The peaks of θ and A ± are at r ≈ r M , as well as that of the contrast η. θ var and A var ± slightly overestimate θ GP and A GP ± , respectively: both θ GP /θ var and A GP ± /A var ± are between 0.90-0.93. This is attributed to the radial kinetic energy that is neglected in the variational calculations. The GP ground state has smaller (θ GP , A GP ± ) than (θ var , A var ± ) where the smaller radial spin gradient corresponds to a smaller radial kinetic energy, and thus the lowest overall energy.
After studying the dependence of the peak contrast η GP on (Δ , r M , R TF ), we vary μ and thus the interaction strengthng 1 /E L at r = r peak ≈ r M , whereng 1 /E L = μ 1 − (r peak /R TF ) 2 /E L (r peak ). We fix Δ = 20, R TF = 50 μm for r M = 5 and 15 μm, respectively. We study μ = h × 21, 46, 93 Hz for both r M , and additionally μ = h × 600 Hz for r M = 5 μm. Figure 5 shows η GP weakly depends onng 1 /E L .  Besides the density contrast, we compare the critical coupling Ω c of the GP results to those given by the variational methods. Using the two-quasiangular-momentum ansatz, equation (4), the critical coupling Ω var0 c is given by dr2πrn(r)Δε var0 = 0, Δε var0 is the energy difference between the stripe phase and the separated phase, θ var0 is the variational solution of θ for the stripe phase, and ε sep is the energy of the separated phase with θ = θ sep and β = 0. When the integral is < (>) 0 at Ω M < (>)Ω c , the ground state is the stripe (separated) phase. Similarly, by using the four-quasiangular-momentum ansatz, equation (9), the critical coupling Ω var c is given by dr2πrn(r)Δε var = 0, In figure 6, we plot Ω c of the GP and from the solutions of equations (28) and (29) versus (r M , R TF ), where Ω c from GP have good agreements with Ω var c . We can understand that Ω c increases with increasing R TF and decreasing r M from a geometric argument. Such dependence on (r M , R TF ) is crucial since a larger Ω c /E L leads to larger stripe contrast η var , see equation (26). We find numerically that the energy difference between the stripe phase and separated phase can be approximated as for small θ, where H is a dimensionless function ofng 1 /E L and H → 1/4 asng 1 /E L → 0. The peak value of f(r) is where the number in the parentheses is a 'geometric' area ratio of the box to that of the distribution f(r); the ratio scales as R 2 box r −2 M Δ 1/2 . By equating equations (31) and (32), assuming a fixed interaction strength ng 1 /E L and using E L ∝ Δ 2 r −2 M , we obtain Ω M = Ω c ∝ R box r −3 M which increases with increasing R box and decreasing r M ; see figure 6.
We then derive the stripe contrast η var versus (Δ , r M , R box ). First, Ω c /E L is given by equating equations (31) and (32), The contrast η var can be expressed as where Ω c /E L is approximately the non-interacting contrast sin 2θ 0 in section 3.2.1, and η var /( Ω c /E L ) < 1 is the contrast-reduction ratio due to interactions; see equation (26) (with expansion to first order in interaction). η var /( Ω c /E L ) decreases with increasing interactionng 1 /E L , which is proportional to Δ −2 r 2 M withng 1 approximately fixed. Therefore, η var /( Ω c /E L ) increases with increasing Δ and decreasing r M . In equation (33), H decreases withng 1 /E L , and thus H −1/2 decreases with increasing Δ and decreasing r M . The overall Ω c /E L decreases with increasing Δ and increases with decreasing r M . When combining equations (33) and (34), η var increases with increasing Δ given fixed (r M , R box ,ng 1 ), and increases with decreasing r M and increasing R box . Such dependence on Δ and r M is shown in figures 3(b) and 4(b), respectively, where the system has harmonic traps with R TF replacing R box .

Ring traps
We show the variational results for the ring trap versus the dimensionless interaction strength g 1 =ng 1 /E L (r 0 ). These are the solutions of the critical coupling Ω c and the contrast at Ω c . We solve the dimensionless critical coupling Ω c = Ω c /E L by using Δε var0 (r 0 ) = 0 and Δε var (r 0 ) = 0, respectively, and plot Ω c vs g 1 in figure 7(a). Ω var c exceeds Ω var0 c for g 1 > 0, and Ω var0 c ≈ Ω var c → 2g 2 /g 1 E L ≈ 0.05 E L as g 1 → 0. we can understand Ω var c > Ω var0 c as the following: at a given Ω, the stripe phase energy is ε var < ε var0 since the smaller contrast η var [see equations (22) and (26)] corresponds to smaller interaction energy. Therefore, the critical coupling determined by Δε = 0, where Δε increases with Ω, is shifted to a larger value for Ω var c . In the g 1 → 0 limit, Δε var0 (r 0 ) ≈ Δε var (r 0 ), where the stripe and separated phases have approximately the same variational solution, θ var0 ≈ θ sep ≈ θ 0 , and thus Δε var0 = (1/4)n g 1 sin 2 2θ 0 − 2g 2 cos 2 2θ 0 = 0 from equation (18) leads to Ω c = 2g 2 /g 1 . The stripe contrasts η var0 and η var at respective Ω c are displayed in figure 7(b); they are ≈5% as g 1 → 0, both decreasing with increasing g 1 and η var < η var0 . To compare with SLMC systems, Ω var0 c /E L for g 1 < 1 agrees well with that of SLMC (see figure 7(a)), which is Ω SLMC c /4E r = 2g 2 /g 1 as g 1 → 0 [14], i.e., Ω SLMC c = 0.2E r ; 4E r is equivalent to E L for SOAMC. The spin-dependent interaction strength g 2 /g 1 determines the critical coupling Ω c /E L and stripe contrast η as g 1 → 0; larger g 2 /g 1 gives larger Ω c /E L and η, i.e., larger miscibility. The contrast of SLMC also agrees well with η var0 , see figure 7(b). We find that Ω c and η of SOAMC with ε var0 and of SLMC, where both are based on equation (4), are incorrect to first order in g 1 . By including higher order OAM in equation (9), the result is correct to first order, as indicated by equations (22) and (26). The stripe contrast of SOAMC is at most ≈5% in the g 1 → 0 limit, and it is independent of either the ring radius r 0 or Δ . On the other hand, a stripe contrast 30% for harmonic traps is achieved with a relatively large R TF = 50 μm and a relatively small r M = 5 μm, and this can be understood from the geometric analysis as shown in equations (31)-(34).
We then perform GP simulations for a ring trap with r 0 = r M = 10 μm and Δ = 20. The atoms are in an annular box potential within 8 μm < r < 12 μm, andng 1 /E L = 0.63 at r = r 0 . The GP result has good agreement with the variational calculation, where θ GP /θ var = 0.95 and A GP ± /A var ± = 0.91.

Conclusions
In summary, we optimize the density contrast of the ground state stripe phase of 87 Rb SOAMC BECs by tuning experimental parameters. A contrast of nearly 30% is achieved for atoms in harmonic traps, and a larger contrast of about 50% is expected by using a twice larger BEC cloud size based on variational calculations. Such high contrasts are achieved owing to the geometry with two length scales in harmonic traps, the Raman Laguerre-Gaussian beam size and the BEC cloud size. While for ring traps, these two scales are the same, leading to maximal contrast about 5%, which is dictated by the spin-dependent interaction strength and is the same as that of the SLMC systems. For both atoms in harmonic traps and ring traps, we perform GP simulations and variational calculations based on the two-quasiangular-momentum ansatz and four-quasiangular-momentum ansatz. We find the results from the simple two-quasiangular-momentum ansatz, which is used in previous papers [14,18,20,25], are consistent with the GP results only in the non-interacting limit. With small interactions, high order OAM components must be included as the four-quasiangular-momentum ansatz; this then leads to correct results to first order in interaction and good agreements with the GP simulations.
We point out that one can improve the stability by using the synthetic clock states instead of bare spin states. The clock states are immune to detuning variations arising from the bias field variations. A 0.1-1 Hz stability of the clock transition frequency is achieved as shown in reference [26]. Thus, the ground state stripe phase within a narrow detuning window of about 1 Hz may be observed for 87 Rb atoms with mean field energy about 1 kHz. The spin-dependent interaction strength in the clock state basis is close to the g 2 /g 1 for bare spin states (see appendix A), leading to similar magnitude of stripe contrast to our simulations using bares spin states. We envision our work to pave the way towards a direct observation of high-contrast stripe phases in spin-orbital-angular-momentum coupled Bose-Einstein condensates, achieving a long-standing goal in quantum gases.

A.1. Spinor wave function ansatz
We show that the spinor wave function ansatz equation (9) is valid for small θ (given small Raman coupling Ω/E L ) and small interactionng 1 /E L . The GPE for |↑ , |↓ with δ = 0 is Here we set V(r) = 0 to simplify the discussion. For the spatially-mixed stripe phase ground state, μ ↑ = μ ↓ . Next we show the nonlinear interaction leads to multiple OAM components in the spinor wave function ansatz. With = Δ /2 = 10, we plug , neglect radial gradients and keep the expansion terms up to θ 2 . The nonlinear interaction terms for ψ ↑ are Besides ↑ = 0, 20, additional OAM terms with ↑ = −20, 40 appear due to the nonlinear interaction, which are of order of θ, θ 2 , respectively, and are not included in equation (4). By keeping up to order θ 2 , the spinor wave function has additional variational parameters A + , A − , B + , B − , given by Next we derive A − for small interactionsng 1 /E L using first order perturbation. We plug into equation (35a) for ψ ↑ , and focus on the coefficient of the e −i20φ term, which is Besides reading out from equation (36), the coefficient of the e −i20φ term in the nonlinear interaction g|ψ ↑ | 2 ψ ↑ + g ↑↓ |ψ ↓ | 2 ψ ↑ can be readily found from the Fourier components of |ψ ↑ | 2 , |ψ ↓ | 2 in equation (5), both of which have OAM = 0, ±20. Then, using μ ↑ ≈ μ ↑ (Ω = 0, g = g ↑↓ = 0) = 0 for small Ω/E L and ng 1 /E L , along with C ± ≈ ∓1/ √ 2, it gives From the GP stripe phase wave function, in figure 8(a) we plot the ratio of the peak values A GP ± /θ GP vs n(r peak )g 1 /E L (r peak ), along with the rationg 1 / √ 2E L given by equation (40), which agrees with A GP ± /θ GP at smalln(r peak )g 1 /E L (r peak ). The dimensionless interactionng 1 /E L evaluated at r peak ≈ r M vs r M for R TF = 12.5, 25, 50 μm is shown in figure 8(b), all with μ = h × 93 Hz. For a fixed μ,ng 1 = μ[1 − (r/R TF ) 2 ] weakly depends on r for r/R TF 0.5.n(r peak )g 1 /E L (r peak ) increases with increasing r M , which is dominated by E L ∝ r −2 M . Similarly, we derive B + by plugging into equation (35a). The coefficient of the e i40φ term is With θ ≈ Ω/2E L (1 −ng 1 /E L ) from θ var , and A + = A − for the ground state, it leads to In our GP data, we have small θ GP , 0.05 < θ GP < 0.16, and small interactions,n(r peak )g 1 /E L (r peak ) < 2 except for the data with the two smallest Δ = 2, 4 in figure 3. The peak values of B GP ± at r peak r M are small, 0.02 < B GP ± /A GP ± < 0.04 forn(r peak )g 1 /E L (r peak ) < 2. Thus it is valid to neglect B ± by using the wave function ansatz equation (9).

A.2. Methods for GP ground state simulations
We run the GP simulations with both the open-source GPELab toolbox [30] and Crank-Nicolson method. The grid size is between 0.093-0.55 μm depending on the spatial resolution we need.
To do analysis of the GP wave function in the cylindrical coordinate, we first make interpolations of the raw data in the Cartesian coordinate. The annular Fourier transform is performed as where q = ↑ , ↓ is the OAM and m =↑, ↓ is the spin label. For the stripe phase withn ↑ (r) =n ↓ (r) =n(r)/2, q |ψ q,m (r)| 2 =n(r)/2.

A.3. Variational calculations
We consider the SOAMC ground state as either the stripe phase with |C + C − | > 0 or the separated phase with |C + C − | = 0, i.e., C + = 1, C − = 0 or C − = 1, C + = 0. The former corresponds to a density stripe and the latter to no density stripe. In the variational calculation using two-quasiangular-momentum ansatz where A ± is absent in the wave function, 0 β = |C + | 2 |C − | 2 1/4 and β = 1/4 corresponds to |C + | 2 = |C − | 2 = |C + C − | = 1/2. We compare the energy of β = 1/4 and of β = 0, and take the lower one as the ground state. This is valid because the lowest energy is at either β = 1/4 or β = 0, i.e., no energy minimum within 0 β 1/4. We find this condition holds by numerically checking the second order derivative of ε var0 , which is negative for 0 β 1/4 andng 1 /E L < 2. As for the calculation using four-quasiangular-momentum ansatz for the stripe phase, we compare the energy ε var of the stripe phase, , and of the separated phase with C + = 1, C − = 0, A ± = 0. It is valid to take for the stripe phase since the GP results confirm that this condition holds, which has the time reversal symmetry, equation (11).

A.4. Dependence of stripe contrast on system parameters
We study the peak stripe contrast versus (Δ , r M , R TF ) and find the contrast increases with increasing Δ , decreasing r M and increasing R TF . We show the contrast η var of the variational calculations has such dependence at the critical Raman coupling Ω c . We employ a simplified model with a cylindrical box trap of radius R box , and obtain equations (33) and (34), showing the dependence of Ω c /E L and η var (Ω = Ω c ) on (Δ , r M , R box ,ng 1 /E L ); here the dimensionless interactionng 1 /E L is evaluated at r peak ≈ r M , whereng 1 /E L is proportional to Δ −2 r 2 M withng 1 approximately fixed (also see appendix A.1). H in equation (33) decreases withng 1 /E L as displayed in figure 9(a), and the resulting Ω c /E L decreases with increasing Δ and increases with decreasing r M . The contrast-reduction ratio η var /( Ω c /E L ) is unity for non-interacting atoms, and decreases with increasing interactionng 1 /E L ; see figure 9(b) with simulations of harmonic traps in figures 3 and 4. Therefore, η var /( Ω c /E L ) increases with increasing Δ and decreasing r M . By combining equations (33) and (34), η var increases with increasing Δ , decreasing r M and increasing R box . Note that Ω c /E L decreases while η var /( Ω c /E L ) increases with Δ , and Ω c /E L has a weaker dependence, leading to η var increasing with increasing Δ . In figure 9(b), the two largestng 1 /E L are about 4 and 20, corresponding to Δ = 4, 2, respectively. Although we expect the variational result to be valid at smallng 1 /E L , we find the contrast η GP and η var agree well for Δ = 2, 4; see figure 3(b).

A.6. Comparison to SLMC systems
We list results of SLMC BECs from reference [14]. Two counter-propagating Raman beams along x transfer linear momentum Δk x = 2k r between spin |↑ and |↓ , producing SLMC. The linear momentum transfer 2k r is analogous to the OAM transfer Δ in SOAMC, and thus 4E r is equivalent to E L in SOAMC. A spinor wave function ansatz analogous to our two-quasiangular-momentum ansatz, equation (4), is employed. For a uniform system with no trapping potentials, the critical coupling is Ω SLMC c ≈ 2g 2 g 1 4E r 1 +n g 1 4E r (47) after expanding to first order inng 1 /4E r for small interactionng 1 /4E r . The stripe contrast is and after expanding to first order inng 1 /4E r , which is the same as η var0 in equation (22) based on ansatz equation (4).

A.7. Scheme of using synthetic clock states
We propose to use synthetic clock states in the SOAMC system of 87 Rb atoms. Here the discussions are based on reference [26]. These clock states are |x , |y , |z , each of which is a radio-frequency-dressed state, and thus a superposition of bare spin states |m F = 0, ±1 . The lowest, middle, and highest-energy dressed state corresponds to |z , |x , |y , respectively. By choosing proper rf parameters, the xz transition frequency can be made fourth-order sensitive to rf detuning, and thus to the bias field. We consider a two-level system of Raman-coupled |x and |z . The mean field energy can be expressed in the basis of |x and |z , where the peak contrast is η GP = 0.217. This is very close to that with g ↑↑ = g ↓↓ = g, where Ω c /2π = 730.0 Hz and η GP = 0.210.