Driven-dissipative many-body pairing states for cold fermionic atoms in an optical lattice

We discuss the preparation of many-body states of cold fermionic atoms in an optical lattice via controlled dissipative processes induced by coupling the system to a reservoir. Based on a mechanism combining Pauli blocking and phase locking between adjacent sites, we construct complete sets of jump operators describing coupling to a reservoir that leads to dissipative preparation of pairing states for fermions with various symmetries in the absence of direct inter-particle interactions. We discuss the uniqueness of these states, and demonstrate it with small-scale numerical simulations. In the late time dissipative dynamics, we identify a"dissipative gap"that persists in the thermodynamic limit. This gap implies exponential convergence of all many-body observables to their steady state values. We then investigate how these pairing states can be used as a starting point for the preparation of the ground state of Fermi-Hubbard Hamiltonian via an adiabatic state preparation process also involving the parent Hamiltonian of the pairing state. We also provide a proof-of-principle example for implementing these dissipative processes and the parent Hamiltonians of the pairing states, based on Yb171 atoms in optical lattice potentials.


Introduction
Quantum simulation using cold atoms in an optical lattice typically requires cooling to low temperatures to see interesting quantum phases with strong correlations [1,2,3,4,5]. An important example in this context is the quantum simulation of the 2-D Fermi-Hubbard Model (FHM) using cold fermionic atoms in an optical lattice [6]. As a crucial first step toward this goal, various experimental groups have been successful in realising the 3-D Fermi-Hubbard model in such a system [7,8,9,10]. But cooling of the system below the critical temperature, and thus into the phases of interest, turns out to be difficult with the conventional cooling schemes (T /T F ∼ 0.25 [9,10], where T F is the Fermi temperature). On the other hand, the atomic physics of these systems opens the possibility for a different approach to the production of interesting many-body states, specifically to dissipatively drive the system into steady states with the desired coherence and symmetry properties by careful engineering of a reservoir [11,12,13,14,15,16,17,18,19]. In this paper, we will discuss schemes by which fermions in a 2-D lattice potential can be dissipatively driven into pairing states with non-trivial correlations even in the absence of attractive interactions, extending our previous work [19]. These many-body pairing states can then also be used as a low-entropy starting point for efficient adiabatic passages, through which the true many-body ground state of many-body models such as the Fermi-Hubbard model might be reached [3,4,5].
For Markovian dissipative processes, the coupling with the reservoir can be modeled by a set of dissipative quantum jump operators, which, when chosen appropriately, drive the system towards a pure steady state with zero entropy starting from any given initial state. Similar ideas have been applied recently for the dissipative preparation of manybody states in optical lattices, e.g. a Bose-Einstein condensate for bosons [12], an ηcondensate for fermions [13], and most recently, d-wave pairing states in optical lattices [19] as well as topological phases [20]. Here, we will focus on the dissipative generation of pairing states for fermions in an optical lattice. Starting from the general principles, we will show in detail how relative phases between the atoms can be imposed by engineering the jump operators, which in turn allows us to prepare many-body pairing states with specific spatial symmetries, e.g. p-wave and d-wave pairing symmetries. Note that in contrast to the equilibrium states of typical non-dissipative processes, the dissipative dynamics is described by a master equation, and the final state that we aim to prepare appears as the steady state of the dynamical process, providing a targeted many-body cooling protocol. As a conceptually important result, we observe that the imaginary spectrum of the effective Hamiltonian of the master equation has a gap which persists in the thermodynamic limit and we refer to it as a "dissipative gap", due to the formal analogy to the pairing gap of conventional BCS paired states [21]. Physically, it leads to exponential convergence of many-body observables to their steady state values in the late time dissipative dynamics. This is a unique feature of the dissipative preparation of paired states with fermions.
We then propose an adiabatic process to connect the steady state of the dissipative process with the ground state of a Hamiltonian with similar symmetry properties. In the case of ideal adiabaticity, the system remains in a pure state and evolves into the ground state of the Hamiltonian without going through any Landau-Zener crossings. We consider this in the context of the Fermi-Hubbard model, where a dissipative process can be used to create dissipatively bound fermion pairs which have d-wave symmetry. As we will see later, this BCS-type mean-field d-wave state lacks the strong correlations associated with the repulsive Fermi-Hubbard model, as the doubleoccupancy is not projected out by the dissipative dynamics. We will show that strong correlations can be built up through proper adiabatic passage in the sense that the steady state can be connected with the ground state of the Fermi-Hubbard model efficiently. The central question in the study of the 2-D FHM is whether the ground state exhibits d-wave superconductivity/superfluidity away from half-filling, and thus captures the universal properties shared by the high-Tc superconducting materials [22,23,24,25,26,27,28,29,30]. If this is indeed the case, the coherence and the d-wave symmetry should be conserved during the adiabatic process.
Finally, we note that a key physical ingredient for the jump operators for pairing states of fermions is the Pauli exclusion principle [31]. Based on this understanding, we propose to implement the jump operators stroboscopically using alkaline-earth-like atoms. The existence of long-lived meta-stable states and rich level structures in these atoms provides us with the freedom to engineer various dissipative processes.
The paper is organised as follows: In Sec. II, we discuss the general formalism for jump operator engineering, and derive both the fixed-number and the fixed-phase jump operators for pairing states with different symmetries; in Sec. III, we study the uniqueness of the steady state under these jump operators for various cases both from the symmetry perspective and with small scale numerical calculations; we then derive a mean-field theory in Sec. IV for the master equation describing the dissipative dynamics, where a "dissipative gap" for the pairing states -a minimal damping rate for the manybody relaxation -is shown to emerge; in Sec. V, having established the dissipative preparation of pure steady states, we illustrate how we may start from these initial states to prepare the ground state of Hamiltonians with similar symmetries, for example the FHM, via an efficient adiabatic process; we discuss the implementation of the jump operators using alkaline-earth-like atoms in Sec. VI; finally we summarise and discuss other possible applications of our dissipative state preparation setting.

General principles of reservoir engineering
Closed systems, where the energy and particle number are conserved, are modeled by a Hamiltonian, which determines the ground state and the dynamics via the Schrödinger equation. In this context, engineering states by realisation of a particular Hamiltonian, for which the desired state is the ground state is often discussed. It is then natural to ask similar questions for open quantum systems, where dissipative processes interrupt coherent evolution. In particular, one can consider engineering the dissipative processes so that their back-actions project the system into the desired subspace of the complete Hilbert space. For Markovian dissipative processes, as appropriate for the systems considered below, the dynamics of the density matrix for an open system is described by a master equation, where the non-Hermitian effective Hamiltonian is given by Here, {j } are non-Hermitian Lindblad operators reflecting the system-bath coupling with rate κ [32]. The Hamiltonian H generates unitary evolution, and describes nondissipative processes of the system. Although H does not have to vanish in general, we will assume H = 0 throughout the state preparation discussion, i.e. the final state is prepared via purely dissipative processes. For fermions loaded into an optical lattice, this can be achieved by increasing the lattice depth to freeze out the kinetic motion, while tuning the inter-particle interaction to zero, e.g. via a Feshbach resonance. Note that by considering a purely dissipative process, we can avoid competition between the Hamiltonian and the dissipative dynamics, which is present when the dark state is not an exact eignestate of the Hamiltonian. In particular, the pairing states that we aim to prepare are not exact eigenstates of the FHM in general.
In the quantum trajectory picture, the system wavefunction |ψ(t) of a given trajectory evolves according to the non-Hermitian Hamiltonian |ψ(t) ∝ e −iH eff t |ψ(0) , and is punctuated by the quantum jump |ψ(t) → j |ψ(t) with rate κ j |ψ(t) 2 . The time-dependent density matrix is then determined by ρ(t) = |ψ(t) ψ(t)| stoch [33], where the average runs over all trajectories. In this picture, we see that for any state satisfying j |BCS N = 0 ∀ , the quantum jumps will never project it to other states. These states are therefore "dark states" of the jump operators, and are necessarily steady states of the master equation evolution [12,14]. If no other stationary solutions exist, the system will be driven to this state by the dissipative dynamics regardless of its initial conditions [11]. Similar ideas have been exploited for the dynamical preparation of a BEC for cold bosons loaded into an optical lattice, where the jump operators give rise to quasi-local phase-locking mechanism, which eventually leads to the condensation of the bosons in the lattice [10]. Here, we will focus on the preparation of pairing states for fermions in an optical lattice. For the fermionic case here, in addition to the phase-locking mechanism as in the case of bosons, the physical foundation of the jump operators is the Pauli blocking (see Fig. 1), i.e. for fermions the spontaneous emission to an already occupied state is blocked [34]. This gives rise to a novel non-equilibrium pairing mechanism for fermions that does not require attractive conservative forces. Figure 1. Illustration of the Pauli blocking in the optical pumping process of a three level system. (a) The spontaneous decay from the excited state |e to the ground state |g 2 is blocked due to Pauli exclusion, leaving the system unchanged; (b) The decay channel is not blocked, and the population in |g 1 is transferred to the state |g 2 .
We are primarily interested in states with a homogeneous product of N identical fermion pairs on a two-dimensional (2-D) square lattice: Here, c † a creates a fermion in mode a, where a = (σ, i) labels spin σ and position i on the 2-D lattice. Typically, when we consider delocalised states, the sum extends over the whole lattice in the position space index. For large systems in the thermodynamic limit, we may adopt the grand canonical ensemble, and the state above becomes the BCS-type coherent state for paired fermions: where we have Fourier transformed Eq. (3) into momentum space, with A = {σ, k} labeling spin σ and momentum k, N being the normalization factor. The pair wavefunction in momentum space is f A,B = i,j η a,b e ikr i +ik r j , where the spin degrees of freedom are not changed. The symmetry of the pairing state is encoded in the form of the pair wavefunction η a,b (f A,B ). In general, the jump operator that drives the system into the pairing state of Eq. (4) can involve multi-particle processes. Single-particle jump operators (involving at most one annihilation mode operator), which act directly on one particle at a time, are generally easier to implement experimentally than two-particle jump operators (involving two annihilation mode operators). However, the two-particle jump operators are the most intuitive, as the dissipative dynamics can be viewed as local phaselocking between pairs of fermions. We will therefore first give a brief description on the derivation of two-particle jump operators for the pairing states, which will provide important clues as to how to proceed with the more practical single-particle operators. In both cases, our focus will be on the d-wave pairing state, though the procedures can be extended to other symmetries as well.
The fixed-number d-wave pairing state is given by the symmetric superposition of d-wave pairs, which are spin-singlet fermion pairs on the bonds in x and y direction whose relative phase changes sign under a π/4 rotation x ↔ y, with e x , e y the unit lattice vectors long the x and y direction, respectively. Following Eq. (4), it is easy to write down the general form of a fixed-phase d-wave pairing state: where α = e iθ |α| is a complex number carrying the phase θ of the pairing state. In the following, we will first focus on the fixed-number state and will come back to the fixed-phase state later. Assuming translational invariance, which is the case for an infinitely large homogeneous lattice or for a finite homogeneous lattice with periodic boundary conditions, we may rewrite the pair creation operators: These operators have the advantage that they are already factorised in the spin degrees of freedom, which allows us to write the pairing jump operators in a similar fashion as well. For convenience, we adopt a shorthand convention and write d † j = ν ρ ν c † j+eν ,↑ c j,↓ , where ρ ±x = 1, ρ ±y = −1. We also find it useful to define the singlet pairing state in one dimension (1-D): where c † i{↑,↓} is the single fermion annihilation operator on the ith site with a certain spin. The 1-D state η † N |vac captures the off-site singlet pairing feature of the d-wave state. Similar to the d-wave case in 2-D, we may use translational symmetry to define the 1-D singlet pairing operator η † i : so that η † = j η † j . In the following, we proceed by devising the jump operators for the 1-D singlet pairing state first before considering the 2-D scenarios where additional spatial phase locking is needed.

Two-particle jump operators
We are interested in jump operators that conserve particle number. This implies that they must carry total charge 0 (or with total global phase 0). Thus, we can write it in normal order as a product of a pure creation and a pure annihilation part. The most general form of the two-particle jump operator thus reads Here we also impose the requirement of quasi-locality, i.e. the functions χ ab (i), ξ ab (i) shall be non-zero only in a small vicinity of site i in position space. Thus the jump operator J i is centered around site i.
To uniquely drive the system into the desired pairing state, it is necessary that the state |Ψ be a dark state of a set of jump operators For a given operator ξ i , we can work out its commutation relation with the creation operator of the pairs: While A carries charge 0 and is composed of a constant plus a normal-ordered second order term, B is a superposition of pair creation operators and carries charge 2, which implies that [B, η † ] = 0 (cf. Appendix A).
With these relations, we find that the commutator of ξ i with the homogeneous product η † N is characterised by the commutators A and B only: The first term on the right-hand side of the equation above gives zero when acting on the vacuum. Similarly, the normal-ordered part of A yields zero on the vacuum. Therefore, in order to satisfy the dark state condition, we need to find a set of quasi-local bilinear operators ξ i and χ † i which for a given η † uniquely solve the two equations For the 1-D case, after some derivation following the arguments leading to Eq. (15) (see Appendix A for details), we find a two-particle jump operator of the form: where C † is an arbitrary superposition of single-fermion creation operators with spin-up. Note that as a spin flip operation only changes the overall sign, jump operators similar to J i but with spins flipped also have the singlet pairing state Eq. (9) as a dark state. This is also true for the jump operators of the pairing states that we consider in this work, as similar symmetries in the spin degrees of freedom hold for all the pairing states that we will consider. For the 2-D case, we find the jump operator (see Appendix A for details): where ρ ±x = 1, ρ ±y = −1. This operator has an interesting structure. It may be seen as a conditional dissipative process. The annihilation of a fermion sitting at site i must only take place if a superposition of fermions located at sites i + e x , i + e y is coherently transferred to a superposition on the four sites {i ± e x , i ± e y } centered around site i. The two-particle jump operators above have the desired pairing states as the dark state. However, numerical simulation shows that in most cases they do not guarantee a unique dark state. Furthermore, these jump operators involve correlated dissipation of two particles and are therefore difficult to implement experimentally. Nevertheless the construction scheme above provides clues for the design of jump operators that are easier to implement and with a unique dark state.

Single-particle jump operators
In this subsection, we will show that counter-intuitively, it is possible to design a set of jump operators for which only single particle operations are required. The singleparticle jump operators dissipatively drive fermions into pairs as well as phase lock the pairs into the desired symmetry. More importantly, we will give arguments later that the dark state of these single-particle jump operators should be unique, and they are easier to implement than the two-particle jump operators. Below, we will describe two ways in which the appropriate single-particle jump operators can be derived.
Firstly, consistent with the discussion in the previous section, in the case of a singleparticle jump operator, the annihilation part ξ i contains a single annihilation operator which must be either c i,↑ or c i,↓ . Hence the operator A carries charge 1, and we have B = [A, η † ] = 0. Taking the 1-D singlet pairing state as an example, it is easy to derive Thus to satisfy Eq. (12), we need χ i A = 0, and the simplest choice is Considering the superpositions of these operators, we have the set of jump operators: with two-spinor c i = (c i,↑ , c i,↓ ) T and Pauli matrices σ a with a = ±, z. Alternatively, we start with the physical intuition that d-wave pairing states may be viewed as delocalised antiferromagnetic order away from half-filling. We may consider a unit cell of the Néel state S ± i,j |vac , where the Néel state unit cell operator S ± i,j = c † i σ ± c † j , and j is one of the nearest neighbors of site i which thus creates an adjacent pair of fermions with opposite spin. We notice that the singlet pairing operator in 1-D can be viewed as the superposition of two antiferromagnetic unit cell operators It is thus instructive to first construct the jump operators for an antiferromagnetic Néel state. To annihilate this unit cell Néel state, we simply need the Lindblad operators of the form: This jump operator generates hopping with a spin flip, which is impossible in the case that antiferromagnetic order is already present, due to the Fermi statistics. Generalising the unit cell operator to a 2-D lattice, we notice the Néel state can be written in eight different forms, |N± = i∈A S ± i,i+eν |vac = (−1) M/2 i∈B S ∓ i,i−eν |vac , with M the lattice size, and e ν = {±e x , ±e y }. The 2-D Lindblad operators corresponding to those in Eq. (19): These operators impose the quasi-local constraint on the steady state that any given site should have opposite spin with its nearest neighboring sites, which guarantees antiferromagnetic order at half-filling. However, there is still a two-fold degeneracy where the antiferromagnetic order differs by a total spin flip. However, there is still a two-fold degeneracy where the antiferromagnetic order differs by a total spin flip. As the Néel states cannot be reached by implementing the jump operators in Eq. () on any other states, and as the complete set of jump operators is invariant under a Hermitian conjugation, the dark subspace, (containing the Néel states) is isolated from the rest of the Hilbert space under the dissipative dynamics given by Eq. (). These problems can be solved, in principle, by adding a single jump operator, or more naturally a set of jump operators of the form: for arbitrary i and its nearest neighbor j, and for either spin σ =↑, ↓. As an example, for i ∈ A and σ =↓, the dark state is |N+ . Note that comparing the jump operators in Eq. (19) with the unit cell operators S ± i,j , we find that the jump operators can be obtained from the state generating unit cell operator via a particle-hole transformation c † j → c j on the central site j. Now that we have found the Lindblad operators for the Néel state, we can proceed to generalise the operators to the d-wave case. First, we identify the "d-wave unit cell operators":D where ρ ±x = −ρ ±y = 1 and a = ±. We then perform a particle-hole transformation on the central site as in the case of the Néel state, and find operators: It is easy to verify that [J α i , jD b j ] = 0 (α = ±, z), which is dictated by the Fermi statistics. Note that the d-wave coherence is established via quasi-local phase locking between adjacent cloverleaves of sites, as illustrated in Fig. 2. In this respect, the jump operators act similarly to those that establish phase coherence in bosonic systems [12]. For fermionic pairing states, however, as explained above, Pauli blocking is an additional key ingredient.
We now consider the general case of quasi-local pairing states that can be factorised in the spin degrees of freedom. The unit cell operator of these pairing states can be expressed as superpositions of the Néel state unit cell operators: where the coefficients ρ ν encode the spatial symmetry of the pairs, and σ 1 and σ 2 can be arbitrary combinations of spin configurations, including spinless fermions σ 1 = σ 2 . Performing the particle-hole transformation as above, we get the following jump operators: which gives [J j , i β † i ] = 0, so long as the coefficients satisfy the following equations: where the spin degrees of freedom have been factored out, and we have omitted the spin index in Eq. (27), as there is only one spin species. For the spinful case, it is easy to verify that Eq. (26) always holds regardless of the structure of ρ ν . The spinless case though, is non-trivial in general. A particularly interesting example following Eq. (27) above is a 2-D p-wave state of spinless fermions generated by The p-wave pairing state in 2-D can be prepared in two different chiralities (see Fig. ): p x + ip y and p x − ip y , which shares the spatial symmetry with the pairing states in topological superconductors. However, as we will show later in the mean-field analysis, the p-wave states prepared in this way are still in the strong pairing limit, such that they are topologically trivial. The generation of topological order in 2-D systems will be discussed in a forthcoming publication. For stable dissipatively induced topological order in 1-D, see [20].

Jump operators for fixed-phase state
In the previous discussion, we have focused on the jump operators for dark states in the form of Eq. (5), i.e. states with fixed total particle number. The dissipative processes characterised by these jump operators necessarily conserve the total particle number of the system, as is clear from the commutation relation [J a,z i ,N ] = 0, whereN is the total particle operator. The total particle number of the dark state in this case is given by that of the initial state.
On the other hand, one can show (c.f. Appendix B) that for any given numberconserving pairing state which is a dark state of single-particle jump operators, one can always construct a set of linear jump operators that have the fixed-phase state in the form of Eq. (7) as a dark state. In the case of the d-wave pairing state, we consider the following jump operators where ρ ±x = −ρ ±y = 1 as given by the d-wave symmetry, and P and Q are complex numbers. It is easy to show that j i,σ |ψ d = 0, provided that P/Q = α, with the fixedphase pairing state |ψ d defined in Eq. (7). Note that the jump operators for the fixed-phase state do not conserve the total particle number, while the average particle density in the dark state is given as is the pair wavefunction in momentum space, and the summation over q runs over the first Brillouin zone. Importantly, the average particle density of the dark state here is determined by the parameters of the dissipative process, i.e. |α| can be chosen to fix a desired average density.
It is straightforward to extend the analysis to pairing states with other symmetries. For the spinless p-wave pairing state for example, the jump operators for the fixed-phase dark state |ψ p = N exp(αp † )|vac are with P/Q = 2α, and the factor of 2 is due to the triplet pairing symmetry of the p-wave pairing state. For a more detailed discussion on jump operators for fixed-phase states, see Appendix B.

Completeness of the jump operators
To achieve the dissipative preparation of many-body states, the final steady state of the master equation should be unique. This requires that (i) the dark state be unique; (ii) the non-existence of stationary solutions other than this dark state [13]. Taking the 1-D singlet pairing state as an example, we shall first analyse the uniqueness of the dark state from the symmetry perspective. We will then illustrate the uniqueness of the steady state with different pairing symmetries both in 1-D and 2-D by directly evolving the master equation for small finite size systems.

Uniqueness of dark state from a symmetry perspective
The problem of showing the uniqueness for all Lindblad operators is equivalent to showing the uniqueness of the ground state of the following positive semi-definite Hamiltonian (dimensionless): where J α i (α = ±, z or x, y, z) are the d-wave jump operators in Eq. (23). To see the equivalence, note that for any operator A, the matrix element ψ|A † A|ψ ≥ 0 is nonnegative. Thus, the spectrum of the Hermitian Hamiltonian is real and non-negative. Any zero energy eigenstate must be a ground state of the Hamiltonian. The energy is zero if and only if each term has zero energy. This is equivalent to J α i |Ψ = 0 ∀ i, α. Hamiltonians with these properties often occur in the context of spin models, where they are constructed as "parent Hamiltonians" for given states [35]. We follow this nomenclature here. This Hamiltonian serves for the uniqueness considerations of this section as well as for the adiabatic passage discussed in Sec. 5.
For later convenience, it is useful to collect the α = ± components into a dimensionless "reduced" parent Hamiltonian, with the normal-ordered form where σ =↑, ↓, and the quartic terms describe effective attractive two-body interactions. Note that the "chemical potential" term proportional to the total particle number is unimportant for fixed particle number states. For α = z, the normal-ordered form (dimensionless) is where the first term describes correlated hopping, and we have H p = H r p +H z p . Following the discussion in the previous paragraph, we see that the d-wave state |Ψ d is a ground state of the full parent Hamiltonian as well as of the reduced parent Hamiltonian. We will further demonstrate below that while the d-wave state is not the unique ground state of the reduced parent Hamiltonian, there are strong indications for it to be the unique ground state for the complete parent Hamiltonian based on symmetry arguments. These considerations will play an important role later in designing the implementation schemes.
By construction of the jump operators, the d-wave states are ground states of the parent Hamiltonian. For any symmetry operation, i.e. a unitary transformation T , that leaves the parent Hamiltonian Eq. (31) invariant T H p T −1 = H p , a necessary but not sufficient condition for the d-wave state to be the unique ground state of the parent Hamiltonian is: where φ T is the phase imposed onto the d-wave state by the unitary operation T . The parent Hamiltonian Eq. (31) has two obvious global symmetries: phase rotation invariance associated with particle number conservation, and translational invariance associated with momentum conservation. In the following, we discuss these symmetries: Phase rotation invariance -The symmetry is generated by Since the d-wave state is an eigenstate for both H p andN , with H p |Ψ d = 0 andN |Ψ d = 2N , it is also an eigenstate of T ϕ with eigenvalue exp 2iϕN . Thus for a given fixed particle number no degeneracies occur according to the above criterion.
Translation invariance -The symmetry is generated by the total center-of-mass momentum operator in 2-D: It is easy to check thatˆ P |Ψ d = 0, i.e. the d-wave state is a momentum eigenstate with zero eigenvalue. Therefore, the translational symmetry does not lead to degeneracies. Note the close relation of the terms in the momentum operator to the jump operators, which may be transformed into each other by flipping the spin on the site in the middle and changing from the symmetric to the antisymmetric superposition on sites {i ± e x , i ± e y }.
Discrete symmetries -On the bipartite lattices, and only for those, we find an additional discrete symmetry for the reduced parent Hamiltonian Eq. periodic boundary conditions in two dimensions. Given a certain finite lattice in any dimension, it is straightforward to check bipartiteness.
On the bipartite lattices we find the symmetry: and analogous for This transformation is canonical. A second quantised representation of the symmetry is highly non-local, similar to Shiba transformations for the Fermi-Hubbard model [23], but in this context we only need its action on the Hamiltonian and the state. Applying the transformation to the d-wave state we have where f (i, a) = 1 for a = +, i ∈ B and for a = −, i ∈ A; f (i, a) = −1 for a = +, i ∈ A and for a = −, i ∈ B. Thus the d-wave is not an eigenstate of T d which implies degeneracy under H r p but not H p . The degeneracy emerging from this is two-fold for any lattice size, which we will see in later sections.
In general these symmetries help to classify the lattice configurations (especially for finite lattices) under which we can expect uniqueness. Provided that there are no other symmetries of the full parent Hamiltonian that transform the d-wave state into other distinctive states, the d-wave state should be the unique dark state. Though we cannot rule out constructively the existence of other symmetries which may bring in additional degeneracies, the analysis here provides useful insights on the possible existence of degeneracy.

Uniqueness of the steady state: Numerical simulations
To further verify the uniqueness of the dark state, and more importantly, the uniqueness of the steady state for the dissipative dynamics given by the jump operators, we have performed numerical simulations of the master equation dynamics on finite size systems. Due to the translational symmetry of the pairing states, we impose periodic boundary conditions on the finite 1-D and 2-D systems. We have also taken periodic boundary conditions on the jump operators to reduce finite size effects and to be consistent with the definition of the pairing states on a finite lattice. In sufficiently large systems, the jump operators drive the atoms in the bulk into the desired phase, and the mixing of states at the boundary becomes negligible so long as the system is in the thermodynamic limit.
We have evolved the master equations for antiferromagnetic Néel states and singlet pairing states for finite size systems in 1-D. The results are shown in Fig. 4, which clearly demonstrate the uniqueness of the steady state anticipated above. Both the fidelity and entropy evolution indicates that in both cases the system is driven into the desired state regardless of the initial state. In comparison, when only jump operators of the reduced parent Hamiltonian are applied (dotted curve in Fig. 4(b)), the final fidelity approaches 0.5. For d-wave pairing states in 2-D, we have carried out a quantum trajectory simulation for small plaquettes. The evolution of the fidelity with respect to the d-wave pairing state indicates that the system approaches the final pure steady state exponentially, which implies the existence of a dissipative gap, analogous to the energy gap for the ground state of the BCS pairing state (see Fig. 5(a)). While in a finite system, a gap of the order L −2 is expected due to the finite linear dimension L, in the thermodynamic limit, this gap vanishes. In the next section, we will derive a meanfield theory of the master equation in the thermodynamic limit, which shows that a dissipative gap appears naturally from the mean-field expansion of the master equation for pairing states, demonstrating the dissipative gap in our numerical simulations is not a mere finite size effect. Similar results can be obtained for the paired states of p-wave symmetry for spinless fermions (see Fig. 5(b)).

Mean-field expansion of the master equation and the dissipative gap
In this section, we will develop a mean-field theory of the master equation for drivendissipative pairing states in the thermodynamic limit which is valid at late times, where the system is close to the BCS-type pairing state. For this purpose, we switch from the fixed-number state representation to the fixed-phase (coherent state) representation. The justification for this procedure builds on two properties. First, the exactly known fixed-number pairing dark states discussed above exhibit phase locking among different fermion pairs. Such a property is reproduced in the coherent state representation |ψ d ∝ n α n (D † ) n n! |vac = exp αD † |vac (α = |α|e iθ ), with a fixed global phase θ [c.f. Eq. (7)]. Second, for the BCS wavefunctions under consideration, these representations The insets indicate the existence of dissipative gaps in both cases, which render the convergence to the steady states exponentially fast. This result is robust in the thermodynamic limit as revealed by our mean-field theory. The fidelity (solid) is calculated by averaging over 1000 trajectories. These trajectories are then bunched into 100-trajectory groups, whose standard deviations are then calculated to show the sampling errors (dashed). are equivalent in the thermodynamic limit, which can be verified explicitly from a consideration of the relative number fluctuations in the fixed-phase state. Indeed, one finds for the variance ∆N = N , whereN is the total particle number operator, the average is taken with respect to the fixed-phase BCS state, and N is the number of degrees of freedom. The general normalised wavefunction of the BCS pairing state then reads [36]: where ϕ q is the pair wavefunction in momentum space, θ is the overall phase of the pairing state, |α| is a real number fixing the average particle number density, and the product over q runs over the first Brillouin zone. Note that for 1-D singlet state, ϕ q = cos q, while for 2-D d-wave state, ϕ q = cos(q x ) − cos(q y ). The fixed-number d-wave state can be obtained from the coherent state representation via a number projection, where N is the total particle number corresponding to the density In the following, we will take α = 1 for simplicity, which corresponds to an initial state with given total particle density n = dq (2π) d 2|ϕq| 2 1+|ϕq| 2 . For initial states with different total particle densities, one just needs to substitute ϕ q in the following discussions with αϕ q , where |α| is determined from the number equation Eq. (40). Finally, we note that in the grand canonical ensemble, the overall phase θ of the pairing state can take any value. Hence we will take θ = 0 in the following discussion and write |D ≡ |D(0) . Indeed, the phase θ is not determined by the microscopic dynamics, since the jump operators are particle number conserving (charge 0). Therefore, the phase will be chosen spontaneously in the thermodynamics limit. This is the analog of the spontaneous symmetry breaking in the dissipative context.
The d-wave state has three non-zero bilinear expectation values for each momentum mode relevant to the corresponding quantities in the subsequent mean-field theory, i.e. particle number and order parameter: where σ =↑, ↓. All other expectation values vanish on the above state. A mean-field theory analogous to the BCS approximation for superconductivity can be set up based on the proximity of the density matrix to the d-wave state, giving rise to an ordering principle for devising a controlled mean-field approximation in the late time evolution, which is thus useful to study the final stages of the master equation evolution.
Starting from the fact that the coherent representation of the pairing state is a product in momentum space, we will require this property also for an approximate ansatz for the density matrix, ρ = q ρ q for the solution of the master equation at late times. This may be viewed as a Gutzwiller factorisation approach in momentum space, which has been used previously for a mean-field decoupling of bosonic master equations [17]. This ansatz will enable us to derive a late time master equation quadratic in the fermion operators which contains information of the complex excitation spectrum, i.e. the damping of the lowest fermionic single particle excitations.
To implement the approximation, we Fourier transform the jump operators to momentum space: where ϕ q is the pair wave function and reflects the similarity in construction of the jump operator and the corresponding pairing state. This gives rise to a Liouvillian in momentum space The mean-field product ansatz for the density operator is where each ρ q spans the subspace for {q ↑, −q ↓} [37]. We thus allow for a residual entanglement of the momentum modes {q, −q}, as in the BCS treatment, which is necessary to describe pairing. The second equation is the projection prescription to obtain the density operator for each mode which can be used to derive the equation of motion ∂ t ρ q . We take the partial trace tr =q on both sides of the master equation, for which we trace over all degrees of freedom outside the subspace {q ↑, −q ↓}. Then, in the spirit of BCS theory, we choose the relevant mean fields, i.e. macroscopically occupied expectation values in the dark state. As mean fields of linear and trilinear correlations vanish due to the Fermi statistics, and that of the quartic and higher order correlations connect momentum modes different from ±q, and thus are small compared to the macroscopic expectation values close to the steady state in the thermodynamic limit, we keep only mean fields of quadratic correlations, i.e. density mean fields or condensate with zero center of mass momentum [cf. Eqs. (41) (42) (43)]. We thus see how the proximity to the steady state can be used as an ordering principle for a mean-field theory at late times, which is based on the exact knowledge of the fixedphase steady state density matrix. Note that we use the commutation (as opposed to anti-commutation) properties during the process which is equivalent to the assumption that there is an even number of fermions in this mode. We then have the following results: where A ≡ dq (2π) d |ϕq| 2 1+|ϕq| 2 ≥ 0, and the integration runs over the first Brillouin zone.
Of particular interest and importance are the {γ q,σ } operators, which coincide with the definition of the Bogoliubov quasiparticles: Indeed, the mean -ield BCS pairing state is the vacuum state for the Bogoliubov quasiparticles, Furthermore, when properly normalised as in Eq. (49), these quasiparticle operators obey the closed fermion algebra: i.e. they are related to the original fermion operators by a canonical transformation. We note that the operators J a,z k do not exhibit a closed algebra structure. From the expressions in Eq. (48), we may identify a damping spectrum in the dissipative part of the master equation: with A = (1 − 1/ √ 2) in 1-D, and A ∼ 0.36 in 2-D. It is important to note that the damping spectrum is gapped, i.e. κ q ≥ 2Aκ is bounded from below. This behavior exhibits strong parallels to the equilibrium problem of paired fermions, where pairing is protected by an energy gap. Furthermore, it has the important implication that the approach to the d-wave dark state will be exponentially fast, in contrast to a bosonic system where wavelengths of arbitrary length make the approach to the dark states with long range order polynomially fast only [12]. This result is reflected in the quantum trajectory simulations in the previous section.
Similarly, for the parent Hamiltonian of the pairing state, it is straightforward to derive Clearly, the dissipative mean-field theory developed here can be applied to pairing states with other spatial symmetries. As an important example, we discuss the result for complex p-wave pairing states. The p-wave pairing state for spinless fermions can be written as: where q runs over half of first Brillouin zone, e.g. q x > 0. The pairing wavefunction ϕ * q = ϕ −q = −ϕ q , which is required for p-wave symmetry: in 1-D, ϕ q = 2i sin q; in 2-D, ϕ q = 2i(sin q x ± i sin q y ).
Following the previous derivations, we define the momentum space jump operators: where the summation over q is over the first Brillouin zone. Taking the partial trace over the degrees of freedom outside the subspace {q, −q} (q spans half the momentum space here, e.g. q x > 0), and identify the mean fields as before, we may arrive at the effective Hamiltonian where the summation of q is over half the first Brillouin zone (q x > 0). The Bogoliubov quasi-particles are given as The dissipative coefficient is given as Finally, we note that the fermionic quasi-particle operators in Eq. (57) formally correspond to the Bogoliubov operators of a p-wave Hamiltonian in the limit µ → −∞ (µ is the chemical potential), which means that the state is in the strong pairing limit describing a state of delocalised tightly bound molecular pairs and is topologically trivial [25]. The generation of stable topological order is however possible in a modified setting, as has been established recently [20].

Adiabatic passage to the ground state of the Hubbard model
As argued above, the pure mean-field state with the correct symmetry (antiferromagnetic at half-filling and d-wave otherwise) is a convenient initial state for the quantum simulation of the ground state of the Fermi-Hubbard model (FHM). With a suitable adiabatic passage, it would be possible to connect the pure dark state with the ground state of the FHM. A first guess as to how to reach this ground state based on the experience with purely Hamiltonian systems might be a simple adiabatic passage, in which the Liouvillian is switched off while the FH Hamiltonian is ramped up. Small scale numerical simulations of the time evolution suggest that this procedure does not work for our combined system with both unitary and dissipative evolution. This is understood from the fact that the mean-field state is not an exact eigenstate of the FH Hamiltonian, such that unitary and dissipative evolutions compete. As a result of this competition, the steady state density matrix in general describes a mixed state, instead of a pure, zero entropy state [12,18]. We thus observe that a dissipative gap cannot play the role of an energy gap in standard adiabatic passage schemes. We therefore propose a modified adiabatic passage, which uses the parent Hamiltonian of the mean-field state, constructed from the complete set of jump operators as in Eq. (31). By construction, this parent Hamiltonian has the mean-field state as a gapped ground state, and therefore provides for an energetic stabilization. The passage now proceeds by first turning off the dissipation while the parent Hamiltonian is applied, then simultaneously ramping down the parameters of the parent Hamiltonian while ramping up the parameters of the FHM. In this way, as long as the symmetry patterns of the mean-field target state and the ground state of the FHM Hamiltonian are the same and no phase transition is crossed, an energy gap persists through the whole passage. Indeed, we show numerically that this modified adiabatic passage ensures an efficient transfer into the desired ground state.
For the antiferromagnetic Néel state at half-filling, the parent Hamiltonian (dimensionless) reads where the jump operators j ± i,j and j i are defined in Eqs. (19,21). We have performed numerical simulations of such an adiabatic passage for a 2 × 2 plaquette with 4 atoms. The result is shown in Fig. 6(a). Indeed the initial Néel state can be adiabatically connected with the ground state of FHM at half-filling with high fidelity.
For the case of d-wave state, the parent Hamiltonian in Eq. (31) by construction has the initial d-wave state as an exact eigenstate and thus supports the d-wave state obtained from the dissipative evolution. From Eq. (53), it is clear that the single fermion excitations on the d-wave state are gapped if the system is sufficiently far away from half-filling. As a consequence, all requirements for an efficient adiabatic passage are met. Note that single fermion excitations above the ground state manifold of the reduced parent Hamiltonian are also suppressed by an energy gap, as H r eff = 1 2 H eff on the mean-field level. We will make use of this important fact in Sec. 6.3 to design a modified adiabatic passage.
The time-dependent Hamiltonian describing the adiabatic passages is given as where the time dependent coefficients h(t), U (t), J(t) give the precise path of the adiabatic passage. In practice, the time dependence of these coefficients is given by the rate at which the lattice potential giving rise to the FHM is ramped up, as well as by the rate of the effective interaction given by the parent Hamiltonian H p . Here, for simplicity, we have chosen linear ramps for these coefficients (see insets of Fig. 6), which already give a clear physical picture of the adiabatic passage. In practice these ramps could be further optimised, so that higher fidelities can be achieved in shorter ramp times. Consistent with the previous discussion, the role of the parent Hamiltonian is to provide an energy gap, and hence energetically stabilise the adiabatic passage. We have performed numerical simulations of the adiabatic passage with various finite size systems. To avoid degeneracies of the ground state of FHM due to finite size effects, we have taken open boundary conditions for the Hubbard Hamiltonian during the adiabatic process, while we retain the periodic boundary conditions for the definition of the initial pairing state and for the jump operators. We expect that the mean-field d-wave state should be efficiently connected to the ground state of the Fermi-Hubbard Hamiltonian so long as the d-wave symmetry of the ground state of the Fermi-Hubbard model is present and not completely destroyed by the finite size effect. We find that this is the case for ladder systems. A typical result of our simulation on a finite ladder is shown in Fig. 6(b). For systems in the thermodynamic limit, the symmetry property of the ground state is not affected by the boundary effect, and we expect an efficient adiabatic passage so long as the symmetries of the ground state are the same as the dissipatively driven initial state [38,39,40].

Physical implementation and modified adiabatic passage
As an illustrative example, we now discuss a proof-of-principle implementation of the single-particle jump operators. The scheme we describe in this section is stroboscopic, and involves realising the action of the jump operators in a series of steps. Though non-trivial to implement in present experiments, this example is made up of elements that are presently accessible in experiments. The example illustrates how the properties of the operators appearing in the previous sections, specifically that they are quasilocal, conserve particle number, and can be implemented based on single-particle 171 Yb 3 P 0 Figure 7. Level scheme using 171 Yb atoms. The physical spin state is encoded in the nuclear spin sublevels of the 1 S 0 manifold. The spin flip operation is implemented via off-resonant coherent coupling to the 3 P 0 manifold with circularly polarised light (red arrows). The long lived 3 P 0 states are coupled to the 1 P 1 level in a two-photon process, from which spontaneous emission into a cavity is induced, leading back to the 1 S 0 manifold. operations, make them favourable for experimental implementation. For an alternative non-stroboscopic, i.e. "always-on" continuous implementation, which is applicable in the case of spinless (spin-polarised) fermions such as the p-wave case discussed above, see [20].
Our example takes advantage of the properties of alkaline-earth-like atoms [41,42,43,44,45,46]. With two valence electrons, these atoms possess metastable triplet levels, and fermionic isotopes have non-zero nuclear spin (e.g., I = 1/2 for 171 Yb, which we will choose here). This nuclear spin acts as an independent degree of freedom in the ground 1 S 0 and lowest excited 3 P 0 manifolds. Here, the nuclear spin will play the role of the physical fermionic spin degree of freedom, and the 3 P 0 manifold will be used as an intermediate state in the dissipative process. These states are depicted in Fig. 7. Note that as 1 S 0 and 3 P 0 are optically separated, they can be trapped in independent lattices using dipole traps at different wavelengths [47].
As a simple example, we will first discuss the implementation of jump operators for driving the system into the antiferromagnetic Néel state at half-filling. We will then move on to the more complicated cases of pairing states.

Antiferromagnetic Néel state
Since the Lindblad jump operators for the Néel states act on unit cells of two sites, we will focus on operations on two adjacent sites i and j. These will be carried out in parallel on pairs of adjacent sites along the lattice. As illustrated in Fig. 7, spins states are encoded in the nuclear spin sublevels of the 1 S 0 manifold. We further assume that atoms in the 1 S 0 and 3 P 0 manifolds are trapped in independent optical lattices, and that 3 P 0 is trapped in a superlattice with a period of two sites, defining pairs of sites, where we label the left well i and the right well j. Initially, the superlattice potential is arranged in such a way that the potential well at site i is much deeper than at site j (with an energy difference of the lowest state in each well of the order of several kHz). The action of the jump operator c † j,↑ c i,↓ can then be realised by performing the following operations (see Fig. 8): (1) apply a circularly polarised π-pulse selectively on site i coupling the 1 S 0 and 3 P 0 manifolds so that any atom originally in the state | ↓, 1 S 0 will end up in | ↑, 3 P 0 ; (2) adiabatically manipulate the superlattice potential so that the population at site i is transferred to site j; (3) couple the 3 P 0 and 1 P 1 manifolds off-resonantly, so that the population in | ↑, 3 P 0 should decay to | ↑, 1 S 0 on site j if and only if the state on site j is empty; (4) repeat the steps (2) and then (1) (in reverse order) to bring any remaining population in 3 P 0 back to the ground state manifold.
Before moving on to extend the scheme to jump operators associated with pairing states, several comments are in order: (a) the jump operator is implemented stroboscopically, which places requirements on the time scale of each step of operations listed above, such that the total time of evolution should be much longer than the time scale of operations; (b) the jump operators are implemented in parallel for each pair of lattice sites along the lattice; (c) during the excitation of the population from 1 S 0 to 3 P 0 , as the line-width of the metastable 3 P 0 is on the order of 10mHz for 171 Yb, the bias between different subwells in the superlattice potential ensures site selectivity; (d) Figure 9. Implementation scheme for jump operators of 1-D d-wave state. (a) Only one of the decay channels is blocked; (b) both decay channels are Pauli blocked, therefore the jump operator does not change the system configuration. The state is a local dark state for this local jump operator. the nuclear spin is conserved during the decay process, which is guaranteed by the large detuning from the 1 P 1 manifold. The nuclear spin conservation can also be realised in this case by applying a large magnetic field so that electronic spin and nuclear spin are decoupled [44].

d-wave pairing state
We now extend the ideas above to the implementation of jump operators for driving the system into d-wave pairing states. For the d-wave jump operators, an additional constraint on the dissipative process is that atoms on quasi-local sites, e.g. site i + e x and i − e x , should decay coherently. To satisfy this requirement, we couple the system to a cavity with a finite linewidth. An atom (or atoms) at the sites i + e x and i − e x will then be coupled collectively to the cavity mode, ensuring that the decay is coherent. For clarity, we first describe the implementation procedures in 1-D, and choose the example The step-by-step implementation scheme is shown in Fig.  9: (1) We first assume that the 3 P 0 state is initially trapped in a lattice of three times the period as that for the 1 S 0 state, defining blocks of three sites in the original lattice. Using this, we excite any spin-down atom in 1 S 0 on central site to the spin-up state of the 3 P 0 manifold, using σ + light. (2) We then add an additional potential, splitting this site into two, and separate these sites so that the mode of atoms confined in them overlap the right and left sites of the original three-site block. (3) We induce decay by coupling atoms in the 3 P 0 state off-resonantly to the 1 P 1 state, with coupling strength Ω, and detuning ∆. If we couple the 1 S 0 -1 P 1 transition to a cavity mode with linewidth Γ and vacuum Rabi frequency g, then the decay will be coherent over the triple of sites.
In the limit ∆ Ω, Γ and Γ max( Ω 2 ∆ , g 2 ∆ , Ωg ∆ ), we adiabatically eliminate the cavity mode and the intermediate far off-resonant state 1 P 1 , and obtain an effective decay rate Γ eff = Ω 2 g 2 ∆ 2 Γ ∼ 9kHz for typical parameters (see Appendix C). Note that Fermi statistics will be observed in this process, and that we assume that the atoms remain in the lowest band, as all parameters are smaller than the trapping frequency in the lattice (see Ref. [31] for more details of Pauli-blocking of spontaneous emissions in this sense).
Other jump operators, J − i and J z i can be implemented by applying rotations in the nuclear spin before and after the three steps above. For J − i , one exchanges the spins with a π-pulse, whereas for J z i , one must apply a π/4 rotation in the nuclear spin basis before and after the operation. In addition, for J z i , both spin states should be excited, and coherence of nuclear spins is maintained throughout the operation. This can be achieved by either going far off-resonant for the field coupling 3 P 0 and 1 P 1 manifold, or by applying a large magnetic field as described in Ref. [44].
This scheme can be generalised to 2-D by considering 3-by-3 plaquettes defined by the appropriate superlattice potential for the 3 P 0 level. As in the 1-D case, we require an adiabatic manipulation of this potential in step (ii), although here the depths of the wells must be adjusted to ensure that the correct relative phases are obtained for atoms "transported" in different directions (see Fig. 10 and its caption). Figure 11. Implementation of (J ± i ) † J ± i in the parent Hamiltonian for 1-D case. Firstly, the population of the outer sites are excited to the upper lattice. The super lattice is then adiabatically tuned from a double well structure to a single well, so that the state |L + |R is projected to the lowest level of the single well potential. The interaction is then induced via a Feshbach resonance for instance after exciting the population of the opposite spin in the central site to the superlattice.

Implementing the reduced parent Hamiltonian and modified adiabatic passage
Here, we extend the scheme above to implement the reduced parent Hamiltonian for the d-wave state stroboscopically. We see from the discussion in Sec. 3 that the mean-field d-wave state is in the ground state manifold of the reduced parent Hamiltonian. As we have discussed in Sec. 5, this degenerate ground state manifold (two-fold in 1-D, four-fold in 2-D) is protected by an energy gap from single fermion excitations under the reduced parent Hamiltonian. Furthermore, we will also show below that an adiabatic passage with high fidelity can be achieved by a modified adiabatic passage scheme.
The implementation of the reduced parent Hamiltonian is similar to that of the jump operators, except that the dissipative part is replaced by an induced phase shift. As shown in Eq. (32), the reduced parent Hamiltonian contains an effective interaction term and a term proportional to the total particle number that is not important for states with fixed particle number. To implement the effective interaction term stroboscopically, as illustrated in Fig. 11 with the example of 1-D Hamiltonian (J + i ) † J + i , the following steps are required: (1) any spin down atoms in the left and right well of the ground state lattice potential are transferred to the superlattice potential of the 3 P 0 manifold; (2) the double-well in the superlattice potential is merged into a single well, during which process the symmetric state in the double-well potential is mapped to the lowest motional state of the final single well potential; (3) a phase shift is then induced to generate an interaction only if the spin-down state in the 1 S 0 manifold and the spin-up state in the 3 P 0 manifold are simultaneously occupied. This can be achieved, e.g., by applying a π pulse between the spin-down state in the 1 S 0 state and the spin-down state in the lowest motional state in the superlattice potential, and inducing an interaction between the different spin states in the 3 P 0 manifold via an optical Feshbach resonance. Finally, to implement (J − i ) † J − i , the spins should be exchanged while the above procedure is carried out. With only the reduced parent Hamiltonian, we find that given an optimised ramping scheme, the mean-field d-wave state can still be adiabatically connected with the ground state of FHM on small lattices. Fig. 12 shows such an example, where ramps with the complete parent Hamiltonian and with the reduced parent Hamiltonian are numerically simulated for 4 atoms on a 2×6 ladder. For the adiabatic passage with the reduced parent Hamiltonian, we ramp up J(t) and U (t) separately. In both cases, we have very high fidelity at the end of the ramp. For ramps with the reduced parent Hamiltonian, the high fidelity is due to the large overlap (∼ 0.95 for most ladder systems) between the mean-field d-wave state and the ground state of the time dependent Hamiltonian at the beginning of the ramping process (when there is no on-site interaction). For the numerical simulations that we considered here, this overlap also sets the upper bound for the final fidelity of the ramps.

Conclusions
We have proposed an approach for the preparation of many-body pairing states of given symmetry for fermionic atoms in an optical lattice via driven dissipative processes based on suitable reservoir engineering. We have discussed in detail the strategy of designing the jump operators making use of the Fermi statistics, which gives rise to the dissipative preparation of the initial pairing state. This process is in general efficient, due to the existence of a dissipative gap. We then argued for the uniqueness of the pairing state as the steady state of the dissipative dynamics, both from symmetry considerations, and via small scale numerical simulations. Note that for realistic finite size systems such as plaquette geometries [48], it is also possible to design jump operators for specific many-body states defined on the finite system, in which case one may need to design special "boundary" jump operators to make the state unique.
We then discussed the adiabatic passage process that could be used to connect the driven-dissipative mean-field state with the ground state of the FHM. As our dwave state is not an eigenstate of the FHM, directly ramping down the dissipation rate while ramping up the FHM leads to competition between the coherent and dissipative dynamics which would not drive the system into the ground state. We therefore introduced the parent Hamiltonian of the d-wave state, a semi-positive Hermitian Hamiltonian constructed from the jump operators. By construction, the parent Hamiltonian has the dark state of the dissipative process as its ground state. We illustrated via small scale numerical simulations that the ground state of the FHM can be adiabatically connected with the mean-field state of the relevant symmetry via optimised adiabatic paths. This is in similar spirit to the recent experimental demonstration of antiferromagnetic order in an optical lattice [5], where the desired eigenstate of the Ising model is prepared via adiabatic passage from a starting state that has low entropy and sufficient overlap with the final state. We note that it is possible to extend these small scale numerical calculations by applying time dependent density matrix renormalisation group (t-DMRG) methods. In fact, quantum trajectories methods could be combined with t-DMRG methods [49] in order to perform larger-scale simulations of the dissipative preparation process and the adiabatic ramp together.
Finally, we discussed a proof-of-principle physical implementation of both the jump operators and the parent Hamiltonian using alkaline-earth-like atoms, which illustrated that the properties of the jump operators discussed here are favourable for implementation. We mainly focused on the implementation of d-wave pairing state, but similar implementations can be readily found for pairing states of other symmetries, so long as the jump operators are quasi-local and involve operations manipulating only single particles.
Appendix A. Two-particle jump operators for d-wave pairing state In this appendix, we derive in detail the two-particle jump operators for pairing states with d-wave symmetry as appeared in Eqs. (16,17).
Appendix A.1. One-dimensional case As an example, we choose The commutation relations then give It is straightforward to show that if we define The symmetry in the expressions above suggest that we may simplify these operators by choosing ξ i = c i+1↓ c i↑ . The commutation relations then have the form: The most straightforward choice of χ i would be χ i = B i , as B 2 i = 0. More generally, χ i B i = 0 is satisfied so long as the pair operators in χ i can be factored out to contain . This actually allows some freedom in choosing the remaining part of the χ i operator.
However, in this second scenario, the existence of a constant term in Eq. (A.5) renders Eq. (14) not equal to zero even if Eq. (15) is satisfied. The resulting jump operator would then not give the desired dark state. To solve this problem, one needs to introduce appropriate symmetry into the design of the jump operator. Notice that assuming the translational symmetry, the creation operator of the state can also be written as: Correspondingly, we examine the following factorised ξ i : Note that the choice of the negative sign here is to ensure that no constant terms appear in the expression for A i . For the commutation relations, we now have This implies that we may satisfy the dark state requirement by choosing a jump operator of the form as given in Eq. (16).

Appendix A.2. Two-dimensional case
We define whose commutation relations are: Following the same derivation as in the previous section, we find where C † is an arbitrary superposition of single-fermion creation operators. Note that this is the most straightforward choice to satisfy χ i B i = 0, other solutions may still exist. Finally, we see that in the case of a d-wave state on a 2-D lattice, the jump operator takes the form: 16) where ρ ±x = 1, ρ ±y = −1, as given in Eq. (17).

Appendix B. Jump operators for the fixed-phase state
In this Appendix, we discuss the general formalism for the construction of jump operators for the fixed-phase state, starting from number-conserving jump operators with known unique dark state with fixe particle number. These jump operators describe dissipative processes for which the total particle number is not exactly conserved, whereas the average particle number approaches the steady state value determined by the parameters of the dissipative process. In the following, we will first discuss the pairing states of spinful fermions, before extending the formalism to spinless fermions.

Appendix B.1. Pairing states with spins
We only consider separable pairing states, i.e. pairing states whose spin degrees of freedom can be factorised. Then the general number-conserving pairing state for spinful fermions can be written as: where C † i = ν ρ ν c † i+eν ,↑ and A † i = µ λ µ c † i+eµ,↓ are translation invariant. Fourier transform the pairing operators into the momentum space, where f k = ν ρ ν e −ik·eν and g k = µ λ µ e −ik·eµ . With these, the fixed-phase correspondence of Eq. (B.1) can be written in the form of a coherent state: where α is a complex number carrying the phase of the pairing state, N is the normalization factor, and q runs over the first Brillouin zone in the product. Without loss of generality, this coherent state can be re-arranged into the standard BCS-type mean-field wave function: where the coefficients u k = 1+|αf k g −k | 2 . Apparently, the coherent state Eq. (B.5) is the vacuum for the Bogoliubov quasiparticle operators: as it is easy to verify the following relations, γ k,↑ |ψ c = γ k,↓ |ψ c = γ −k,↑ |ψ c = γ −k,↓ |ψ c = 0. Therefore, these Bogoliubov quasiparticle operators are the momentum space jump operators for the fixed-phase state Eq. (B.5). Based on Eqs. (B.6, B.7), it is easy to find a more general form of the momentum space jump operators where ϕ ± (k) are arbitrary functions of k. Fourier transforming Eqs. (B.8, B.9) back to the coordinate space, we immediately get the quasi-local jump operators that we look for.
As an illustrating example, let us investigate the simple case with g k = 1, ϕ ± (k) = 1, which implies the structure of the pairing state should be completely encoded in the spin-up degrees of freedom. The coherent state in this case becomes We then Fourier transform the momentum space jump operators to the coordinate space, For the d-wave pairing state, ρ ±x = −ρ ±y = 1, and we recover the fixed-phase jump operators for the d-wave pairing state in Sec. 2.3. As a consistency check, we demonstrate below that the coherent state in Eq. (B.4) is a dark state of the jump operator in Eq. (B.14): Finally, we note that the fixed-phase jump operators presented here are superpositions of a creation operator and an annihilation operator, and hence do not conserve the total particle number. The average particle number on the other hand, is driven to the final steady state value during the dissipative process. The average particle number in the final steady state is given as: (B.16)

Appendix B.2. Spinless pairing states
One can easily extend the formalism above to the pairing states of spinless fermions, with modifications to the pairing parameters due to the triplet pairing symmetry. The pairing state with fixed phase can be written as where indicates that k only runs over half of the Brillouin zone, e.g. k x > 0. The coefficients here are similar to the spinful case, with u k = 1+|2αf k g −k | 2 , f k = ν ρ ν e −ik·eν and g k = µ λ µ e −ik·eµ . As a simple example, we consider the case with g −k = 1. The coherent state then becomes where due to the triplet pairing symmetry, f −k = f * k = −f k . Following the previous approach, the general form of the momentum space jump operators are where ϕ(k) is an arbitrary function of k. For simplicity, we take ϕ(k) = 1, and the Fourier transform of Eq. (B.19) gives the quasi-local jump operators: Similar to Eq. (B.15), it is straightforward to check that γ i |ψ p = 0. Finally, the average particle number is given by where the summation runs over the entire first Brillouin zone. With ρ x = −ρ −x = −iρ y = iρ −y = 1, we recover the results for p-wave pairing states in Sec. 2.3.