Exact matrix product solutions in the Heisenberg picture of an open quantum spin chain

In recent work Hartmann et al [Phys. Rev. Lett. 102, 057202 (2009)] demonstrated that the classical simulation of the dynamics of open 1D quantum systems with matrix product algorithms can often be dramatically improved by performing time evolution in the Heisenberg picture. For a closed system this was exemplified by an exact matrix product operator solution of the time-evolved creation operator of a quadratic fermi chain with a matrix dimension of just two. In this work we show that this exact solution can be significantly generalized to include the case of an open quadratic fermi chain subjected to master equation evolution with Lindblad operators that are linear in the fermionic operators. Remarkably even in this open system the time-evolution of operators continues to be described by matrix product operators with the same fixed dimension as that required by the solution of a coherent quadratic fermi chain for all times. Through the use of matrix product algorithms the dynamical behaviour of operators in this non-equilibrium open quantum system can be computed with a cost that is linear in the system size. We present some simple numerical examples which highlight how useful this might be for the more detailed study of open system dynamics. Given that Heisenberg picture simulations have been demonstrated to offer significant accuracy improvements for other open systems that are not exactly solvable our work also provides further insight into how and why this advantage arises.

2 times. Through the use of matrix product algorithms the dynamical behaviour of operators in this non-equilibrium open quantum system can be computed with a cost that is linear in the system size. We present some simple numerical examples that highlight how useful this might be for the more detailed study of open system dynamics. Given that Heisenberg picture simulations have been demonstrated to offer significant accuracy improvements for other open systems that are not exactly solvable, our work also provides further insight into how and why this advantage arises.

Introduction
Developing a more detailed understanding of the numerous intriguing phenomena displayed by strongly correlated quantum systems is one of the major theoretical challenges in physics today. To meet this challenge a formidable arsenal of non-perturbative, renormalization and numerical techniques have been devised. The success of these approaches has for the most part been routed in situations at or close to equilibrium, while comparatively little is known about the physics of strongly correlated systems far-from-equilibrium. Yet non-equilibrium systems are both ubiquitous and of significant practical interest in physics. A typical example is where a finite sized strongly correlated quantum system is driven far-from-equilibrium by introducing couplings to several different macroscopic reservoirs, forming an open quantum system [1]. Under these circumstances both analytical and numerical descriptions of the behaviour of the system become highly nontrivial.
An immediate need to study open quantum systems is given by the inevitable decoherence and dissipation present in any realistic experimental realization of a strongly correlated quantum system. Numerous examples of such experiments now exist ranging from arrays of Josephson junctions [2], ultra-cold atoms in optical lattice [3,4], ion traps [5]- [7] and arrays of coupled microcavities [8]- [10]. Beyond this, however, open systems are becoming increasingly relevant in themselves as efforts are made to both understand, and also potentially exploit, the 3 beautiful and subtle interplay of the coherent many-body dynamics and incoherent quantum processes they possess. One example can be found in quantum information processing [11] where the suppression of noise is typically considered a prerequisite. Despite this it is found that certain dissipative processes can in fact assist in the preparation of highly entangled quantum states [12]- [14]. More generally some of the most common occurrences of nonequilibrium physics are in transport problems [15] relevant to numerous systems including quantum contacts [16], molecular motors [17], molecular junctions [18] and other lowdimensional heat conducting quantum systems [19]. In addition to revealing a wealth of non-equilibrium phenomena, including non-diffuse heat transfer [20] and negative differential conductance [21], the presence of noise has been found, contrary to expectations, to enhance transmission efficiency through a dissipative quantum network [22,23]. There, noise has a beneficial influence thanks to its interplay with destructive quantum interference and energy mismatches [24]. It is therefore of technological relevance to better understand quantum mechanical effects in driven dissipative strongly correlated systems in order to exploit them in achieving more robust and efficient energy transfer in artificial structures [24] and nanomaterials [21]. Finally, open quantum systems present a virtually unexplored landscape of non-equilibrium phases transitions whose properties are likely to differ considerably from conventional equilibrium transitions [25].
In this work, we shall adopt a master equation [1] description of an open quantum system. Our attention is focused on a specific class of open quantum systems, described in detail in section 2, which are governed by a quadratic spinless fermionic Hamiltonian and coupled to baths described by linear 9 Lindblad operators. Only very recently was this class of open system solved semi-analytically [26] with this solution later providing strong evidence [27] identifying, somewhat unexpectedly, a phase transition in the far-from-equilibrium open XY spin chain with boundary pumping, but no losses otherwise. Despite being a specialized type of system its relevance is elevated by the fact that only a very limited number of exactly solvable master equation models are known, namely those involving a single particle, harmonic oscillator or spin. Indeed the presence of a non-equilibrium transition in this solvable model suggests that it may well come to represent a paradigm for such phenomena analogous to the Ising model for quantum phase transitions.
In this work, we present an entirely complementary exact solution of this system for the Heisenberg picture evolution of commonly required observables by employing a matrix product ansatz [28]- [30] for the operators, often called a matrix product operator (MPO). This result is a significant extension of the exact MPO solution presented in [31] for the purely coherent limit of the same system. By using this approach our work provides, through the formal structure of the MPO solution, further physical insight into why this model is exactly solvable. The utility of our solution, however, is only truly revealed when it is combined with powerful matrix productbased numerical methods, of which density matrix renormalization group (DMRG) [32,33], and more recently its quantum information inspired extension to time evolution [34]- [36], are leading examples. These numerical methods enable the exact Heisenberg picture MPO solution of the dynamical evolution of operators to be computed under very general situations, which are not easily accessible otherwise. Given the rich properties of the model [27] this in itself may be very useful. However, perhaps of even greater importance, it also provides a significant 4 nontrivial example of where Heisenberg picture MPO numerics is exact for an open system. It was shown [31] recently that in some cases it is much more efficient and accurate to simulate open quantum systems in the Heisenberg picture and thus the underlying numerical method used here can be readily applied to more general interacting systems. For the equivalent Schrödinger picture MPO numerics which solve for the density operator of an open system [37,38] the study of entanglement has provided a crucial understanding of its strengths and limitation [39]- [42]. In contrast the merits of performing the numerics in the Heisenberg picture for general situations is far less clear. Our result provides further evidence in understanding when rigourously exact or good approximate Heisenberg picture MPO solutions may exist. We note also that promising results have been found very recently in combining solutions of the Heisenberg equations of motion with matrix product representations of states for bosonic systems [43].
The structure of this paper is as follows. In section 2, we describe the master equation and system we shall solve exactly and introduce a particular spin-chain model that our later numerical calculations will focus on. Our solution exploits the MPO formalism and so section 3 describes all the necessary details. We then show in section 4 that the Heisenberg picture solution of many operators for a closed coherent system possesses an MPO form with a finite dimension, a fact which is crucial for the exact solution to be numerically accessible. In section 5, we introduce an ancilla construction which reproduces the underlying master equation introduced in section 2. A key component of this construction is the tracing out of ancillae and the effect of this on an MPO is described in section 6. We then combine these observations in section 7 to demonstrate that an MPO solution with a bounded dimension, identical to that of the purely coherent case, exists for the open system considered. The versatility of this result is highlighted in section 8 where we numerically determine the MPO solutions for several situations, including the approach to stationarity and a sudden quench of the transverse field. Finally, in section 9 we conclude and comment on future work.

Model
In its most general form the physical system considered in the work is a one-dimensional (1D) system of spinless fermions governed by a quadratic Hamiltonian which reads where c † j is fermionic creation operator for site j. By demanding H s to be Hermitian we can choose a to be real symmetric and b to be real antisymmetric matrices. To describe an open fermi lattice we adopt the quantum master equation approach leading to a Heisenberg picture evolution of an operator O(t) that is governed by the Lindblad master equation of the form [1] (usingh = 1 throughout) where the Hamiltonian and Lindblad superoperators are time independent and defined as

5
... The coherent evolution of the spin chain is described by an X Y -type Hamiltonian H x y which maps to an effective quadratic fermionic Hamiltonian. Boundaries of the chain are subject to couplings to baths, which are described by Lindblad operators L γ each of which map to linear fermionic operators.
respectively. Here L γ are Lindblad operators specifying the coupling of the system to a set of Markovian baths. We place a restriction on the operators L γ that they are linear in the fermionic creation and annihilation operators with the form where γ j and l γ j are complex coefficients. A final constraint, which shall been seen in section 5 to be essential to our result, is that O(t) must be an even ordered operator, so PO(t)P = O(t) where P = N j=1 (1 − 2c † j c j ) is the parity operator. Since parity is conserved by equation (2) with quadratic H and linear L γ 's we require only that the initial operator O(0) is even.
Very recently this class of open systems was solved semi-analytically [26] by an entirely different approach to that which will be described here. In [26] a sophisticated method of constructing a Fock space of operators was employed which maps the Liouvillian into a form which can be diagonalized by a procedure analogous to the famous solution of the XY Hamiltonian [44]. This solution gives access to a range of properties of the system including expectation values of observables for the non-equilibrium stationary state and excitations, as well as the spectrum of so-called rapidities [26]. We shall exploit this solution later in section 8 for testing the approach to stationarity in a dynamical setting.
The fermionic model outlined has considerable freedom in the non-locality of the terms in H s and the operators L γ . We shall consider a concrete example within this class of open fermi systems composed of an XY spin chain with boundary pumping, as depicted in figure 1. As is well known the Jordan-Wigner transformation which relates spin ladder operators σ ± to fermionic creation and annihilation operators maps the XY spin-chain Hamiltonian directly to a spinless fermionic Hamiltonian [44] of the type in equation (1). Here J is the strength of the nearest-neighbour spin coupling, γ is the anisotropy, B is a transverse magnetic field, and σ α j is the α = {x, y, z} Pauli spin operator on the jth spin. The boundary pumping is described by the set of Lindblad operators give the temperature of the thermal state that the baths drive the boundary spins to. This spin-chain setup is not only of importance to heat and spin transport problems [15] in 1D but also strong numerical evidence suggests it possesses a non-equilibrium phase transition as B is varied [27]. Later in section 8 we shall present some exact numerical results for the dynamical behaviour of this system possible only through the solution that we will now describe.

MPOs
The framework in which we cast our exact solution of equation (2) is the matrix product representation of operators. Given a system composed of N sites each with a local ddimensional Hilbert space spanned by the states | j , we define the tensor-product basis states as |j = | j 1 | j 2 . . . | j N where j = ( j 1 , j 2 , . . . , j N ) is a vector of physical indices. An arbitrary operator O acting on this system can then be expanded in the operator basis |j k| as O = j,k o j,k |j k |. A MPO is where the coefficients o j,k of this expansion are expressed in the following form [29,30]: where A [n] j n k n is a matrix, of dimension χ × χ , for each site n, selected by two independent physical indices j n and k n for that site, while L| and |R are χ-dimensional row and column boundary vectors, respectively. A schematic representation of this quantity is given in figure 2. Each expansion coefficient o j,k is therefore encoded as a particular ordered product of A matrices associated with each site, which is contracted to a scalar by the fixed boundary vectors.

7
Given that there are in general exponentially many coefficients o j,k , a matrix product representation in equation (7) yields a highly compact description of an operator if it requires only a small dimension χ. For this reason, and others, matrix product representations for both states and operators have been applied with considerable success in a variety of related numerical methods. The key to their success is that many states or operators, for 1D systems at least, can be very accurately approximated by a matrix product representation of small dimension despite formally requiring a much larger intractable dimension to be exact. In contrast to this the MPO solutions we shall present require only a bounded dimension for the representation to be exact when describing operators evolving according to the open system introduced in section 2. This means that by utilizing one of these matrix product methods, namely the time-evolving-block-decimation (TEBD) algorithm, we can evaluate the exact solution numerically. However, beyond this much is learnt about the nature of the solution by examining the structure of the formal MPO solution itself. For this purpose we utilize an entirely lower triangular form for all A-matrices, introduced in [45]- [47], which permits exact low-dimensional MPO representations for many operators to be constructed easily. The key feature of this approach is that the lower triangular form is preserved under the standard matrix product manipulations such as direct sum or direct product. This means that if an operator O A has a MPO representation with matrices A of dimension χ A , and an operator O B has one with matrices B with dimension χ B , then the operator N over a system of N sites and is an MPO with a χ = 1) can be extended to highly nontrivial operators with an MPO dimension greater than unity.
For our purposes we need only consider the simplest MPO with a general 2 × 2 lower triangular form. For an operator O we assign the following matrices to each site where p, q and r are d × d matrices representing local operators on a site. Note that in this compact form of A the physical indices j and k are subsumed into physical operators p, q and r , while the row and column indices of the A matrix are the internal χ = 2 indices of the MPO representation. To compute the full operator described by assigning A to every site we note that the standard multiplication of two A matrices is equivalent to the tensor product of the physical operators they contain as For a longer string of multiplications A × A × . . . × A this generalizes to yield an operator in the bottom left corner which is the sum of all terms of the form r ⊗ · · · ⊗ r ⊗ q ⊗ p ⊗ · · · ⊗ p with the location of the q operator in the string translated. Finally, for a lower triangular MPO the full operator is extracted via the left and right boundary states L| = (0, 1) and |R = (1, 0) T , which select the bottom left operator 'matrix element' from the matrix product. Using appropriate choices of a, b and c many useful single-particle operators can be formed, for example j σ z j is formed by each site having a matrix [45] A = 1 0 σ z 1 .

8
Notice that an MPO representation is based on a tensor-product structure and therefore implicitly assumes commutativity between local operators appearing in the A matrices for different lattice sites. The local operators cannot therefore be fermionic directly. For products of such operator sums, which we shall consider shortly, this means that MPOs always arrange the resulting local operators in lattice site ordering.

Exact MPO solution for a closed system
Using the lower triangular MPO formalism we reexpress the finite-dimensional MPO solution described in [31] for any fermionic operator governed by purely coherent evolution with a quadratic Hamiltonian H s . To do this we need only consider an arbitrary local sum of creation and annihilation operators C = xc + yc † . The formal solution to the equation of motion of this operator has the standard form

It can be readily shown that the action of H on C for a quadratic H s is
and thus C is transformed into a sum of linear operators spread across the lattice. The linearity of H implies that its repeated application any integer number of times p as H p {C } generates only a linear operator. Now since the formal solution of the equation of motion can be expanded as we establish the well-known fact that the Heisenberg picture unitary time evolution of the operator C governed by a quadratic Hamiltonian is closed. The general solution can then be written as where α j (t) and β j (t) are time-dependent complex coefficients containing all the nontrivial features of the evolution. To recast this solution in MPO form we apply an inverse Jordan-Wigner transformation back to the equivalent spin representation giving The spin operator on the right-hand side can be expressed as a simple 2 × 2 lower triangular MPO, independent on the number of sites N , with site-dependent matrices where X n j (t) = α j (t)σ + j + β j (t)σ − j . From our earlier discussion the bottom left X j (t) operator inserts the necessary site and time-dependent superposition of spin raising and lowering operators into the product, while the bottom right σ z j operator creates the Jordan-Wigner operator string, which establishes the appropriate anti-commutative behaviour.
Since the evolution is unitary the solution to c † j (t) and c j (t) for all sites j automatically provides the time evolution for any string of local sums of creation and annihilation operators, i.e.
Two consequences of this are that the dynamics of a quadratic Hamiltonian H conserves the order of any initial fermionic operator and the MPO solution for the operator string is simply the direct product of the MPO solution for each constituent C-operator given in equation (11). The latter then straightforwardly determines the fixed matrix product dimension χ required for the solution of any given operator string so χ = 2 n for an nth order operator. For example, a general quadratic operator C p (t)C q (t) has a 4 × 4 MPO representation, independent of p and q given by matrices for each site j as where X pj (t) and X q j (t) are the site-dependent X operators associated with C p (t) and C q (t).
The Kronecker product of the MPO solutions gives the appropriately enlarged boundary vectors L| = (0, 0, 0, 1) and |R = (1, 0, 0, 0) T which select the accumulated operators in the bottom left corner as A general feature of such solutions for strings of C operators is that each constituent C operator contributes its own X operator to the representation and shows how highly constrained the evolution of operators is in the space of operators, a fact that has ultimately permitted such a compact representation.
Common spin-chain observables such as σ z j , σ x j σ x j+1 and σ y j σ y j+1 are contained in this class of 4 × 4 MPOs. Long-range correlations like σ z p σ z q involve quartic fermionic operators, independent of p and q and thus require χ = 16. However, the behaviour of some operators can be very different. For local spin observables such as σ x j and σ y j the fermionic representation obtained via an inverse Jordan-Wigner transformation acquires a linearly growing order with the site index j due to the string of (1 − 2c † c) operators which appear. Such an operator could then require an exponentially growing MPO dimension χ = 2 | j| to describe its exact solution. Correlations like σ y p σ y q behave similarly with an exponentially growing dimension χ = 4 | p−q| dependent on their separation. What we shall now show in the remainder of this paper is that the χ required for the MPO solution of even ordered fermionic operators evolving according to the open system described in section 2 is identical to that of the purely coherent system.

Ancilla master equation construction
Open quantum systems typically arise when the system of interest interacts with a large bath or reservoir, often identified as the system's environment. Using this approach Lindblad master equations can be rigorously derived using various microscopic models of the systemenvironment interactions under the Born-Markov approximation and in the limit of extremely large reservoirs [1]. To prove that an exact finite-dimensional MPO representation exists for the open systems introduced in section 2 we shall instead employ a derivation of a master equation similar to that of non-selective continuous measurements [1]. While this construction itself is perhaps less physically motivated it has the advantage for our purposes that it yields a Lindblad master equation exactly with no additional approximations.
A non-selective continuous measurement process involves dividing time into small intervals of length δt with each interval associated with a separate independent ancilla (or probe) forming a time-ordered chain. At the beginning of each interval δt the system evolves coherently and interacts with the associated ancilla, which is subsequently measured at the end of the interval. Depending on the interaction, measurement and ancilla initial state this setup represents a general indirect continuous monitoring of the system [1,48]. In the case where the indirect measurement is ideal the evolution of the system is frozen by the quantum-Zeno effect. For more general imperfect measurements the system evolves according to a master equation with Hermitian Lindblad operators. In order to model the linear fermionic Lindblad operators introduced in equation (5) we modify this construction slightly by considering a different class of system-ancilla coupling and trace out rather than measure the ancilla at every time step. As we shall show below this setup, depicted in figure 3, produces in the continuous limit an effective evolution of the system that is again described exactly by a Markov master equation with the chain of ancillae representing a manifestly delta-correlated environment in time.
The construction begins by augmenting the system of N sites with a chain of τ + 1 ancilla sites described by the fermionic modes a t with t = 0, 1, . . . , τ . Occupation states of the system + ancillae are chosen to be defined by the specific mode ordering where n = (n 1 , . . . , n N ) and m = (m 0 , . . . , m τ ) are binary vectors of occupation numbers over the ancillae and system modes, respectively. By placing the ancillae modes to the right this choice, in conjunction with the Jordan-Wigner transformation defined in equation (6), ensures that any system operator has a spin equivalent of the form 10 where 1 a is the identity over the corresponding ancillae spins. This enables the tracing of spins to be completely equivalent to the tracing of the corresponding fermionic mode. The ancilla mode label t is essentially a time label denoting at which time interval the full Hamiltonian of the system will involve that ancilla mode, as depicted in figure 3. In equation (13), we have also ordered the ancillae amongst themselves with their time label increasing inwards from the right so tracing can proceed iteratively from the boundary. The full Hamiltonian of the system + ancillae is composed of two parts; the time-independent system Hamiltonian H s involving only system modes, and H i (t) which is a time-dependent interaction Hamiltonian between the system and ancillae modes. The time-dependence of H i (t) is taken to be piece-wise constant over intervals δt giving a full Hamiltonian where S is a system operator and (t) is the Heaviside function. Notice that the interaction between the lattice and ancilla in equation (14) depends on δt and is singular in the limit δt → 0. This is physically required in order for the ancilla to have a finite influence on the lattice in the limit of a vanishingly small interaction time [48]. Also the ancillae possess a zero self-Hamiltonian so the only dynamics acting upon them is that generated by the terms in H i (t).
As a final definition for this construction we take the initial time t = 0 state ρ of the system + ancillae to have all ancillae modes unoccupied, but otherwise arbitrary. Let us focus on a particular time t = T δt. Between the time t and t + δt the Hamiltonian H (t) is time independent and only involves the system modes and the ancilla mode a T . For this reason we shall, without loss of generality, restrict our considerations to these modes only 11 . To make the connection to MPOs transparent we perform a Jordan-Wigner transformation and work with a spin representation. The initial density matrix at time t = 0 becomes a spin state in which the lattice and all ancilla spins are uncorrelated ρ = ρ s ⊗ | ↑ ↑ | T with ρ s being an arbitrary spin state of the system. At time t the initial operator O(t) for the system + ancilla spin is composed of system modes only and so it transforms to is the system operator resulting from earlier evolution. Similarly, the relevant interaction terms 10 For notational clarity we do not distinguish symbolically between a fermionic operator and its Jordan-Wigner transformed spin equivalent. It should be clear from the context, which is implied. 11 Ancillae modes T + 1, . . . , τ , which are yet to interact are spectators in the proceeding calculation since neither O(t) nor H (t) contain any of these modes. The ancillae modes 0, . . . , T − 1, which have previously interacted may be contained in O(t). The proceeding calculation is the same regardless of whether these modes are traced out before or after the considered time interval δt. Thus for brevity, we assume that they have been traced out before in the same fashion as we shall trace out ancilla mode T below.
in H (t) transform to the spin operators where P = N j=1 σ z j is the spin equivalent of the parity operator for the system. During the time interval δt the formal solution for the evolution is which can be simplified considerably after the partial expectation value is taken due to the special choice of interaction H i (t) and ancilla initial state | ↑ T . In particular, ↑ |H i (t)| ↑ T = 0 signifying that the ancilla has no direct back action on the system [1,48]. The surviving second-order term involving Using the resulting evolution an equation of motion is formed by taking the continuous limit as An implicit inverse Jordan-Wigner transformation back to spinless fermions can be assumed whereupon we see that the construction has yielded a standard Lindblad master equation [49,50] with Lindblad operators L = PS. The construction can be straightforwardly extended to account for multiple Lindblad operators L γ by introducing more ancillae and additional interactions of the form in equation (14) at each time interval. It also admits the option of having explicitly time-dependent Lindblad operators L γ (t).
The presence of the P operator relating the coupling S to the resulting Lindblad operator L has important consequences. Following our requirements outlined in section 2 our aim is for this construction to model linear L operators. For even parity operators O(0) the P operator plays no role making the coupling S equivalent to the Lindblad operator and therefore linear also. For the same choice of linear coupling S odd parity operators O(0) instead evolve according to a different Lindblad superoperatorL of the form with a sign flip of the first term signifying that the Lindblad operators L are now higher order. Alternatively, for odd parity operators to evolve with linear Lindblad operators L the interaction S must instead include the parity and be higher order.

Tracing out ancilla within an MPO
A crucial step in the master equation construction is the repeated tracing out of an initially uncorrelated ancilla. We now detail the consequences this step has on the resulting MPO representation of the system operator O s (t) given an MPO of the full operator O(t). Specifically, let us take O(t) as being represented by MPO matrices A [ j] for each site j of dimension χ, and, without loss of generality, take the ancilla to be the last site j = N + 1.
Given an initial density matrix ρ s ⊗ ρ a figure 4 shows that the MPO representing the system operator O s (t), satisfying tr sa [O(t)ρ s ⊗ ρ a ] = tr s [O s (t)ρ s ], can be found by contracting in isolation the single-site ancilla density matrix ρ a with the matrices The contribution of the ancilla to the remaining MPO is then reduced to a matrix T = jk A [N +1] jk (ρ a ) jk whose effect is simply to transform the right boundary vector as T|R = |R . The MPO for O s (t), defined only over system sites j = 1, . . . , N , then retains the same set of matrices A [ j] for those sites, but possesses the new right boundary vector |R . Thus, so as long as the initial density matrix between the system and ancilla is uncorrelated, the dimension of the MPO for O s (t) is identical to that of O(t). This conclusion can be readily seen to hold for any number of initially uncorrelated ancilla located at the right edge of the total system. For clarity let us consider what type of operators arise from using an arbitrary left boundary vector with lower triangular MPOs. Using the earlier formal solution B [ j] for a generic quadratic operator C p (t)C q (t) given in equation (12), the introduction of an arbitrary right boundary vector |R = (α, β, γ , δ) T , while keeping L| = (0, 0, 0, 1), gives a weighted sum of the bottom row of operator 'matrix elements'. More precisely it yields an operator which now contains, modulo a parity operator P, linear and zeroth-order operators that are derived from the operators appearing in the quadratic operator string with the admixture determined by the components of |R . In general, an arbitrary right boundary vector |R for an nth order string C p (t)C q (t) . . . C k (t) generates a sum of all operators of order less than or equal to n derived from the constituents of this string. When the parent string operator is an even parity operator then lower order odd parity terms, such as the linear terms in equation (18), acquire an additional factor of P. The opposite occurs when tracing a parent string operator which has an odd parity. The appearance of P operators coincides with terms which do not share the parity of the parent string C p (t)C q (t) . . . C k (t) and thus they violate parity conservation. These terms, however, are never generated by the ancilla construction introduced, since tracing out the ancilla spin in the specific initial state ρ a = | ↑ ↑ | does not generate an arbitrary boundary vector |R . Instead the allowed structure for |R can be readily discerned by again considering the operator O(t) = C p (t)C q (t) starting with the standard |R . Tracing out the ancilla spin j = N + 1 results where the complex number ζ = ↑ |X p (t)X q (t)| ↑ is not necessarily zero. Absorbing this matrix into the boundary gives |R = (1, 0, 0, ζ ) T signifying that only zeroth and quadratic terms can appear. Depending on the nature of the system-ancilla interaction tracing out a single ancilla spin for a general nth order operator string can in principle generate a boundary vector |R corresponding to a sum of operators of order n, n − 2, n − 4, . . . , mod(n, 2). This is then consistent with the resulting incoherent evolution preserving the parity of the initial fermionic operators.
This shows that a given lower triangular MPOs already possess the capacity to describe a very specific class of mixed order operators simply by varying one of the boundary vectors. As described in section 4 the coherent evolution according to a quadratic Hamiltonian H s of any one C operator in a string is described by a specific time-dependent X operator in its MPO representation. This is true regardless of the boundary vector, and so the same coherent quadratic evolution for this type of mixed order operator is automatically captured by this MPO solution.

Building an MPO solution
The results of the preceding sections can be readily combined to demonstrate that the bounded dimension of MPOs seen for coherent quadratic Heisenberg picture evolution also applies to even parity operators evolving according to the specific open system introduced in section 2. As mentioned in section 5 for even parity operators the requirement for linear Lindblad operators is met by using a linear coupling operator S between the system and ancilla. This, along with a quadratic system Hamiltonian H s , makes the full time-dependent system + ancillae Hamiltonian H (t) quadratic. If all the ancillae modes are retained, as depicted in figure 5(a), then the subsequent evolution would be entirely coherent and would represent a purification of the open dynamics of the system alone.
Since we have a coherent quadratic evolution, following the discussion in section 4, an exact MPO solution of fixed dimension therefore exists for the full operator O(t). To extract the reduced operator O s (t) for the system for any time 0 t τ δt the entire ancillae chain is traced out, of which only those labelled up to t have any relevance. Given that the ancillae and system are uncorrelated initially the tracing out of the ancillae has no effect on the resulting MPO dimension for O s (t). The tracing out of the ancilla sites yields a product of timeordered T matrices 12 as shown in figure 5(b). The incoherent effects induced by the ancilla are entirely captured by a time-dependent boundary vector |R(t) = T t . . . T δt T 0 |R . This therefore establishes that for any even-ordered initial system operator O s (0) there is an equality of the required MPO dimension for its coherent evolution with a quadratic Hamiltonian and its incoherent evolution with this special type of open system. Since mixed order operators arising from tracing out any single ancilla can also be coherently evolved, with no change in their MPO dimension, this conclusion is independent of when the tracing is performed. In particular ancilla may be traced out immediately after they interact, as done explicitly in section 5, and thus the bounded MPO dimension applies to the continuum limit as well. We demonstrate this with some numerical examples in section 8.

Properties of the MPO solution
Here, we make some additional comments on the MPO solution found. Firstly, the existence of an exact solution for this open system is not simply a consequence of the closure of the equations of motion for the lowest order as it is for quadratic coherent evolution in section 4. The lack of unitarity of the evolution means that (C p C q . . . C k )(t) = C p (t)C q (t) . . . C k (t) in general so knowledge of the evolution of lower order operators does not furnish us with knowledge of the evolution of higher order products. Secondly, as shown in section 6 the formal structure of the MPO solution with a time-varying boundary vector with |R(t) = (1, 0, . . . , 0) T implies that the evolution of an initial nth-order operator will involve a special type of mixed order operator composed of equal or lower order operators only. This behaviour for the operator order, revealed by the formal MPO solution, is entirely consistent with what is seen for other related open systems [1].
A classic example of this is provided by a modified version of the damped harmonic oscillator. Here the coherent part of the master equation in equation (3) has H = ωb † b where ω is the oscillator frequency and b is its corresponding bosonic annihilation operator. We then take the Lindblad superoperator L in equation (4) as being described by a single linear Lindblad operator L = αb + βb † analogous to the fermionic model introduced in section 2. Since any initial operator of the system can be expanded as O(t) = nm o nm (t)(b † ) n (b) m with n, m 0 and coefficients o mn (t), we need only consider the effect of the right-hand side of equation (2) on a general term within this expansion. The action of the coherent part is simply 12 For ancilla related to later times which have yet to interact the transformation matrix T = 1.
Thus, whenever |β| > 0 the action of (H + L) is to leave the order of a constituent term (b † ) n (b) m unchanged as (n + m) or reduced by two. The formal solution O(t) = exp(H + L)t{O(0)} for an initial operator O(0) of order n (odd or even) will in general include contributions only from orders n, n − 2, n − 4, . . . , mod(n, 2). An identical type of analysis can be performed for an open bosonic lattice defined by annihilation operators b j for each site j and again governed by a quadratic Hamiltonian whereā andb are real symmetric matrices, along with linear Lindblad operators This readily confirms that the lack of growth of an operators order seen for a single oscillator also applies to the fully bosonic version of the model introduced in section 2.
Applying a similar analysis on the fermionic lattice model itself reveals that for the specific Lindblad superoperator L defined in equation (4) only initially even ordered operators display this closure property. In contrast odd-ordered operators can be shown to acquire a proliferating order under the repeated application of L. When such growth in the order occurs the link between operator order and MPO dimension seen for the coherent solution in section 4 suggest that the dimension will not be bounded 13 . Notice that our ancilla construction applied to odd parity operators does not model L, but ratherL. Incoherent evolution according toL reverses the situation with odd parity operators now displaying no growth in their order. Thus the ancilla construction presented in section 5 models an incoherent evolution where all operators O(t) have a bounded order, which in turn permits the bounded dimension MPO solution. It is only for even parity operators, however, that this evolution corresponds to the precise open system defined in section 2.
Finally, the MPO solution offers a more efficient representation than the closure of the operator order can provide on its own. In particular, once the exact χ is used for the MPO its description only grows linearly with the number of sites N . This is also an improved scaling compared to alterative approaches to this open fermi system exploiting fermionic Gaussian states [52]. They display a N 2 scaling and moreover are restricted to considering initial states which are of Gaussian form. The Heisenberg picture MPO approach used here can compute properties for any initial state which can itself be well approximated by a matrix product ansatz. Importantly, this can be either a matrix product state or, as figure 4 illustrates, it can be a mixed state so long as its density operator has an accurate MPO representation with a tractable dimension. Indeed it has already been shown that such MPO representations of thermal states for 1D systems can be obtained numerically [37,38] so finite temperatures are also accessible within this framework.

Numerical examples
Having shown that there exists a formal MPO solution with a specific fixed dimension χ for a given system operator we now show that this solution can be determined numerically via 13 Numerical evidence following calculations like those to be presented in section 8 confirms this.
Heisenberg picture evolution with the TEBD algorithm [31,34,35]. While the formal solution presented has no restrictions regarding the locality of the terms in H s and the Lindblad operators L γ , efficient integration of the equation of motion via the TEBD algorithm requires that terms are nearest neighbour. Moreover since linear Lindblad operators involving fermionic creation and annihilation operators away from the boundaries acquire a many-spin Jordan-Wigner σ z string the numerical solution is restricted to noise terms on or one site in from the boundary. This means that specific open XY spin-chain model introduced in section 2 can be solved with this numerical method.
The numerical solution determined by TEBD unmistakably demonstrates the existence of the bounded MPO dimension proven above for this open system, just as it does for the coherent limit [53]. When normalized according to the Frobenius norm, where j,k |o j,k | 2 = 1, the MPO solution produced by TEBD 14 is in canonical form [34,51] where the Schmidt decomposition [11] of the operator for any contiguous bipartition is explicitly contained in the representation [37]. Here [n] j n k n are sets of matrices, different to the lower triangular matrices A [n] j n k n used earlier but performing the same function. The new important addition to this representation are the diagonal matrices [n] for the bulk (i.e. n = 2, . . . , N − 1) with diagonal elements equal to the Schmidt coefficients λ ν of a bipartition after site n. The appropriately sized boundary vectors are now [1] | = (1, 0, . . . , 0) and | [N +1] = (1, 0, . . . , 0) T , representing the single unit Schmidt coefficient before site 1 and after site N . Once in this canonical form Schmidt coefficients allow the effective MPO dimension of the operators to be identified by counting the number of significant Schmidt coefficients λ ν 1, where is some small threshold.
For an XY chain of length N = 50 (see figure 6 for the parameters used) we have calculated the evolution of the operators σ z 25 , σ z 1 σ z 50 and σ x 24 σ y 27 . In figure 6(a), the central Schmidt coefficients for the chosen operators after a time J t = 50 of evolution are shown. A clear cut-off in the λ ν 's is seen where their value drops in excess of 11 orders of magnitude. This cut-off is robust to time-evolution and the insignificant λ ν 's are numerical noise that may be safely truncated away. The effective MPO dimension given by this cut-off coincides with the dimension expected from the formal solution. We may therefore rigidly enforce the exact MPO dimension required and given the lack of truncation the only error in the time integration comes from the customary Trotter expansion used in TEBD. In figure 6(b), the resulting time evolution of the central z-magnetization σ z 25 and boundary-boundary z correlation σ z 1 σ z 50 is shown for an initial spin-polarized state | ↓↓ · · · ↓ . The transient evolution displays plateaus caused by the time it takes for the influence of the boundary pumping to propagate across the chain. This is better illustrated in figure 7 where the evolution of the z-magnetization profile of the entire chain is plotted up to a time J t = 500. Being plotted with a logarithmic timescale it is apparent that the majority of the z-magnetization in the bulk is eroded rapidly by the dynamics from its initial value. However, in figure 8(a) a more detailed comparison of the z-magnetization profile at a time J t = 500 and the stationary profile [26] reveals that for N = 50 spins this time is only sufficient to drive the boundary z-magnetization to their stationary values and that the bulk  The z-magnetization profile for the entire chain over time plotted against τ = log 10 (1 + J t) up to a time of J t = 500. For visual convenience the stationary t → ∞ magnetization profile determined from the exact solution given in [26] is plotted at τ = 3. See also figure 8(a) for a comparison of the stationary profile with that attained after a time J t = 500. The parameters used are identical to those stated in figure 6.
is still far from stationary. Tests reveal that a significantly longer evolution time is needed to achieve convergence of the bulk z-magnetizations.
To demonstrate a time-dependent dynamical scenario 15 we consider the simplest case of an abrupt quench of the transverse field. Specifically we evolve the initial spin-polarized state with B/J = 10 for a time J t = 500, analogous to the previous example. Then at the time J t = 500 the transverse field is switched instantaneously to B/J = 1 and the evolution is continued. In figure 8(b) the evolution of the central z-magnetization σ z 25 is shown as a function of time around the quench point. For times J t < 500 we see that there is a slow change in σ z 25 and it is still quite far from its stationary value of σ z 25 s = −0.0161. In the same region of time figure 8(b) shows (dashed line) the evolution of σ z 25 from the previous example where B/J = 1 throughout which is slightly closer to its stationary value of σ z 25 s = −0.0391 but displays a similar rate of convergence. For time J t > 500, after the quench, there is initially a rapid change in σ z 25 which after a time of approximately 15/J then settles down with small oscillations around a new value, which again differs from the stationary value of the new transverse field. Instead this newly acquired z-magnetization is very close to the non-stationary value obtained via constant evolution with B/J = 1. This shows that even after a comparatively long evolution time the system has retained a significant memory of its initial spin-polarized state.

Conclusions
We have presented a detailed study of the MPO description of a specific class of open quantum systems governed by a master equation with a quadratic spinless fermionic Hamiltonian and linear fermionic Lindblad operators. By mapping this master equation to an entirely coherent quadratic evolution involving additional ancillae we have shown that the MPO representation for the evolution of operators with even parity possesses a finite and fixed dimension. This has revealed the quadratic nature of the evolution underlying this class of master equations and our ancilla construction gives decisive insight into why it is exactly solvable. The formal structure of the MPO representation also indicates how a given initial operator can evolve into a specific type of mixed order operator, consistent with behaviour seen in other simpler open systems. Exploiting the fixed MPO dimension the TEBD algorithm allows the dynamical evolution of operators in this non-equilibrium open quantum system to be computed with a cost that is linear in the system size. The dynamical behaviour accessible via the MPO solution presented therefore complements the existing exact solution for this models stationary states and spectral properties [26]. We have exemplified this by computing some examples involving the approach to stationarity and the response of the z-magnetization to a sudden quench in the transverse field.
An interesting calculation, beyond the scope of the current work, is to perform a dynamical quenching through the non-equilibrium quantum phase transition. Such a dynamical calculation appears to be very demanding with the Schrödinger picture [27]. The non-equilibrium transition manifests itself as a discontinuous change in the σ z i σ z j − σ z i σ z j correlations, but not in other local observables such as energy and magnetization. From the MPO perspective of this work the behaviour of this transition appears to be very reminiscent of the matrix-producttype equilibrium quantum phase transitions [54]. Computing a dynamical crossing of this non-equilibrium transition could help determine the realistic adiabacity requirements for its observation.
Beyond this our work has provided an important and nontrivial class of open systems with an exact Heisenberg picture MPO representation. This may yet aid in determining other models where such solutions exist. For instance, it remains to be seen whether Heisenberg picture simulability is readily related to the integrability of the underlying model [53]. For example, a finite-sized XXZ chain can be made integrable with appropriate boundary fields; however, it is not clear that an efficient representation exists for commonly required local observables like σ z . This raises the question as to whether finite-sized MPO representations of certain types of operators are possible for systems possessing a Bethe-ansatz solution. This is an interesting open problem and would reveal if the MPO formalism can aid in evaluating otherwise very complicated quantities from these solutions. For the presently studied XY model the noninteracting nature of the effective fermi system for both the open and closed system appears to be a crucial property permitting simulability, which is more constraining than integrability alone.
Finally, the MPO solution introduced may allow a better understanding of the trade-off between efficiencies possible by changing pictures. Future work [55] will look at how quickly the accuracy of Heisenberg picture simulations breakdown when they are applied to models which are only weakly perturbed from the exact solution presented here. In the context of spin chains the most obvious extensions outside the exact solution would be additional σ z j σ z j+1 interaction terms and/or dephasing noise.