Optical lattice with spin-dependent sub-wavelength barriers

We analyze a tripod atom light coupling scheme characterized by two dark states playing the role of quasi-spin states. It is demonstrated that by properly configuring the coupling laser fields, one can create a lattice with spin-dependent sub-wavelength barriers. This allows to flexibly alter the atomic motion ranging from atomic dynamics in the effective brick-wall type lattice to free motion of atoms in one dark state and a tight binding lattice with a twice smaller periodicity for atoms in the other dark state. Between the two regimes, the spectrum undergoes significant changes controlled by the laser fields. The tripod lattice can be produced using current experimental techniques. The use of the tripod scheme to create a lattice of degenerate dark states opens new possibilities for spin ordering and symmetry breaking.


Introduction
Traditionally, optical lattices are created by interfering two or more light beams, so that atoms are trapped at minima or maxima of the emerging interference pattern depending on the sign of the atomic polarizability [1]. Optical lattices are highly tunable and play an essential role in manipulation of ultracold atoms [2][3][4]. The characteristic distances over which optical lattice potentials change are limited by diffraction and thus cannot be smaller than half of the optical wavelength λ. Yet the diffraction limit does not necessarily apply to optical lattices [5][6][7][8][9] or other sub-wavelength structures [10][11][12][13][14][15][16][17] relying on coherent coupling between atomic internal states.
It was demonstrated theoretically [5,6,8,9,18] and experimentally [19,20] that a periodic array of sub-wavelength barriers can be formed for atoms populating a long lived dark state of the Λ-type atom-light coupling scheme. These barriers appear in regions of steep change in the dark state due to the geometric scalar (Born-Huang) potential [18,21,22]. The Λ scheme has a single dark state, so no spin (or quasi-spin) degree of freedom is involved for the atomic motion in the dark state manifold affected by the sub-wavelength barriers.
Here we consider a way of creating a periodic array of spin-dependent sub-wavelength barriers that act differently on atoms depending on their internal state. For this we analyze a tripod atom light coupling scheme [23][24][25][26] characterized by two dark states playing the role of quasi-spin states. Previously, the tripod scheme was studied for creating a monopole field [24,25,27] and homogenous spin-orbit coupling [28][29][30][31][32], as well as for cooling atoms and ions [33,34], but not for producing a lattice characterized by spin-dependent sub-wavelength barriers. Inclusion of the spinor degree of freedom provides new possibilities for controlling the spectral and kinetic properties of atoms in the lattice. This paper is structured as follows. In the following Sec. 2 we introduce the tripod scheme of atom-light coupling characterized by two degenerate dark states, and assume the spatially periodic atom-light coupling. In Sec. 3 we consider the adiabatic motion of atoms in the manifold of two dark states playing the role of quasi-spin states. The atomic motion is then affected by a periodic array of sub-wavelength potential barriers acting mostly on the second dark state, as well as by a matrix-valued vector potential inducing transitions between the dark states. In Sec. 4 the spectrum of the tripod lattice is analyzed using the exact diagonalization of the Hamiltonian including all four atomic states, as well as exactly solving the eigenvalue problem for the adiabatic atomic motion projected onto the dark state manifold. The latter spectrum is also compared with the one relying on the scattering approach within the dark state manifold. The concluding Sec. 5   We consider atom-light interaction in the tripod scheme in which three laser beams characterized by the Rabi frequencies Ω j = Ω j (x), with j = 1, 2, 3, couple three atomic ground states |j to an excited state |0 detuned by ∆ and characterized by the decay rate Γ, as illustrated in Fig. 1. The operator describing such coupling is given by in the rotating frame for the atomic internal stateŝ where we have set = 1, so the energy is measured in units of frequency.
To create spin-dependent sub-wavelength barriers along the x axis, we consider the following configuration of Rabi frequencies where α is the spatial phase difference between the standing waves Ω 2 (x) and Ω 3 (x), and a is the lattice spacing. The amplitudes of the first two laser fields Ω 1 (x) and Ω 2 (x) (to be referred to as the probe fields) are taken to be much smaller than the amplitude of the control field Ω 3 (x): = Ω p /Ω c 1 .
In that case the control field is dominant everywhere except around its zeros at x = na/2, with n being an integer. This results in a periodic array of state-dependent sub-wavelength barriers at x = na/2, as will be shown in Sec. 3.2.2. At each barrier the ratio between the probe beams Ω 2 /Ω 1 = (−1) n cos α depends on the angle α, so by changing α one can significantly alter the spin-dependence of the sub-wavelength lattice. It is convenient to express the real valued Rabi frequencies Ω j ≡ Ω j (x) entering the atomlight interaction operator in terms of spherical angles 0 ≤ θ ≤ π and 0 ≤ φ < 2π, given by (see Fig. 2) where Ω = Ω 2 1 + Ω 2 2 + Ω 2 3 is the total Rabi frequency. The zeros of the control field correspond to θ = π 2 . By including the kinetic energy along the x direction, the full atomic Hamiltonian readŝ where m is the atomic mass. The corresponding state-vector describes a combined atomic internal and center of mass motion, where ψ j (x) is a wavefunction for the atomic center of mass motion in the jth internal state.

Spatial shift by half the lattice constant
The Hamiltonian (7) is invariant with respect to the spatial shift by the lattice constant a, such that Ĥ , e ipxa = 0. On the other hand, since the Rabi frequencies Ω 1 (x) and Ω 2,3 (x) given by Eqs.
(2)-(4) are, respectively, symmetric and anti-symmetric with respect to the spatial shift by half the lattice constant a/2, the Hamiltonian also commutes with the combined shift operatorT a/2 =Û e ipxa/2 , which involves both a spatial a/2 shift and a change of atomic internal state, the latter described by the operator withÛ 2 =Î. Therefore the Hamiltonian and the translation operatorT a/2 have a common set of eigenstates where g (q) (x) is invariant with respect to the combined shift operator In that case the components of the expanded state-vector (8) read j (x)e iqx , with g (q) where the upper and lower signs correspond to j = 2, 3 and j = 0, 1, respectively. The solution given by Eqs. (10)- (12) is characterized by the quasi-momentum q covering the first Brillouin zone (1BZ) (−2π/a, 2π/a] which is twice larger than the 1BZ of the original Bloch solution for a lattice with periodicity a. Neglecting the decay rate Γ, the system is governed by time-reversal symmetry, so the complex conjugation of the components g (q) j (x), accompanied by reversion of the quasi-momentum q → −q, does not change the eigenvalue equation giving Such a condition is well maintained for atomic dynamics in the dark state manifold in which atomic decay is suppressed, as one can see in the energy dispersions presented in Figs 4 and 8 in Sec. 4.

Dark states
The atom-light operatorV (x) couples the excited state |0 to a special superposition of the ground states |B ≡ |B (x) , known as the bright state: Two other superpositions of the ground states |D 1 ≡ |D 1 (x) and |D 2 ≡ |D 2 (x) are referred to as the dark states [23][24][25]. They are orthogonal between each other D 1 |D 2 = 0 and to the bright state B |D l = 0, and thus are immune to the atom-light coupling, i.e.V (x) |D l = 0 with l = 1, 2.
The dark states are not uniquely defined. One of them, such as |D 1 , can be chosen to be an arbitrary state vector orthogonal to |B . In that case the second dark state |D 2 should be orthogonal to both |B and |D 1 . A common choice for the dark states is the following [23][24][25]: the first state |D 1 is taken to be the dark state of the Λ system composed of the bare states |1 and |2 : so the second dark state is then given by The first dark state (14) is characterized by the angle φ and thus depends exclusively on the ratio between the Rabi frequencies of the first two laser beams Ω 1 /Ω 2 . The second dark state (15) additionally depends on the Rabi frequency of the third laser beam Ω 3 through the angle θ. As we will see in Sec.3.2, this provides steep potential barriers for atoms in the second dark state at the zero points of the dominant Rabi frequency Ω 3 , where the derivative θ becomes much larger than the inverse lattice constant 1/a.

Projection onto the dark state manifold
If the total Rabi frequency Ω is large enough compared to the characteristic energy of the atomic center of mass motion, one can employ the adiabatic approximation by restricting the atomic motion to the manifold of degenerate dark states described by a projected state vector where ψ D l (x) is the wave-function for the atomic center of mass motion in the l-th dark state, and I D = 2 l=1 |D l D l | is the unit operator for projection onto the dark state manifold. Thus the atomic motion is described by a two-component wave-function The corresponding 2 × 2 Hamiltonian governing the evolution of such a spinor wave-function ψ D (x) reads where the 2 × 2 matrices A D (x) and V D (x) are the geometric vector and scalar potentials presented in Ref. [25] for the general tripod atom-light coupling scheme. For the real valued Rabi frequencies Ω j considered here, the geometric vector potential A D (x) is proportional to the Pauli matrix σ y : On the other hand, the geometric scalar potentialV D (x) can be expressed as a product of the row vector c 1 , c 2 = φ sin θ , −θ and its corresponding column vector: Here the primes label the spatial derivatives of the angles θ and φ that parametrize the Rabi frequencies of the laser fields through Eq. (6). For the Rabi frequencies given by Eqs.
(2)-(4), both dark states are invariant with respect to the combined a/2 shift:T a/2 |D l (x) = |D l (x) . Consequently the periodicity of the geometric vector and scalar potentials is half of the original lattice constant a: so the eigenstates of the 2 × 2 matrix Hamiltonian H D are the two component Bloch spinors

Spin-dependent sub-wavelength barriers
The state-dependent scalar potential (20) affects the atoms in a special superposition |D ≡ |D (x) of the original dark states |D 1 ≡ |D 1 (x) and |D 2 ≡ |D 2 (x) proportional to The coefficient c 2 = −θ is associated with changes in the second dark state, and its absolute value is much larger than the inverse lattice constant 1/a at the zeros of the dominant Rabi frequency Ω 3 (x) (appearing at x = na/2 with integer n), as one can see in Fig. 3. On the other hand, the angle φ is defined by the ratio of the first two Rabi frequencies Ω 2 /Ω 1 and thus changes smoothly, so the absolute value of the coefficient c 1 = φ sin θ is of the order of 1/a or smaller. This means that at x = na/2, the factor c 1 associated with the first dark state plays a considerably less important role in the superposition (22) than the factor c 2 associated with the second dark state, as shown in Fig. 3. Furthermore, one can neglect c 1 also beyond the points where x = na/2. Using these arguments, one can simplify Eq. (20) by omitting c 1 , leading to the following approximate scalar potential Calling on the relation θ = −(cos θ) / sin θ, together with Eq. (6) for cos θ = Ω 3 /Ω and Eqs.
(2)-(4) for Ω 1,2,3 , one finds where k = 2π/a. Since = Ω p /Ω c 1, Eq. (24) represents a periodic set of barriers positioned at x = na/2, where neighboring barriers are separated by a distance of a/2, much larger than the characteristic width of each barrier. Thus one can represent Eq. (24) as a sum of narrow potential barriers is the effective ratio of the Rabi frequencies.
In this way, only atoms in the second dark state are scattered by the periodic array of sub-wavelength potential barriers located at x = na/2. The barriers appear due to rapid changes in the mixing angle θ featured in Eq. (15) for the second dark state. On the other hand, between the barriers the vector potential A D (x) ∝ σ y describes transitions between the dark states (see Fig. 3). These two processes determine the shape of the energy dispersion. Approaching the barriers at x = na/2, the second dark state reduces to for |x − na/2| a. Consequently, if cos α = 0, the dark state |D 2 ≡ |D 2 (x) represents a different physical state in the vicinity of the barriers at x = na/2 for even and odd values of n, leading to scattering of different internal states at even and odd barriers.
Applying the methods presented in the Supplementary material of Ref. [5], the atom in the second dark state with the kinetic energy E tunnels over an individual potential barrier v (x) given by Eq. (25) with the following transition and reflection amplitudes t = t(E) and r = r(E) for aQ /π 2 1: where Q = √ 2mE is the momentum corresponding to the energy E. We will use these amplitudes t and r to get an approximate dispersion represented by the 4 × 4 eigenvalue equation (30) for a periodic array of spin-dependent scatterers.

Dispersion for a periodic array of spin-dependent scatterers
To gain deeper insight into the tripod lattice, we present an effective approach based on scattering by individual barriers. Such an analysis relies on the fact that the distance between neighboring potential barriers positioned at x = na/2 is much larger than the characteristic width of each barrier ζ =˜ a/2π. Further away from the potential barriers, one can neglect the potential V D (x) in the Hamiltonian (18). In this spatial region, the atomic motion is affected only by the vector potential A D = φ cos θσ y , which induces transitions between the dark states due to the Pauli matrix σ y . (The spatial dependence of φ cos θ is plotted in Fig. 3 for different values of α.) Therefore the Bloch eigensolution characterized by the quasi-momentum q reads for na/2 − ζ < x < (n + 1) a/2 + ζ, i.e. further away from the barriers: with γ n (x) = x na/2 dxφ x cos θ, where σ y is a Pauli matrix acting on the two component spinors B (q) and C (q) .
Connecting the solution (29) at different sides of the barriers via the spin-dependent transition matrix W Q given by Eq. (31), one arrives at the eigenvalue equation The 4 × 4 matrix W Q has the following block structure with where I Q = diag e −iQa/2 , e iQa/2 is a diagonal 2 × 2 matrix reflecting the phases acquired during the forward and backward propagation of atoms between the barriers. Here also is a matrix for scattering of atoms in the second dark state by a single barrier, t(E) and r(E) being the corresponding transmission and reflection coefficients, given by Eq. (28). The atoms in the first dark state are not affected by the potential barriers, hence the free propagation matrix I Q that is featured in the first line of the block diagonal matrix diag (I Q , I Q M ) in Eq. (31). Finally, the operator e −iγ 0 σy entering W Q describes a transition between the dark states during propagation between the barriers. The 4 × 4 eigenvalue equation (30) provides a relation between the quasi-momentum q and the momentum Q = Q q defining the dispersion E q = Q 2 q /2m. In Eq. (32) for γ 0 , the function cos θ is close to unity everywhere in the integrand except for a narrow range of distances close to the boundaries at x = 0 and x = a/2, where the control field Ω 3 is no longer dominant. Thus, in the zero order of approximation, cos θ can be omitted in Eq. (32), giving γ 0 ≈ φ(0) − φ(a/2). In particular, for α = 0, one has φ(0) = π/4 and φ(a/2) = −π/4, so that γ 0 ≈ π/2. In that case the eigenvalue equation given by Eqs.
Thus one arrives at the following approximate eigenvalue equation where the 2 × 2 matrix −I 2 Q M describes the scattering of atoms in the second dark state at even or odd barriers, separated by the distance a and positioned at na or (n + 1/2) a. Hence the scattering of atoms in the second dark state takes place at every second barrier. This is because propagation of atoms between neighboring barriers separated by a/2 is associated with the matrix e iγ 0 σy , converting the second atomic dark state into the first one, which is not affected by the barrier.
Equation (35) has the same form as the corresponding scattering equation for the Λ subwavelength lattice [5], but with a twice larger periodicity. Furthermore, in comparison to the Λ-type coupling scheme, the quasi-momentum q is shifted by π/a due to the minus sign in the scattering matrix −I 2 Q M = −I 2Q M , the sign change appearing due to the interchange between the dark states during propagation between the barriers. Thus in contrast to the Λ setup [5,6], the quasi-momentum q = 0 does not correspond to the energy minimum in the lowest energy band, as one can see in the exactly calculated energy dispersion for α = 0 in Fig. 4(c).
Since we are considering the expanded 1BZ covering the range −2π/a < q ≤ 2π/a, the dispersion described by Eq. (35) is twice degenerate, the eigenenergies being the same for q = q and q = q + 2π/a. In a more accurate analysis, one should take into account that γ 0 = π/2. This means that the energy dispersion no longer repeats itself by shifting the quasimomentum q → q + 2π/a in the expanded 1BZ, which is apparent in the exactly calculated dispersion in Figs. 4(c)-(d).
In the limit where α → π/2, the angle γ 0 defined by the integral in Eq. (32) goes to zero, as in that case the vector potential shown in Fig. 3 averages to zero after integration, giving γ 0 = 0. As such, for α = π/2, the eigenvalue equation given by Eqs. (30)- (31) shows that atoms move freely in the first dark state, while atoms in the second dark state are scattered at every barrier (in contrast to scattering at every second barrier for α = 0). This is confirmed by the exact calculations of the spectrum in this limit presented in Sec. 4.3.

Analysis of the spectrum
The spectrum of the tripod lattice is analyzed using the exact diagonalization of the Hamiltonian including all four atomic states using numerical algorithms of Ref. [36][37][38][39][40], as well as exactly solving the eigenvalue problem for the adiabatic atomic motion projected onto the dark state manifold. The latter spectrum is also compared with solutions of the 4 × 4 eigenvalue equation (30) relying on the scattering approach within the dark state manifold.

Energy dispersion for small and moderate phase α
Parts (a) and (b) of Fig. 4 display exact calculations of the three lowest dispersion branches for small values of the angle α characterizing the phase of the Rabi frequency of the second probe beam Ω 2 (x) with respect to that of the control beam Ω 3 (x) in Eq. (3). Black lines correspond to the real part E q,s of the complex eigenvalues E q ≡ E q,s = E q,s + iE q,s , and red areas indicate the widths of the dispersion lines represented by the imaginary part E q,s , where the index s = 1, 2, . . . labels the dispersion bands. The minima of the dispersion branches E q,s are seen to fit well with the discrete spectrum of a square well with a width equal to a, which is twice larger than the spacing a/2 between the neighboring potential barriers:  Figures 4(c)-(d) display the enlarged dispersion of the first energy band s = 1, showing a characteristic double well shape with a smaller height at the center q = 0 than at the edges qa = ±2π of the expanded 1BZ. The linewidths are much smaller than the width of the energy bands. Losses are negligible in the lower part of the dispersion curves E q,s , corresponding to the minima of the double well dispersion at qa ≈ ±π for odd bands s = 1, 3, . . ., as well as at the center or edges of the expanded 1BZ at qa = 0, ±2π for even bands s = 2, 4, . . .. As discussed in Sec. 3.2.3, this is opposite to the Λ scheme in which both the real part of the dispersion and losses are minimum at q = 0 for odd bands [5]. The dispersion shown in Fig. 4 can be well described in terms of the dissipative tight binding model including the nearest neighbor (NN) and the next nearest neighbor (NNN) couplings: where J n = J n +iJ n with |J n | |J n |. Interestingly, the NNN coupling is much more significant than the NN coupling, |J 2 | |J 1 | and |J 2 | |J 1 | for small α. Furthermore, J 1 and J 2 have opposite signs, as one can see in Figs. 9 and 10 presented in Sec. 4.3. In particular, for the lowest band (s=1) one has J 1 < 0 and J 2 > 0, which explains the double well type dispersion featured in Fig. 4(c)-(d). Note also that the real and imaginary parts of the NNN matrix element J 2 and J 2 have opposite signs, with J 2 being negative (positive) for even (odd) bands s. Therefore for even bands losses vanish at the dispersion minimum at q = 0, whereas losses for odd bands are negligible in the vicinity of the minima of the double well dispersion at qa ≈ ±π. The solid dispersion lines in Figs. 4 and 5 agree well with the dispersions obtained for the adiabatic motion within the dark state manifold (dotted lines). More generally, we have checked that the effective dark state Hamiltonian (18) recreates the exact band structure for any α, not necessarily small, as shown in Fig. 8 of Sec. 4.3, provided the detuning ∆ is considerably smaller than the amplitude of the control field Ω c , and the ratio between the amplitudes of the probe and control fields is not too small, = Ω p /Ω c 0.05. For smaller values of , non-adiabatic losses become important due to significant changes in the composition of the atomic dark states at the potential barriers. Note that in Fig. 4, deviations of the adiabatic dispersion curves (dashed lines) from the exact dispersion (solid lines) are mainly due to non-zero values of ∆ and Γ, the difference becoming minimum for zero detuning ∆ = 0 shown in Fig. 5.
The effective 4 × 4 eigenvalue equation (30) relying on the scattering approach within the dark state manifold reproduces the band structure of the odd energy bands accurately for α π/4, ∆ Ω c and 0.05, as one can see in parts b) and d) of Fig. 5. While the overall band structure and band gaps shown in parts a) and b) of Fig. 5 are maintained for any α, the fine details of the even band energy dispersion curves are not precisely reproduced. Additionally, we note that the solutions of the scattering approach become less accurate for the higher energy bands.

Wannier functions
To gain more insight into the tight binding model, we have considered the multi-component Wannier function localized around x = na/2, with n being an integer: where the alternating factor (−1) n reflects the fact that the bright and excited states alter their sign after applying a combined a/2 shift operatorT a/2 , defined in Sec. 2.2: As discussed in more detail in Appendix A, the Wannier state vector |W n (x) maximally localized around x = na/2 is obtained by superimposing the Bloch states of an individual band in the extended Brillouin zone with the phase factors e −iqna/2 . Thus the Wannier functions are strictly orthogonal for different values of n, despite the Hamiltonian being non-Hermitian due to the imaginary part in the excited state energy. The multicomponent Wannier functions W n (x) ≡ W s,n (x) with the same index n but belonging to different Bloch bands s are orthogonal to great accuracy, as we are considering the atomic motion in the dark state manifold with only a tiny contribution from the lossy excited state. Figure 6 shows the Wannier functions W D 1,2 (x), W B (x) and W 0 (x) centered at x = 0 for the first and second dispersion bands (s = 1, 2) and α = 0. Despite the reduction of the unit cells to a/2 due to expansion of the 1BZ, the Wannier functions extend over a distance equal to the lattice constant a, as the atoms are confined between the barriers located at x = ±a/2. These barriers scatter atoms in the second dark state |D 2 (x) , so the Wannier function W D 2 (x) dominates over W D 1 (x) in the vicinity of the barriers. Between the barriers there is an exchange among the two dark states due to the vector potential term A D (x) ∝ σ y . The second dark state is thus converted to the first one at the center x = 0, where the Wannier function W D 2 (x) vanishes, whereas W D 1 (x) reaches its maximum magnitude for the first band, as shown in Fig. 6(a). Since only the second dark state is affected by the potential barriers, the barrier at x = 0 has little influence over the multi-component Wannier function centered at x = 0. The Wannier functions of the bright and excited states (W B (x) and W 0 (x)) emerge due to non-adiabatic transitions. Thus they are localized close to the barriers at x = ±a/2, where the Rabi frequency of the control field goes to zero (see Figs. 6(b),(d)). A tiny residual non-adiabatic effect due to the barrier at the center of the Wannier function appears as a small bump at x = 0 in Fig. 6(d) for the Wannier functions of the bright and excited states.
Atoms populating the multi-component Wannier function W n (x) centered at x = na/2 can tunnel to the next neighboring Wannier functions W n±2 (x) centered at x = (n ± 2) a/2 through the barriers located at x = (n ± 1) a/2. These processes are described in terms of the matrix elements J 2 of the Hamiltonian between such next neighboring multi-component Wannier functions. A small imaginary part J 2 emerges due to a small non-adiabatic contribution by the lossy excited state shown in Figs. 6(b),(d). On the other hand, the real part J 2 is mostly due to the Wannier function of the second dark state W D 2 (x), which is anti-symmetric with respect to the center x = 0 for the first band (see Fig. 6(a)), leading to the dispersion minima at non-zero quasi-momentum q in the lowest band shown in Fig. 4(c)-(d).

Brick-wall lattice
The system can be visualized in terms of a two layer brick-wall lattice shown in Fig. 7 (a), in which the bricks in the lower (upper) layer represent the multi-component Wannier function W n (x) with odd (even) values of n. NNN coupling between the Wannier functions described by the matrix element J 2 corresponds to horizontal tunneling between the adjacent bricks of individual layers. Additionally, there is an inter-layer tunneling described by the NN coupling element J 1 , which appears due to a residual effect of the barrier at the center of the Wannier function. For the odd bands, the effect of NN coupling becomes more important as α grows, leading to a more shallow double well dispersion in Fig. 4(d) in comparison to Fig. 4(c). In fact, for larger values of α, the vector potential term converts the second dark state into the first one to a lesser extent at the center of the Wannier function, which increases scattering of atoms at x = na/2, and thus NN coupling becomes more significant.
By further increasing the phase α, the band structure undergoes significant modifications, as illustrated in Fig. 8. For the first band, NN tunneling grows quickly with α and eventually overcomes NNN tunneling. The odd bands (s = 1, 3, . . .) are very sensitive to α and eventually merge to form a free particle energy dispersion branch as α → π/2. In contrast, the even bands (s = 2, 4, . . .) remain relatively intact with increasing α; as α approaches π/2, they transform into energy bands of Λ-like atom-light coupling [5,6], characterized by a twice smaller distance a/2 between the barriers. This can be understood by observing that, for α = 0, the Wannier functions of both dark states vanish at the center x = na/2 for even bands (see Fig. 6(c)), so these Wannier functions are not as sensitive to the increasing effect of the barrier at x = na/2 as α grows. In the limit where α → π/2, tunneling to the left or right takes place at every barrier located at na/2 for atoms in the second dark state, as shown in the upper part of Fig. 7(b), whereas atoms move freely in the first dark state (lower part). This is because for α = π/2 there is no longer a vector potential type term converting the first dark state to the second one (and vice versa) during atomic propagation between the barriers. Note that the localization centers of the Wannier functions become shifted by a/4 in the upper part of Fig. 7(b) as a consequence of scattering at every barrier.  Figure 9(a) shows the α dependence of the tight binding tunneling parameters J v between the first four neighboring lattice sites (v=1,2,3,4) for the first energy band (s = 1). As α becomes non-zero, one can see a steep increase in the tunneling matrix element J v for odd values of v, whereas for even v the parameter J v start increasing at much larger values of α. When α approaches π/2, all the tunnelling parameters J v are a few orders of magnitude larger than in the case where α = 0, as the atoms loaded into the bands with odd s behave like free particles for α → π/2. We also note that the 4 × 4 eigenvalue equation (30) relying on the scattering approach accurately predicts the free particle behavior of bands with odd s for α → π/2. For the second energy band (s=2) shown in Fig. 9(b), the phase α has a very different effect on the tunnelling parameters. As α approached π/2, all the tunnelling parameters J v vanish, with the exception of the NN tunnelling corresponding to v = 1. As such, atoms loaded into the even bands are effectively confined by steep potential barriers at x = na/2 separated by a/2.
Up to now we have investigated the situation where the detuning ∆ is considerably smaller than the amplitude of the control field Ω c . Going beyond this assumption, in Fig. 10 we observe that the NN tunneling matrix element J 1 does not change with ∆. On the other hand, the NNN tunneling element J 2 decreases with ∆ > 0 and undergoes sign reversal after reaching the zero point, in which NNN tunneling is suppressed due to destructive interference between different components of the multi-component Wannier functions. When the NNN tunneling element J 2 vanishes, only the small NN tunneling parameter J 1 contributes to the energy dispersion, which results in an almost flat dispersion band. A similar effect can be observed also for the sub-wavelength lattice obtained using the Λ coupling scheme, where the NN tunnelling element also goes to zero at a certain value of the detuning. For the tripod scheme, the detuning causes the NNN tunnelling element to vanish in the same fashion, as illustrated in Fig. 10, while the NN element is largely unaffected because the Wannier function is spread over two unit cells. Additionally, the tight binding parameters can be controlled with the Rabi frequency ratio = Ω p /Ω c , as in the Λ sub-wavelength lattice [19]. 0

Concluding remarks
We have analyzed the tripod atom light coupling scheme characterized by two dark states playing the role of quasi-spin states. It is demonstrated that by properly choosing the laser fields which couple the atomic internal states, one can create a lattice with spin-dependent sub-wavelength barriers acting differently on atoms in different atomic dark states. Inclusion of the spinor degree of freedom allows to flexibly alter the atomic motion ranging from the tight-binding type atomic dynamics in the effective brick-wall type lattice shown in Fig. 7(a) to the free motion of atoms in one dark state and the tight binding lattice with a twice smaller periodicity for atoms in another internal state, Fig. 7 (b). Between the two regimes, the spectrum undergoes significant changes controlled by the laser fields, as illustrated in Fig. 8.
The tripod scheme has already been experimentally realized by coupling three metastable atomic ground states to a common excited state [41,42]. In particular, it was shown that the tripod setup can provide a robust and variable beam splitter for a beam of metastable neon atoms [41]. On the other hand, the tripod scheme was used to implement non-Abelian adiabatic geometric transformations for cold fermionic gas of strontium-87 atoms by laser coupling of three hyperfine ground states |F g = 9/2; m g with m g = 9/2 , 7/2 , 5/2 to their common excited state |e = |F e = 9/2; m e = 7/2 using two co-propagating beams with opposite circular polarizations and a third linearly polarized beam orthogonal to the first two beams [42]. The tripod-type scheme has also been implemented for photonic systems [43]. The lattice with spin-dependent sub-wavelength barriers considered here can be implemented using a setup similar to that in Ref. [42], with the co-propagating circularly polarized waves replaced by the standing waves describing the Rabi frequencies Ω 2 (x) and Ω 3 (x) in Eqs. (3)-(4). In this way, it is feasible to produce the sub-wavelength tripod lattice using currently available experimental techniques.
The lattice with sub-wavelength barriers has already been implemented for the Λ scheme by adiabatically loading the ultracold atoms to the lowest band [19]. A similar adiabatic loading can be applied for the lattice with spin-dependent sub-wavelength barriers involving the tripod scheme. The atoms loaded to the lowest Bloch band will be maintained in it, as the difference in energies between the ground and first excited states is of the order of the recoil energy corresponding to temperature which greatly exceeds the nK temperature range of the ultracold atomic gas. Since the tripod scheme has two degenerate dark states, there is an extra possibility to populate specific dark states during the adiabatic loading. This can be done by transferring the ultracold atoms to a desired dark state by the tripod STIRAP (stimulated adiabatic Raman passage) [23,24] prior to the adiabatic loading.
The use of the tripod scheme to create a lattice of degenerate dark states opens new possibilities for spin ordering and symmetry breaking. In conventional spinor systems degeneracy is usually controlled by an external magnetic field that typically has noise on the scale of the spin-dependent interactions set either by exchange interactions or the native spin-dependent collisions. One might speculate that the present system would support new forms of magnetic ordering, or new routes to create established forms of magnetic order.
Note added. Shortly after submitting this manuscript to arXiv and SciPost, a related preprint appeared in arXiv [44] which also analyzes the spin-dependent sub-wavelength barriers using the tripod atom-light coupling scheme. The latter work complements the present study by covering some other aspects of the problem, including the case where the amplitudes of the probe fields Ω 1,2 (x) are not the same.
where N is the number of q values in the extended Brillouin zone −2π/a < q ≤ 2π/a, and |g (q) (x) is the periodic part of the Bloch state vector. Each term in the summation over q is demanded to be real [45], which concludes the optimization procedure. For the second energy band, the Wannier function Wn(x) vanishes at the center point x = na/2 representing a node (Fig. 12). As such, the optimal phase factors are determined by investigating the projected Wannier function at a point x = na/2 ± a/4 shifted from the center x = na/2 by a/4 to the right or the left. This gives Wn(na/2 ± a/4) = 1 √ N q e ±iqa/4 n|g (q) (na/2 ± a/4) , where, once again, all the components in the sum over q are optimized by demanding them to be real. As α approaches π/2, one must optimize the Wannier function of the second band in a similar fashion, but now centered at x = na/2 + a/4 and without any additional shifting. In that case the Wannier function of the second energy band localizes between the neighboring potential barriers situated at x = na/2 and x = (n + 1)a/2. For small α, the Wannier function Wn(x) = n |W n (x) shown in Figs. 11-12 represents the dominant contribution to the full Wannier state vector |W n (x) . Furthermore, the Wannier functions Wn(x) localized around the neighboring sites n = 0 and n = 1 are seen to have opposite signs due to theT a/2 symmetry of the Hamiltonian discussed in Sec. 2.2.
The methods described above work well for small or moderate values of α. Figure 13 shows the Wannier functions Wn(x) = n |W n (x) and Wn⊥(x) = n ⊥ |W n (x) of the first is a state vector orthogonal to |n . One can see in Fig. 13 that the contribution from Wn⊥(x) is tiny for α = 0 and increases with α, leading to the growth of the NN tunneling matrix element J 1 shown in Fig. 9(a). This is because for small angle α the state vector |n ⊥ almost completely overlaps with the state vector |n at the neighboring sites n = n±1, the latter |n providing a dominant contribution to the neighboring Wannier state vectors |W n±1 (x) . An additional increase in the angle α leads to further growth of Wn⊥(x) and subsequent spread of both Wn(x) and Wn⊥(x). This enables tunneling between more distant sites, as illustrated in Fig. 9(a). When α approaches 90 • , the description of the odd energy bands in terms of the Wannier functions becomes problematic as the particle dynamics in the odd bands reduce to free atomic motion, as one can see in Fig. 8(d). On the other hand, for α → 90 • the Wannier functions of even bands become localized between the neighboring barriers separated by a/2, like in the Λ scheme [5,6]. The Wannier functions of the dark, bright and excited states W n,D 1,2 (x), W n,B (x) and W 0 (x) are displayed in Fig. 6 of the main text. They are found by projecting the full Wannier state vector |W n (x) (obtained using the above methods) to the corresponding atomic internal states |D l ≡ |D l (x) , |B ≡ |B (x) and |0 . The approximate Wannier functions of the dark states W D 1,2 (x) calculated in terms of the adiabatic dark state Hamiltonian (18) are presented by the dotted lines in Fig. 6. They are optimized in the following way: 1. For the first band, we demand W D 1 (x) to be real at x = 0.
2. For the second band and, we take the superposition W D 1 (x) + W D 2 (x) at x = a/4 and demand it to be real.
The methods discussed here are implemented in the numerical codes that are available in the Supplementary Material [35].