{'fraction_non_alphanumeric': 0.038483369281416734, 'fraction_numerical': 0.010239527755704393, 'mean_word_length': 4.54735516372796, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dissipative topological superconductors in number-conserving systems

We discuss the dissipative preparation of p-wave superconductors in number-conserving one-dimensional fermionic systems. We focus on two setups: the first one entails a single wire coupled to a bath, whereas in the second one the environment is connected to a two-leg ladder. Both settings lead to stationary states which feature the bulk properties of a p-wave superconductor, identified in this number-conserving setting through the long-distance behavior of the proper p-wave correlations. The two schemes differ in the fact that the steady state of the single wire is not characterized by topological order, whereas the two-leg ladder hosts Majorana zero modes, which are decoupled from damping and exponentially localized at the edges. Our analytical results are complemented by an extensive numerical study of the steady-state properties, of the asymptotic decay rate and of the robustness of the protocols.


I. INTRODUCTION
Topological quantum computation has recently emerged as one of the most intriguing paradigms for the storage and manipulation of quantum information [1,2]. The defining features of topological order, namely the existence of degenerate ground states which (i) share the same thermodynamic properties and (ii) can only be distinguished by a global measurement, portend for a true many-body protection of quantum information. Additionally, the non-Abelian anyons which typically appear in these models are crucial for the active manipulation of the information, to be accomplished through their adiabatic braiding [3,4].
Among the several systems featuring topological order, free p-wave superconducting systems with symmetry protected topological properties have lately attracted a significant amount of attention [5][6][7]. On the one hand, they are exactly-solvable fermionic models which help building a clear physical intuition of some aspects of topological order [8,9]. On the other one, they are physically relevant, and several articles have recently reported experimental evidences to be linked to p-wave-like superconductors featuring zero-energy Majorana modes [10][11][12][13][14].
Whereas up to now these experimental results have been obtained in solid-state setups, it is natural to ask whether such physics might as well be observed in cold atomic gases [15], which owing to their well-controlled microscopic physics should allow for a more thorough understanding of these peculiar phases of matter. Important theoretical efforts have thus proposed a variety of schemes which exploit in different ways several properties of such setups [16][17][18][19][20][21][22].
Among these ideas, that of a dissipative preparation of interesting many-body quantum states [23,24] is particu-larly appealing: rather than suffering from some unavoidable open-system dynamics, such as three-body losses or spontaneous emission, one tries to take advantage of it (see Refs. [19,[25][26][27][28] for the case of states with topological order, such as p-wave superconductors). The key point is the engineering of an environment that in the long-time limit drives the system into the desired quantum state. This approach has the remarkable advantage of being a workaround to the ultra-low temperatures necessary for the observation of important quantum phenomena which constitute a particularly severe obstacle in fermionic systems. The trust is thus that the mentioned "non-equilibrium cooling" may open the path towards the experimental investigation of currently unattainable states, e.g. characterized by p-wave superconductivity.
In this article we discuss the dissipative engineering of a p-wave superconductor with a fixed number of particles, an important constraint in cold-atom experiments. We consider two different setups: (i) A single quantum wire, introduced in Ref. [19]; this system displays the typical features of a p-wave superconductor but it is not topological in its number conserving variant. (ii) A twoleg ladder [21,[29][30][31][32][33][34], supporting a dissipative dynamics which entails a two-dimensional steady state space characterized by p-wave superconducting order with boundary Majorana modes for every fixed particle number.
We identify the p-wave superconducting nature of the steady states by studying the proper correlators, which saturate to a finite value in the long distance limit. Their topological properties are best discussed using a mathematical connection between dark states of the Markovian dynamics and ground states of a suitable parent Hamiltonian. In both setups we demonstrate that the dissipative gap closes at least polynomially in the system size and thus that the typical decay time to the steady state di-verges in the thermodynamic limit. This contrasts with the case where number conservation is not enforced. In this case typically the decay time is finite in the thermodynamic limit [19,26], and reflects the presence of dynamical slow modes related to the particle-number conservation [35,36], which also exist in non-equilibrium systems (see also Ref. [37][38][39]).
Our exact analytical findings are complemented by a numerical study based on a matrix-product-operator representation of the density matrix [40,41], one of the techniques for open quantum systems which are recently attracting an increasing attention [42][43][44][45][46][47][48][49][50]. These methods are employed to test the robustness of these setups to perturbations, which is thoroughly discussed.
The article is organized as follows: in Sec. II we review the key facts behind the idea of dissipative state preparation using the dark states of a many body problem, and exemplify them recalling the problem studied in Ref. [19]. A simple criterion for signalling the divergence of the decay-time with the system size is also introduced. In Sec. III we present the exact analytical study of the single-wire protocol, and in Sec. IV a numerical analysis complements the previous discussion with the characterization of the robustness to perturbations of these setups. In Sec. V we discuss the protocol based on the ladder geometry. Finally, in Sec. VI we present our conclusions. The dissipative dynamics considered in this article is Markovian and, in the absence of a coherent part, can be cast in the following Lindblad form: where L is the so-called Lindbladian super-operator and theL j are the (local) Lindblad operators. We now discuss a fact which will be extensively used in the following. Let us assume that a pure state |Ψ exists, with the property: A simple inspection of Eq. (1) shows that |Ψ is a steady state of the dynamics, and it is usually referred to as dark state. Although the existence of a state satisfying Eq. (2) is usually not guaranteed, in this article we will mainly consider master equations which enjoy this property. A remarkable feature of dark states is that they can be searched through the minimization of a parent Hamiltonian. Let us first observe that Eq. (2) implies that Ψ|L † jL j |Ψ = 0 and since every operatorL † jL j is positive semi-definite, |Ψ minimizes it. Consequently, |Ψ is a ground state of the parent Hamiltonian: Conversely, every zero-energy ground state |Φ of Hamiltonian (3) is a steady state of the dynamics (1). Indeed,Ĥ p |Φ = 0 implies that Φ|L † jL j |Φ = 0 for all j = 1, . . . , m. The last relation means that the norm of the statesL j |Φ is zero, and thus that the states themselves are zero:L j |Φ = 0. As we have already shown, this is sufficient to imply that |Φ is a steady-state of the dynamics.
In order to quantify the typical time-scale of the convergence to the steady state, it is customary to consider the right eigenvalues of the super-operator L, which are defined through the secular equation L[θ λ ] = λθ λ . The asymptotic decay rate for a finite system is defined as The minus sign in the previous equation follows from the fact that the real part of the eigenvalues of a Lindbladian super-operator satisfy the following inequality: (λ) ≤ 0.
Remarkably, for every eigenvalue ξ ofĤ p there is an eigenvalue λ = −ξ/2 of L which is at least twofold degenerate. Indeed, given the state |ψ ξ such that H p |ψ ξ = ξ|ψ ξ , the operators made up of the dark state |Ψ and of |ψ ξ satisfy the appropriate secular equation. This has an important consequence: ifĤ p is gapless, then λ ADR L→∞ − −−− → 0 in the thermodynamic limit, where L is the size of the system. Indeed: for every eigenvalue ξ ofĤ p ; if ξ closes as L −α (α > 0), then the dissipative gap closes at least polynomially in the system size. Note that this argument also implies that if L is gapped, then the parent Hamiltonian is gapped as well.
It is important to stress that the spectral properties of the parent HamiltonianĤ p do not contain all the information concerning the long-time dissipative dynamics. As an example, let us assume that the Markovian dynamics in Eq. (1) (i) supports at least one dark state and (ii) has an associated parent Hamiltonian which is gapped. If the Lindblad operators are Hermitian, then the fullymixed state is a steady state of the master equation too. The presence of such stationary state is not signaled by the parent Hamiltonian, which is gapped and only detects the pure steady states of the dynamics.
Whereas the some of the above relations have been often pointed out in the literature [23,24], to the best of our knowledge the remarks on the relation between the spectral properties of L andĤ p are original.

B. The Kitaev chain and the dissipative preparation of its ground states
Let us now briefly review the results in Ref. [19] and use them to exemplify how property (2) can be used as a guideline for dissipative state preparation in the number non-conserving case. This will be valuable for our detailed studies of its number conserving variant below.
The simplest model displaying zero-energy unpaired Majorana modes is the one-dimensional Kitaev model at the so-called "sweet point" [8]: where the fermionic operatorsâ The ground state has energy E 0 and is two-fold degenerate: there are two linearly independent states |ψ e and |ψ o which satisfy: The quantum number distinguishing the two states is the parity of the number of fermions,P = (−1) â † jâ j , which is a symmetry of the model (the subscripts e and o stand for even and odd ). Both states |ψ σ are p-wave superconductors, as it can be explicitly proven by computing the expectation value of the corresponding order parameter: It is thus relevant to develop a master equation which features |ψ e and |ψ o as steady states of the dynamics [19,26]. Property (11) provides the catch: upon identification of theˆ j operators with the Lindblad operators of a Markovian dynamics, Eq. (2) ensures that the states |ψ σ are steady states of the dynamics and that in the long-time limit the system evolves into a subspace described in terms of p-wave superconducting states. This becomes particularly clear once it is noticed that the parent Hamiltonian of this Markov process coincides witĥ H K in Eq. (8) apart from an additive constant.
Let us conclude mentioning that the obtained dynamics satisfies an important physical requirement, namely locality. The Lindblad operatorsˆ j only act on two neighboring fermionic modes; this fact makes the dynamics both physical and experimentally feasible. On the other hand, they do not conserve the number of particles, thus making their engineering quite challenging with coldatom experiments. The goal of this article is to provide dissipative schemes with Lindblad operators which commute with the number operator and feature the typical properties of a p-wave superconductor.

III. SINGLE WIRE: ANALYTICAL RESULTS
The simplest way to generalize the previous results to systems where the number of particles is conserved is to consider the master equation induced by the Lindblad operators [19,26]: for a chain with hard-wall boundaries and spinless fermions: {L † jL j ,ρ} ; γ > 0; (14) where γ is the damping rate. This Markovian dynamics has already been considered in Refs. [19,26]. Using the results presented in Ref. [33], where the parent Hamiltonian related to the dynamics in Eq. (14) is considered, it is possible to conclude that for a chain with periodic boundary conditions (i) there is a unique dark state for every particle number density ν = N/L, and (ii) this state is a p-wave superconductor. A remarkable point is that theL j are local and do not change the number of particles: their experimental engineering is discussed in Ref. [19], see also [51].
Here we clarify that for the master equation for a single wire with hard-wall boundaries, the steady state is not topological and does not feature Majorana edge physics, although they still display the bulk properties of a p-wave superconductor (instead, the two-wire version studied below has topological properties associated to dissipative Majorana zero modes). The asymptotic decay rate of the master equation is also characterized. An extensive numerical study of the stability of this protocol is postponed to Sec. IV.

A. Steady states
In order to characterize the stationary states of the dynamics, let us first observe that Eq. (11) implies [33] so that:L Thus, |ψ σ are steady states of the dynamics. Let us define the states whereΠ N is the projector onto the subspace of the global Hilbert (Fock) space with N fermions (Π N |ψ σ = 0 when the parity of N differs from σ and thus we avoid the redundant notation |ψ σ,N ). Since [L j ,N ] = 0, wherê N = jâ † jâ j is the particle-number operator, it holds thatL j |ψ N = 0 for all j = 1, . . . , L − 1 and thus the |ψ N are dark states. Let us show that there is only one dark state |ψ N once the value of N is fixed. To this end, we consider the parent Hamiltonian (3) associated to the Lindblad operators (13): wheren j ≡â † jâ j and J > 0 is a typical energy scale setting the units of measurement. Upon application of the Jordan-Wigner transformation, the modelĤ p is unitarily equivalent to the following spin-1/2 chain model: whereσ α j are Pauli matrices. Apart from a constant proportional to L−1,Ĥ p,spin is the ferromagnetic Heisenberg model. The particle-number conservation corresponds to the conservation of the total magnetization along theẑ direction. It is a well-known fact that this model has a highly degenerate ground state but that there is only one ground state for each magnetization sector, both for finite and infinite lattices. Thus, this state corresponds to the state |ψ N identified above; therefore, the possibility that the ground state ofĤ p is two-fold degenerate (as would be required for the existence of Majorana modes) for fixed number of fermions and hard-wall boundary conditions is ruled out.
Summarizing, the dynamics induced by the Lindblad operators in (13) conserves the number of particles and drives the system into a quantum state with the properties of a p-wave superconductor (in the thermodynamic limit |ψ e andΠ N |ψ e have the same bulk properties, as it is explicitly checked in Ref. [19,26], but see also the discussion below). Since the steady states of the system for open boundary conditions are unique, they do not display any topological edge property.

B. P-wave superconductivity
Let us explicitly check that the states |ψ N have the properties of a p-wave superconductor. Since each state has a definite number of fermions, the order parameter defined in Eq. (12) is zero by symmetry arguments. In a number-conserving setting, we thus rely on the p-wave pairing correlations: (20) If in the long-distance limit, |l − j| → ∞, the expectation value saturates to a finite value or shows a powerlaw behavior, the system displays p-wave superconducting (quasi-)long-range order. If the decay is faster, e.g. exponential, the system is disordered.
In this specific case, the explicit calculation shows a saturation at large distance (see also Ref. [33]): in the thermodynamic limit. The saturation to a finite value captures the p-wave superconducting nature of the states. Note that the breaking of a continuous symmetry in a one-dimensional system signaled by Eq. (21) is a nongeneric feature: a perturbation of HamiltonianĤ p would turn that relation into a power-law decay to zero as a function of |j − l| (see Ref. [33] for an explicit example).

C. Dissipative gap
An interesting feature ofĤ p,spin is that it is gapless; the gap closes as L −2 due to the fact that the low-energy excitations have energy-momentum relation ω q ∼ q 2 , as follows from well-known properties of the ferromagnetic Heisenberg model. The Jordan-Wigner transformation conserves the spectral properties and thusĤ p is also gapless. Thus, according to the discussion in Sec. II A, the asymptotic decay rate λ ADR associated to the Lindbladian L closes in the thermodynamic limit. This is true both for periodic and hard-wall boundary conditions. This fact has two important consequences. The first is that the dissipative preparation of a fixed-number p-wave superconductor through this method requires at least a typical time τ that scales like L 2 . In Sec. IV we numerically confirm this polynomial scaling. Although this requires an effort which is polynomial in the system size, and which is thus efficient, it is a slower dynamical scenario than that of the non-number-conserving dynamics considered in Refs. [19,26] and summarized in Sec. II B, where τ does not scale with L (the super-operator L in that case is gapped), and thus the approach to stationarity is exponential in time. The difference can be traced to the presence of dynamical slow modes related to exact particle number conservation, a property which is abandoned in the mean field approximation of Refs. [19,26].
The second consequence is that a gapless Lindbladian L does not ensure an a priori stability of the dissipative quantum state preparation. Roughly speaking, even a small perturbation M ( 1) to the Lindbladian L such that the dynamics is ruled by L + M has the potential to qualitatively change the physics of the steady-state (see Refs. [54][55][56] for some examples where the presence of a gap is exploited for a perturbative analysis of the steady states). This concerns, in particular, the long-distance behavior of correlation functions. To further understand this last point, in Sec. IV we have analyzed the effect of several perturbations through numerical simulations. In the case in which the steady state has topological properties, they may still be robust. We further elaborate on this point in Sec. V, where we study the ladder setup.
Notwithstanding the gapless nature of the Lindbladian L , we can show that waiting for longer times is beneficial to the quantum state preparation. If we define p 0 (t) = tr P 0ρ (t) , whereP 0 is the projector onto the ground space of the parent HamiltonianĤ p , then the following monotonicity property holds: which is a non-negative operator because for any state |φ it holds that: where {|ψ α } are a basis of the ground space of the parent HamiltonianĤ p . If we consider the spectral decomposition ofρ(t) = β p β |φ β φ β |, with p β > 0, we obtain Eq. (22).

IV. SINGLE WIRE: NUMERICAL RESULTS
Although the previous analysis, based on the study of the dark states of the dynamics, has already identified many distinguishing properties of the system, there are several features which lie outside its prediction range. Let us list for instance the exact size scaling of the asymptotic decay rate or the resilience of the scheme to perturbations. In order to complement the analysis of the dissipative dynamics with these data, we now rely on a numerical approach.
The numerical analysis that we are going to present is restricted to systems with hard-wall boundary conditions. In order to characterize the time evolution described by the master equation (14), we use two different numerical methods. The first is a Runge-Kutta (RK) integration for systems of small size (up to L = 10) [57]. This method entails an error due to inaccuracies in the numerical integration, but the density matrix is represented without any approximation.
On the contrary, the second method, based on a Matrix-Product-Density-Operator (MPDO) representation of the density matrix, allows the study of longer systems through an efficient approximation ofρ [40,41,43]. The time evolution is performed through the Time-Evolving Block Decimation (TEBD) algorithm, which is essentially based on the Trotter decomposition of the Liouville super-operator e tL . Although this method has been shown to be able to reliably describe problems with up to ∼ 100 sites [47], in this case we are not able to consider lengths beyond L = 22 because of the highlyentangled structure of the states encountered during the dynamics. It is an interesting perspective to investigate whether algorithms based on an MPDO representation of the density matrix, which compute the steady state through maximization of the Lindbladian super-operator L , might prove more fruitful in this context [46,49].
Finally, we have also performed Exact-Diagonalization (ED) studies of system sizes up to L = 5 in order to access properties of L , such as its spectrum, which cannot be observed with the time-evolution.

A. Asymptotic decay rate
Let us first assess that the asymptotic decay rate of the system closes polynomially with the system size (from the previous analysis we know that it closes at least polynomially. As we discuss in Appendix A, in the asymptotic limit, it is possible to represent the expectation value of any observableÂ as: where Â (t) = tr Âρ (t) , Â ss = lim t→∞ Â (t) and κ is a non-universal constant. The notation −λ ADR is due to the fact that λ ADR is positive, being defined through the additive inverse of the real part of the eigenvalues, see Eq. (4). It is possible to envision situations where κ = 0 and thus the long-time decay is dominated by eigenvalues of L with smaller real part.
The study of the long-time dependence of any observable can be used to extract the value of λ ADR ; among all the possible choices, we employ the pairing correlator G (20)] because of its special physical significance. In Fig. 1(top), we consider L = 10 and plot the time evolution of G and N = (L + 1)/2 for L odd).
In order to benchmark the reliability of the RK integration for getting the steady state, we compare the expectation value of several observables (in particular of pairing correlators) with the exactly-known results (Sec. III provides the exact wavefunction of the steady state, from which several observables can be computed). In all cases the absolute differences are below 10 −6 . Similar results are obtained for smaller system sizes, where it is even possible to compute the trace-distance of the RK steadystate from the λ = 0 eigenstate of the Liouvillian computed with ED.
In the long-time limit, the observable (20) displays a clear stationary behavior, G (p) j,l ss = lim τ →∞ G (p) j,l (τ ), consistently with Eq. (24). Once such stationary value is subtracted, it is possible to fit λ ADR from the exponential decay of The subtraction is possible to high precision because the value of G (p) j,l ss is known from the previous analytical considerations. Moreover, as we have already pointed out, the evolution continues up to times such that G (p) j,l (t) differs in absolute terms from the analytical value for 10 −6 , which makes the whole procedure reliable. In Fig. 1(bottom) we display the quantity in (25) for various lattice sizes L. It is clear that the convergence of the observable requires an amount of time which increases with L. A systematic fit of λ ADR for several chain lengths allows for an estimate of its dependence on L [see Fig. 1(bottom)]: the finite-size dissipative gap scales as The exact diagonalization (ED) of the Liouvillian up to L = 5 allows a number of further considerations. First, the Liouvillian eigenvalues with largest real part ( (λ) 0) are independent of the number of particles (the check has been performed for every value of N = 1, ..., 5). Second, comparing the ED with the previous analysis, we observe that the λ ADR in Eq. (26) coincides with the second eigenvalue of the Liouvillian, rather than with the first [here the generalized eigenvalues are ordered according to the additive inverse of their real part − (λ)]. Numerical inspection of small systems (up to L = 5) shows that the first excited eigenvalue of L is two-fold degenerate and takes the value −ξ/2, where ξ is the energy of the first excited state ofĤ p (see the discussion in Sec. II A). Our numerics suggests that it does not play any role in this particular dissipative evolution, hinting at the fact that the chosenρ(0) does not overlap with the eigensubspace relative to −ξ/2. In this case, the value of κ in Eq. (24) is zero.

B. Perturbations
In order to test the robustness of the dissipative scheme for the preparation of a p-wave superconductor, we now consider several perturbations of the Lindbladian L of both dissipative and Hamiltonian form. The robustness of the dissipative state preparation of the p-wave superconductor is probed through the behavior of the correlations G (p) j,l (t), which define such phase.

Perturbations of the Lindblad operators
Let us define the following perturbed Lindblad operator: which allows for slight asymmetries in the action of the dissipation between sites j and j+1. The continuity equation associated to the dynamics, ∂ tni = −(ĵ i −ĵ i−1 ), is characterized by the following current operator:ĵ i = n i − (1 − ) 2n j + ( 2 − 2 )n ini+1 . When = 0,ĵ i is not anymore odd under space reflection around the link between sites i and i + 1, so that in the stationary state a non-zero current can flow even if the density profile is homogeneous (and even under the previous space-inversion transformation), which is quite intuitive given the explicit breaking of inversion symmetry in this problem.
We employ the MPDO method to analyze the steadystate properties of a system with size L = 22 initialized in the ground state of the free HamiltonianĤ 0 for N = 11 and subject to such dissipation. The results in the inset of Fig. 2 show that the steady state is not homogeneous and that a relatively high degree of inhomogeneity n L − n1 n L/2 ≈ 1 is found also for small perturbations = 0.05. This is not to be confused with the phaseseparation instability which characterizes the ferromagnetic parent HamiltonianĤ p,spin . Indeed, if PBC are considered, the system becomes homogeneous and a current starts flowing in it (not shown here). P-wave superconducting correlations are affected by such inhomogeneity. Whereas for = 0 the correlations G j /( n j ss n j+1 ss ). An exponential decay behavior appears as a function of |j − l|, which becomes more pronounced when is increased. Even if the simulation is performed on a finite short system, for significant perturbations, = 0.1, the value of G (p) j,l ss decays of almost two decades, so that the exponential behavior is identified with reasonable certainty.
In Appendix B we discuss some interesting analogies of these results with the properties of the ground state of the parent HamiltonianĤ p, = J jL † j, L j, . It should be stressed that, sinceĤ p, does not have a zero-energy ground state, there is no exact correspondence between its ground state and the steady states of L .
Concluding, we mention that a similar analysis can be done introducing an analogous perturbation in the operatorĈ † j ; our study did not observe any qualitative difference (not shown).

Perturbations due to unitary dynamics
An alternative way of perturbing the dynamics of L in Eq. (14) is to introduce a Hamiltonian into the system, chosen for simplicity to be the already-introduced free HamiltonianĤ 0 : Using the MPDO method to characterize the steady state of the dynamics, we analyze the spatial decay of the pairing correlations for L = 22 and at half-filling (N = 11); the initial state is set in the same way as in the previous section. In Fig. 3 (top) we display the results: even for very small perturbations the pairing correlator G (p) 4,j ss decays rapidly in space. The long-distance saturation observed in the absence of perturbations is lost and qualitatively different from this result. In Fig. 3 (bottom) we highlight that the decay is exponential. Summarizing, in all the cases that we have considered, the p-wave pairing correlations of the stationary state G (p) j,l ss are observed to decay as a function of |j −l|. Due to the interplay between the targeted dissipative dynamics and the perturbations, which do not support a p-wave ordered dark state, the steady state is mixed, similar to a finite temperature state. From this intuition, the result is easily rationalized: Any (quasi) long range order is destroyed in one-dimensional systems at finite temperature. We note that the true long range order found in the unperturbed case (correlators saturating at large distance; opposed to the more generic quasi-long range order defined with algebraic decay) is non-generic in onedimensional systems and a special feature of our model, see [33] for a thorough discussion. However, the destruction of any such order via effective finite temperature effects must be expected on general grounds. The absence of quasi-long-range p-wave superconducting order, which in one-dimension only occurs at zero-temperature for pure state, is likely to be in connection with this fact.

Perturbation strength
Finally, we perform a quantitative investigation of the dependence of the pairing correlations on the perturbation strength, .
Lindblad perturbation -In Fig. 4 we plot the p-wave superconducting correlation G , is also included). Our data confirm that correlations undergo a clear suppression in the presence of = 0, which in one case is exponential in and in the other in 2 . The calculation is performed through RK integration of the dynamics.
Hamiltonian perturbation -We begin with the two cases:Ĥ 0 andĤ nn = −J jn jnj+1 . Fig. 5 shows, in both cases, an exponential decay to zero of G (p) 2,L−2 ss when is increased. On the contrary, a Hamiltonian which introduces p-wave correlations in the system, such asĤ Concluding, we have shown that in all the considered cases, perturbations of both dissipative and Hamiltonian form are detrimental to the creation of a p-wave superconductor. This is rationalized by the mixedness of the stationary state in that case, and parallels a finite temperature situation. In any generic, algebraically ordered system at T=0, one has gapless modes.

V. TWO WIRES
An intuitive explanation of why the dissipative setup discussed in Sec. III does not show topological dark states with fixed number of particles is the fact that this constraint fixes the parity of the state, and thus no topological degeneracy can occur. It has already been realized in several works that a setup with two parallel wires can overcome this issue [21,[29][30][31][32][33][34]. In this case it is possible to envision a number-conserving p-wave superconducting Hamiltonian which conserves the parity of the number of fermions in each wire: such symmetry can play the role of the parity of the number of fermions forĤ K in Eq. (8). Several equilibrium models have already been discussed in this context; here we consider the novel possibility of engineering a topological number-conserving p-wave superconductor with Markovian dynamics.

A. Steady states
Let us study a system composed of two wires with spinless fermions described by the canonical fermionic operatorsâ j . For this model we consider three kinds of Lindblad operators: We now characterize the dark states of the Markovian dynamics induced by these operators for a two-leg ladder of length L with hard-wall boundary conditions: (32) In particular, we will show that, for every fermionic density different from the completely empty and filled cases, there are always two steady states.
It is easy to identify the linear space S N of states which are annihilated by theL a,j andL b,j and have a total number of particles N : 0 |ψ b,N , |ψ a,1 |ψ b,N −1 , . . . , |ψ a,N |ψ  From the conditionĈ † j |ψ σ = −Â j |ψ σ we obtain: and when we impose the conditionL I,j |ψ = 0: The result is α m = α m+2 , so that two linearly independent states can be constructed which are annihilated by all the Lindblad operators in (31): The subscripts ee and oo refer to the fermionic parities in the first and second wire assuming that N is even; N N,ee and N N,oo are normalization constants [33]. For N odd one can similarly construct the states |ψ N,eo and |ψ N,oe . By construction, the states that we have just identified are the only dark states of the dynamics.
It is an interesting fact that at least two parent Hamiltonians are known for the states in (37), as discussed in Refs. [33,34]. We refer the reader interested in the full characterization of the topological properties of these steady-states to those articles.
Finally, let us mention that the form of the Lindblad operators in (31) is not uniquely defined. For example one could replaceL I,j in Eq. (31c) with the following: without affecting the results [33]. The latter operator is most realistic for an experimental implementation, as we point out below.

B. P-wave superconductivity
Let us now check that the obtained states are p-wave superconductors. Similarly to the single-wire protocol discussed in Eq. (21), the explicit calculation [33] shows that p-wave correlations saturate to a final value at large distances in the thermodynamic limit [for the two-leg ladder we consider ν = N/(2L)] This relation clearly highlights the p-wave superconducting nature of the states.

C. Dissipative gap
In order to demonstrate that the asymptotic decay rate λ ADR associated to L tends to 0 in the thermodynamic limit, we consider the parent Hamiltonian of the model: where J > 0 is a typical energy scale setting the units of measurement. This Hamiltonian has been extensively analyzed in Ref. [33]. Numerical simulations performed with the density-matrix renormalization-group algorithm assess thatĤ p is gapless and that the gap is closing as 1/L 2 . According to the discussion in Sec. II A, the asymptotic decay rate λ ADR associated to the Lindbladian L closes in the thermodynamic limit with a scaling which is equal to ∼ L −2 or faster. This is true both for periodic and hard-wall boundary conditions.

D. Experimental implementation
The Lindblad operators in Eqs. (31a), (31b) and (38) lend themselves to a natural experimental implementation. The engineering of terms likeL a,j andL b,j has been extensively discussed in Ref. [19] starting from ideas originally presented in Ref. [23]. As we will see, the Lindblad operatorL I,j in Eq. (38) is just a simple generalization.
The idea is as follows: a superlattice is imposed which introduces in the system additional higher-energy auxiliary sites located in the middle of each square of the lower sites target lattice. Driving lasers are then applied to the system, whose phases are chosen such that the excitation to the auxiliary sites happens only for states |ϕ such that (Â a,j +Â b,j )|ϕ = 0. If the whole system is immersed into, e.g., a Bose-Einstein condensate reservoir, atoms located in the auxiliary sites can decay to the original wire by emission of a Bogoliubov phonon of the condensate. This process is isotropic and, for a wavelength of the emitted phonons comparable to the lattice spacing, gives rise to the four-site creation part with relative plus sign:Ĉ † a,j +Ĉ † b,j .

E. Perturbations
An important property of topological Hamiltonians is the robustness of their edge physics to local perturbations. Similar features have been highlighted in the case of topological superconductors where the setup is not number conserving [19,26]. The goal of this section is to probe the resilience of the twofold-degenerate steady states of L . A conclusive analysis is beyond our current numerical possibilities; here we present some preliminary results obtained via exact diagonalization methods.
We consider the natural choice of Lindblad operators Eqs. (31a,31b,38), subject to perturbations: Those define a perturbed Lindbladian L . They are a simple generalization of those defined in Eq. (27) for the single-wire setup. Let us begin our analysis by showing that for small sizes L ∼ 6 the degeneracy of the steady space for = 0 is broken. Let us first remark that for = 0 the steady space is four-fold degenerate; a possible parameterization is: |ψ N,oo ψ N,ee |, |ψ N,oo ψ N,oo |}.
A direct inspection of the eigenvalues of L shows that this degeneracy is broken once = 0. Results, shown in Fig. 6 for a fixed lattice size L = 6 and N = 6, display a quadratic splitting of the steady steady degeneracy with the perturbation strength.
Let us now check the behavior with the system size of the first eigenvalues of the system for longer system sizes. In order to obtain a reasonable number of data, the extreme choice of setting N = 2 in all simulations has been taken, which allows us to analyze system sizes up to L = 20. Results shown in Fig. 7 (top) show that the Liouvillian eigenvalues related to the steady-state degeneracy display an algebraic scaling λ ADR ∼ L −1 in the accessible regime of system sizes for small perturbations ( = 10 −2 ), while they are gapped for larger perturbations ( = 10 −1 ). Note that, for the system sizes which could be accessed, larger eigenvalues clear display an algebraic decay, as shown in Fig. 7 (bottom), also for = 0.1. The scaling of the eigenvalues related to the steady state degeneracy is not exponential and thus in principle should not be connected to the topological properties of the system. However, these preliminary considerations suffer from two significant biases: (i) the small considered sizes, (ii) the fact that they are not performed at exactly fixed density, and (iii) the very low filling. A more thorough analysis is left for future work.

VI. CONCLUSIONS
In this article we have discussed the dissipative quantum state preparation of a p-wave superconductor in one-dimensional fermionic systems with fixed number of particles. In particular, we have presented two protocols which have been fully characterized in the presence of hard-wall boundaries. Whereas the former does not display topological property, the latter features a twodimensional steady space to be understood in terms of boundary Majorana modes for any number of fermions. Through the analysis of a related parent Hamiltonian, we are able to make precise statements about the gapless nature of the Lindbladian super-operators associated to both dynamics.
The peculiar form of the master equations considered in this article allows for the exact characterization of several properties of the system, and in particular of the steady state, even if the dynamics is not solvable with the methods of fermionic linear optics [59,60] exploited in Refs. [19,26]. This result is very interesting per se, as such examples are usually rare but can drive physical intuition into regimes inaccessible without approximations. It is a remarkable challenge to investigate which of the properties presented so far are general and survive to modifications of the environment, and which ones are peculiar of this setup.
Using several numerical methods for the study of dissipative many-body systems, we have presented a detailed analysis of the robustness to perturbations of these setups. Through the calculation of the proper p-wave correlations we have discussed how external perturbations can modify the nature of the steady state. In the ladder setup, where the steady states are topological, we have presented preliminary results on the stability of the de-generate steady-space of the system.
The analysis presented here has greatly benefited from exact mathematical relations between the properties of the Lindbladian and of a related parent Hamiltonian. Since the study of closed systems is much more developed than that of open systems both from the analytical and from the numerical points of view, a more detailed understanding of the relations between Lindbladians and associated parent Hamiltonian operators stands as a priority research program.
In order to discuss the long-time properties of the dissipative dynamics, it is convenient to start from the spectral decomposition of the Lindbladian. Since L is in general a non-Hermitian operator, its eigenvalues are related to its Jordan canonical form [58]. Let us briefly review these results. The Hilbert space of linear operators on the fermionic Fock space, H, can be decomposed into the direct sum of linear spaces M j (usually not orthogonal) such that if we denote with P j the projectors onto such subspaces (usually not orthogonal) and with N j a nilpotent super-operator acting on M j , the following is true: The {λ j } are the generalized complex eigenvalues of the super-operator L and, for the case of a Lindbladian, have non-positive real part; the N j can also be equal to zero. By this explicit construction it is possible to observe that the {P j } and {N j } are all mutually commuting (P j P k = δ j,k P j , P j N k = N k P j = δ j,k N j and N j N k = δ j,k N 2 j ). Using these properties, the time evolution can be writ- which highlights that at a given time t only the terms of the sum such that | (λ j ) t| 1 play a role. In the long-time limit, it is possible to represent the expectation value of any observableÂ as: (A3) Eq. (A3) is the mathematical formula motivating Eq. (24) in the text, defining also the meaning of κ.
Let us mention that in the example discussed in the text N ADR = 0: this is observed by explicit inspection via exact diagonalization of small systems (L = 5). Since the presence of a non-zero nilpotent super-operator is a fine-tuned property, it is reasonable to assume that the situation remains similar for longer systems.

Appendix B: Analogies with the parent Hamiltonian
In this Appendix we discuss some interesting analogies between the steady stateρ ss of the dissipative dynamics for the perturbed Lindblad operatorL j, in Eq. (27) with the ground state |g of its parent HamiltonianĤ p, = J jL † j, L j, . It should be stressed that, sinceĤ p, does not have a zero-energy ground state, there is no exact correspondence between both states.
We first study a small lattice with L = 8 sites at half-filling, performing a Runge-Kutta integration of the master equation. The initial state of the evolution is the ground state ofĤ 0 . In Fig. 8 it is shown that both the purity of the steady state P(ρ ss ) = tr ρ 2 ss and its fidelity with the ground state of the parent Hamiltonian decrease with the perturbation strength. Notice, however, that for small perturbations the fidelity F(ρ ss , |g ) = g|ρ ss |g remains close to one, thus revealing the similarity of the states in such regime.
Such feature is also observed for larger lattices. Using the MPDO method forρ ss and an algorithm based on matrix product states for |g , we analyze a lattice with L = 22 sites at half-filling. We compare the pairing correlations and density profiles for both states, which differ only for O(10 −2 ), when the perturbation strength is 0.05 (not shown). Let us explicitly show the results for the Hamiltonian case. In Fig. 9 we show that, for a lattice with L = 40 sites at half-filling, even a small perturbation ( ∼ 10 −3 ) produces a non-negligible inhomogeneity. Moreover, the pairing correlations decay, indicating that such perturbation breaks the p-wave ordered nature of the purely dissipative dark state.
This similarity encourages the possibility of accessing some steady-state properties for large lattices through the study of the ground states of the corresponding parent Hamiltonians, even if no mathematical connection is present and the mixedness of the state is expected to act like a finite temperature, washing out several groundstate properties.