Exact solution of a quantum asymmetric exclusion process with particle creation and annihilation

We consider a Lindblad equation that for particular initial conditions reduces to an asymmetric simple exclusion process with additional loss and gain terms. The resulting Lindbladian exhibits operator-space fragmentation and each block is Yang-Baxter integrable. For particular loss/gain rates the model can be mapped to free fermions. We determine the full quantum dynamics for an initial product state in this case.


Introduction
Whilst most standard tools for many-body quantum mechanics only apply to closed systems, real systems are invariably influenced by their environment. Under a Markovian approximation, an effective description of this interaction can be given in terms of the Lindblad equation [1,2,3] for the evolution of the density matrix. The main approaches that have been used in the literature to study Lindblad equations for many body systems are either perturbative [4,5] or numerical [6,7,8,9,10]. Given that solvable models have provided deep insights into the non-equilibrium dynamics of closed many body systems [11,12,13,14] it is natural to ask if there are any exact results that can be obtained for many particle Lindblad equations.
The first step in this direction was the realisation that certain Lindblad equations can be cast in the form of imaginary-time Schrödinger equations with non-Hermitian "Hamiltonians" that are quadratic in fermionic or bosonic field operators [15], which then can be analyzed by standard methods for free theories to extract physical properties [16,17,18,19,20,21,22,23]. A characteristic feature of these models is the fundamental boson or fermion operators fulfil linear equations of motion and concomitantly so do the Green's functions of interest. Another step towards obtaining exact solutions of many particle Lindblad equations was the discovery that there exist classes of models in which some or all local correlation functions satisfy closed hierarchies of equations of motion [24,25,26,27,28,29,30]. This permits one to obtain some exact results on the dynamics although full solutions typically remain out of reach. Another class of solvable Lindblad equations are "triangular" models which add particle loss and dephasing terms to otherwise number conserving integrable models [31,32,33]. Recently a new direction for constructing solvable many particle Lindblad equations was identified through the discovery of Lindblad equations that can be related to interacting Yang-Baxter integrable models [34,35,36,37,38,32,33,39,40,41,42]. The approach of Refs [34,38] is based on a superoperator representation of the Lindblad equations, which gives rise to solvable "two-leg ladder" quantum spin chain models. Importantly the equations of motion for correlation functions do not generally close in these models but form an infinite hierarchy of coupled nonlinear equations. More recently a method for constructing Yang-Baxter integrable Lindblad systems was developed [42].
A related but different route of constructing Yang-Baxter integrable Lindblad equations was discovered in Ref. [43]. It is based on a "fragmentation" of the space of operators into an exponential (in system size) number of subspaces that are left invariant under the dissipative evolution. Importantly, this mechanism applies to the quantum version of the simple asymmetric exclusion process (ASEP) [44,45,46,47]. The corresponding Lindblad equation can be obtained [48] as the averaged dynamics of a stochastic quantum model of particles hopping with random amplitudes first introduced in its symmetric form in Ref. [49] and further analyzed in Refs. [50,51,52,53]. In this case, it was shown in [43] that the Lindbladian restricted to each invariant subspace can be mapped onto an XXZ Heisenberg Hamiltonian with integrable boundary conditions. In particular, in the subspace of diagonal density matrices the model reduces to the classical ASEP, which is exactly solvable [54,55] and for which many exact results have been derived using integrability methods [56,57,58,59,60,61,62]. While [43] established the integrability of the quantum ASEP in each fragmented sector, the full solution of the dissipative quantum dynamics remains an open problem except in the special case of the quantum symmetric simple exclusion process. In order to show how the operator-space fragmentation can be exploited in practice to obtain a full solution of the dissipative dynamics we here consider a generalization of the quantum ASEP. As we will show, the Lindbladian of this model exhibits operator-space fragmentation and in each sector can be mapped onto a Lindbladian that is quadratic in fermions. The resulting dynamics can then be solved exactly.
The rest of this paper is organised as follows: in Section 2 we introduce the model of interest, which can be viewed as an ASEP with additional loss/gain terms, and the model exhibits operator-space fragmentation, with each subspace labelled by a sequence of "defects". In Section 3 we analyse the Lindbladian's projection on to each of these subspaces. We then focus on a particular line in parameter space, on which the Lindbladian in each sector can be mapped onto a bilinear form in auxiliary fermions. We show that the subspace of diagonal density matrices is invariant under time evolution and reduces to a classical stochastic process similar to ones that have been previously studied in the literature [63,64,65,66]. We employ Jordan-Wigner and Bogoliubov transforms to solve the dynamics in this sector and show that it has an infinite temperature steady state. In Section 5 we consider the defect problem and outline how to efficiently find the spectrum of the Lindbladian. In section 6 we consider evolution out of an initial product state and compute the transverse spin-spin correlation function. Lastly, we relegate some technical calculations necessary for the conclusions in the main text to two appendices.

Lindblad equation
For a system interacting with its environment, the Lindblad equation for the time evolution of the reduced density operator of the system ρ is given by where the jump operators L a describe the interactions of the system with the environment, J a are the corresponding rates and {·, ·} denotes an anti-commutator. The Lindblad equation (1) describes the time evolution of the system degrees of freedom after averaging over Markovian bath degrees of freedom [3]. In order to study the fluctuations of system degrees of freedom that are induced by coupling to the bath -a question that has been extensively studied for classical systems (see e.g. [67,68,69,70,71,72,73]) -it is necessary to go beyond this description, see e.g. [49,48,52], but these fluctuations can still be described in terms of a quantum master equation of Lindblad form [74]. In contrast, quantum measurement noise [75,76,77,78,79,80,81] is captured by the description (1). Since (1) is manifestly linear in ρ it can be recast in terms of a Lindblad superoperator that generates time evolution in the same way the Hamiltonian does in closed quantum systems, with the major difference being the time evolution need no longer be unitary. That is, there exists a (super)operator L acting on the vector space of linear operators on the Hilbert space such that Here we have written |ρ to stress that we are considering ρ as a vector in a larger vector space whose dimension is the square of that for the original Hilbert space. In this work we consider an open spin 1/2 chain with periodic boundary conditions and no coherent dynamics (H = 0) described by four jump operators [82] L (1) In terms of Jordan-Wigner fermions the first two of these correspond to hopping left and right, whilst the latter two represent pair creation and annihilation on neighbouring sites respectively. In general the rates of these may all be different and one obtains a four parameter family of models [82]. The case J 3 = J 4 = 0 reduces to the quantum ASEP [49,48,43]. In contrast to the latter case the additional jump operators describe processes that violate spin rotational invariance around the z-direction (or equivalently particle number conservation at the level of Jordan-Wigner fermions) so that the magnetization is no longer conserved. As we will see this leads to interesting new effects compared to the ASEP. The Lindblad equation (1) with jump operators (3) can be obtained by coupling our quantum spins across each bond of our chain to an environment modelled by appropriate quantum Brownian motions as in [48] and then averaging over the bath degrees of freedom. Our choice of model is not motivated by any particular experimental setup, but aims to address a problem in mathematical physics, namely to obtain a many-particle Lindblad equation exhibiting operator-space fragmentation that can be solved exactly in practice. Having said this, in a particular parameter regime and for diagonal initial density matrices our model reduces to a classical master equation that has been argued to describe the kinetics of excitons in certain polymers [63] and it would be interesting to investigate whether quantum effects could be relevant to this system. In order to recast the Lindblad equation (1) with jump operators (3) in the superoperator formalism we note that the density operator is expressed as Then right multiplication by an operator L must turn into left multiplication by some superoperator L R such that which implies that L R = 1⊗ L β γ |β γ| . Note that this has indices swapped compared to 4, indicating that the right multiplication action is implemented via the transpose of the original operator, along with acting on bras instead of kets. Left multiplication is simply implemented via the operator acting on kets. This can be summarised as L L = L ⊗ 1 and L R = 1 ⊗ L T . In order to obtain an explicit expression for the Lindblad superoperator L we pick the following basis of the local Hilbert space of operators acting on site j A basis of local superoperators acting on these states in then given by For convenience, we split the Lindbladian up as where L Diag leaves invariant the subspace of diagonal density matrices These diagonal density matrices correspond to classical probability distributions. We have If we initialize the system in a purely diagonal density matrix the Lindblad equation (1) reduces to a classical master equation with transition matrix L Diag . This describes generalizations of the asymmetric simple exclusion process [44,45,46,47] similar to the diffusion-annihilation models studied in [63,64,65,66]. If we set J 3 = J 4 = 0 we recover the ASEP with periodic boundary conditions.

Operator-space fragmentation
The origin of operator-space fragmentation in the model (8), (10), (11) is the presence of strictly local conservation laws These conservation laws imply that particles of species 2, 3 are left invariant by the dynamics and we therefore refer to these as "defects". The Hilbert space of operators thus breaks up into exponentially many invariant subspaces with fixed occupancies of defects. This is somewhat reminiscent of the Hilbert space fragmentation found in certain fractonic circuits [83,84,85]. The fragmentation of operator-space does not rely on the fact our model is one dimensional. Indeed, operator-space fragmentation occurs if we consider a square lattice and jump operators defined on all nearest neighbour bonds In this case the 2L 2 operators E 22 (i,j) , E 33 (i,j) are then strictly conserved. By focusing on one dimensional models however we allow for the possibility that the Lindbladian's action on each subspace can be mapped to an integrable model. However, the ocurrence of fragmentation will have implications for the dynamics in higher dimensions as well.
This operator-space fragmentation then allows observables to be computed by analyzing each sector separately. In the case of the ASEP (J 3 = J 4 = 0), the key result is that restricted to each defect-subspace the Lindbladian can be mapped to a collection of disjoint finite XXZ chains with diagonal boundary fields and is thus integrable on every subspace. Integrability techniques can be similarly applied to (8), (10), (11) [66,82] but we do not pursue this line of enquiry here and instead impose a particular constraint on the rates J 1 , . . . , J 4 which will allow us to employ mappings to free fermion systems (see below).
It should be stressed that for a particular observable, it may not be necessary to deal with very many invariant subspaces. This is illustrated by the transverse spin-spin correlation function This depends only on the subspace with a type 3 defect at site 0 (equivalently site L) and a type 2 defect at site + 1. To see this, note that in the superoperator formalism traces are replaced by inner products with the state An immediate consequence of the fact that the time evolution operator e Lt preserves traces is that 1 cl | is a left eigenvector of the time evolution operator with eigenvalue 1.
If there is a unique steady state of the system then it is also the only left eigenvector with this property, a fact that we will use later. The spin operators act by left multiplication so in the superoperator formalism they are mapped to Since 1 cl | only contains states 1, 4 the only terms that survive in the trace are then where we have introduced Eqn (16) shows that the correlation function depends only on the projection of |ρ(t) on the single subspace described above. This means the correlation function can be written in terms of propagators defined on open chain segments: In the ASEP case these propagators involve computing the overlap of a time evolved state in the finite length XXZ model (with diagonal boundary fields) with the state 1 cl |. The rest of this paper will consider a different subspace of the full four-parameter model which reduces to free fermions, thus allowing the calculation of G [a,b] for some initial states, although its calculation for general states is still difficult .

"Classical" sector
As we noted earlier, the subspace of diagonal density matrices (9) is invariant under the dynamics. The "classical" part L Diag of the Lindbladian acts on this 2 L dimensional subspace of diagonal density matrices and can be expressed in terms of Pauli matrices τ j defined by We find We now observe that under the constraint the model (21) can be mapped to a free fermionic theory by means of a Jordan-Wigner transformation. In the periodic case we use that whereN is the total fermion number operator. Since each term in the Lindbladian preserves fermion parity, the operator (−1)N is conserved and we can work in definite parity sectors where it equals +1 (periodic, or Ramond, boundary conditions) or −1 (anti-periodic, or Neveu-Schwarz, boundary conditions). It will furthermore be convenient in the following to define in terms of which the Lindbladian can be written (defining c L+1 = (−1)N c 1 ) as We largely focus on the special case µ = 0 in the following but do discuss the steady state in the imbalanced case in Section 4.1. Crucially the constraints (24) enforce that J 3 + J 4 = 0, which takes us away from the ASEP limit J 3 = J 4 = 0. Hence the exact solutions presented below cannot be related to known results for the ASEP.

Two defect sector
We now consider the case where there are two defects that without loss of generality can be taken to be located at positions + 1 and L. Inspection of (8), (10), (11) shows that on the corresponding subspace the Lindbladian takes the form where and the constant c is given by

q defect sector
The Lindbladian for the entire chain restricted to the invariant subspace with q defects at locations 1 , . . . q is simply a sum of Lindbladians for the q disjoint finite chains obtained by removing the sites j from the original chain Here 0 = q so that for instance in the 1 defect sector the corresponding Lindbladian L [ +1, −1] corresponds to the original ring with a single site removed. If the defects j , j+1 are not immediate neighbours then these are exactly as given in Eq (28). If there are two neighbouring defects then the only term in the full Lindbladian that acts on them is which contributes c = −2J + if the neighbouring defects are different species and 0 if they are the same.

Dynamics in the classical subspace
As a first step to understanding this model we solve it exactly in the diagonal subspace. We focus initially on the balanced (µ = 0) case. Our system has periodic boundary conditions in terms of the original spins and (anti)-periodic boundary conditions for Jordan-Wigner fermions in sectors of (even) odd fermion parity. We therefore go to Fourier space where δ = 0, 1/2 for states with odd or even fermion parity respectively. We then carry out a Bogoliubov transformation to diagonalize the Lindbladian Despite the fact that L is non-Hermitian, this transformation is still unitary. We have where the non-Hermitian nature of the Lindbladian presents through the complex eigenvalues The time-evolved operators are We can now immediately conclude that the stationary state is unique and given simply by the Bogoliubov vacuum b k |0 = 0.
This implies that 0|L = 0 and exploiting uniqueness we therefore have This is turn shows that the stationary state |0 is the completely mixed (infinite temperature) state, which we now demonstrate in more detail. An important question is what operators of the original spin-chain problem can have finite expectation values within the defect-free subspace. To answer this we project the original Pauli matrices σ j on to the diagonal subspace and write the result in terms of the τ j operators. Defining projection operators by we have This shows that the only physical operators with non-zero expectation in the stationary state are O j1,...,jn = n j1 . . . n jn .
The expectation value of O j1,...,jn can be obtained using Wick's theorem with the help of the elementary two-point functions Here we have replaced the state 1 cl | used to compute traces with the left Bogoliubov vacuum 0| following the discussion above. As a result, we find that all such expectations factorise We now make use of the fact that a density operator is fully determined by the expectation values of a complete set of operators to conclude that in terms of the original problem the stationary state is the infinite temperature state

Imbalanced loss and gain
We now briefly discuss the nature of the steady state with imbalanced loss and gain. When µ = 0 the Lindbladian in the defect-free sector is given by where the appropriate Ramond or Neveu-Schwarz boundary conditions are assumed. We make use of the translational symmetry in the defect-free problem to Fourier transform this to give where c k = c k c † −k T the matrix A k (µ) is 2 × 2 and non-Hermitian Here we have introduced the dimensionless parameters The parameter ν satisfies −1 ≤ ν ≤ 1 where the extreme case of ν = −1 corresponds to only particle loss and ν = 1 to only gain. In terms of ν the eigenvalues of A k (µ) are ± k (µ) = 2J + (i∆ sin k ± (1 + ν cos k)) .
For ν = ±1 these are always distinct. Degerate eigenvalues only occur for ν = 1, k = π and ν = −1, k = 0 which both yield A k (µ) = 0. A k (µ) is thus always diagonalisable, however it is not unitarily diagonalisable if µ = 0. In this case we cannot perform a canonical transformation as for the balanced case. We can however perform an almost canonical transformation by defining the matrix chosen such that S −1 AS is diagonal and det(S) = 1. We then define These are almost canonical fermions in that they satisfy the relations We note that b +,−k = −b −,k due to the choice of normalisation in the definition of S k , which allows us to consistently define This then allows us to write the Lindbladian in terms of the almost canonical fermion operators as where we have used that + k = − − −k . The constant can be seen to be 0 by carefully keeping track of the constants discarded throughout this argument. We can now define left and right vacua by Since L|L = 0 has only one solution, we conclude that where |0 0 is the fermionic vacuum state The expression for L| in terms of the original fermions in (56) is easily verified by acting with b k and using (53). The right eigenstate can be expressed as a squeezed state via where N is chosen such that L|R = 1. This can be verified by acting with b k and using (53). As before we may consider the expectation values of all operators in the classical subspace -the operators σ x , σ y project to zero and all physical operators are given in terms of fermions by products of densities O j1...jr = n j1 . . . n jr .
The simplest such expectation value is As we have seen above, for µ = 0 the steady state corresponds to a completely mixed state and we have n j ∞ = 1/2. When µ = 0 and assuming that the steady state is uncorrelated we have the following relation expressing the balance between particle gain and loss If this relation holds we may solve for the particle density where we have used (47). We now verify (62) by direct calculation. Due to translation invariance n j ∞ is the same on each site and we can instead calculate the average occupation of the k modes In the thermodynamic limit this turns into an integral This indeed agrees with (62). We note that the stationary state (58) has a simple product form in terms of the spin states |1 , |4 on which the spin operators τ α j act, cf. (20). where

Time dependence
We return to considering only the balanced case of µ = 0 and now consider the time dependent problem. As we have seen above, on the diagonal subspace we have This allows us to identify where j k are the positions of down spins ordered such that j 1 < j 2 · · · < j n and |0 0 is the fermionic vacuum, which is related to the Bogoliubov vacuum state by Using this and an initial density matrix in our subspace we can calculate We can thus compute the expectations of any observables in this subspace at arbitrary times using free-fermion techniques. As an example we now compute n j (t) for a system initially in the classical Néel state In terms of fermions this can be written as In practice it will be more useful to work with the original fermion operators than the Bogoliubov ones. Solving their equations of motion gives where f (k, t) = cos 2 (k/2)e (−k)t + sin 2 (k/2)e − (k)t , h(k, t) = cos 2 (k/2)e − (k)t + sin 2 (k/2)e (−k)t .
This then allows us to write where we have definedf A straightforward calculation then gives The time evolution of the particle density (78) is shown in Fig. 1. As we are working at balanced particle creation and annihilation (J 3 = J 4 ), n j (t) relaxes to 1/2 at late times for all values of j.

Two defect sector
We now turn to the defect sector problem. The Lindbladian in the sector with q defects can be written as a sum over quadratic open chain Lindbladians of the form As these are not Hermitian the standard analysis of Lieb, Schultz and Mattis [86] for diagonalizing Hamiltonians quadratic in fermionic creation/annihilation operators does not apply. We therefore proceed as in Section 4.1, but find it advantageous to switch to Majorana fermions [15] In terms of the Majorana operators L M is expressed as where here and elsewhere (·) represents the dot product with no complex conjugation- Assuming A to be diagonalizable, anti-symmetry ensures its eigenvalues come in pairs ±β j which we order as β 1 , −β 1 . . . . We then normalize the eigenvectors according to v r · v s = (σ x ⊗ 1) rs .
In fact, the complex eigenvalues also come in complex conjugate pairs. This can be seen by noting that one obtains A * from A by conjugating C by σ z . In particular this means that if A v = β v then also Finally we define new fermion operators by These fulfil simple anticommutation relations due to (83) and diagonalise the Lindbladian For the matrix A in our problem it is not a simple matter to find a closed form analytic expression for the spectrum, but we can gain insight into what the solutions look like by deforming our Lindbladian by adding a boundary term J − (n 1 − n L ). We stress that the resulting Lindbladian is unphysical. Then the matrix A is modified to A It is now straightforward to obtain the eigenvectors of the matrix A . We make an For this to be an eigenvector we require A z = z 2 and z to satisfy The associated eigenvalues are then given by This only gives rise to M linearly independent eigenvectors, all with non-negative eigenvalues. We however get the full spectrum using this ansatz by reflecting in the imaginary axis. Thus in this case we find that the positive real part eigenvalues consist of M − 1 values of z that are roots of unity z = e −ik which recovers the periodic boundary condition result. There are also two eigenvalues that are exactly 0 -for our actual boundary conditions these become two small real eigenvalues ±λ 0 (they cannot be complex as the requirement that β * is an eigenvalue would then give four nearly zero eigenvalues which is too many). We plot the eigenvalues in the complex plane in Figure 2 for both A and A to highlight the impact of removing the boundary potential.

Transverse correlation function
We now turn to observables that involve defects. We focus on the particular example of an initial product state with ferromagnetic order along some direction in spin space Our aim is to determine As we showed above in (17) this involves only the projection of |ρ(t) onto the subspace with two defects Figure 2. Eigenvalues of A' (red diamonds), A (blue triangles) for L = 12, J + = 1.0, J − = 0.9. Central inset is magnified such that both axes run from ±3 × 10 −5 .
Here we have defined γ = |α| 2 since all quantities depend only on γ and have separated out an overall factor γ(1 + γ) N for convenience. The propagatorsG are defined on the finite chains discussed in Section 5 and can be expressed in terms of fermions as As shown in Appendix A we can rewrite this as where Using fermion parity conservation this simplifies tõ The two terms above can be written in the form where e Y ,e L N t , e γ 2 Y † are all manifestly Gaussian as are theρ (α) since they are the ground states of the quadratic Hamiltonians Thus (103) is now in the form of the trace of a product of Gaussian operators and can be evaluated. The procedure for this is given in detail in Appendix A. Here we outline the two key steps to the evaluation. The first step is to realise that since a product of Gaussian operators is Gaussian, we have Here, Z (α) and Z(W (α) ) are two different normalisation factors defined each defined such that Tr ρ (a) = 1. The Z (a) are calculated in Appendix B and given by B.8. Writing the time evolution operator in the form we then obtain the following expression for the propagators (cf. Appendix A ) The second step is to use the fact that a Gaussian is determined by its second moments to change from working with the density matrix e a·W ·a/4 itself to instead working with its correlation matrix Γ mn = Tr[ρa n a m ] − δ mn . We calculate the latter in Appendix B by rewriting the trace as an inner product which can be evaluated in terms of Jordan-Wigner spins. Once Γ is found W is obtained through Γ = tanh W 2 , or equivalently This then leads to an apparent difficulty since the correlation matrices corresponding to ρ (α) (which are fixed through our choice of initial condition) satisfy (Γ (α) ) 2 = 1, implying that they have only eigenvalues equal to ±1 and the method set out above appears to break down. This issue can be dealt with by noting that where Π (α) ± are projectors onto the two N dimensional subspaces corresponding to eigenvalues 1 and −1 of Γ (α) respectively. These would correspond to eigenvalues ±∞ in W , which we regulate by setting them equal to ±Λ and taking the limit Λ → ∞ in the end of the calculation. That is, we put: which simplifies (107) to read This yields a simple expression for the propagator Substituting this into (98) then gives the transverse correlation function The determinants d (α) N can now be straightforwardly computed numerically. In Fig. 3 we plot the transverse correlator at separation 2, S +− 0,2 (t), as a function of time. We observe that the correlator decays quite quickly, and monotonically, from its initial value. We note however that this is the full correlation function and that more physically interesting is the connected correlator Here we have use the translation invariance of our initial condition to express the correlation function in terms of only the distance between the defects (note that this is + 1 and not ). The one point functions depend on the same propagators as the 2-point functions since which gives Where we have expressed this in terms of G N = (1 + γ) −NG N as this is more natural. In particular, since our initial state was a product state we have G N (0) = 1 for all N and so the connected correlation function is initially 0, indicating no correlations. We then expect that the Lindblad evolution will correlate neighbouring sites. This is countered by the fact that the steady state values of observables are all governed by the diagonal subspace values and so the connected correlations must go to zero at long times. In Figure 4(a) we plot the connected correlation between sites 1 and 3. We are able to observe that the dissipative dynamics does produce some correlations although they are small. Given that they also exponentially decay, the correlation generation would most likely not be visible had we started in an initially correlated state. In Figure 4(b) we plot the corresponding values for varying site separations and note an approximately exponential decrease with distance. We perform these calculations for total chain lengths of L = 30, one might wonder if this is large enough to be essentially in the thermodynamic limit (in the sense that finite size effects are small enough to neglect). In fact, we find that the numerical values of the connected correlator vary very little as we increase L so long as it is larger than twice the separation + 1. To show this we plot the connected correlator for = 3 for L = 8, 9, 10 in Figure 5(a). Since the difference between the result for 9, 10 is too small to be visible, we plot the residual (along with the corresponding residual for L = 10, 11) in Figure 5

Conclusions
We have considered a dissipative many-particle quantum system described by a Lindblad equation that for particular initial conditions reduces to an asymmetric simple exclusion process with additional pair creation and annihilation terms. The Lindbladian exhibits operator-space fragmentation and for particular pair creation/annihilation rates the model can be mapped to free fermions. The model thus extends the class of solvable Lindblad systems and in particular provides a concrete example of a setting where operator-space fragmentation can be used to compute correlation functions exactly. We have restricted attention to initial product states in order to make calculations simpler as well as to allow us to see the generation of correlations through dissipation in our model. Even though the initial states we have considered here are quite simple the analysis is not straightforward. It would be interesting to attempt to generalize our analysis to the case of entangled initial states. It would also be interesting to study the time evolution of entanglement measures such as entanglement negativity within this model.
To arrive at Eq (100) the key identity is that for any collection of mutually anticommuting variables {ζ i } N i=1 the following holds This identity immediately provides a convenient decomposition into even and odd fermion parity parts. It can be proven by focussing on the odd and even components and using induction. To do so note that the even terms have the form The counterpart for the odd terms is completely analogous. When multiplying this by N j=1 ζ j the result will be a sum of (N − k) N k non-zero terms, each of which contain k + 1 distinct ζ's. Moreover, in each term k of the ζ's will be in ascending order by construction, with the final one appearing in each possible position. Thus k/2 pairs will cancel and the remaining (N −k) k+1 N k = (N − k − 1) N k+1 terms are precisely those in the corresponding expansion of the odd part. We thus need only prove by induction the statement about the even terms To this end we assume the induction hypothesis up to N − 1 and notice that for N sites we can rewrite the product of quadratic terms as We then use the induction hypothesis on the first factor on the right hand side. When multiplied by the second factor two things can happen: (i) it gets multiplied by 1, thus generating all possible even terms not including ζ N , (ii) it gets multiplied by ( For the latter note that multiplying by ζ ω generates all possible odd expressions without ζ N and multiplying by ζ N at the end then gives the desired result. Along with the observation that the base case of N = 0 is trivial this completes the proof of (A.1).
In the context of Eq (100) we set ζ i = γc † i so that we have This is the desired simplification upon defining X, Y in the main text and applying the standard result that e Y † = m<n (1 + c † m c † n ) where Y † = m<n c † m c † n .

Appendix A.2. Trace of Gaussian operators
The main result required in the derivation of (107) is the identity [87] Tr[e a·W ·a/4 ] = det(e W/2 + e −W/2 ). (A.6) Here a j are Majorana fermions and W an antisymmetric matrix. We note that (A.6) is easy to establish if W is diagonalizable. In that case its eigenvalues come in pairs ±β k and the left hand side becomes Because W is anti-symmetric, its eigenvalues come in ± pairs and so the determinant contains exactly two copies of each factor on the right hand side of (A.7). This then establishes (A.6). If W is not diagonalizable then we define

Appendix B. Correlation matrices
The two correlation matrices we need are given by the inner products The denominators are equal to the normalization factors Z (1) , Z (2) that appear in the final result (113). Both numerators and denominators can be found by making use of (A.1) in the form A simple calculation then gives that g ++ mn = g −− mn and g +− mn = g −+ mn and likewise for Z σσ . Explicit expressions for g σσ mn and Z σσ are readily obtained by reverting to their respective representations in terms of spins (i.e. undoing the Jordan-Wigner transformation). We find The 2 × 2 blocks are given by where we have defined One can verify that (Γ (α) ) 2 = 1. Since Γ (α) is anti-symmetric it therefore has equal numbers of eigenvalues ±1, which we used in the main text. The eigenvectors depend on x, y and we find these numerically to determine the correct projectors to use.