Topological lattice using multi-frequency radiation

For the remainder of the manuscript we will use dimensionless units. All energies will be expressed in units of ${\hslash }\omega$, derived from the Floquet frequency ω; time will be expressed in units of inverse driving frequency ω⁻¹, denoted by τ = ωt; spatial coordinates will be expressed in units of inverse recoil momentum K₀⁻¹, denoted by lowercase letters (x, y) = K₀(X, Y). In these units, the Hamiltonian (1) takes the form

$$\begin{eqnarray}&&h(\tau )=\displaystyle \frac{H(\tau /\omega )}{{\hslash }\omega }={E}_{{\rm{r}}}{{\boldsymbol{k}}}^{2}+\displaystyle \frac{1}{2}{\boldsymbol{\Omega }}(\tau )\cdot {\boldsymbol{\sigma }},\end{eqnarray} \tag{ 4 }$$

where ${E}_{{\rm{r}}}={{\hslash }}^{2}{K}_{0}^{2}/(2M{\hslash }\omega )$ is the dimensionless recoil energy associated with the recoil wavenumber K₀; ${\boldsymbol{k}}={\boldsymbol{K}}/{K}_{0}$ is the dimensionless wavenumber. The dimensionless coupling

$$\begin{eqnarray}&&{\boldsymbol{\Omega }}(x,y,\tau )=(2\mathrm{Re}\,u(y,\tau ),2\mathrm{Im}\,u(y,\tau ),\beta x)\end{eqnarray} \tag{ 5 }$$

includes a combination of position-dependent detuning and Raman coupling. Here $\beta ={{\rm{\Delta }}}^{{\prime} }/({\hslash }\omega {k}_{0})$ describes the linearly varying detuning in dimensionless units; the function $u(y,\tau )={v}_{0}{\sum }_{n}\{\exp [{\rm{i}}(y-2n\tau )]$ $+\,\exp [{\rm{i}}(-y-(2n+1)\tau )]\}$ is a dimensionless version of the sum in equation (3) with ${v}_{0}={V}_{0}/({\hslash }\omega )$.

In the time domain the coupling given by equation (5) is

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{\boldsymbol{\Omega }}(\tau )\cdot {\boldsymbol{\sigma }}=\displaystyle \frac{1}{2}\beta x{\sigma }_{3}+\displaystyle \sum _{l}{v}_{l}(y)\delta (\tau -\pi l),\end{eqnarray} \tag{ 6 }$$

with

$$\begin{eqnarray}&&{v}_{l}(y)=\pi {v}_{0}[{{\rm{e}}}^{{\rm{i}}y}+{(-1)}^{l}{{\rm{e}}}^{-{\rm{i}}y}]| \downarrow \rangle \langle \uparrow | +{\rm{H}}.{\rm{c}}.\end{eqnarray} \tag{ 7 }$$

In this way we have separated the spatial and temporal dependencies in the coupling (6).

We continue our analysis by deriving an approximate Hamiltonian that describes the complete time evolution over a single period from τ = 0 −  to $\tau =2\pi -\epsilon$ with $\epsilon \to 0$. This evolution includes a kick v₀ at the beginning of the period τ₊ = 0 and a second kick v₁ in the middle of the period ${\tau }_{-}=\pi ;$ between the kicks the evolution includes the kinetic and gradient energies. In the full time period, the complete evolution operator is a product of four terms:

$$\begin{eqnarray}&&U(2\pi ,0)\equiv \mathop{\mathrm{lim}}\limits_{\epsilon \to 0}U(2\pi -\epsilon ,0-\epsilon )={U}_{0}{U}_{\mathrm{kick}}^{(1)}{U}_{0}{U}_{\mathrm{kick}}^{(0)}.\end{eqnarray} \tag{ 8 }$$

Here

$$\begin{eqnarray}&&{U}_{0}=\exp \left\{-{\rm{i}}\pi \left[{E}_{{\rm{r}}}{{\boldsymbol{k}}}^{2}+\displaystyle \frac{1}{2}{\sigma }_{3}\beta x\right]\right\}\end{eqnarray} \tag{ 9 }$$

is the evolution operator over a half period, generated by kinetic energy and gradient. The operator

$$\begin{eqnarray}&&{U}_{\mathrm{kick}}^{(l)}=\exp [-{\rm{i}}{v}_{l}(y)]\end{eqnarray} \tag{ 10 }$$

describes a kick at $\tau =l\pi$.

We obtain an effective Hamiltonian by assuming that the Floquet frequency ω greatly exceeds the recoil frequency, $1\gg {E}_{{\rm{r}}}$, allowing us to ignore the commutators between the kinetic energy and functions of coordinates in equation (8). We then rearrange terms in the full time evolution operator (8) and obtain (see appendix)

$$\begin{eqnarray}&&{U}_{{\rm{eff}}}=\exp [\{-{\rm{i}}2\pi ([{E}_{{\rm{r}}}{{\boldsymbol{k}}}^{2}+{v}_{{\rm{eff}}})]]\},\end{eqnarray} \tag{ 11 }$$

where v_eff is an effective coupling defined by

$$\begin{eqnarray}&&\exp (-{\rm{i}}2\pi {v}_{\mathrm{eff}})={{\rm{e}}}^{-{\rm{i}}\pi {\sigma }_{3}\beta x/2}{U}_{\mathrm{kick}}^{(1)}{{\rm{e}}}^{-{\rm{i}}\pi {\sigma }_{3}\beta x/2}{U}_{\mathrm{kick}}^{(0)}.\end{eqnarray} \tag{ 12 }$$

The function v_l(y) entering the kick operators ${U}_{\mathrm{kick}}^{(l)}$ is spatially periodic along the y direction with a period 2π. This period can be halved to π by virtue of a gauge transformation $U=\exp (-{\rm{i}}y{\sigma }_{3}/2)$. Subsequently, when exploring energy bands and their topological properties, this prevents problems arising from using a twice larger elementary cell. Following this transformation the evolution operator becomes

$$\begin{eqnarray}&&{U}_{\mathrm{eff}}=\exp \{-{\rm{i}}2\pi [{E}_{{\rm{r}}}{({\boldsymbol{k}}+{\sigma }_{3}{{\boldsymbol{e}}}_{y}/2)}^{2}+{v}_{\mathrm{eff}}]\},\end{eqnarray} \tag{ 13 }$$

where v_l(y) featured in the kick operators ${U}_{\mathrm{kick}}^{(l)}$ has now the spatial periodicity π along the y direction, i.e. it should be replace to

$$\begin{eqnarray}&&{v}_{l}(y)=\pi {v}_{0}[{{\rm{e}}}^{{\rm{i}}2y}+{(-1)}^{l}]| \downarrow \rangle \langle \uparrow | +{\rm{H}}.{\rm{c}}.\end{eqnarray} \tag{ 14 }$$

The algebra of Pauli matrices allows us to write the effective coupling ${v}_{\mathrm{eff}}({\boldsymbol{r}})$ featured in the evolution equations (12) and (13) as:

$$\begin{eqnarray}&&{v}_{\mathrm{eff}}({\boldsymbol{r}})=\displaystyle \frac{1}{2}{{\boldsymbol{\Omega }}}_{\mathrm{eff}}({\boldsymbol{r}})\cdot {\boldsymbol{\sigma }},\end{eqnarray} \tag{ 15 }$$

where ${{\boldsymbol{\Omega }}}_{\mathrm{eff}}=({{\rm{\Omega }}}_{\mathrm{eff},1},{{\rm{\Omega }}}_{\mathrm{eff},2},{{\rm{\Omega }}}_{\mathrm{eff},3})$ is a position-dependent effective Zeeman field which takes the analytic form

$$\begin{eqnarray}&&\exp (-{\rm{i}}2\pi {v}_{\mathrm{eff}})={q}_{0}-{\rm{i}}{q}_{1}{\sigma }_{1}-{\rm{i}}{q}_{2}{\sigma }_{2}-{\rm{i}}{q}_{3}{\sigma }_{3}.\end{eqnarray} \tag{ 16 }$$

Here q₀, q₁, q₂ and q₃ are real functions of the coordinates (x, y), allowing to express the effective Zeeman field as

$$\begin{eqnarray}&&{{\boldsymbol{\Omega }}}_{\mathrm{eff}}={\pi }^{-1}\displaystyle \frac{{\boldsymbol{q}}}{| | {\boldsymbol{q}}| | }\arccos \,{q}_{0},\end{eqnarray} \tag{ 17 }$$

where ${\boldsymbol{q}}$ is a shorthand of a three dimensional vector (q₁, q_2, q₃). In general the equation (16) gives multiple solutions that correspond for different Floquet bands. Our choice (17) picks only to the two bands that lie in the energy window from $-1/2$ to 1/2 covering a single Floquet period.

Comparing (12) and (16) and multiplying four matrix exponents give explicit expressions

$$\begin{eqnarray}&&{q}_{0}=\cos {f}_{1}\cos {f}_{2}\cos (\pi \beta x),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{q}_{1}=\sin {f}_{1}\cos {f}_{2}\cos (y+\pi \beta x)-\cos {f}_{1}\sin {f}_{2}\sin (y),\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{q}_{2}=\sin {f}_{1}\cos {f}_{2}\sin (y+\pi \beta x)+\cos {f}_{1}\sin {f}_{2}\cos (y),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{q}_{3}=\cos {f}_{1}\cos {f}_{2}\sin (\pi \beta x)-\sin {f}_{1}\sin {f}_{2}\end{eqnarray} \tag{ 21 }$$

with

$$\begin{eqnarray}&&{f}_{1}(y)=2\pi {v}_{0}\cos (y),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{f}_{2}(y)=2\pi {v}_{0}\sin (y).\end{eqnarray} \tag{ 23 }$$

These explicit expressions show that the resulting effective Zeeman field (17) and the associated effective coupling (15) are periodic along both ${{\boldsymbol{e}}}_{x}$ and ${{\boldsymbol{e}}}_{y}$, with spatial periods ${a}_{x}=2/\beta$ and ${a}_{y}=\pi$ respectively. Therefore, although the original Hamiltonian containing the spin-dependent potential slope $\propto x{\sigma }_{3}$ is not periodic along the x direction, the effective Floquet Hamiltonian is. The spatial dependence of the Zeeman field components ${{\rm{\Omega }}}_{\mathrm{eff},1}$, ${{\rm{\Omega }}}_{\mathrm{eff},2}$ and ${{\rm{\Omega }}}_{\mathrm{eff},3}$ is presented in the figure 2 for β = 0.6 giving an approximately square unit cell. In figure 2 we select v₀ = 0.25 where the absolute value of the Zeeman field ${{\rm{\Omega }}}_{\mathrm{eff}}$ is almost uniform, as is apparent from the nearly flat adiabatic bands shown in figure 3 below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before moving further to an explicit numerical analysis of the band structure, we develop an intuitive understanding by performing an adiabatic analysis of motion governed by effective Hamiltonian

$$\begin{eqnarray}&&{h}_{\mathrm{eff}}({\boldsymbol{r}})={E}_{{\rm{r}}}{({\boldsymbol{k}}+{\sigma }_{3}{{\boldsymbol{e}}}_{y}/2)}^{2}+\displaystyle \frac{1}{2}{{\boldsymbol{\Omega }}}_{\mathrm{eff}}\cdot {\boldsymbol{\sigma }}\,\end{eqnarray} \tag{ 24 }$$

featured in the evolution operator U_eff, equation (13). The coupling field ${{\boldsymbol{\Omega }}}_{\mathrm{eff}}({\boldsymbol{r}})$ is parametrized by the spherical angles $\theta ({\boldsymbol{r}})$ and $\phi ({\boldsymbol{r}})$ defined by

$$\begin{eqnarray}&&\cos \theta =\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{eff},3}}{{{\rm{\Omega }}}_{\mathrm{eff}}},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&\tan \phi =\displaystyle \frac{{{\rm{\Omega }}}_{\mathrm{eff},2}}{{{\rm{\Omega }}}_{\mathrm{eff},1}}.\end{eqnarray} \tag{ 26 }$$

This gives the effective coupling [12]

$$\begin{eqnarray}\displaystyle \frac{1}{2}{{\boldsymbol{\Omega }}}_{\mathrm{eff}}\cdot {\boldsymbol{\sigma }}=\displaystyle \frac{1}{2}{{\rm{\Omega }}}_{\mathrm{eff}}\left[\begin{array}{cc}\cos \theta & {{\rm{e}}}^{-{\rm{i}}\phi }\sin \theta \\ {{\rm{e}}}^{{\rm{i}}\phi }\sin \theta & -\cos \theta \end{array}\right],\end{eqnarray} \tag{ 27 }$$

characterized by the position-dependent eigenstates

$$\begin{eqnarray}&&| + \rangle =\left(\begin{array}{c}\cos (\theta /2)\\ {{\rm{e}}}^{{\rm{i}}\phi }\sin (\theta /2)\end{array}\right),\,| - \rangle =\left(\begin{array}{c}-{{\rm{e}}}^{-{\rm{i}}\phi }\sin (\theta /2)\\ \cos (\theta /2)\end{array}\right).\end{eqnarray} \tag{ 28 }$$

The corresponding eigenvalues

$$\begin{eqnarray}&&{v}_{\pm }({\boldsymbol{r}})=\pm \displaystyle \frac{1}{2}{{\rm{\Omega }}}_{\mathrm{eff}},\end{eqnarray} \tag{ 29 }$$

are shown in figure 3 for various value of the Raman coupling v₀. As one can see in figure 3, for v₀ = 0.25 the resulting bands ${v}_{\pm }({\boldsymbol{r}})$ (adiabatic potentials) are flat and have a considerable gap $\approx \omega /2$, a regime suitable for a description in terms of an adiabatic motion in selected bands [27].

As in [28], we consider the adiabatic motion of the atom in one of these flat adiabatic bands with the projection Schrödinger equation that includes a geometric vector potential

$$\begin{eqnarray}&&{{\boldsymbol{A}}}_{\pm }({\boldsymbol{r}})=\pm \displaystyle \frac{1}{2}(\cos \theta -1){\rm{\nabla }}\phi .\end{eqnarray} \tag{ 30 }$$

This provides a synthetic magnetic flux density ${{\boldsymbol{B}}}_{\pm }({\boldsymbol{r}})={\rm{\nabla }}\times {{\boldsymbol{A}}}_{\pm }({\boldsymbol{r}})$. The geometric vector potential ${{\boldsymbol{A}}}_{\pm }({\boldsymbol{r}})$ may contain Aharonov–Bohm type singularities, that give rise to a synthetic magnetic flux over an elementary cell

$$\begin{eqnarray}&&{\alpha }_{\pm }=-\sum {\oint }_{\mathrm{singul}}{\rm{d}}{\boldsymbol{r}}\cdot {{\boldsymbol{A}}}_{\pm }({\boldsymbol{r}}).\end{eqnarray} \tag{ 31 }$$

The singularities appear at points where $\theta =\pi$, where the angle ϕ and its gradient ${\rm{\nabla }}\phi$ are undefined and $\cos \theta =-1$. The term $\cos \theta -1$ in (30) is nonzero and does not remove the undefined phase ${\rm{\nabla }}\phi$. Our unit cell contains two such singularities located at ${\boldsymbol{r}}=({a}_{x},3{a}_{y})/4$ and ${\boldsymbol{r}}=(3{a}_{x},{a}_{y})/4$, containing the same flux, so that they do not compensate each other, giving the synthetic magnetic flux $\pm 2\pi$ in each unit cell. Note that usually the optical lattices are sufficiently deep, and the $\pm 2\pi$ flux per elementary is topologically trivial. In that case the tight binding model can be applied, with the tunneling taking place only between the nearest-neighboring sites of the square plaquette. The $\pm 2\pi$ flux over the square plaquette can then be eliminated by a gauge transformation. Yet if the lattice is shallow enough, the tight binding model is not applicable and the above arguments do not work. In the present situation, the most interesting topological lattice appears for a flat adiabatic trapping potential shown by a solid red curve in figure 3. In such a situation there are no well defined lattice sites, and the $\pm 2\pi$ flux per elementary cell results in topologically non-trivial bands explored in the next subsection.

For a weak coupling (such as v = 0.05) the geometric flux density ${\boldsymbol{B}}({\boldsymbol{r}})\equiv {{\boldsymbol{B}}}_{\pm }({\boldsymbol{r}})$ is concentrated around the intersection points of the gradient slopes shown in in figure 3 and has a very weak y dependence. With increasing the coupling v, the flux extends beyond the intersection areas and acquires a y dependence. Figure 4 shows the geometric flux density ${\boldsymbol{B}}({\boldsymbol{r}})\equiv {{\boldsymbol{B}}}_{+}({\boldsymbol{r}})$ for the strong coupling (v₀ = 0.25) corresponding to the most flat adiabatic bands. In this regime the flux develops stripes in the x direction and has a strong y dependence. For the whole range of coupling strengths $0\leqslant {v}_{0}\leqslant 1/2$ the total synthetic magnetic flux per unit cell is 2π and is independent of the Floquet frequency ω and the gradient β.

Zoom In Zoom Out Reset image size

Now let us discuss the effect of an extra spin-independent trapping potential. The present scheme requires a large spin-dependent energy gradient which would have a huge influence on the relative trapping for the two spin states without the Raman coupling or for a weak Raman coupling. In that case one would expect that the stable positions for any trapped sample of the two spin states would live at entirely distinct locations, possibly with no overlap. Yet we are interested mostly in a sufficiently strong Raman coupling where the two spin states get mixed, and the atomic motion takes place in almost flat adiabatic potentials shown in red in figure 3. Therefore the atoms are no longer affected by the steep spin-dependent potential slopes, and the spin-independent trapping potential would not cause separation of different spin states. Instead, the extra spin-independent parabolic trapping potential would simply make the flat adiabatic trapping potentials parabolic. Of course, one needs to be all the time in the regime where the Raman coupling is strong compared to the characteristic energy of the spin-dependent potential slope. That is why we propose to introduce the spin-dependent potential gradient only at the final stage of the adiabatic protocol discussed in section 3.5.

We analyze the topological properties of this Floquet flux lattice by explicitly numerically computing the band structure and associated Chern number using the effective Hamiltonian (24) without making the adiabatic approximation introduced in section 3.3. Again the gradient of the original magnetic field is such that we approximately get a square lattice, β = 0.6. Furthermore, we choose the Floquet frequency to be ten times larger than the recoil energy, E_r = 0.1. Note that one can alter the length of the plaquette along the x direction (and thus the flux density) by changing β representing the potential gradient along the x axis.

First, let us consider the case where v₀ = 0.25 corresponding to the most flat adiabatic potential. In this situation the Chern numbers of the first five bands appear to be equal to the unity, as one can see in the left part of figure 5. Thus the Hall current should monotonically increase when filling these bands. This resembles the QHE involving the Landau levels. Second, we check what happens when we leave the regime v₀ = 0.25 where the adiabatic potential is flat, and consider lower and higher values of the coupling strength v₀. Near v₀ = 0.175 we find a topological phase transition where the lowest two energy bands touch and their Chern numbers change to c₁ = 0 and c₂ = 2, while the Chern numbers of the higher bands remain unchanged, illustrated in figure 6. In a vicinity of v₀ = 0.3 there is another phase transition, where the second and third bands touch, leading to a new distribution of Chern numbers: c₁ = 1, c₂ = − 1, c₃ = 3, c₄ = 1. Interestingly the Chern numbers of the second and the third bands jump by two units during such a transition.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, we explore the robustness of the topological bands. The right part of figure 5 shows the dependence of the band gap Δ₁₂ between the first and second bands on the coupling strength v₀ and the potential gradient β. One can see that the band gap is maximum for v₀ = 0.25 when the adiabatic potential is the most flat. The gap increases by increasing the gradient β, simultaneously extending the range of the v₀ values where the band gap is nonzero. Therefore to observe the topological bands, one needs to take a proper value of the Raman coupling v₀ ≈ 0.25 and a sufficiently large gradient β, such as β = 0.6.

We now make some numerical estimates to confirm that this scheme is reasonable. We consider an ensemble of ${}^{87}\mathrm{Rb}$ atoms, with $| \uparrow \rangle =| f\,=\,2,{m}_{F}=2 \rangle$ and $| \downarrow \rangle =| f\,=\,1,{m}_{F}=1 \rangle ;$ the relative magnetic moment of these hyperfine states is $\approx 2.1\ \mathrm{MHz}\,{{\rm{G}}}^{-1}$, where 1 G = 10⁻⁴ T. For a reasonable magnetic field gradient of 300 G/cm, this leads to the ${{\rm{\Delta }}}^{{\prime} }/{\hslash }\approx 2\pi \times 600\ \mathrm{MHz}\,{\mathrm{cm}}^{-1}=2\pi \times 60\ \mathrm{kHz}\,{\mu {\rm{m}}}^{-1}$ detuning gradient. For ${}^{87}\mathrm{Rb}$ with λ = 790 nm laser fields the recoil frequency is ${\omega }_{r}/2\pi =3.5\ \mathrm{kHz}$. Along with the driving frequency ω = 10ω_r, this provides the dimensionless energy gradient $\beta ={{\rm{\Delta }}}^{{\prime} }/({\hslash }\omega {k}_{0})\approx 1.3$, allowing easy access to the topological bands displayed in figure 5.

Adiabatic loading into this lattice can be achieved by extending the techniques already applied to loading into Raman dressed states [29]. The loading technique begins with a Bose–Einstein condensate (BEC) in the lower energy $| \downarrow \rangle$ state in a uniform magnetic field B₀. Subsequently one slowly ramps on a single off resonance RF coupling field and the adiabatically ramp the RF field to resonance (at frequency δω). This RF dressed state can be transformed into a resonant Raman dressed by ramping on the Raman lasers (with only the ω₀ + δω frequency on the k⁻ laser beam) while ramping off the RF field. The loading procedure then continues by slowly ramping on the remaining frequency components on the k⁻ beam, and finally by ramping on the magnetic field gradient (essentially according in the lattice sites from infinity). This procedure leaves the BEC in the q = 0 crystal momentum state in a single Floquet band.