Dissipative preparation of antiferromagnetic order in the Fermi-Hubbard model

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 temperature $T\sim U$. If the temperature is decreased further and reaches the so-called 'Néel temperature,' ${T}_{{\rm{N}}}\sim 4{t}^{2}/U$, 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}_{{\rm{N}}}\approx 1.42$ 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–28] to preparation of quantum phases in optical lattices.

In this work, we propose an efficient scheme for the preparation of AF order in the Fermi-Hubbard model, based on the ideas of dissipative state engineering, which have recently emerged in the context of many-particle systems [19, 22, 25, 26, 29–45] and have been implemented experimentally [46–52]. In such scenarios, a many-body state of interest (here: states exhibiting AF order) is prepared as a steady state of the quantum master equation governing the open system dynamics, as opposed to the ground state of the Hamiltonian. Such steady state can undergo quantum phase transitions to an ordered state of matter, which can be classified in close analogy to equilibrium systems [22, 45, 53–59].

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 ⁴⁰K, currently used in several laboratories [7, 8, 15, 60–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–68]. In our analysis of the dissipative Fermi-Hubbard model we use two complementary techniques: the Monte Carlo wave function (MCWF) [69–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.

We start with the Fermi-Hubbard Hamiltonian

$$\begin{eqnarray}&&\hat{H}=\displaystyle \sum _{{\bf{i}},{\bf{j}},\sigma }{t}_{{\bf{ij}}}{\hat{c}}_{{\bf{i}}\sigma }^{\dagger }{\hat{c}}_{{\bf{j}}\sigma }+U\displaystyle \sum _{{\bf{i}}}{\hat{n}}_{{\bf{i}}\uparrow }{\hat{n}}_{{\bf{i}}\downarrow },\end{eqnarray} \tag{ 1 }$$

which has been experimentally realized in a range of systems such as ⁶Li [14] and ⁴⁰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.

The dynamics of an open quantum system is governed by the master equation for the system's density matrix

$$\begin{eqnarray}&&\dot{\rho }=-{\rm{i}}[\hat{H},\rho ]+\mathop{\displaystyle \sum ^{\prime} }\limits_{{\bf{i}},{\bf{j}},\sigma ,\alpha }\left({\hat{j}}_{{\bf{ij}},\sigma }^{(\alpha )}\,\,\rho \,\,{\hat{j}}_{{\bf{ij}},\sigma }^{{(\alpha )}^{\dagger }}-\displaystyle \frac{1}{2}\{{\hat{j}}_{{\bf{ij}},\sigma }^{{(\alpha )}^{\dagger }}\,\,{\hat{j}}_{{\bf{ij}},\sigma }^{(\alpha )},\rho \}\right),\end{eqnarray} \tag{ 2 }$$

where we set ${\hslash }\equiv 1$ and the primed sum runs over nearest-neighbor sites. Since we start with a disordered sample, all possible nearest-neighbor configurations, including $| \uparrow ;\uparrow \rangle$, $| \downarrow ;\downarrow \rangle$, and $| \downarrow \uparrow ;0\rangle$ will be present. The jump operators, therefore, need to convert the latter into those with the local AF order, $| \uparrow ;\downarrow \rangle$.

We choose the jump operators to be as follows:

$$\begin{eqnarray}&&{\hat{j}}_{{\bf{ij}},\uparrow }^{(1)}=\sqrt{{\gamma }_{1}}\,{\hat{n}}_{{\bf{i}}\uparrow }(1-{\hat{n}}_{{\bf{j}}\downarrow }){\hat{c}}_{{\bf{i}}\downarrow }^{\dagger }{\hat{c}}_{{\bf{j}}\uparrow },\qquad {\hat{j}}_{{\bf{ij}},\uparrow }^{(2)}={({\hat{j}}_{{\bf{ij}},\uparrow }^{(1)})}^{\dagger },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\hat{j}}_{{\bf{ij}},\downarrow }^{(1)}=\sqrt{{\gamma }_{1}}\,{\hat{n}}_{{\bf{i}}\downarrow }(1-{\hat{n}}_{{\bf{j}}\uparrow }){\hat{c}}_{{\bf{i}}\uparrow }^{\dagger }{\hat{c}}_{{\bf{j}}\downarrow },\qquad {\hat{j}}_{{\bf{ij}},\downarrow }^{(2)}={({\hat{j}}_{{\bf{ij}},\downarrow }^{(1)})}^{\dagger },\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\hat{j}}_{{\bf{ij}}}^{(3)}=\sqrt{{\gamma }_{2}}\,(1-{\hat{n}}_{{\bf{i}}\uparrow }){\hat{n}}_{{\bf{j}}\uparrow }{\hat{c}}_{{\bf{i}}\downarrow }^{\dagger }{\hat{c}}_{{\bf{j}}\downarrow },\end{eqnarray} \tag{ 5 }$$

where ${\gamma }_{1}$, ${\gamma }_{2}$ are the dissipation rates. The resulting dissipative dynamics is visualized in figure 1(a): the jump operators ${\hat{j}}_{{\bf{ij}},\sigma }^{(1)}$ (with $\sigma =\uparrow ,\downarrow$) turn the configurations on neighboring sites from $| \sigma ;\sigma \rangle$ to $| \uparrow \downarrow ;0\rangle$, whereas ${\hat{j}}_{{\bf{ij}},\sigma }^{(2)}$ act in the opposite direction. These processes are labelled by the amplitude ${\gamma }_{1}$ in the figure. Finally, the jump operators ${\hat{j}}_{{\bf{ij}}}^{(3)}$ turn the configurations $| \uparrow \downarrow ;0\rangle$, into those with the AF order, $| \uparrow ;\downarrow \rangle$ (as labelled by ${\gamma }_{2}$). As a result, the dissipative dynamics drives the system towards the AF phase.

Zoom In Zoom Out Reset image size

We note that the jump operators ${\hat{j}}_{{\bf{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 $| \uparrow \rangle$-state atom away from the double-occupancy configuration.

As we discuss in the following section, such choice of the jump operators is straighforward to realize in experiment using incoherent laser scattering.

Simulating the Fermi-Hubbard model requires mapping the two spin states, $| \uparrow \rangle$ and $| \downarrow \rangle$, 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 ${}^{40}{\rm{K}}$ [73], with the levels $| \uparrow \rangle \equiv \left|F=\tfrac{7}{2};{m}_{F}=-\tfrac{7}{2}\right.\rangle$ and $| \downarrow \rangle \equiv \left|\tfrac{9}{2};-\tfrac{7}{2}\right.\rangle$. Furthermore, in order to realize the dissipative part of the dynamics, we introduce an auxiliary state, $| X\rangle \equiv \left|\tfrac{9}{2};-\tfrac{9}{2}\right.\rangle$, belonging to the ${}^{2}{S}_{1/2}$ manifold, as well as an electronically excited ${}^{2}{P}_{3/2}$ state, $| e\rangle \equiv \left|\tfrac{11}{2};-\tfrac{9}{2}\right.\rangle$.

The first-stage jump operators, ${\hat{j}}_{{\bf{ij}},\sigma }^{(1)}$ and ${\hat{j}}_{{\bf{ij}},\sigma }^{(2)}$, can be implemented using Raman-assisted hopping, as illustrated in figure 1(b). In such a process, transitions between two quantum states are induced using two laser beams, which are detuned from some excited state (here the ${}^{2}{P}_{3/2}$ state). For example, the Raman beams ${{\rm{\Omega }}}_{r1}$ and ${{\rm{\Omega }}}_{r3}$ realize the Raman-assisted hopping between the $| \uparrow ;\uparrow \rangle$ and $| 0;\uparrow \downarrow \rangle$ states. If the Raman beams, ${{\rm{\Omega }}}_{r1,\ldots ,r3}$, 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 $| \uparrow ;\downarrow \rangle$ state at this stage. Otherwise, the jump operators would also lead from the $| \uparrow ;\downarrow \rangle$ state back to the $| \uparrow ;\uparrow \rangle$ and $| \downarrow ;\downarrow \rangle$ 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, ${\hat{j}}_{{\bf{ij}}}^{(3)}$, in a one-directional fashion, we use Raman-assisted hopping (lasers ${{\rm{\Omega }}}_{r4}$ and ${{\rm{\Omega }}}_{r5}$) from the $| \uparrow \downarrow ;0\rangle$ state to the $| \uparrow ;X\rangle$ configuration with an auxiliary $| X\rangle$ state. This $| X\rangle$ state is then pumped (laser Ω) to an excited $| e\rangle$ state, which can decay to the $| \downarrow \rangle$ state completing the process, as illustrated in figure 1(c). The $| e\rangle$ state cannot decay to the $| \uparrow \rangle$ state because of selection rules on the F quantum number. The Zeeman splitting, ${{\rm{\Delta }}}_{X\downarrow }$, and the energy, $U+{{\rm{\Delta }}}_{X\downarrow }$, 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 $| \uparrow ;X\rangle$ and $| \uparrow X;0\rangle$ states as the final states of the Raman-assisted hopping process we need nonzero on-site interaction between the $| \uparrow \rangle$ and $| X\rangle$ states. The typical values of the background scattering lengths, $a=105\,{a}_{0}$, in units of the Bohr radius a₀ [73, 74], should be sufficient for this purpose. Spontaneous emission from the $| e\rangle$ state ensures that the resulting dissipative processes are unidirectional and take place from the $| \uparrow \downarrow ;0\rangle$ to the $| \uparrow ;\downarrow \rangle$ state with AF ordering.

In total our dissipative preparation scheme requires six lasers (${{\rm{\Omega }}}_{r1}$, ..., ${{\rm{\Omega }}}_{r5}$, and Ω) in order to induce the dissipative transitions, in addition to the lasers used to trap the atoms and prepare them in the $| \uparrow \rangle$ and $| \downarrow \rangle$ states. Due to the s-wave character of the ${}^{2}{S}_{1/2}$ manifold, the optical lattice potential for the $| \uparrow \rangle$, $| \downarrow \rangle$, and $| X\rangle$ states will be the same, provided the optical-trapping lasers are sufficiently far-detuned from the ${}^{2}{P}_{3/2}$ states (see the discussion in [75], section 2.B).

The values of ${\gamma }_{1}$ and ${\gamma }_{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, ${\gamma }_{1}$, and ${\gamma }_{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 $| \uparrow \rangle$, $| \downarrow \rangle$, $| X\rangle$, and $| e\rangle$. In particular, it should be possible to use the states of ⁴⁰K, for which the Feschbach resonance is already known, i.e. $| \downarrow \rangle =\left|\tfrac{9}{2};-\tfrac{9}{2}\right.\rangle$, $| \uparrow \rangle =\left|\tfrac{9}{2};-\tfrac{7}{2}\right.\rangle$ or $| \downarrow \rangle =\left|\tfrac{9}{2};-\tfrac{9}{2}\right.\rangle$, $| \uparrow \rangle =\left|\tfrac{9}{2};-\tfrac{5}{2}\right.\rangle$. In the first case one could use $| X\rangle =\left|\tfrac{7}{2},-\tfrac{7}{2}\right.\rangle$ and $| e\rangle =\left|\tfrac{11}{2},-\tfrac{11}{2}\right.\rangle$, which requires a two-photon process from $| X\rangle$ to $| e\rangle$. This could be realized with Raman beams detuned from the ${}^{2}{D}_{5/2}$ state or directly as a two-photon excitation. In the second case one could use $| X\rangle =\left|\tfrac{9}{2},-\tfrac{7}{2}\right.\rangle$ and $| e\rangle =\left|\tfrac{9}{2},-\tfrac{9}{2}\right.\rangle$ ($| e\rangle$ is from the ${}^{2}{P}_{3/2}$ manifold), which would require resolving the $| X\rangle \leftrightarrow | e\rangle$ transition from $| \uparrow \rangle \leftrightarrow \left|\tfrac{9}{2},-\tfrac{7}{2}\right.\rangle$.

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–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.

4.1. 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}_{{\bf{ij}}}\equiv -t$ (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 ${\gamma }_{1}$, ${\gamma }_{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 (${\hat{c}}_{{\bf{i}}\sigma }^{\dagger }{\hat{c}}_{{\bf{j}}\overline{\sigma }}$), 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 $\tau \times t\approx 50\mbox{--}100$, which corresponds to $\tau =1\mbox{--}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\approx 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.

Zoom In Zoom Out Reset image size

To quantify the AF ordering of the system, we evaluate the spin-structure factor, as defined by

$$\begin{eqnarray}&&I({\bf{Q}})\equiv \displaystyle \frac{4}{{N}^{2}}\displaystyle \sum _{{\bf{i}},{\bf{j}}}{{\rm{e}}}^{{\rm{i}}({{\bf{R}}}_{{\rm{i}}}-{{\bf{R}}}_{{\rm{j}}}){\bf{Q}}}\langle {\sigma }_{{\bf{i}}}^{z}{\sigma }_{{\bf{j}}}^{z}\rangle .\end{eqnarray} \tag{ 6 }$$

Here, ${\bf{Q}}$ is the ordering vector, which for the case of the AF phase is equal to ${{\bf{Q}}}_{\mathrm{AF}}=[\pi ,\pi ]$, ${{\bf{R}}}_{{\rm{i}}}$ and ${{\bf{R}}}_{{\rm{j}}}$ are the lattice site vectors, N is the number of particles (N = 9 in our system), and ${\sigma }_{{\bf{i}}({\bf{j}})}^{z}=\pm \tfrac{1}{2}$ is the z component of the particles' spin. In the case of the AF ordering, the relevant spin-structure factor is given by ${I}_{\mathrm{AF}}\equiv I({{\bf{Q}}}_{\mathrm{AF}})$. For the steady-state it can be as large as ${I}_{\mathrm{AF}}\approx 0.8$ (see also figure 3), quite close to the fully polarized Néel state, for which ${I}_{\mathrm{AF}}=1$. Note that the Néel state is an 'ideal' AF state, to which the dissipative processes would drive the system in the limit of ${\gamma }_{1},{\gamma }_{2},U\gg t$. 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.

Zoom In Zoom Out Reset image size

Figure 2(b) illustrates the decrease of the system's entropy per particle, $S\equiv {S}_{\mathrm{tot}}/N=-1/N\,\mathrm{Tr}\rho \mathrm{log}\rho$, with time. Due to the large size of the system's density matrix ρ, in our numerical calculations we use the equivalent formula, ${S}_{\mathrm{tot}}=-\mathrm{Tr}A\mathrm{log}A$. Here ${A}_{{ij}}=\langle {\psi }_{i}| {\psi }_{j}\rangle$ 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\approx 0.1-0.3$ (see also figure 3). For the Néel phase the entropy per particle is equal to $\tfrac{1}{9}\mathrm{log}2\,\approx \,0.077$ 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}_{\uparrow }$, and spin-down, ${N}_{\downarrow }$, atoms in the steady state is related to breaking of the SU(2) symmetry by the jump operators ${\hat{j}}_{{\bf{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}_{\uparrow }\gt {N}_{\downarrow }$. For a larger system, as well as for a system with even number of sites we would have ${N}_{\uparrow }\approx {N}_{\downarrow }$. 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 ${\gamma }_{2}$ and U, see figures 3(b) and (c), however, only weakly dependent on ${\gamma }_{1}$, see figure 3(a). Furthermore, while ${I}_{\mathrm{AF}}$ grows substantially with increasing U and ${\gamma }_{2}$, for ${\gamma }_{1}$ a saturation effect is observed and increasing the magnitude above ${\gamma }_{1}\approx 0.5$ 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 $| \uparrow ;\downarrow \rangle$ type. The processes that drive the system away from the ordered state are related to coherent hopping from the $| \uparrow ;\downarrow \rangle$ state to the $| \uparrow \downarrow ;0\rangle$ state. In the large-U limit, the timescale of such processes is given by $4{t}^{2}/U$. Therefore, increasing U reduces the contribution of the processes that destroy AF ordering. Increase of ${I}_{\mathrm{AF}}$ with ${\gamma }_{2}$ is expected, as the related dissipative processes drive the system directly into the AF ordered state. The saturation effect for ${\gamma }_{1}$ can be due to the bidirectional character of the related dissipative processes. For sufficiently large ${\gamma }_{1}$, the value of the spin-structure factor is determined by an interplay between the dissipative processes related to ${\gamma }_{2}$ and coherent hopping processes with a time scale governed by $4{t}^{2}/U$.

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].

4.2. 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³ of the master equation (2)

$$\begin{eqnarray}&&| | \dot{\rho }| | =| | -{\rm{i}}[H,\rho ]+{ \mathcal D }(\rho )| | .\end{eqnarray} \tag{ 7 }$$

Here ${ \mathcal D }$ 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–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 so-called 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, $\rho ={\rho }_{p}={\prod }_{{\bf{i}}}{\rho }_{{\bf{i}}}$ and minimize the upper bound of the norm

$$\begin{eqnarray}&&| | \dot{\rho }| | \leqslant \displaystyle \sum _{\langle {\bf{ij}}\rangle }\mathrm{Tr}| {\dot{\rho }}_{{\bf{ij}}}| \to \min ,\end{eqnarray} \tag{ 8 }$$

where the reduced two-site operators are defined as ${\dot{\rho }}_{{\bf{ij}}}={\mathrm{Tr}}_{\not\hspace{-2pt}{{\bf{i}}}\not\hspace{-2pt}{{\bf{j}}}}\,\dot{\rho }$. It is sufficient to minimize the norm $| | {\dot{\rho }}_{{\bf{ij}}}| |$ of a single bond, which for the case of an AF order can be expressed as

$$\begin{eqnarray}&&| | {\mathop{\dot{\rho }}\limits^{}}_{{\bf{ij}}}| | \equiv | | {\mathop{\dot{\rho }}\limits^{}}_{{AB}}| | =\mathrm{Tr}| {\mathop{\dot{\rho }}\limits^{}}_{{AB}}| ,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\mathop{\dot{\rho }}\limits^{}}_{{AB}} & = & -{\rm{i}}[{H}_{{AB}},{\rho }_{{AB}}]+{{ \mathcal D }}_{{AB}}({\rho }_{{AB}})\\ & & +\displaystyle \sum _{{A}^{\prime} }{\mathrm{Tr}}_{{A}^{\prime} }\{-{\rm{i}}[{H}_{{BA}^{\prime} },{\rho }_{{ABA}^{\prime} }]+{{ \mathcal D }}_{{BA}^{\prime} }({\rho }_{{ABA}^{\prime} })\}\\ & & +\displaystyle \sum _{{B}^{\prime} }{\mathrm{Tr}}_{{B}^{\prime} }\{-{\rm{i}}[{H}_{{B}^{\prime} {A}},{\rho }_{{B}^{\prime} {AB}}]+{{ \mathcal D }}_{{B}^{\prime} {A}}({\rho }_{{B}^{\prime} {AB}})\}.\end{array}\end{eqnarray} \tag{ 10 }$$

Here A and B label the two sublattices and, e.g., ${\rho }_{{AB}}\equiv {\rho }_{{A}}\otimes {\rho }_{{B}}$, ${\rho }_{{ABA}^{\prime} }\equiv {\rho }_{{A}}\otimes {\rho }_{{B}}\otimes {\rho }_{{A}^{\prime} }$, while ${{ \mathcal D }}_{{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).

Zoom In Zoom Out Reset image size

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 spin-structure factor larger than 0.5) when ${\gamma }_{1}={\gamma }_{2}\gtrsim 1$. 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.

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 Raman-assisted 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 ⁴⁰K [6], the method is general and can be also applied to other fermionic atoms currently available in laboratory, such as ⁶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–62, 88], opens the door to the preparation of more sophisticated many-particles states.

We acknowledge stimulating discussions with Ken Brown, Tommaso Calarco, Andrew Daley, Suzanne McEndoo, Tobias Osborne, Cindy Regal, Luis Santos, Michał Tomza, and Martin Zwierlein. The work was supported by the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement no. [291734], by the Volkswagen Foundation, and by DFG within SFB 1227 (DQ-mat).

The Jordan–Wigner transformation between fermionic operators ${\hat{c}}_{{\bf{i}}\sigma }^{\dagger }$, ${\hat{c}}_{{\bf{i}}\sigma }$ and spin operators (${\hat{\sigma }}^{+},{\hat{\sigma }}^{-},{\hat{\sigma }}^{z}$) for spin-1/2 particles is defined as

$$\begin{eqnarray}&&{\hat{c}}_{{\bf{i}}\sigma }^{\dagger }={\hat{\sigma }}_{k({\bf{i}},\sigma )}^{+}{{\rm{e}}}^{+{\rm{i}}\pi {\displaystyle \sum }_{l\lt k({\bf{i}},\sigma )}{\hat{n}}_{l}},\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{\hat{c}}_{{\bf{i}}\sigma }={\hat{\sigma }}_{k({\bf{i}},\sigma )}^{-}{{\rm{e}}}^{-{\rm{i}}\pi {\displaystyle \sum }_{l\lt k({\bf{i}},\sigma )}{\hat{n}}_{l}},\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{\hat{c}}_{{\bf{i}}\sigma }^{\dagger }{\hat{c}}_{{\bf{i}}\sigma }=\displaystyle \frac{1}{2}({\hat{\sigma }}_{k({\bf{i}},\sigma )}^{z}+1)\equiv {\hat{n}}_{k({\bf{i}},\sigma )},\end{eqnarray} \tag{ A.3 }$$

where the function $k({\bf{i}},\sigma )$ enumerates the fermionic particles for a given lattice site, ${\bf{i}}=({i}_{x},{i}_{y})$, and particle's spin, $\sigma =\uparrow ,\downarrow$. In this manner, fermions are positioned on a chain, where the position along the chain is given by $k({\bf{i}},\sigma )$. The factor ${{\rm{e}}}^{\pm {\rm{i}}\pi {\sum }_{l\lt k({\bf{i}},\sigma )}{\hat{n}}_{l}}$, where the summation runs over all fermions 'before' the one at site ${\bf{i}}$ with spin $\sigma$, evaluates to $+1(-1)$ for even (odd) number of fermions 'before' the given one. Thereby, the anticommutation rules for the fermionic operators ${\hat{c}}_{{\bf{i}}\sigma }^{\dagger }$, ${\hat{c}}_{{\bf{i}}\sigma }$ are fulfilled. We can also use the relation ${{\rm{e}}}^{\pm {\rm{i}}\pi {\sum }_{l\lt k({\bf{i}},\sigma )}{\hat{n}}_{l}}={{\rm{\Pi }}}_{l\lt k({\bf{i}},\sigma )}(-{\hat{\sigma }}_{l}^{z})$, where $-{\hat{\sigma }}_{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 $k({\bf{i}},\sigma )$ can be chosen simply as $k({\bf{i}},\sigma )=2| {\bf{i}}| +{\delta }_{\sigma ,\downarrow }$. In two dimensions, however, there are a few possibilities to perform the Jordan–Wigner transformation [76–79] (see [76] for a review). Here, we consider an ${N}_{x}\times {N}_{y}$ system with lattice site coordinates ${i}_{x(y)}\in \{1,2,...,{N}_{x(y)}\}$ and use the following function to enumerate particles (with $m\in {\mathbb{N}}$)

$$\begin{eqnarray}k({\bf{i}},\sigma )=\left\{\begin{array}{cc}2[{N}_{x}({i}_{y}-1)+{i}_{x}-1]+{\delta }_{\sigma ,\downarrow }, & {i}_{y}=2m+1,\\ 2[{N}_{x}({i}_{y}-1)+{N}_{x}-{i}_{x}]+{\delta }_{\sigma ,\downarrow }, & {i}_{y}=2m.\end{array}\right.\end{eqnarray} \tag{ A.4 }$$

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)).

Zoom In Zoom Out Reset image size

By performing the Jordan–Wigner transformation as described in appendix A, we can transform all terms of the master equation (2), which we exemplify here by considering the hopping, ${\hat{c}}_{{\bf{i}}\sigma }^{\dagger }{\hat{c}}_{{\bf{j}}\sigma }$. Using the transformation (A.1)–(A.3) the hopping is expressed as

$$\begin{eqnarray}&&{\hat{c}}_{{\bf{i}}\sigma }^{\dagger }{\hat{c}}_{{\bf{j}}\sigma }={\hat{\sigma }}_{k({\bf{i}},\sigma )}^{+}[{{\rm{\Pi }}}_{k({\bf{i}},\sigma )\lt l\lt k({\bf{j}},\sigma )}(-{\hat{\sigma }}_{l}^{z})]{\hat{\sigma }}_{k({\bf{j}},\sigma )}^{-}.\end{eqnarray} \tag{ B.1 }$$

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 ${\bf{j}}={\bf{i}}+(1,0)$ and $k({\bf{i}},\sigma )\lt k({\bf{j}},\sigma )$, we have

$$\begin{eqnarray}&&{\hat{c}}_{{\bf{i}}\uparrow }^{\dagger }{\hat{c}}_{{\bf{j}}\uparrow }={\hat{\sigma }}_{k({\bf{i}},\uparrow )}^{+}[-{\hat{\sigma }}_{k({\bf{i}},\downarrow )}^{z}]{\hat{\sigma }}_{k({\bf{j}},\uparrow )}^{-},\end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray}&&{\hat{c}}_{{\bf{i}}\downarrow }^{\dagger }{\hat{c}}_{{\bf{j}}\downarrow }={\hat{\sigma }}_{k({\bf{i}},\downarrow )}^{+}[-{\hat{\sigma }}_{k({\bf{j}},\uparrow )}^{z}]{\hat{\sigma }}_{k({\bf{j}},\downarrow )}^{-}.\end{eqnarray} \tag{ B.3 }$$

If the hopping takes place between two sites in different rows, the Jordan–Wigner string includes also sites between the ${\bf{i}},{\bf{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 $k({\bf{i}},\sigma )$. Consequently, at half-filling, i.e. with one particle per site on average, there is an even number of particles, $2{N}_{1}$, along the 'non-local' part of the Jordan–Wigner string. In such case, the string evaluates to the factor ${(-1)}^{2{N}_{1}}=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, $| | {\mathop{\dot{\rho }}\limits^{}}_{{AB}}| | =| | {\mathrm{Tr}}_{\not\hspace{-2pt}{A}\not\hspace{-2pt}{B}}\,\mathop{\dot{\rho }}\limits^{}| |$, on the example of the first term in equation (10), is given by

$$\begin{eqnarray}\begin{array}{rcl}-{\rm{i}}[{H}_{{AB}},{\rho }_{{AB}}] & = & -{\rm{i}}t[{\hat{\sigma }}_{k({A},\uparrow )}^{+}{\hat{\sigma }}_{k({A},\downarrow )}^{z}{\hat{\sigma }}_{k({B},\uparrow )}^{-}\\ & & +{\hat{\sigma }}_{k({A},\downarrow )}^{+}{\hat{\sigma }}_{k({B},\uparrow )}^{z}{\hat{\sigma }}_{k({B},\downarrow )}^{-}+{\rm{h.c.}},{\rho }_{{AB}}],\end{array}\end{eqnarray} \tag{ B.4 }$$

where we assumed $k({A},\sigma )\lt k({B},\sigma )$. The treatment of other terms in the norm is analogous and, similarly, results in the appearance of parity factors ($-{\hat{\sigma }}_{k({A},\sigma )}^{z}$ and $-{\hat{\sigma }}_{k({B},\sigma )}^{z}$), 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}\times {N}_{y}$ (however, they make sense only for ${N}_{x},{N}_{y}\geqslant 4$) and hence correspond to the thermodynamic limit.