Dissipative Preparation of Antiferromagnetic Order in the Fermi-Hubbard Model

The Fermi-Hubbard model is one of the key models of condensed matter physics, which holds a potential for explaining the mystery of high-temperature superconductivity. Recent progress in ultracold atoms in optical lattices has paved the way to studying the model's phase diagram using the tools of quantum simulation, which emerged as a promising alternative to the numerical calculations plagued by the infamous sign problem. However, the temperatures achieved using elaborate laser cooling protocols so far have been too high to show the appearance of antiferromagnetic and superconducting quantum phases directly. In this work, we demonstrate that using the machinery of dissipative quantum state engineering, one can efficiently prepare antiferromagnetic order in present-day experiments with ultracold fermions. The core of the approach is to add incoherent laser scattering in such a way that the antiferromagnetic state emerges as the dark state of the driven-dissipative dynamics. In order to elucidate the development of the antiferromagnetic order we employ two complementary techniques: Monte Carlo wave function simulations for small systems and a recently proposed variational method for open quantum systems, operating in the thermodynamic limit. The controlled dissipation channels described in this work are straightforward to add to already existing experimental setups.


Introduction
Experimental progress with ultracold fermions in optical lattices [1,2] leads the way to achieving one of the key goals of quantum simulation [3]-mimicking realistic condensed matter systems. To date, the experiments covered a broad range of systems and interaction regimes, from probing the BEC-BCS crossover in lattices [4], to the observation of a fermionic Mott insulator [5,6], to studying short range magnetism [7] and multiflavor spin dynamics [8], to realizing topological Haldane model [9] and artificial graphene sheets [10]. These discoveries pave the way to use ultracold atoms to reveal the properties of the repulsive Fermi-Hubbard model [11,12]. The latter is of particular importance since it represents a playground to get insight into the physics of high-temperature superconductivity and related phenomena observed in the cuprates [13].
In the case of one particle per site and large on-site interaction, U, the Fermi-Hubbard model exhibits the transition to the Mott-insulating state [5,6] around the temperatureT U. If the temperature is decreased further and reaches the so-called 'Néel temperature,'T t U 4 N 2 , where t gives the hopping rate between neighboring sites, the transition to the antiferromagnetic (AF) phase is expected [14,15]. Currently the temperatures achievable in experiment are slightly above the Néel temperature where AF correlations can already be observed, for instance, » T T 1.42 N has been reached in [14]. Ultimately, in order to study the superconducting phase or other phenomena related to pairing in high-temperature superconductors, the temperature needs to be substantially lower. Therefore, due to the experimental limitations inherent to the standard laser cooling techniques, it is crucial to develop alternative approaches [16][17][18][19][20][21][22][23][24][25][26][27][28] to preparation of quantum phases in optical lattices.
We start with a system of fermions in an optical lattice as described by the Fermi-Hubbard model. The parameters of the Hamiltonian are left intact, instead we introduce dissipative channels on top of the unitary evolution. As a result, fermions remain mobile in the optical lattice during the entire dissipative preparation stage, which should help with retaining coherence after the dissipation channels are switched off. Furthermore, the dissipation channels of our scheme are implemented using the level structure of fermionic 40 K, currently used in several laboratories [7,8,15,[60][61][62]. Consequently, the presented scheme can be readily implemented into already existing experimental setups.
Theoretical description of open many-body quantum systems represents a challenging task and is currently an active field of research [45,[63][64][65][66][67][68]. In our analysis of the dissipative Fermi-Hubbard model we use two complementary techniques: the Monte Carlo wave function (MCWF) [69][70][71] and the variational method [45,64], which is generalized here to the description of fermionic systems at half-filling. By using these two methods we demonstrate that a substantial AF magnetization is present in the system both for an exact solution on a 3×3 lattice, as well as in the thermodynamic limit.

The dissipative Fermi-Hubbard model
We start with the Fermi-Hubbard Hamiltonian which has been experimentally realized in a range of systems such as 6 Li [14] and 40 K [6]. Our goal is to design dissipative processes in such a way that the state with an AF order is the dark state of the dissipative dynamics and the time evolution of the open system will drive it towards such a dark state.


where we set º  1 and the primed sum runs over nearest-neighbor sites. Since we start with a disordered sample, all possible nearest-neighbor configurations, including   ñ | ; ,   ñ | ; , and   ñ | ; 0 will be present. The jump operators, therefore, need to convert the latter into those with the local AF order,   ñ | ; . We choose the jump operators to be as follows: (as labelled by g 2 ). As a result, the dissipative dynamics drives the system towards the AF phase.
We note that the jump operatorsˆ( ) j ij 3 break the SU(2) symmetry, as the down-spin atom becomes more mobile. This, however, does not constitute a limitation of our scheme. Moreover, it is possible to reestablish the symmetry by using additional auxiliary states to induce hopping of the  ñ | -state atom away from the doubleoccupancy configuration.
As we discuss in the following section, such choice of the jump operators is straighforward to realize in experiment using incoherent laser scattering.

Experimental implementation of the jump operators
Simulating the Fermi-Hubbard model requires mapping the two spin states,  ñ | and  ñ | , onto the fine or hyperfine components of the ground electronic state manifold of an ultracold atom. The on-site interaction between the spin components, U, can be tuned using a Feshbach resonance [72]. We exemplify the scheme using the atomic level structure of fermionic K 40 [73], with the levels  ñ º = = -| F m ; F W r3 realize the Raman-assisted hopping between the   ñ | ; and   ñ |0; states. If the Raman beams, W ¼ r r 1, , 3 , are not phase-locked such hopping processes are dissipative. Since the Raman-assisted hopping takes place directly between the initial and final states of the jump operators, the related dissipative processes are bidirectional. Therefore, we need to avoid populating the   ñ | ; state at this stage. Otherwise, the jump operators would also lead from the   ñ | ; state back to the   ñ | ; and   ñ | ; states. The on-site interaction energy U (which we assume to be on the order of a few kHz) can be used for this purpose. We note that this step can be implemented in a coherent way as well, however, the required phase-locking of the lasers would introduce an additional complication into the experimental setup.
In order to implement the second-stage jump operators,ˆ( ) j ij 3 , in a one-directional fashion, we use Ramanassisted hopping (lasers W r4 and W r5 ) from the   ñ | ; 0 state to the  ñ | X ; configuration with an auxiliary ñ |X state. This ñ |X state is then pumped (laser Ω) to an excited ñ |e state, which can decay to the  ñ | state completing the process, as illustrated in figure 1(c). The ñ |e state cannot decay to the  ñ | state because of selection rules on the F quantum number. The Zeeman splitting, D  X , and the energy, + D  U X , differentiate among the three states from the lower band in figure 1(c). In order to resolve between these three states it is sufficient to use selection rules. To resolve between the  ñ | X ; and  ñ | X; 0 states as the final states of the Raman-assisted hopping process we need nonzero on-site interaction between the  ñ | and ñ |X states. The typical values of the background scattering lengths, = a a 105 0 , in units of the Bohr radius a 0 [73,74], should be sufficient for this purpose. Spontaneous emission from the ñ |e state ensures that the resulting dissipative processes are unidirectional and take place from the   ñ | ; 0 to the   ñ | ; state with AF ordering. In total our dissipative preparation scheme requires six lasers (W r1 , ..., W r5 , and Ω) in order to induce the dissipative transitions, in addition to the lasers used to trap the atoms and prepare them in the  ñ | and  ñ | states. Due to the s-wave character of the S 2 1 2 manifold, the optical lattice potential for the  ñ | ,  ñ | , and ñ |X 5 , the pumping beam by Ω, and the decay rate is given by γ; D  gives the hyperfine splitting between the  ñ | and  ñ | states, D  X gives the Zeeman splitting between the ñ |X and  ñ | states, while U denotes the on-site interaction between the  ñ | and  ñ | states.
states will be the same, provided the optical-trapping lasers are sufficiently far-detuned from the P 2 3 2 states (see the discussion in [75], section 2.B).
The values of g 1 and g 2 are unrelated to the nearest-neighbor hopping integral t, however, all three of them are proportional to certain integrals involving two Wannier functions on the nearest-neighboring sites. As a safe estimate we consider t, g 1 , and g 2 to be on the same order of magnitude.
Please note that the presented dissipative-preparation scheme is general and other states can be used as  ñ | ,  ñ | , ñ |X , and ñ |e . In particular, it should be possible to use the states of 40 K, for which the Feschbach resonance is already known, i.e.  ñ = -| ; 2 , which requires a two-photon process from ñ |X to ñ |e . This could be realized with Raman beams detuned from the D 2 5 2 state or directly as a two-photon excitation. In the second case one could use ñ = -|X , 9 2 7 2 and ñ = -|e , 2 .

Time evolution and steady-state properties
In order to reveal the properties of the system we use two complementary techniques: the MCWF technique [69-71] on a 3×3 lattice and the variational method [45,64] in the thermodynamic limit. While the latter was originally formulated for bosons, here we extend it to fermionic systems. In both methods we start from the Jordan-Wigner transformation in two spatial dimensions [76][77][78][79]. The related Jordan-Wigner strings restrict the applicability of our variational scheme to the situation with one particle per site (half-filling).
Experimentally, this is the most interesting regime as it corresponds to the maximal Néel temperature.

Monte Carlo wave function
We study the dynamics of the driven-dissipative system governed by equation (2) using the MCWF method implemented in the QuTiP numerical library [80,81]. We consider only the nearest-neighbor hopping ºt t ij (we use t as the unit of energy hereafter) and study the time evolution of the half-filled 3×3 lattice as a function of the parameters g 1 , g 2 , and U. Since the lattice dimensions are odd numbers, we use the antiperiodic boundary conditions. This is required in the presence of the AF ordering because a particle hopping to its nearest neighbor across a boundary does not change the sublattice index, whereas for the same process within the boundaries such a change occurs. Therefore, when the hopping process takes place across the boundary, we introduce an additional spin-flip ( s ŝ † c c i j ), which ensures that hopping processes within and across the boundaries are equivalent. We also introduce analogous corrections to the jump operators in equations (3)- (5).
The initial states for the MCWF realizations were chosen with randomly-positioned spin-up or spin-down particles (also allowing for double occupancies), however, the steady-state properties were found to be independent on the initial conditions.
The time evolution of the systemʼs properties is shown in figure 2. One can see that the steady state is reached for t´»t 50 100, which corresponds to t = -1 2 s for t=50 Hz. Even for a large value of the on-site interaction, U=100, a small number of double occupancies is still present in the system, » D 0.04, see figure 2(a). These states are involved in the dissipation processes as an intermediate step towards preparation of the AF ordered phase, see figure 1(a). Non-zero double occupancies can lead to inelastic losses of atoms [11], which however are more problematic for the attractive [82], than for the repulsive [6] potassium gas. In the latter case the inelastic decay time for atoms on doubly occupied sites was reported [6] to exceed 850ms. Consequently, such inelastic losses should not constitute a limitation of our scheme.
To quantify the AF ordering of the system, we evaluate the spin-structure factor, as defined by . Note that the Néel state is an 'ideal' AF state, to which the dissipative processes would drive the system in the limit of g g U t , , 1 2  . This is because the jump operators suppress intersite coherence in the system. However, the fact that fermions remain mobile in the optical lattice during the entire dissipative preparation stage should help with retaining coherence after the dissipation channels are switched off. Figure 2(b) illustrates the decrease of the systemʼs entropy per particle, Tr log tot , with time. Due to the large size of the systemʼs density matrix ρ, in our numerical calculations we use the  Here y y = á ñ | A ij i j is the matrix of overlaps of wave functions obtained from single realizations of the Monte Carlo algorithm. The resulting steady-state entropy per particle can be as low as » -S 0.1 0.3 (see also figure 3). For the Néel phase the entropy per particle is equal to » log 2 0.077 1 9 due to the two-fold degeneracy corresponding to flipping of all spins. The relaxation time (usually below 1 s) increases with the increasing on-site interaction, U. For larger systems the relaxation times can be longer, due to possible formation of domains, as it is the case for coherent preparation strategies. A slightly different number of spin-up,  N , and spin-down,  N , atoms in the steady state is related to breaking of the SU(2) symmetry by the jump operatorsˆ( ) j ij 3 . As a result, the spin-down atom becomes more 'mobile'. Additionally, in the 3×3 lattice the number of sites is odd. Therefore, inherently in the steady state there is a spin-direction imbalance with >   N N . For a larger system, as well as for a system with even number of sites we would have »   N N . This, however, does not preclude the formation of the AF order. Figure 3 shows the steady-state properties: the entropy per particle and spin-structure factor as a function of the parameters. In the employed range of parameters these features turn out to be strongly dependent on g 2 and U, see figures 3(b) and (c), however, only weakly dependent on g 1 , see figure 3(a). Furthermore, while I AF grows substantially with increasing U and g 2 , for g 1 a saturation effect is observed and increasing the magnitude above g » 0.5 1 does not improve the efficiency of the scheme significantly. These observations can be qualitatively understood from figure 1(a). Namely, when the system is close to the AF phase most of the nearest-neighbor configurations are of the   ñ | ; type. The processes that drive the system away from the ordered state are related to coherent hopping from the   ñ | ; state to the   ñ | ; 0 state. In the large-U limit, the timescale of such processes is given by t U 4 2 . Therefore, increasing U reduces the contribution of the processes that destroy AF ordering. Increase of I AF with g 2 is expected, as the related dissipative processes drive the system directly into the AF ordered state. The saturation effect for g 1 can be due to the bidirectional character of the related dissipative processes. For sufficiently large g 1 , the value of the spin-structure factor is determined by an interplay between the dissipative processes related to g 2 and coherent hopping processes with a time scale governed by t U 4 2 . While the values of U required for an efficient preparation of the AF order are quite large, they are within experimental reach, e.g. = U t 180 was reported in [6]. Increasing U even further might lead to appearance of non-standard terms on top of the Fermi-Hubbard model [83].

Variational scheme
In order to describe the steady-state properties in the thermodynamic limit, we use a recently introduced variational principle [45,64]. In this method, one has to minimize a suitable variational norm 3 of the master equation (2) Here  is the dissipative part as given by equation (2). The variational method was originally formulated for bosonic systems where a local ansatz on the density matrix can be used. For fermionic systems such a procedure is not possible directly. However, the Jordan-Wigner transformation [76][77][78][79] can be used to map the fermionic creation and annihilation operators to spin operators (see appendix A for details). The resulting system of spin-1/2 particles is equivalent to a system of hard-core bosons, however, with the fermionic statistics included. Namely, the equivalence with the starting fermionic system is ensured by long-range interaction terms-the socalled Jordan-Wigner strings-manifesting the anticommutation relations of the original fermionic particles. To apply the variational method we first need to approximate the 'non-local' parts of the Jordan-Wigner strings and, thereby, transform the master equation (2) to the form with at most two-site interactions (see appendix B for details of this procedure).
We next consider variational states of the product-state type, r r r = =   A , while  AB gives the dissipative part with the jump operators acting on the sites A and B. The first two terms of equation (10) correspond to an exact treatment of a single bond, which already goes beyond the mean-field description, whereas the next ones describe interaction with the surrounding sites treated on the mean-field level, as visualized by the dashed lines in figure 4(a).
In figure 4(b) we compare the spin-structure factor obtained from our variational scheme and the MCWF method as a function of magnitude of the jump operators (which we set equal here without the loss of generality). According to the results of both methods, the system exhibits substantial ordering (e.g. with spinstructure factor larger than 0.5) when g g =  1 1 2 . Although the variational method contains terms going beyond mean field, its results for AF ordering do not depend on the value of U, as opposed to the exact approach. This happens due to restriction of the density matrices to the form of product states. Moreover, the variational method overestimates the AF ordering due to the absence of fluctuations in our variational manifold. Therefore, the employed approaches are complementary to each other, and both indicate the formation of an AF order of substantial magnitude in the steady state.

Conclusions
In this paper, we proposed a scheme for dissipative preparation of AF order in ultracold fermions trapped in an optical lattice. We demonstrated that by using a combination of two dissipative processes based on Ramanassisted hopping it is possible to engineer the dissipative dynamics in such a way that the AF phase emerges as its dark state. By using a combination of an exact and variational approaches, we observed the formation of a strong AF order on the timescales achievable in present-day experiments.
We note that the technique presented here can be readily implemented in the setups already used to search for the AF order [15], and thereby paves the way to an experimental realization of the AF phase in the Fermi-Hubbard model. While we exemplified the approach using the atomic level structure of 40 K [6], the method is general and can be also applied to other fermionic atoms currently available in laboratory, such as 6 Li [14], Er [84], Dy [85], Yb [86], and Cr [87]. After preparation of the AF phase with low entropy it should be possible to explore the phase diagram of the Hubbard model, including the pseudogap regime, by coherently removing a fraction of the atoms from the trap thereby introducing hole carriers into the system. Finally, extending these ideas to single-site addressable lattices as offered by the fermionic quantum gas microscopes [60][61][62]88], opens the door to the preparation of more sophisticated many-particles states. , where s -ˆl z evaluates to −1 when the lth spin is up (correspondingly, when there is a fermionic particle in the mode l) and to 1 in the opposite case.
For a one-dimensional system the function s ( ) k i, can be chosen simply as In two dimensions, however, there are a few possibilities to perform the Jordan-Wigner transformation [76][77][78][79] (see [76] for a review). Here, we consider anŃ N x y system with lattice site coordinates x y x y and use the following function to enumerate particles (with Î  m ) A.4 x y x y The chain of fermions and the related Jordan-Wigner string goes from left to right in the first row of the system, then in the second row goes to the left and forms a zig-zag (see figure 5(a)).
If the hopping takes place between two sites in the same row, the Jordan-Wigner string is local (as illustrated in figure 5(b)). For example, if = + ( ) , , we have  If the hopping takes place between two sites in different rows, the Jordan-Wigner string includes also sites between the i j , pair and the (left or right) edge of the system (as illustrated in figure 5(c)). Note that the number of these 'non-local' sites is always even due to our choice of the enumerating function s ( ) k i, . Consequently, at half-filling, i.e. with one particle per site on average, there is an even number of particles, N 2 1 , along the 'nonlocal' part of the Jordan-Wigner string. In such case, the string evaluates to the factor -= ( ) 1 1 N 2 1 . In our variational method we use this value as an approximation for the 'non-local' part of the string, whereas the local terms are retained. Such a procedure is similar to the 'mean-field' treatment of the strings applied, e.g. in [76,78,89]. As a result, the hopping between rows is expressed in the same way as in equations (B.2) and (B.3), which is visualized in figure 5(d).
For the jump operators we use the same approximation. Thereby, the master equation (2) acquires the form with at most two-site interaction terms (site here refers to the original fermionic site), for which the variational method can be readily applied. Note that in the MCWF calculations the 'non-local' part of the strings needs to be preserved.
Explicitly, the contribution from the hopping term to the norm of a single bond, r r = || || || || Tr

AB
A B  , on the example of the first term in equation (10) where we assumed , . The treatment of other terms in the norm is analogous and, similarly, results in the appearance of parity factors ( s - , which originate from the fermionic anticommutation rules.
Finally, let us note that the results of our variational method do not depend on the system size,Ń N x y (however, they make sense only for  N N , 4 x y ) and hence correspond to the thermodynamic limit.