A dark state of Chern bands: Designing flat bands with higher Chern number

We introduce a scheme by which flat bands with higher Chern number $\vert C\vert>1$ can be designed in ultracold gases through a coherent manipulation of Bloch bands. Inspired by quantum-optics methods, our approach consists in creating a"dark Bloch band"by coupling a set of source bands through resonant processes. Considering a $\Lambda$ system of three bands, the Chern number of the dark band is found to follow a simple sum rule in terms of the Chern numbers of the source bands: $C_D\!=\!C_1+C_2-C_3$. Altogether, our dark-state scheme realizes a nearly flat Bloch band with predictable and tunable Chern number $C_D$. We illustrate our method based on a $\Lambda$ system, formed of the bands of the Harper-Hofstadter model, which leads to a nearly flat Chern band with $C_D\!=\!2$. We explore a realistic sequence to load atoms into the dark Chern band, as well as a probing scheme based on Hall drift measurements. Dark Chern bands offer a practical platform where exotic fractional quantum Hall states could be realized in ultracold gases.

: (a): The Λ scheme involving three discrete levels. The Rabi couplings Ω 1,2 and frequency detunings ∆ and δ are indicated; see Eq. (2). A stable dark state |D is created when δ ≈ 0; see Eq. (3). (b) The Λ scheme involving three Bloch bands E 1,2,3 (k). The detunings and Rabi frequencies depend on the quasimomentum k of the Bloch states. The Chern number of the resulting "dark band" C D follows a simple sum rule in terms of the Chern numbers C 1,2,3 of the bare bands.

Λ system of Bloch bands and the sum rule
For the sake of concreteness, we consider the Bloch bands of the emblematic Harper-Hofstadter (HH) model [61], which can be engineered in optical-lattice experiments [7][8][9][10]. The corresponding Hamiltonian readŝ H HH = −J n=(n,m) a † nâ n+1y +â † nâ n+1x e i2πmφ + h.c., (5) where n = (n, m) enumerates lattice sites on a 2D square lattice with lattice constant a, J is the hopping amplitude, and where a n /a † n annihilate/create a particle at site n. The complex phase factor, which accompanies hopping along the x direction and generates a uniform magnetic flux 2πφ per plaquette, can be tuned in experiments [7][8][9][10].
Considering a flux of the form φ = 1/q, with q an odd integer, and applying periodic boundary conditions (PBC), the spectrum displays q non-degenerate Bloch bands ε ν (k), where ν = 1, . . . , q. We define the first Brillouin zone (BZ) as k ∈ [0, 2π/aq) × [0, 2π/a). Each band is associated with a topological Chern number [62] where ψ ν (k) denotes an eigenstate in the νth band. For q a generic odd integer, the Chern number of the central band [ν = (q + 1)/2] reads C ν = (−q + 1), while C ν = 1 for all the other bands. Hence, except from the central band, the bands of the HH model are reminiscent of Landau levels (LL). In particular, these LL-like bands become flat in the limit q 1. We now discuss how a coherent coupling of these LL-like bands allow for the generation of flat bands with higher Chern number |C| > 1 within the HH model. This is achieved by building a Λ system [Eq. (2)] from three selected bands E 1,2,3 (k) of the HH spectrum {ε ν }; we will denote the mean energy of each band asĒ 1,2,3 . In order to couple these bands quasiresonantly, we consider a two-frequency drive with frequencies set tohω 1,2 =Ē 1,2 −Ē 3 ; see Fig. 1(b). Assuming that the coupling field has a spatial periodicity that is compatible with the (magnetic) unit cell, the system forms a collection of decoupled Λ systems,Ĥ Λ (k), one at each quasi-momentum k, and involving the three relevant states where |ψ n (k) denotes the Bloch state in band n and quasimomentum k; see Fig. 1(b). Applying the rotating wave approximation (RWA) [37], as detailed in Appendix B, and neglecting the other bands, this setting is indeed well described by Eq. (2), with momentum-dependent Rabi frequencies and detunings δ(k), ∆(k). In this setting, one can define a dark state |D(k) and two bright states |B ± (k) at every quasi-momentum k; the corresponding bands will be referred to as the "dark band" and the two "bright bands". As a central result of this work, we find that the Chern number of the dark band is given by a simple sum rule: where C 1,2,3 denote the Chern numbers of the selected bands E 1,2,3 . This simple rule is independent of: (i) the flux φ = 1/q, (ii) the operator associated with the coupling field, and (iii) the chosen bands in the HH spectrum. Altogether, this provides a unique way to produce bands with predictable higher Chern number. We now provide the proof of Eq. (9) and discuss its regime of applicability. The formula (9) is valid whenever an effective Λ (three-band) configuration is achieved and the dark-state band is well separated from the two bright-state bands. To prove it, let us first remark that it is sufficient to show that the Chern number of the well-isolated bright bands is equal to C 3 . Indeed, the couplings Ω 1,2 only mix the three bands E 1,2,3 (k), so that this implies the sum rule in Eq. (9). The hypothesis in Eq. (10) is demonstrated via the existence of a smooth deformation |h(k, t) that connects the bright states |B ± = |h(k, 1) to the bare state |ψ 3 (k) = |h(k, 0) . This deformation is provided by the expression where N (k, t) is a normalization factor, and where we used the definition of the Rabi frequencies Ω 1,2 in Eq. (2). The deformation |h(k, t) in Eq. (11) (12)). The Chern number of a band is preserved under such a homotopy, hence leading to the announced result in Eq. (10); this closes the proof of the sum rule in Eq. (9). Importantly, the assumptions underlying the sum rule (9) impose a series of constraints on the chosen Bloch bands and system parameters, as we now explain.
First, the finite bandwidth of the bare bands E 1,2,3 (k) produces detunings δ(k) and ∆(k) in Eq. (2), which affects the flatness of the dark band and reduces the gap to the bright bands; see Appendix C.1. This effect can be limited by noting that the width of the HH bands decreases exponentially with q (except for the central band); see Appendix C.2.
Then, the Rabi frequencies Ω 1,2 (k) should satisfy the inequalities δ(k) Ω 1,2 (k) |E ± (k)| to allow for a good separation of the dark band and optimize its flatness; this condition can already be reached for moderate q ∼ 7 − 10. Besides, one should avoid pathological zeros Ω(k) = 0, which would invalidate the dark state construction; see Appendix C.4.
Furthermore, for large q 1, the HH spectrum forms a ladder of quasi-equally spaced Landau levels, which indicates that only special choices of bare bands E 1,2,3 can lead to a genuine Λ configuration. In addition, the RWA is only valid when the resonant frequencies ω 1,2 (and their differences ∆ω) are much larger than all other frequency scales, which also sets an important constraint on the chosen bands.
Finally, the unique bare band with C ν = 1 is the central band [ν = (q + 1)/2], for which the bandwidth to bandgap ratio is O (1). This unfortunately rules out involving this band in our construction, which eventually implies the disappointing result C D = 1 [Eq. (9)].
In order to overcome these limitations and constraints, we slightly generalize our scheme by introducing different atomic species (i.e. "spins"), as we describe in the next section.
Before doing so, we note that similar sum rules were numerically investigated in the multilayer configuration of Ref. [63]. The latter work considered the stacking of trivial (boronnitride-type) and non-trivial (Haldane-type) 2D lattices, and it explored the total Chern number of the system at half-filling as a function of the inter-layer coupling strength. In that context, the calculated Chern number reflects the topology of a three-fold degenerate band, which is associated with the three underlying layers; using an open-system approach, this total Chern number was decomposed as a sum of sub-system indices [64]. We point out that the approach developed in the present work is radically different, as it involves the Chern number of a non-degenerate and well-isolated dark band; in particular, the sum rule in Eq. (9) directly relates the dark band's Chern number to the Chern numbers of the individual bare bands, without requiring the use of an open-system formalism.

The multi-species configuration
We now describe the multi-species configuration, which allows to produce dark bands with higher Chern number. Specifically, we propose to generate a flat and non-degenerate dark band with C D = 2, by constructing a Λ system made of two HH bands [C 1,2 = 1] and a trivial band [C 3 = 0]. The latter is provided by the lowest band of a square lattice without flux [Eq. (5) with φ = 0, denotedĤ 0 ], whose unit cell area A cell = qa × a ensures a common BZ withĤ HH ; the corresponding sites are located atñ = (qn, m). This setting could be realized using different atomic (internal) states trapped in state-dependent potentials [66][67][68]; see Refs. [7,69,70] for schemes realizing state-dependent synthetic flux. In order to guarantee the validity of the RWA, we henceforth consider that each of the three selected bands E 1,2,3 is populated by a specific internal state σ = {1, 2, 3}. The Λ coupling between the three bands is then performed by properly coupling the internal states with microwave fields. For convenience, we take the trivial band (E 3 ) to be extremely flat (i.e. the corresponding statedependent lattice is assumed to be very deep, with hopping amplitudeJ J). We note that using only two internal states could also be envisaged in practice. Hofstadter bands associated with two uncoupled atomic internal states (σ = 1, 2), for a flux φ = 1/11. The black dispersions correspond to edge-state branches [65]; the system is diagonalized on a cylindrical geometry. (b) Adjusting the detunings δ σ inĤ Λ such that two Hofstadter bands (populated by states σ = 1, 2) become degenerate with a flat trivial band (populated by σ = 3) in the uncoupled limit A 1,2 = 0. (c) Activating the coupling (A 1,2 = 0, max |Ω 1,2 (k)| = J/5) splits the three overlapping bands into a flat dark band at zero energy (red shaded), and two bright bands (green/purple shaded). The flatness of the dark band is emphasized by a logarithmic scale. The edge-state branches, associated with the two edges (I,II) of the cylinder, indicate that the Chern number of the dark band is C D = 2.
In a rotating frame, and upon applying the RWA, the Hamiltonian of this setting readŝ where a † n,σ creates an atom at siteñ in internal state σ;Ĥ 1,2 =Ĥ HH ⊗P 1,2 andĤ 3 =Ĥ 0 ⊗P 3 , whereP σ projects onto the σ component; A 1,2 denote the amplitudes of the coupling fields, and the detunings δ σ are controlled by tuning the driving frequencies out of resonance. In particular, these detunings can be chosen such that two bands of the HH spectrum (populated by states σ = 1, 2) become degenerate with the trivial band (populated by σ = 3) in the decoupled limit A 1,2 = 0; see Figs. 2(a)-(b). Upon activating the coupling (A 1,2 = 0), these three bands split into the dark band and the two bright bands shown in Fig. 2(c). While the bright bands acquire a small dispersion, due to the k-dependence of the Rabi frequencies Ω 1,2 (k) in Eq. (8), the dark band at zero energy remains almost perfectly flat; see Appendix C.4.
One confirms that the Chern number of the dark band is C D = 2, as dictated by the sum rule (9); this is readily obtained by analyzing the edge-state branches of opposite chirality that enter and leave the dark band in Fig. 2(c); see Ref. [65]. One also finds that the Chern number of the bright bands are zero [Eq. (10)], such that the total Chern number of the three coupled bands is indeed conserved. We have verified that these results are generic, in the sense that they do not depend on the specific form of the coupling operator.

Higher Chern number from centre-of-mass responses
The Chern number C D of the constructed dark band could be measured through different probes in ultracold atoms, such as center-of-mass responses [9], edge-state spectroscopy [71,72], and circular dichroism [20]. Here, we validate our approach by simulating the center-ofmass displacement of an atomic cloud, loaded in the dark band, and perturbed by a linear potential gradient [9,[38][39][40].
As a first step, we numerically simulate the following protocol [38]: We initially confine the system, using sharp rectangular walls; we generate a dark band [as illustrated in Fig. 2(c)] in this geometry and completely fill it with non-interacting fermions; we remove the confining walls, and act on the particles with a weak linear potential gradient aligned along the y direction,F = F a s=1,2,3 n=(n,m) mn n,s ; finally, we calculate the time evolution of the particle density ρ(x, y), and evaluate the center-of-mass displacement ∆x(t). The latter observable is related to the Chern number of the populated band (C D ) through the relation [9,[38][39][40] We show the simulated time-evolved density in Figs. 3(a)-(b), which demonstrate a clear transverse drift of the cloud (along the x direction). We note that the density modulation along the x direction reflects the coupling to the "trivial" lattice, of area A cell = qa × a, which supports the component σ = 3 (which is absent in the dark state). In Fig. 3(a),(b), the applied force is weak compared to the gap ∆ db separating the dark band from the bright bands, F a ∆ db , such that the particles' motion adiabatically follows the dark band. In this linear-response regime, Eq. (14) is applicable [38,39], and we extract an "experimental" value  (14)]. (c) Density at the end of the proposed loading sequence, using a total ramp duration T = 3 · 10 3h /J and system size 242 × 22. (d) Timeevolved density of the realistically prepared state, when exposed to the linear gradient along y for the same duration as in (a)-(b), t = 10(h/F a); here we set F = 5 · 10 −3 (J/a). (e) Center of mass drift for the manually-loaded dark band (blue line) as in (b), the bare Hofstadter band (red line) as in (a), and for the realistically-prepared state (black crosses) as in (c); in the latter case, residual bright-state contributions were filtered out (main text), hence leading to the Chern number measurement of C exp D ≈ 1.92.
for the Chern number of the dark band C exp D = 1.99 from the center-of-mass drift [ Fig. 3(b)]. In contrast, the force is stronger and comparable to ∆ db in the case depicted in Fig. 3(d), which leads to a dynamical repopulation of the bright states. Importantly, these undesired excitations are clearly identified by narrow stripes of highly populated sites, separated by qa along the x direction. For comparison, we also calculated the time-evolved density upon populating a bare Hofstadter band (with C = 1) instead of the dark band: the resulting Fig. 3(a), when compared with Fig. 3(b),(d), reveals that the center-of-mass velocity indeed differs by the predicted factor of two; see also Fig. 3(e). Finally, we explore a realistic sequence for loading atoms into the designed dark band of Fig. 2. We start the sequence by completely filling the lowest Hofstadter band with a single species (σ = 1), in the absence of coupling to the other species (σ = 2, 3). This initial state is a Chern insulator with Chern number C 1 = 1. Considering the Λ scheme of Eq. (13), one then progressively activates the couplings (A 1,2 = 0) so as to adiabatically transfer the atoms from the bare Hofstadter band to the dark band with C D = 2. Because of the change in the Chern number (which requires a gap-closing event), non-adiabatic effects are unavoidable, so that this final state cannot be reached with 100% fidelity. Moreover, the edge modes that cross the bulk gaps also provide a natural channel for non-adiabatic processes. For small systems, where the number of edge modes O(L x + L y ) is not marginal compared to the number of bulk modes O(L x L y ), we expect the dominant loss source to originate from the edge modes contribution.
An optimal population in the target dark band can be obtained by slowly ramping up the coupling A 1 to its final value [9,73]. We perform a simulation of this sequence, by calculating the time evolution of single-particle states according to the time-dependent Hamiltonian H Λ [A 1 (t)], where A 1 (t) = A 1 (t/T ) α . In order to reduce the initial energy-degeneracy of three eigenstates at t = 0, we add a small and time-dependent detuning δ(t) = −δ 0 (1 − t/T ) α+1 to δ 1 , which indeed allows us to increase the efficiency of the process [73]. The particle density in the final state is shown in Fig. 3(c), where the thin vertical stripes separated by qa again highlight the residual population in the bright states.
We illustrate the efficiency of the simulated loading scheme through Table 1. For instance, for a lattice of size 363 × 33, using a linear ramp of duration T ≈ 10 4 (h/J) and a final value of the coupling strength A s = 2J, we obtain a residual population fraction in the component σ = 3 of 3.2%; we remind that this population is a direct measure of the loss in our procedure. As visible in Table 1, the ramp duration is a crucial parameter: going from 2 · 10 3 (h/J) to 10 4 (h/J) allows for a nearly threefold decrease of loss. One also observes that the detuning parameter δ 0 introduced above significantly increases the loading efficiency (Table 1 compares δ 0 = J/50 with δ 0 = 0). Besides, we note that increasing the coupling strengths A 1,2 increases the dark to bright bands gap; accordingly, when considering the small value A 1,2 = 0.5J, the losses significantly increase (12.5%) for the same lattice configuration.
We note that this undesired population can be significantly reduced by optimizing the ramp function (as in [74]), or by increasing the system size and ramp duration; other statepreparation protocols, based on optimal-control theory [75][76][77] and nonunitary dynamics [78], could also be designed to further maximize the fidelity. As far as the system size is concerned, we observe that the loading fidelity saturates as one increases the system size (close to 242×22 in our simulations); this indicates that the remaining losses are no longer dominated by transfer through the edge modes, but rather stem from the bulk gap closing during the ramp.
Finally, we calculate the time-evolution of the realistically prepared state under the action of a linear gradient along the y direction, and show the resulting density in Fig. 3(d) for a drift duration t = 10(h/aF ). After filtering out the bright-state contribution to the density, by removing the aforementioned stripes and restricting the measurement to well populated sites only, one obtains the center-of-mass motion depicted by crosses in Fig. 3(e), from which one extracts the "experimental" value C exp D = 1.92. These simulations demonstrate that a dark band with C D = 2 can be loaded and probed accurately using a realistic Λ system [Eq. (13)].
The time-scales discussed in this work are expressed in terms ofh/J. The values of J achieved in modern experiments correspond to several hundreds of Hz (e.g. 880 Hz×h in [79]). Considering such a value of J, the loading ramp duration discussed in the main text and in Fig. 3 would correspond to T = 3 · 10 3h /J = 540ms. Besides, the total evolution time in the tilted lattice (Chern-number measurement) is T = 2 · 10 3h /J = 360ms. Those durations are compatible with the coherence times of cold-atom experiments [79].

Conclusion
We proposed a realistic scheme by which flat bands with predictable higher Chern number can be constructed through a coherent manipulation of Bloch bands. While this work focused on the simplest Λ configuration, more bands could be involved in the construction in view of building N -pod settings [80,81]; this strategy suggests a promising route to reach even higher Chern numbers |C D | 2. The dark Chern bands deriving from our scheme offer a platform where exotic fractional quantum Hall states could be explored with cold atoms, including generalized Moore-Read and Read-Rezayi states [24,26], topological nematic states [23] and "genon" defects [23,29,82]. Interesting perspectives concern the fate of dark Chern bands in the presence of dissipation, and the coupling of Bloch bands belonging to other topological classes [83].

A Dark state formula
The dark state formula [Eq. (3)] for the three level Λ system in Eq. (2) is valid in the case where δ = 0; no simple analytical formula for the dark state is known for the general Hamiltonian in Eq. (2). As moderate values of δ, ∆ are used in practice, we can use δ/Ω 1,2 and ∆/Ω 1,2 as small parameters from which a perturbative expansion can be constructed. For finite δ, admixtures lead to corrections to the dark state's energy, and to the state itself, and it is then referred to as the "gray" state. The above equation reveals the lowest-order correction to the dark state, which only involves the excited state |3 . Next terms that are O(δ 2 ) and O(δ∆) involve all three base states |1 , |2 , |3 . Importantly, when setting δ = 0, the detuning ∆ does not affect the E 0 eigenvalue and the state |D(0) remains to be an accurate dark state. As a result, in a Λ scheme made of Bloch bands [ Fig. 1(b) in the main text], the flatness requirements on the band associated with the "third" state |ψ 3 (k) ≡ |3 are substantially relaxed with respect to the other two bands.
In the standard situation where Ω 1 and Ω 2 couple different atomic states, the excited state |3 undergoes spontaneous emission. However, the spontaneous emission rate (typically in MHz range) dominates over single kHz energy scales for ultracold atom dynamics, and consequently even a small admixture to the dark state is detrimental.
In the setting considered in this work, where the excited state |3 is chosen as a stable optical lattice band, the admixture implied by (16) does not imply larger losses, and the Chern number is invariant with respect to perturbations of the topological dark band. This holds as long as the band gap separating dark and bright bands does not close. In fact, as seen in the following section, this requires δ/Ω 1,2 1.

B Rotating Wave Approximation
We now consider three Bloch states with fixed quasimomentum k that are coupled by two time-dependent processes of frequencies ω 1,2 . We hereby discuss the validity of Eq. (1) in the main text in describing this setting. We consider three states |ψ 1 (k) , |ψ 2 (k) , |ψ 3 (k) with energies E 1,2,3 . The Hamiltonian is written as where A, B are operators that describe the couplings, and similarly for B. We now transform to the rotating frame with the transformation ψ = The transformation to the co-rotating frame is achieved by: and analogously for B .
In the situation where all the states |ψ i (k) belong to different spin manifolds, the Rotating Wave Approximation (RWA) directly applies and all the time-dependent, rapidly oscillating terms in the above matrix A may be neglected. One then obtains If two states |1 and |2 are taken within the same spin manifold, the above approximation should be replaced with Indeed, the time-dependent terms ±A ( * ) 23 e ∓i(ω 1 −ω 2 )t contain the frequency (ω 1 − ω 2 ), which is orders of magnitude smaller then all others frequencies. In fact, it is typically of the order of J/h, as established by the Bloch band's bandwidth. However, we point out that the approximation in Eq. (19) is still justified in this case, as long as the amplitude verifies A 23 J; this approximation was assumed to be valid in the main text, when describing the Λ system by Eq. (2). For instance, when A is generated by a lattice modulation, the coefficients A 23 , A * 23 may be minimized by simply lowering the amplitude of the modulation. Alternatively, when different spin states are involved, A 23 can vanish due to selection rules. The justification for neglecting the terms oscillating with frequency (ω 1 − ω 2 ) is much more evident when the states |1 and |2 correspond to different hyperfine states; thenh(ω 1 − ω 2 ) is of the order of the hyperfine splitting, which is several orders of magnitudes larger than J. The same reasoning applies to the operator B, with B 12 and B 23 replacing the coefficients A 12 and A 23 .

C Dark state stability and the flatness of bands C.1 Dark state stability
The stability of the dark state is understood as the existence of a finite gap to the remaining two bright eigenstates in Eq. (2). Here we analyze the effects of finite δ on these gaps.
For δ(k) = 0 the bright state eigenenergies are: Since the dark state energy is E 0 = 0, |E ± (k)| also gives the amplitude of the gap between bright and dark states. It is simple to check numerically that if both Ω s (k) = 0, then the gap remains finite (though possibly small) for certain values of δ(k) or ∆(k). When exactly one of the Rabi couplings (Ω 1 (k) or Ω 2 (k)) is zero, one numerically obtains that the gap may close. The other, non-zero Rabi frequency can then be chosen as the energy unit.
We first consider the case when Ω 1 (k) = 0 and Ω 2 (k) = 0. Fig. 4(a) shows the eigenvalues of H Λ as solid, thin lines in that case. For δ(k) = 0, as expected, the dark state energy is precisely 0. We see that for large |δ(k)| the energy of one of the bright states become close to the energy of the dark/gray state, though never crosses it. The gap value is, therefore, substantially lowered which may lead to the depletion of the dark/gray state in the experiment. For ∆(k) = 0 the dark/gray energy level is crossed by [see Fig. (4)(a)] the bright state at a finite value of δ(k) (dashed lines). For ∆(k) = 0 the crossing occurs for finite value of δ(k). For example for ∆(k) ≈ 0.4Ω 2 (k), the crossing takes place when δ(k)/Ω 2 (k) ≈ −0.6, and the smaller the |∆(k)|, the larger δ(k) is necessary for the crossing to occur.
When Ω 2 (k) = 0 and Ω 1 (k) = 0, the dark/gray state energy line is exactly linear with δ(k) as |D(k) = |ψ 1 and |D(k) is decoupled from |ψ 2 (k) and |ψ 3 (k) . In particular, for δ(k) = ±Ω 2 (k)/2 the gap to the bright state closes. For Ω 2 (k) = 0 this energy level crossing becomes avoided [see Fig. 4(b)]. In this case the value of the detuning ∆(k) does not change the location of the crossing, or whether it is avoided or not. The presence of the energy level crossings (either exact or narrowly avoided) would have a detrimental effect on the construction of the dark-state band with a well-defined Chern number. However, in both cases discussed above, the crossings require large values of δ(k): i) Ω 1 (k) = 0 and δ(k) = O(Ω 2 (k)) or ii) Ω 2 (k) = 0 and δ(k) = O(Ω 1 (k)). Nevertheless, when the dependence of Ω 1,2 (k) on the quasimomentum k is considered, one notices that Ω 1,2 (k) may have zeros where Ω 1,2 (k) is also small relative to its maximum value. When the values of Ω 1,2 (k) rapidly oscillate across the BZ, a strong upper bound on δ(k) is therefore implied.
The first two requirements can be summarized in the following inequality (for each quasimomentum k): Here∆E i is the minimum distance of the band associated with bare state |i to its nearest neighbors. In the example discussed in the main text,∆E 1 = O(J) simply corresponds to a gap between the lowest two Hofstadter bands; the gap∆E 2 is between the top and secondtop bands [in the specific choice of bands described in the main text,∆E 2 =∆E 1 due to the symmetry of the spectrum; see Fig. 2(a)]; the gap∆E 3 is the standard separation between s and p bands, in an optical lattice potential V (x, y) = V x cos 2 (kx/q) + V y cos 2 (ky). The latter where E R is the recoil energy; sufficiently large optical-potential amplitudes V x , V y can easily ensure sufficiently large∆E 3 . Altogether, onlȳ ∆E 1,2 give an upper bound on practical Rabi frequencies Ω 1,2 (k). It is useful to consider the ratio measuring the total variation of the Rabi frequency across the BZ: The value of this quantity provides an indication on how flat the bands associated with the states |ψ 1 (k) and |ψ 2 (k) have to be in view of tuning the overall amplitude of the Rabi frequencies so as to satisfy inequality (21). The flatness of the state |ψ 3 (k) is not so crucial [see Appendix A]. The definition of g depends implicitly on the choice of the bands b 1,2 ; we will use the subscripted notation g b 1 b 2 below, where b 1 , b 2 indicate the bands chosen from the Hofstadter spectrum for the states |1 , |2 in the Λ system.
For a particular choice of the model and bands, we consider the following flatness f -factor: where ∆E i indicates the bandwidth of a given band. The inequality (21) can be satisfied when The following two subsections will discuss the values of f and g in different cases, with an emphasis on cases where the former is maximized and the latter minimized.

C.3 Band flatness -flatness factor f
We consider the Hofstadter-Harper model under PBC in the thermodynamical limit [for numerics, system sizes O(100 × 100) suffice]. In Fig. 5(a), (b) and (c) we show the bandgap∆E i , bandwidth ∆E i and the flatness ratio f for different q, for the lowest three bands and the central band (with C = −q + 1) of the HH model. We notice that∆E i /J is always O(1), except for the band with C = 1 where it quickly drops, and that the bandwidth ∆E i /J drops exponentially with q already for the considered moderate values of q. As a result, the factor f exponentially increases with q for bands with C = 1.
The behavior for C = −q +1 is starkly different. The bandgap decreases as ∼ 1/q 2 and the bandwidth also decreases like ∼ 1/q with q. In the end, the maximal flatness factor f ≈ 6.5 for the middle band is for q = 6 and it drops down to zero as ∼ 1/q.
In light of this analysis, the use of the middle band (with C = 1), to provide states |1 or |2 in our Λ system, is questionable.  − d) show k-dependence across the entire BZ of log 10 |Ω(k)|, the Rabi frequency coupling a trivial flat s-band to lowest, first-excited, top, and middle band of a HH model, for a flux φ = 1/11. Here, max k |Ω(k)| = 1 sets the normalization. Panel (e) shows factor g 1,11 (g 1,2 ) as a function of m = max k∈BZ |Ω 1 (k)| with solid, black (dashed red) lines. Panel (f ) shows log 10 g b 1 b 2 for m maximizing its value for different (b 1 , b 2 ). Red squares correspond to (b 1 , b 2 ) where g b 1 b 2 = 0.

C.4 Rabi frequencies total variation -factor g
Here we detail the dependence of the Rabi frequencies Ω 1,2 on the quasimomentum k for the HH model. This aspect is important in view of satisfying Eq. (21).
First let us remark that the Rabi frequency coupling two bands with different Chern number has to be zero somewhere in the BZ. Indeed, if we assume a contrario, that Ω s (k) = 0, for all k ∈ BZ, we can construct the following system: H = δ|ψ s ψ s |+Ω s (k) (|ψ s ψ 3 | + |ψ 3 ψ s |) for all k ∈ BZ. If the detuning δ is swept from a large positive value δ |Ω s (k)| to a large, negative value: δ −|Ω s (k)| then the eigenstate of this model changes adiabatically between |ψ 3 and |ψ s , which would imply C 3 = C s , violating our assumptions. This argument says nothing about Ω s (k) connecting bands with same Chern number.
We numerically find that in the example studied in the main text, where the two bands b 1 and b 2 come from the HH model, and the band b 3 is a s band of a (topologically trivial) 2D optical potential, Ω 1 (k) and Ω 2 (k) have coinciding zeros for roughly half of the choices of the two Hofstadter bands [see Fig. 6(f)] .
The value of g b 1 b 2 does not change if both Ω 1 (k) and Ω 2 (k) are multiplied by a common factor. We assume max k |Ω 2 (k)| = 1 and max k |Ω 1 (k)| = m. Figure 6(e) shows g 1,11 as a function of m [compare panels (a) and (c)]. For some m = m * the g b 1 b 2 attains the maximum value. It is also evident that the band choice can alter the factor g b 1 b 2 by several orders of magnitude. For the choice of (b 1 , b 2 ) = (1, 11) for q = 11 the maximal value of g 1,11 is g 1,11 ≈ 0.3. Figure 6 (f ) shows log 10 |g b 1 b 2 | for different choice of b 1 , b 2 bands as a color array plot (always for m = m * ). We notice that for some values (indicated by deep red square) of b 1 , b 2 (for example b 1 = b 2 or b 1 = 1, b 2 = 3) we have g b 1 b 2 = 0. This is due to coinciding zeros of Ω 1 (k) and Ω 2 (k).

C.5 Numerical evaluation of Chern numbers
The Chern number is defined as the integral of the Berry curvature over the BZ [Eq. (6)]. The integral can be directly computed using standard methods for numerical integration. In this work we used a regular discretization of BZ into N × N evenly-sized pieces, and integration using trapezoid prescription. When N is sufficiently large, the integral for C gives an error E ∼ N −2 as expected. In particular, applying this method, the approximation for C is not an integer number.
An alternate method has been proposed by Fukui et. al. [84]. There the integration is performed also by summation over rectangular plaquettes, but the approximant for the Chern number integral is manifestly gauge invariant, and the Chern number approximation is guaranteed to be an integer. Figure 7(a) shows the comparison of the two methods. There we have computed the Chern numbers of the lowest band (C = 1) for the standard HH model for q ∈ {3, 11} using brute-force discretization and Fukui's method. In this case the method proposed by Fukui offers an obvious numerical advantage, and returns correct Chern numbers even for smallest considered discretizations. In Panel (b) of the same Figure we compute the Chern number of the dark state (C D = 2) as discussed in the main text. Panel (c) shows the computation of the Chern number of a dark state band where |1 , |2 correspond to the lowest and middle bands of the HH model (C D = −q + 2). In the latter case the brute force approach converges up to the correct result much faster than the Fukui's method.