Dissipative preparation of fractional Chern insulators

We report on the numerically exact simulation of the dissipative dynamics governed by quantum master equations that feature fractional quantum Hall states as unique steady states. In particular, for the paradigmatic Hofstadter model, we show how Laughlin states can be to good approximation prepared in a dissipative fashion from arbitrary initial states by simply pumping strongly interacting bosons into the lowest Chern band of the corresponding single-particle spectrum. While pure (up to topological degeneracy) steady states are only reached in the low-flux limit or for extended hopping range, we observe a certain robustness regarding the overlap of the steady state with fractional quantum Hall states for experimentally well-controlled flux densities. This may be seen as an encouraging step towards addressing the long-standing challenge of preparing strongly correlated topological phases in quantum simulators.


I. INTRODUCTION
Strongly correlated topological phases that conceptually elude both the independent-particle approximation and a classification in terms of local order parameters exhibit some of the most fascinating phenomena known in nature [1][2][3]. The influence of dissipation on such systems, while being to some extent inevitable, may be seen as a mixed blessing. On one hand, thermal fluctuations and strong environmental couplings may challenge the desirable topological quantization of observables as well as the coherence of quantum information encoded in topological states. On the other hand, from a viewpoint of state preparation, some form of dissipation is essential to reach topological phases, as is clear from their very definition as equivalence classes under local unitary transformations [2,4]. More specifically, local coherent physical processes, mathematically described precisely by such local unitary transformations [4], would never allow a system to dynamically enter a topological phase.
For electronic materials, electron-phonon coupling is typically capable of establishing thermal equilibrium with their surroundings, and the practical challenge in preparing a topological state thus amounts to reaching cryostat temperatures below the scale set by the energy gap protecting the targeted topological phase in a given material [1]. By contrast, a main feature of synthetic materials, e.g., based on ultracold atoms in optical lattices [5], is their high degree of quantum coherence and experimental control over microscopic processes [6][7][8][9][10]. While detrimental effects of dissipation may thereby be strongly contained, a generic analog of the aforementioned cryostat cooling is not naturally present in such quantum simulators [5,11]. Instead, engineered dissipation [12][13][14][15] has been proposed as a means to prepare complex manybody states in a nonequilibrium fashion as steady states of a quantum master equation [16,17] from basically arbitrary initial states. This approach may be understood as a tailor-made cooling or entropy reduction protocol for FIG. 1. The band structure for (a) the KM model [35] with q = 2 and (b) the Hofstadter model [39] with q = 3. The red arrows indicate the dissipative pumping of particles from higher bands to the lowest band, as described by the Lindblad operators Lm in Eq. (4). a specific target state rather than a generic thermalization process, and has been discussed for a range of topological states [18][19][20][21][22][23][24][25] within the realm of non-interacting topological insulator phases [26]. By contrast, there is relatively little corresponding effort [27,28] on the engineered dissipation of strongly correlated states.
Here, turning to strongly correlated dissipative systems, we microscopically simulate the dissipative state preparation of paradigmatic fractional quantum Hall (FQH) states [1,[29][30][31][32][33] in the hard-core boson limit [34][35][36] by numerically exact methods. In particular, we construct and study a quantum master equation from which lattice FQH states can be shown to emerge as exact and unique steady states (see Fig. 1 for an illustration). Considering various deviations from this analytically amenable model, we study the robustness of the resulting steady states and their topological properties. While the notion of adiabatic continuity can be readily generalized from Hamiltonian systems to Gaussian states evolving under bilinear Liouvillians [18,37,38], the absence of such general continuity arguments for our present setting of master equations involving quartic terms serves as a main motivation for our in-depth computational study.
Notably, for the celebrated Hofstadter model [39] of strongly interacting bosons subject to a magnetic flux [40], we provide numerical evidence that the Laughlin FQH state [30,34] can be prepared by a dissipation simply pumping particles to the lowest band of the corresponding noninteracting model. While a pure (up to topological degeneracy) FQH steady state is only reached in the low-flux limit, our data indicate that a satisfactory approximation corresponding to an FQH phase at low temperature may already be reached at experimentally well-studied finite flux [21,41]. This result is encouraging regarding the long-standing goal of preparing strongly correlated topological phases in synthetic materials. Complementary to our present study, protocols for dissipatively stabilizing FQH states via the quantum Zeno effect [42], self-stabilization of initial FQH states [43], and preparing few-body FQH states [44] have been reported, and a nonequilibrium topological field theory has been constructed for free Chern insulators from quartic Liouvillians [24]. Furthermore, open systems in an FQH regime have been recently studied [45] employing an effective non-Hermitian Hamiltonian approach [46,47], as well as considering their stability against quantum jumps [48].

A. Static Hamiltonian
We study N on-site interacting bosons in a periodic two-dimensional (2D) square lattice in the xy plane with unit lattice spacing. Each lattice site j is labeled by its position (x j , y j ), where x j = 0, 1, · · · , L x − 1 and y j = 0, 1, · · · , L y − 1. The lattice is pierced by a uniform magnetic field, such that the number of flux quanta in an elementary plaquette is a rational number φ = p/q with coprime integers p and q, leading to q magnetic Bloch bands. We choose q sites in the x direction as a magnetic unit cell, so the numbers of unit cells are N x = L x /q and N y = L y in the x and y directions, respectively. In this paper, we assume L x to be divisible by q to ensure an integer number of unit cells and focus on φ = 1/q. Based on the scenario described above, we identify desirable target states for our subsequent dissipative state preparation analysis by first considering the static physics governed by the tight-binding Hamiltonian where a † j (a j ) creates (annihilates) a boson on site j, U > 0 is the strength of the onsite repulsion, and t jk is the hopping coefficient between sites j and k.
The single-particle physics is determined by t jk . If we choose where with x = x j −x k and y = y j −y k , the setup corresponds to the Kapit-Mueller (KM) model [35] in the Landau gauge. Note that in this case the hopping strength decays exponentially with the hopping range, and this decay is quicker for smaller φ. At φ = 1/q, the lowest band of the KM model has an elegant analytical property: It is spanned by the discretized version of the continuum lowest Landau level (LLL) wavefunctions and is exactly flat carrying Chern number C = 1 [35]. This analytical similarity to the LLL also allows for the generalization of a number of many-body results found in the continuum LLL. In particular, the ground state of the static Hamiltonian Eq. (1) at filling ν ≡ N/(N x N y ) = 1/2 is the exact Laughlin state discretized to the lattice and possesses zero energy [35]. These exact Laughlin states on the lattice completely reside in the lowest band of the KM model, and their vanishing property (i.e., two bosons cannot be located at the same position) guarantees that the occupation of bosons forming these states is no more than one on each lattice site. These properties survive even when the interaction strength U → +∞. In this limit, we can see from Eq. (1) that the exact Laughlin state is the unique state at ν = 1/2 that completely resides in the lowest KM band.
Besides the KM model, we will also consider other lattice models with t jk truncated to the nearest-neighbor (NN) and next-nearest-neighbor (NNN) sites. When only the NN hopping exists, we have the usual Hofstadter model [39]. Note that the capability of realizing the exact ν = 1/2 Laughlin state on the KM lattice originates from the analytical similarity between the lowest KM band and the LLL, and is thus in general impossible in other lattice models.

B. Dissipative dynamics
Our aim is to theoretically model and analyze the preparation of the ν = 1/2 bosonic Laughlin phase through dissipative dynamics. Therefore we fix the filling at ν = 1/2 throughout this paper. At t = 0, we prepare the system in a state described by a density matrix ρ(0). Under the assumption of a Markovian process, the dissipative dynamics is governed by the master equation where ρ ≡ ρ(t) is the density matrix of the system at time t and L m 's are the Lindblad operators of strength κ m accounting for dissipative processes within the Born-Markov approximation. Motivated by the fact that for the KM model the exact ν = 1/2 Laughlin state is the unique state completely residing in the lowest band when U → +∞, we choose the dissipative processes that pump bosons from excited bands to the lowest band and meanwhile assume the hard-core limit for bosons. Such processes may, for example, be induced by collisional coupling of the bosons to atoms in a low-temperature Bose-Einstein-condensate (BEC) [49]. The generic setting of dissipative pumping to the lowest band will be considered not only for the KM model but also for models with truncated hopping. In this scenario, the Lindblad operator is of the form k , which moves a particle from an excited band s > 1 to the lowest band and meanwhile scatters its crystal momentum from k to k (see Fig. 1 for an illustration). Here, γ s k † creates a particle with momentum k in band s, and k is restricted in the first Brillouin zone (1BZ). There are (q − 1)(N x N y ) 2 Lindblad operators, where q−1 is the number of higher bands and N x N y is the number of k points in the 1BZ. For simplicity, we set κ m = κ for all Lindblad operators.
To impose the hard-core boson condition, we Fouriertransform these Lindblad operators and the master equation Eq. (4) from momentum space to real space (see details in Appendix A). Then we construct the real-space Fock basis in which the number of bosons on each lattice site is no more than 1. The Lindblad operators, the many-body Hamiltonian, and the density matrix are all represented as matrices under this basis. In the hardcoreboson limit, only the hopping term exists in the manybody Hamiltonian Eq. (1), while the hardcore constraint is encoded in the restricted Hilbert space as compared with bosons with finite interaction.
When simulating the dynamics, we divide time into many steps of small interval ∆t and discretize Eq. (4). Given the state ρ(t) at time t, the state at time t + ∆t is evaluated by For our specific setting, our examination shows that it is typically sufficient to choose κ∆t = 0.1 and further reduction of ∆t gives very similar results. Within the limit of our numerical resources, we can only deal with at most four bosons at a numerically exact level. In this context, it is worth noting that the dimension of the considered vector space in which the density matrix is defined grows quadratically with the Hilbert space dimension of pure state vectors. We will focus on the pure dissipation case, i.e., we neglect the effect of the − i [H, ρ] term. This is reasonable when the Hamiltonian and the Lindblad operators come from the same lattice model such that they are compatible with each other, or in a scenario where the target system may be viewed as a (to good approximation) flat band of a Hamiltonian. We have checked that bringing the − i [H, ρ] term with a compatible Hamiltonian back leads to even quantitatively quite similar results and identical conclusions.
For a specific lattice model (KM or truncated model), we characterize the dynamics of ρ(t) by two quantities. First, its overlap with the twofold degenerate topologically ordered ground states |Ψ i of the corresponding static Hamiltonian on the torus, which are numerically obtained by diagonalizing the Hamiltonian shown in Eq. (1). For the KM model in the hardcore limit, the Hamiltonian ground states are simply the exact Laughlin states, and otherwise FQH states in the same topological phase as the ν = 1/2 Laughlin state. Second, its weight in the lowest band. The calculation of O(t) is straightforward because both |Ψ i and ρ(t) are expressed in the real-space basis corresponding to the square lattice geometry. When evaluating W(t), we need to transform γ 1 k † γ 1 k to the real space first (see details in Appendix A).

III. DYNAMICS OF THE KM MODEL
As for the KM model the exact Laughlin state is the unique state at ν = 1/2 that completely resides in the lowest band, we expect that our state-preparation protocol will select the exact ν = 1/2 Laughlin state as the steady state in this case for any initial state with the right particle number. To numerically confirm this, we consider a pure initial state, which can be either a product state randomly chosen from the Fock basis, or a random superposition of all basis states. We find similar results in both cases, so we demonstrate the dynamics starting from a product state in the following.
The dynamics at flux densities φ = 1/2 and φ = 1/3 is shown in Fig. 2. For all system sizes that we have studied, both O(t) and W(t) tend to 100% in the long-time limit. This unambiguously indicates that for the KM model the dissipation governed by particles' jumping from higher bands to the lowest band indeed drives the system to a steady state which is completely located in the subspace of exact Laughlin states, i.e., a unique steady state up to topological degeneracy. The purity of ρ(t), defined as Trρ 2 (t), is 0.5 in the long-time limit, meaning that the steady state is a mixed state with 50% probability on each of the two exact Laughlin states on the torus.

IV. BEYOND THE KM MODEL
Having confirmed that we can perfectly prepare the exact ν = 1/2 Laughlin state in the KM model by the dissipative dynamics, we find it interesting to investigate to what extent the strategy of pumping all particles to the lowest band of a single-particle picture still works when deviating from the KM model. The answer to this question is unclear due to the absence of analytical similarity to the LLL in generic lattice models. In particular, it is not obvious whether a state close to the exact ν = 1/2 Laughlin state exists completely in the lowest band of the pertinent model in the hard-core limit. Furthermore, in general there is no strict adiabatic continuity in the steady state of non-Gaussian Liouvillians towards changes in the Lindblad operators.
To be concrete, let us first focus on the truncated   Fig. 3(b)] for the corresponding system sizes.
The worse performance of preparing the Laughlin phase in the Hofstadter model (especially at high flux density) is due to the absence of a state thereof at ν = 1/2 that is close to the exact Laughlin state and completely resides in the lowest Hofstadter band in the hard-core limit. In Fig. 4, we plot the weights in the lowest Hofstadter band for all eigenstates of the Hofstadter model at ν = 1/2 in the hard-core limit. While the ground states have higher weights than excited states, they do not completely reside in the lowest Hofstadter band, which is more obvious for larger φ (Fig. 4). Therefore the strategy of pumping particles to the lowest Hofstadter band does not necessarily drive the system to the ground-state subspace of the Hofstadter model, which has a high overlap with the exact Laughlin state. Instead, some excited eigenstates far from the exact Laughlin states may be involved, leading to the relatively low overlap between the steady state and the exact Laughlin state. This is also supported by our numerical observation that the weight of the steady state in the lowest Hofstadter band is lower than that of the static ground state (95.7 versus 99% for three bosons on the 2 × 3 Hofstadter lattice at φ = 1/3), suggesting the contribution from static excited eigenstates to the steady state [ Fig. 3(b) and Fig. 4]. By contrast, for the KM model at ν = 1/2, all particles can jump to the lowest KM band and meanwhile satisfy the hard-core condition by forming the exact Laughlin state, such that our dissipation mechanism naturally selects the Laughlin state as the steady state.
There are two ways to gain better preparation of the Laughlin state in the Hofstadter model. One is to decrease the flux density, because in this case the Hofstadter model is closer to the KM model due to the fast decaying of longer-range KM hopping. Indeed, for three bosons on the 2 × 3 Hofstadter lattice, the overlap between the steady state and the static Hofstadter ground state increases from 57.0 to 88.8% (the overlap with the Alternatively, we can add longer hopping in the Hofstadter model. Remarkably, we observe a significant improvement of the results once we add the next-nearestneighbor (NNN) hopping. In Fig. 5, we show the results obtained by keeping the NN and NNN hopping in the original KM model for three bosons on the 2 × 3 lattice. Compared with the data in Fig. 3, the long-time overlap with the static ground state of the model dramatically grows from 57.0 to 99.4% (the overlap with the exact Laughlin state dramatically grows from 60 to 99.4%) for φ = 1/3 [ Fig. 5(a)], accompanied by an increasing of the lowest-band weight from 95.7 to 99.9% [ Fig. 5(b)]. A similar improvement can also be seen for φ = 1/4 (Fig. 5).
In summary, the proposed state-preparation protocol of band pumping still works excellently at low flux densities for realistic models with only the NN hopping. Further improvement can be obtained by adding slightly longer range hopping.

V. CONCLUDING DISCUSSION
Motivated by the important open challenge of preparing strongly correlated topological phases in quantum simulators, we have numerically studied the Liouvillian dynamics of microscopic lattice models towards fractional quantum Hall steady states. For the experimentally wellstudied Hofstadter model in the low-flux limit, as well as for the Kapit-Mueller model characterized by an extended hopping range, bilinear Lindblad jump operators (leading to a quartic Liouvillian) that yield a pure FQH steady state (up to topological degeneracy) can readily be constructed guided by the physical picture of pumping particles into the lowest Chern band. Beyond these exact limits, when generalizing the aforementioned band pumping picture to experimentally studied finite flux densities, we find an encouraging robustness in the sense of moderate and continuous reduction of the steady state overlaps with FQH model states in a wide parameter range. On a more general note, compared with adiabatic preparation protocols of topological phases the paradigm of dissipative state preparation does not suffer from an unavoidable critical slow-down when approaching the thermodynamic limit.
While this paper provides a first fully microscopic study on the dissipative preparation of FQH states from arbitrary initial states, several issues clearly remain interesting subjects of future research. First, the pumping rates of the Lindblad jump operators are assumed constant for simplicity. Estimating the momentum dependence of such rates in an experimentally realistic setting of cold atoms with collisional coupling to a BEC is one direction devising a feasible protocol based on our general analysis. Second, in the interest of computational feasibility, our numerically exact study is based on systems of hard-core bosons, which naturally raises the question of the influence of finite contact interactions. To address this issue qualitatively, we argue that taking the hard-core limit is by no means expected to be pathological in the present context. Specifically, an experimental protocol could start from a Mott state without initial double occupations. The dynamical creation of such energetically very costly double occupations (that would kick the system out of the considered hard-core Hilbert space) may then be practically prohibited already at large but finite interaction strength. This is because only energetically resonant processes (in the combined system and bath setting) contribute to the Lindblad operators due to the underlying Born-Markov approximation. Finally, although generalizations of the KM model to non-Abelian and higher-Chern-number states provide a natural starting point [50], a similarly simple mechanism for the preparation of more complex FQH states than Laughlin phases remains to be identified. where R m is the position of the mth unit cell, a † m,α creates a boson in the αth site of the mth unit cell, and [v s 1 (k), · · · , v s q (k)] is the eigenvector of band s. The master equation with pure dissipation can then be written in real space as ×e ik·(Rm−R m ) e −ik ·(Rn−R n ) a † m,α a n,β |i j|a † n ,β a m ,α , ×e ik·(Rm−R m ) e −ik ·(Rn−R n ) a † n ,β a m ,α a † m,α a n,β |i j| + |i j|a † n ,β a m ,α a † m,α a n,β , where |i is the Fock basis state of the many-body Hilbert space of dimension D and ρ i,j is the matrix element under this basis. Once we express the master equation in real space as Eqs. (A4) and (A5), we can easily impose the hardcore constraint by working in the Fock basis where the boson occupation per site is either 0 or 1. The weight in the lowest band can be evaluated by β * (k)e ik·(Rm−Rn) j|a † m,α a n,β |i . (A6)