Floquet engineering of optical lattices with spatial features and periodicity below the diffraction limit

Floquet engineering or coherent time-periodic driving of quantum systems has been successfully used to synthesize Hamiltonians with novel properties. In ultracold atomic systems, this has led to experimental realizations of artificial gauge fields, topological bandstructures, and observation of dynamical localization, to name a few. Here we present a Floquet-based framework to stroboscopically engineer Hamiltonians with spatial features and periodicity below the diffraction limit of light used to create them by time-averaging over various configurations of a 1D optical Kronig–Penney (KP) lattice. The KP potential is a lattice of narrow subwavelength barriers spaced by half the optical wavelength (λ/2) and arises from the nonlinear optical response of the atomic dark state. Stroboscopic control over the strength and position of this lattice requires time-dependent adiabatic manipulation of the dark-state spin composition. We investigate adiabaticity requirements and shape our time-dependent light fields to respect the requirements. We apply this framework to show that a λ/4-spaced lattice can be synthesized using realistic experimental parameters as an example, discuss mechanisms that limit lifetimes in these lattices, explore candidate systems and their limitations, and treat adiabatic loading into the ground band of these lattices.


Introduction
Time-periodic driving of quantum systems is ubiquitous in quantum mechanics. Small amplitude driving of a quantum system probes its linear response [1], while strong driving allows for Hamiltonian engineering [2][3][4][5][6]. Optical potentials and in particular optical lattices have proven to be a powerful tool for manipulating ultracold atomic systems and are used in a wide range of experiments [7][8][9]. However, the spatial features and periodicity of these potentials (generally arising from the second order ac-Stark shift) in the far field are constrained by the diffraction limit to be of order the wavelength of light used to create them. In particular, the Fourier decomposition of these far-field optical potentials cannot have components with wavelength less than λ/2, and thus the minimum lattice spacing is λ/2. As the lattice spacing determines many of the energy scales in cold-atom lattice systems, it has been of interest to produce optical lattices with smaller spacings in order to increase relevant energy scales [10,11]. Approaches to making subwavelength-spaced optical lattices have been proposed [12] and realized [13,14] based on multiphoton effects, and on adiabatic dressing of different spin dependent lattices [15,16].
Recently, optical lattices based on the nonlinear optical response of dark states [17,18] were realized [19] with λ/2 periodicity but strongly subwavelength structure within a unit cell, consisting of a Kronig-Penney-like (KP) lattice of narrow repulsive barriers of width ;λ/50. Time averaging a stroboscopically applied lattice potential with high spatial frequency Fourier components can give rise to an average potential with periodicity and spatial features smaller than λ/2 [10]. Since the dark-state KP lattice has high spatial frequency Fourier components, it is a candidate progenitor lattice with which to realize such a time-averaged, subwavelength-featured lattice. Here, we explore the implementation of a time-averaged dark-state KP lattice, taking into account realistic imperfections in the dark-state system. After careful consideration of the adiabaticity requirements, we show that lattices with λ/4 period can be realized as an example, and discuss the prospects for lattices with smaller spacings and features. Reference [20] explores related ideas about painting arbitrary subwavelength optical potentials.
In the time-averaged approach, a time-periodic progenitor potential W 0 (x, t) is applied such that the atoms experience the time-averaged potential W avg (x): where T=2π/ω T is the period of W 0 (x, t) and ω T is the Floquet frequency. In order to successfully realize W avg (x) while avoiding heating, ω T must be much faster than the timescale associated with the motional degree of freedom in the lattice, which is set by the energy gaps between bands in the lattice [10,21]. This requirement suggests that ω T be as large as possible. As we discuss below, the particular realization of W 0 (x, t) using a darkstate lattice [19] has an additional requirement of spin adiabaticity that limits the maximum allowable ω T . The dark-state lattice is an artificial scalar gauge potential [17-19, 22, 23] experienced by an atom in the dark-state eigenfunction of a three-level Λ-system with a spatially dependent spin composition. Dynamically manipulating the height, barrier width, and position of the lattice requires time-dependent manipulation of the spin composition of the dark-state eigenfunction. This spin manipulation can be seen as a stimulated Raman adiabatic passage (STIRAP) process [24] and adiabaticity requirements set an upper bound on the window for usable ω T within which the atoms are simultaneously motionally diabatic and spin adiabatic. Understanding the practical limits of these constraints requires a detailed consideration of the system dynamics, which we apply to the specific 171 Yb system previously used to demonstrate the dark-state lattice [19].

Time-dependent dark-state potentials
We consider the creation of time-periodic potentials for the dark-state channel, W DS (x, t) (which serves as W 0 (x, t) in equation (1)), by coupling the three atomic levels in a Λ-system with a spatially homogeneous probe light field Ω p (t), and a spatially inhomogeneous control light field. The inhomogeneous control light field is composed of two counter propagating fields with equal magnitudes driven simultaneously, Ω c (x, t)=Ω c (t) cos(kx+f(t)) where k=2π/λ, as shown in figure 1(a). Working in the spatially and temporally local dressed state basis of the Λ-system determined by the coupling fields Ω p (t) and Ω c (x, t), the Hamiltonian is given by (equation (A.5)) where W DS (x, t) and W ± (x, t) are the dark-state and bright-state potentials in the three Born-Oppenheimer (BO) channels and H x p t , , od ( )represents the off-diagonal couplings between these channels (see appendix A). W DS (x, t) and W ± (x, t) include the BO potentials as well as the non-adiabatic corrections to these potentials. The dressed state coupling induced by H x p t , , od ( )is detrimental, since it mixes bare excited state 3ñ | into the darkstate channel through the bright-state channels, inducing photon scattering in the otherwise lossless dark state.
The spin composition of the dark-state eigenfunction for the Λ-system in figure 1 The non-adiabatic correction to the dark-state BO potential that gives rise to W DS (x, t) is determined by the spatial gradient of the spin composition [17,18] which for the light-field configuration considered here is a lattice of narrow repulsive barriers with temporally modulated strength and position. We take here a stroboscopic approach, where W DS (x, t) is repeatedly pulsed on and off in magnitude at N different positions for time T i with the position of W DS (x, t) being shifted in between the lattice pulses (here T T i = å ). In addition, W DS (x, t) can be held on and off for t on,i and t off,i (figure 1(d)). Time averaging over the N different pulsed KP lattice potentials with arbitrary strength and position can produce an arbitrary time-averaged potential W avg (x) [20].
The ability to paint potentials requires real-time control over the position, strength and width of the barriers (equation (3)). The strength of the barriers can be controlled via the Rabi frequencies Ω p (t) and Ω c (t) (figures 1(b), (c)) with the height and width of the barriers being proportional to 1/ò 2 (t) and ò(t) respectively [17][18][19] where ò(t)=Ω p (t)/Ω c (t) (for ò(t)=1). The barriers are located at the nodes/minimums of Ω c (x, t) (figure 1(c)), and their positions can be controlled by the control beam phases f 1 (t)andf 2 (t) (figures 1(b), (c)). Stitching N different sub-Floquet periods together (while ensuring continuity in the Rabi pulses between the sub-Floquet periods) into one Floquet period allows for the versatility in the time-averaged potential W avg (x) that can be generated. Each sub-Floquet period of duration T i pulses a KP potential at a different position x 0i (determined by the phase f 0i ) with a strength and width determined by ò i . Figure 1(d) shows the pulses Ω p(i) (t), Ω c(i) (t) and f (i) (t) for the ith sub-Floquet period −T i /2tT i /2.

Adiabaticity considerations
Without explicit time-dependence, H x p t , , od ( )has static, off-diagonal terms depending on the spatial gradient of the dark-state spin composition that couple the dark-state channel to the lossy bright-state channels. This loss mechanism was theoretically [17,18] and experimentally [19] shown to limit lifetimes in the KP lattice. Large energy gaps between the BO channels via large Rabi frequencies Ω c and Ω p generally aid in suppressing this loss [19]. When explicit time dependence to the Rabi frequencies is included i.e. Ω c (x, t) and Ω p (t), the spatially dependent loss mechanism has a trivial time dependence due to the temporally periodic nature of the changing dressed states, and the loss is quantified by averaging over one Floquet period. There is, however, an additional loss mechanism mediated via an explicitly time-dependent term in H x p t , , od ( )(see appendix A). This term mediates non-adiabatic couplings between the dark-state channel and the bright-state channels, but can be suppressed by careful pulse shaping.
Our goal is to design Ω p (t) and Ω c (x, t) to be simultaneously motionally diabatic and spin adiabatic. In order to design pulses that are spin adiabatic, we consider the three inequalities that quantify the sufficiency requirements for adiabaticity [25] defined at single photon resonance Δ=0 (see appendices A.2, B.2): , , 4 An ideal Λ-system with inverse lifetime Γ and single-photon detuning Δ. One leg of the Λ-system is coupled by a spatially homogeneous and temporally varying probe light field Ω p (t) and the other leg by a spatially inhomogeneous and temporally varying control light field Ω c (x, t). (b) The geometry of the light fields with arbitrary control over the envelope, Ω c1 (t), Ω c2 (t), Ω p (t) and phase, , , called the local adiabatic criterion [24], states that to ensure adiabaticity during pulsing, the energy gap between the dark and bright eigenstates (set by Ω rms (x, t)) must be much greater than the off-diagonal couplings between them ( forces the pulses to be smooth while both equations (4b) and (4c) set bounds on their rise time and fall times.
To design pulse shapes that satisfy equations (4a)-(4c), we parameterize the condition equation 4(a) through a parameter r(t): where T i =t Si +t off,i +t on,i and t Si /2 is the rise or fall time. Generally, it is easiest to be spin adiabatic for large energy separation between the dark and bright-state channels. However at the nodes of Ω c (x, t), this energy gap is the smallest with a value of ÿΩ p (t)/2 for Δ=0. Therefore, we consider pulse schemes that change the positions of the nodes only when the energy gap at the nodes is large and the spin composition is essentially homogeneous (ò(t)?1). We consider two ways to achieve the homogeneous condition in between pulses: (1) Ω p (t)?Ω c (x, t) achieved by turning up Ω p while turning off both control beams, (2) Ω c1 ?Ω p ?Ω c2 (t) achieved by turning off Ω c2 (t) while Ω p and Ω c1 are kept constant.
We note that the pulsing schemes considered here are not unique. Control over Ω c1 (t), Ω c2 (t), Ω p (t), f 1 (t), and f 2 (t) allows for multiple ways by which arbitrary potentials can be painted, and we refer the reader to [20] for other variants.
For pulse scheme (1), the position where the local adiabatic criterion is the hardest to satisfy, x h , occurs between the nodes. (Pulse scheme (2), for which only one of the two control beams is driven, is treated in the appendix). We choose is provided by the hold times t off,i and t on,i . the rms average of the Rabi frequencies to be constant at Solving equations (5) and (6) simultaneously, the expressions for t Si , Ω c(i) (t) are as follows (figure 2): As a specific example, we explore creation of λ/(2N)-spaced lattices where N=2, 3, 4 K. These lattices are created by time-averaging N λ/2-spaced progenitor KP lattice potentials, each shifted in position by iλ/(2N) for (i−1)T/NtiT/N and pulsed for a period of T i =T/N [10], where i=0, 1, K, N−1 (figure 3(a)). More flexibility is possible by pulsing the progenitor lattice with different strengths and relative positions, realizing for example the Rice-Mele model [26,27] as shown in figure 3 The goal to create λ /(2N)-spaced lattices that significantly confines the ground band sets constraints on the lattice parameters. Without requirements of spin adiabaticity, time-averaging the KP potential creates λ/(2N)spaced lattices with barriers of maximum average height of (1/N)E R /ò 2 . Due to the reduction in the size of the unit cell by N, the characteristic energy increases to N 2 E R , which is also approximately the energy of the lowest band in a KP lattice. Hence for the λ/(2N)-spaced lattice to provide significant confinement, which can be controlled by choosing ò (limited by requirements on spin adiabaticity) and t off,i and t on,i . Of course, non-zero values for t off,i and t on,i decrease the Floquet frequency (ω T ) as on, Reducing ω T makes it more difficult to be fully motionally diabatic, so that the operational window between the two constraints rapidly decreases with increasing N.

Solving for the Bloch-Floquet bandstructure
We solve the Bloch-Floquet bandstructure for our Hamiltonian is the quasienergy operator defined in an extended Hilbert space where time is treated as a coordinate with periodic boundary conditions [2,3,5]. The extension of the Hilbert space is symbolically represented by the double ket notation for the Bloch-Floquet mode u x t , We solve the eigenvalue problem in equation (12) to calculate the Bloch-Floquet bandstructure (appendix C.1).

Results
The loss due to the off-diagonal coupling terms in H x t , od ( )that arises from the spatial gradient of the dark-state spin composition increases with smaller ò [17][18][19]. This suggests that in order to generate potentials that have reasonable lifetimes with realistic values for Rabi frequencies, it is desirable to work at large ò, as allowed by equation (9). In figure 4, we investigate the creation of a λ/4-spaced lattice potential that significantly confines the ground band. With a choice of ò i =0.2, r 0i =0.01, ci 0 W =600Γ, and t Si =0.36T i (equation (10)) we use the pulse shape in equation (7) to create this potential ( figure 4(a)). The green trace is the pulse profile for Ω c (t) (and therefore t t while the blue trace is the phase profile f(t) during one Floquet period. The numbered red dots enumerate the different time slices during the pulse. One Floquet period of pulsing involves stitching together two sub-Floquet periods that have a relative phase f (i) (t) differing by π/2. Note that the phase is suddenly switched during an off period when there is no spatial variation to the dark state. The sub-Floquet periods are color coded and labeled as I and II. This pulse yields an effective λ/4-spaced W avg (x) potential with ∼8E R barriers as shown in figure 4 In figure 4(c), we search for the window of operational ω T within which the bands of an effective λ/4-spaced lattice Hamiltonian are clearly defined by monitoring the Floquet spectrum E q at q=0 as a function of ω T /ω R . As ω T is increased, the pulsing becomes more motionally diabatic, but at the cost of increased r 0i (equation (8)). The increased r 0i results in stronger admixing of the dark-state channel with the bright-state channels. The loss rate −Im(E q ) is encoded in the color of the points in figure 4(c). The gray dots have loss (−Im(E q )) much larger than the highest value in the color bar.
In ignore H x t , od ( )and therefore exclude loss due to non-adiabatic couplings with the bright states. In figure 4(f), we show the Bloch-Floquet bandstructure of H x t , ( )(equation (11)), which includes the non-adiabatic bright-state couplings. The avoided crossings exist in the Bloch-Floquet bandstructure at the same place (q, E q ) for the same parameters in figures 4(e) and (f), suggesting that these crossings arise from couplings with high-lying dark-state eigenfunctions. The ground Bloch-Floquet band for H x t , ( )has the same shape as the static λ/4-spaced lattice (except near the avoided crossings). The calculated average lifetime in the time-averaged potential for the ground band in figure 4(f) is 1 ms, which can be substantially improved with a lower ω T . In general, lifetimes can be increased and the avoided crossings can be removed by operating at larger Rabi frequencies.
We also calculate the Bloch-Floquet bandstructures at a lower Floquet frequency of ω T =50ω R (pink vertical line in figure 4(c)). The Bloch-Floquet bandstructure of the dark-state channel Hamiltonian H x t , DŜ ( )is shown in figure 4(g), and of H x t , ( )in figure 4(h). The average lifetime of the ground band is 32 ms (the colored regions in figure 4(h)) and is much longer than the lifetime of 1 ms for the Bloch-Floquet bandstructure obtained at ω T =150ω R in figure 4(f). The avoided crossings due to coupling to high-lying dark states, however, are larger for the same Rabi frequencies.
In figure 5, we show the dynamics of the spatial probability densities of the dark-state Bloch-Floquet mode and its spin composition at time slices 1-5 ( figure 4(a)) in the (q, E q ) configuration labeled by the brown star in figure 4(f). The purple trace is the scaled probability density of the Bloch-Floquet dark-state mode and it has λ/4 periodicity. It is roughly stationary except for the small wiggles that correspond to micromotion. Meanwhile, the spin composition (the black and red traces) of the Bloch-Floquet mode changes dramatically as a function of time during the Floquet period. The population in 3ñ | (yellow trace) remains quite small. To create subwavelength-spaced lattices with larger N, the window of operational ω T will be smaller by N (equation (10)). In addition, since larger N requires working with smaller ò (equation (9)), the loss mechanisms limiting lifetime in these lattices become more significant. However, higher Rabi frequencies can combat both these limitations. Using larger Rabi frequencies, we show an example demonstrating the feasibility of a λ/6spaced lattice.
With increased Rabi frequency of ci 0 W =1800Γ, we use the pulse shape in equation (7) to design a composite pulsing profile (figure 6(a)) to create a time-averaged λ/6-spaced lattice potential W avg (x) that has ∼9E R tall barriers ( figure 6(b)). Here ò i =0.135, r 0i =0.01, and t Si =0.59T i . One Floquet period of pulsing involves stitching together three sub-Floquet periods with the relative phase f (i) (t) differing by π/3 between adjacent sub-Floquet periods. For ω T =280ω R (yellow vertical line in figure 6(c)), we plot the bandstructure of the timeaveraged Hamiltonian p m W x 2 2 avg + ( ) ( )in figure 6(d), the Bloch-Floquet bandstructure of the dark-state The folded bandstructures are indicative of a λ/6-spaced lattice. The particular phase sequence in figure 6(a) used to create the λ/6-spaced lattice breaks time-reversal symmetry: the Bloch-Floquet bandstructures are therefore asymmetric about q=0. The avoided crossings enclosed in the red circles exist in both Bloch-Floquet bandstructures at the same place (q, E q ) for the same parameters in figures 6(e) and (f), suggesting that these particular crossings arise from couplings with high-lying dark-state eigenfunctions. The avoided crossing enclosed in the green circle in figure 6(f) arises due to couplings with bright-state eigenfunctions.

Experimental considerations and limitations
While working with large Rabi frequencies reduces losses, a potential disadvantage is that the Λ-system approximation may break down. Perfect Λ-systems are rare in nature, and Ω p (t) and Ω c (t) can couple offresonantly to states outside of the Λ-system. These off-resonant couplings manifest as effective two-photon detunings for the approximate Λ-system. Non-zero two-photon detunings are detrimental to STIRAP [24,32], although spatially homogeneous detuning could in principle be compensated with time-dependent laser detuning. Two-photon detunings originating from Ω c (x, t), however, are temporally and spatially modulated  and may not be completely compensated without the significant experimental overhead of adding more spatially dependent compensating laser fields. In addition to added two-photon detuning, the lifetime in the timeaveraged lattices is further limited due to admixing of excited states outside the Λ-system. Hence, there are tradeoffs when increasing the magnitude of the Rabi frequencies: while the dark-state evolution is more adiabatic with less bright-state admixture, the off-resonant scattering from states outside the Λ-system also increases. Appendix D presents calculations for a realistic system consisting of 171 Yb atoms, which was used to create KP lattices [19].
A number of techniques can be used to verify the creation of these subwavelength lattices. For example, nanoresolution microscopy [33] can be used to directly map out the ensemble-averaged probability density of atoms in the ground band of the λ/(2N)-spaced lattices. In addition, Bloch oscillations [20,34] or time-of-flight measurements of the momentum distributions [35] could be used to measure the N times larger Brillouin zones. In fact, λ/4-spaced lattices were recently realized using the techniques explored in this paper with 171 Yb atoms [36].
To adiabatically load into the ground band of the time-averaged λ/(2N)-spaced lattice potential, the stroboscopic lattice should be turned on slower than the motional timescale set by

Summary and outlook
In this paper, we evaluate the idea of stroboscopically generating potentials using the repulsive barriers of a darkstate KP potential. We analyzed the competing requirements of maintaining dark-state spin adiabaticity and simultaneous motional diabaticity during pulsing of the KP potentials in the presence of realistic imperfections. We showed that it is possible to create such potentials in an experimental system of 171 Yb atoms by calculating the Floquet spectrum of atoms in a stroboscopically generated λ/4-spaced lattice. This approach is applicable to any three-level system, although it needs to be well isolated from coupling to other levels in order to ensure good lifetimes. While we have treated 1D systems here, this method can be readily generalized to 2D. Using progenitor lattices of subwavelength attractive traps [37] in conjunction with barriers can allow for flexibility in tailoring arbitrary time-averaged potential landscapes not limited by diffraction.

Acknowledgments
We would like to thank Jean Dalibard

Appendix A. Sufficiency conditions for adiabaticity
We start with the time-dependent Hamiltonian H x t , ( )is non-Hermitian due to the iΓ/2 term and requires a biorthogonal set of eigenvectors to diagonalize it [38]. Due to the non-Hermitian nature of H x t , ( ), the eigenvectors are not guaranteed to be orthogonal to each other, but still form a linearly independent set that spans the Hilbert space [38]. To derive the artificial gauge potentials and for quantifying the sufficiency conditions for adiabaticity, we transform H x t , ( )using a rotation transform R x t , ( )composed of the right eigenvectors [38] of the spin-light field coupling part of H x t , ( ). The expression for R x t , ( )is:   (1) where both control beams are changed simultaneously with equal magnitude: Ω c (x, t)=Ω c (t)sin(kx+f(t)) resulting in θ(x, t)=0. For pulse scheme (2) where only one control beam is pulsed: ˆˆˆˆˆrotates H x t , ( )into the dressed-atom picture of the Λ-system. The effective Hamiltonian after the transformation is [22,24,32,[39][40][41]] . The term B x t , ( )arises from the explicit time dependence of H x t , ( ). Careful pulse shaping can help suppress the terms in B x t , ( )that couple the dark-state channel with the brightstate channels. The coupling terms in H x p t , , od ( )depend only on the ratio of the Rabi frequencies ( , a a ¢  ) and not on their absolute magnitudes, while the energy separation between the channels (E x t , BÔ ( )) depend on absolute magnitudes of the Rabi frequencies. Thus at higher Rabi frequencies the BO channels become increasingly decoupled. In addition, , ensures that the bright-state channels are well separated from the dark-state channel.

A.1. Floquet scalar gauge potentials
Here we derive the expression for the Floquet scalar gauge potential for the dark-state channel. The expression for Â when both control beams are driven simultaneously as in pulse scheme (1) is The expression for the scalar gauge potentials for the BO channels is [23,40]: where the first matrix contains the scalar gauge potentials for each of the BO channels. The scalar gauge potential for the dark-state channel is The expressions for A x t , ( )and W DS (x, t) for pulse scheme (2) (θ(x, t)¹constant), are and

A.2. Formulating the sufficiency conditions for adiabaticity
The general expression for B is analogous to Â, except that the derivatives are with respect to x in Â and with respect to t in B, which for pulse scheme (1) is: The local adiabatic criterion (equation (4a)) for the instantaneous dark state requires that the minimum of the spatially and temporally varying energy gap between the dark and bright states must be much larger than the largest off-diagonal couplings between them, which we quantify as For Δ=0, where fastest STIRAP pulses are guaranteed [24,32,42,43] x t t x t , , , , which in addition to a  and α depends on q  and θ.
For this scheme, x h is at the nodes of Ω c (x, t) since the energy gap between the dark-state and bright-state channels is the smallest at the nodes and the spin at the node must completely flip from 2ñ | to 1ñ | ( figure 1(a)) at the end of the pulse.

B.2. Verifying the spin-adiabaticity requirements and choice for r 0i
We use the off-diagonal coupling terms in equation (A.11) and set Δ=0 to recast the sufficiency conditions in [25] (the inequalities equations (4a)-(4c)). The first condition equation (4a) implies where we have used equations (5) and (6). For r 0i =0.01, this inequality is well satisfied. The stronger version [25] of the second inequality equation (4b) states: For r 0i =0.01, this inequality is also well satisfied. We note that equation (B.5) also enforces that r(t) must be differentiable. The stronger version of the third inequality equation (4c) is [25] x t r t t , where we have substituted equation (8). Again, this inequality is well satisfied for r 0i =0.01.
Appendix C. Bloch-Floquet bandstructure C.1. Bandstructure of H x t , ( ) We evaluate the matrix elements of the quasienergy operator K q derived in section 4. K q is expressed in dimensionless units x and t where x , and the tildes over x and t are dropped for convenience, as follows: We expand the Hilbert space of the Bloch-Floquet modes in a plane wave basis: . We solve equation (C.2) by diagonalizing K q . The matrix elements of the spin-independent components of K q l m j K lmj It is important [17] to appropriately choose the number of plane waves L and M to be large enough to accurately represent the couplings between the dark-state channel and bright-state channels. We solve for the lowest few dozen eigenstates near zero energy of these sparse matrices with dimensions 3(2L+1)(2M+1)×3(2L+1)(2M+1)∼ 10 5 ×10 5 using the Arnoldi algorithm. We find that the solution converges with M as low as 25, however for all our calculations we use M;210.
C.2. Bandstructure of H x t , DŜ ( ) In this subsection, we outline the method used to numerically solve for the Bloch-Floquet bandstructure of the dark-state channel ignoring non-adiabatic couplings to the bright-state channels: