Complexity of Fermionic Dissipative Interactions and Applications to Quantum Computing

Interactions between particles are usually a resource for quantum computing, making quantum many-body systems intractable by any known classical algorithm. In contrast, noise is typically considered as being inimical to quantum many-body correlations, ultimately leading the system to a classically tractable state. This work shows that noise represented by two-body processes, such as pair loss, plays the same role as many-body interactions and makes otherwise classically simulable systems universal for quantum computing. We analyze such processes in detail and establish a complexity transition between simulable and nonsimulable systems as a function of a tuning parameter. We determine important classes of simulable and nonsimulable two-body dissipation. Finally, we show how using resonant dissipation in cold atoms can enhance the performance of two-qubit gates.

Interactions between particles are usually a resource for quantum computing, making quantum many-body systems intractable by any known classical algorithm. In contrast, noise is typically considered as being inimical to quantum many-body correlations, ultimately leading the system to a classically tractable state. This work shows that noise represented by two-body processes, such as pair loss, plays the same role as many-body interactions and makes otherwise classically simulable systems universal for quantum computing. We analyze such processes in detail and establish a complexity transition between simulable and nonsimulable systems as a function of a tuning parameter. We determine important classes of simulable and nonsimulable two-body dissipation. Finally, we show how using resonant dissipation in cold atoms can enhance the performance of two-qubit gates.
Understanding whether a particular quantum system is easy or hard to simulate from the perspective of classical computation is a crucial task serving several goals. The first goal, as a primary step of many numerical studies, is to find efficient classical algorithms describing the desired quantum phenomena. Another goal arises in quantum computing, where finding many-body systems lacking an effective classical description may be worthwhile for constructing quantum computation [1] and simulation [2,3] devices. The versatility of the classical simulability problem can be illustrated by considering the sampling problem for noninteracting and interacting fermions [4][5][6][7]. There are efficient classical algorithms to simulate fermions described by a quadratic Hamiltonian: the amplitudes of time-evolved many-body configurations are expressed by an efficiently computable analytical formula [6,8]. The existence of an efficient algorithm makes the free-fermion approximation a numerically accessible and valuable method with applications to condensed-matter systems. At the same time, simulating interacting fermions is believed to be classically intractable. Indeed, simulating general interacting fermions is as hard as simulating the output of a universal quantum computer [9]. A similar practical differentiation between easy and hard problems can be applied to other systems [10][11][12][13][14].
In this work, we study the fate of classical simulability of fermionic systems in the presence of dissipation, both for computing local observables and for sampling from the manybody output distribution (to be defined shortly). To obtain a classification of the complexity of simulating free fermions with dissipation, we consider a general class of Markovian processes, i.e. dynamics that depends only on the instantaneous system state and is independent of the preceding evolution [15]. In previous studies, it was shown that Markovian single-fermion loss or gain terms keep the noninteracting system classically tractable [16,17]. As a step forward, we consider a much wider class of quadratic-linear Lindblad jump operators. Using the method of stochastic trajectories [18,19], we establish a wide subclass of problems that can be simulated classically. At the same time, we demonstrate that, in general, quadratic Lindblad jump operators are at least as hard to simulate as most unitary interacting systems. More precisely, we establish a connection between dissipative inter-ψ r ψ r'

H(t), A k (t)
FIG. 1. Classical simulability. We look for the existence of an efficient algorithm running on a classical computer and producing (sampling) the many-body configurations with the probability distribution close to the physical system after measurement in some basis. We show that, for fermionic systems with Hamiltonian H(t) and with dissipation described by quadratic Lindblad jump operators A k (t), such an algorithm exists for at least a restricted number of problems, while the worst-case scenario requires a quantum computer in order to be solved efficiently. The three optical lattices illustrate the state of the system at initial, intermediate, and final times.
actions and fault-tolerant universal quantum computation exploiting the quantum Zeno effect [20][21][22][23][24]. Therefore, evolution under quadratic Lindblad jump operators is equivalent in power to quantum computation. This effect can be compared with parity measurements, which can also make free-fermion dynamics universal [25]. The tractability and intractability results together show that simulation of quartic dissipative Liouvillian operators is a problem whose complexity can be changed from hard to easy by varying one or more parameters in the system [13].
We provide a classification of dissipative fermionic processes into easy (efficiently simulable) and hard (not efficiently simulable) classes according to their worst-case computational complexity. The classical simulability problem may be phrased in two ways, either in terms of evaluation of few-body observables or sampling from the full probability distribution on many-body outcomes. In the first task (fewbody observables), a classical computer is required to output the expectation value of an observable supported on k sites, where k does not grow with the system size. In the second task (sampling), a classical computer is tasked with producing samples from the same distribution as the one obtained by measuring the time-evolved state in some canonical basis (see Fig. 1). Both tasks allow for the computer to make a small error , measured appropriately in each case [66]. The task of sampling is computationally harder; an algorithm producing samples in some product-state basis can also be used to obtain expectation values of few-body observables in the same basis. Therefore, in this work, we focus mainly on the easiness of sampling in arbitrary product-state bases as a criterion for overall easiness and on the hardness of computing few-body observables as a criterion for overall hardness. This choice gives the stronger of the two results for both easiness and hardness.
We note that a limited version of classical simulability for some models below was also studied in previous works [67][68][69]. In particular, it was shown that two-point correlation functions in such models can be evaluated via solving a closed set of equations. This result establishes classical simulability of local observables and can be used in various problems such as dissipative transport or optical response. However, this result alone is not sufficient to establish the simulability of sampling. In contrast to local observables, simulated sampling requires the full knowledge of the many-body output probabilities, therefore the sampling complexity of systems with simulable low-order correlations remains unclear. As we revisit below, Gaussian systems are the exception that allow reducing these output probabilities to two-point correlation functions via Wick's theorem; other systems we study below do not have such simple reduction (see Appendix D). To overcome this problem, we develop the easiness proof that does not require applicability of Wick's theorem. In conclusion, sampling is a stronger notion of simulation compared to twopoint correlators in a sense that any local observables can be efficiently obtained using an oracle producing sampling outcomes Let us emphasize the importance of the provided complexity analysis. While established easy dissipation classes are limited to certain fine-tuned processes, such limited simulable models have an important application for quantum computing. For example, classical models can be used in calibration of quantum computers, simulation of the impact of noise on sampling, and analysis of fermionic quantum error-correcting codes [70]. More fundamentally, identifying easy classes is an important first step to analyze easy-hard transitions in open fermionic systems, as we analyze below. At the same time, the hardness result we obtain in this work establishes utilizing dissipative interactions as an alternative path toward building a universal quantum computer. This conclusion is surprising, since dissipative interactions generally produce mixed states. However, dissipative interactions can be used only in a manner utilizing a blockade mechanism induced by the quantum Zeno effect, as we show below. In cold atomic systems, controlling dissipative interactions differs from photonic systems studied before [22,23] and can be achieved using an atomic Feshbach resonance. In this paper, we analyze in detail a scheme for universal quantum computation with 40 K atoms and illustrate that, with realistic experimental parameters, an entangling gate with low error rate of roughly 10 −4 is possible. Existence of both easy and hard classes for two-body dissipation establishes it as a valuable model for physical analysis of noisy intermediate-scale quantum devices.
Model.-We consider dynamics generated by the Lindblad master equation [15,71] where {X, Y } ≡ XY +Y X is the anticommutator, ρ(t) is the density matrix of the system, H(t) is a noninteracting Hamiltonian, and A k (t) ∈ A(t) form a set of k A Lindblad jump operators. We set = 1 throughout the paper unless specified otherwise. Both the Hamiltonian and the Lindblad jump operators may depend explicitly on the time but not on the state itself. The corresponding map ρ(t 2 ) = V(t 2 , t 1 )ρ(t 1 ) between arbitrary times t 1 and This divisibility condition is commonly referred to as the most general definition of Markovian dynamics [72]. The master equation in Eq. (1) is invariant under certain transformations of the set of Lindblad jump operators A(t), such as operator permutations, multiplying any Lindblad operator by a phase factor, or splitting and merging of identical operators as As a physical system of interest, we consider a fermionic many-body problem where N ≤ L spinless fermions initially occupy L available levels. Such systems are commonly described by the second quantization method, which expresses  any operator, including the Hamiltonian and Lindblad jump operators, in terms of fermionic Fock operators c † n and c n , n ∈ {0, 1, . . . L − 1}. Fock operators create and annihilate a single fermion in a particular mode and satisfy the canonical commutation relations {c n , c m } = 0, {c † n , c m } = δ nm . Though the conventional fermion operators are suitable in most physical problems, in the absence of fermion number conservation it is convenient to use the 2L Hermitian Majorana fermion operators γ 2n = c n + c † n and γ 2n+1 = −i(c n − c † n ), due to their simple anticommutation relations {γ i , γ j } = 2δ ij , i, j ∈ {0, 1, . . . , 2L − 1}. We consider the most general form of a noninteracting Hamiltonian [73] where a k (t) are antisymmetric 2L × 2L matrices, b k (t) are 2L vectors, and d k (t) are numbers; all the parameters are complex-valued in general.We assume that the magnitude of all entries of α(t) and β(t) and their time derivatives scale at most polynomially with system size.
In this work, we focus on the classical resources needed to approximately sample from the fermion distribution at time t, where the vectors r and r denote the positions of occupied modes in the initial and final (measured) product-state configurations, respectively, and |ψ r is a product state defined as |ψ r = c † r1 . . . c † r N |0 = γ 2r1 . . . γ 2r N |0 . Here |0 is the vacuum state defined as the state satisfying c n |0 = 0 for all n. Importantly, because the dynamics may not conserve the total fermion number, the final number of fermionsÑ can, in general, be different from the initial number: N =Ñ .
We establish the sufficiency of polynomial resources for classically simulating dynamics due to arbitrary noninteracting Hamiltonians in Eq. (2) and a limited set of Lindblad jump operators A k (t) ∈ A(t) in the worst case. In order to prove polynomial-time simulability (also called easiness) for limited classes of dissipative dynamics, we reduce the problem to that of simulating unitary noninteracting fermionic evolution, an easy problem for a classical computer. In order to prove hardness for more general Lindblad jump operators, we exploit the ability of dissipative dynamics to perform arbitrary quantum computation (i.e. we prove that simulating universal quantum computation reduces to simulating Lindbladian dynamics).
The results of this work are briefly illustrated in Table I. First of all, we define three classically tractable classes of Lindblad jump operators (defined as easy classes 1, 2, and 3). All of these cases allow for polynomial-time sampling of any Hamiltonian and Lindblad jump operators from the given class on a classical computer, with error scaling inversepolynomially with L. Easy Class 1 (EC1) allows for simulation of self-adjoint sets of quadratic Lindblad jump operators: all Lindblad jump operators in the set A(t) come with their Hermitian conjugate. This class includes such widely used examples as dephasing, incoherent particle shuffle, and classical fluctuations of the number of fermions and of the number of fermion pairs. Easy class 2 (EC2) works with unitary quadratic Lindblad jump operators. Finally, easy class 3 (EC3) describes the loss or gain of a single particle in the system and can be used in combination with EC1 and/or EC2. At the same time, there exists a class of Lindblad jump operators with a nonzero measure that is hard to classically simulate. Examples from this class include pair loss or gain and incoherent fermion hopping. Below we explore each class separately.
We focus on quadratic-linear Lindblad jump operators of the form where a k (t) and b k (t) are arbitrary complex-valued 2L × 2L matrices and 2L-vectors respectively, and d k (t) is a number. Notably, the same master equation allows representation with more than one set of jump operators A = {A k }. A set can be reduced to a smaller one if its jump operators are linearly dependent [15]. Therefore, the number k A of Lindblad jump operators in the smallest set does not exceed the number of linearly independent quadratic operators, i.e. L(L + 1). Also, as with the Hamiltonian, we assume that the magnitude of the entries of a k (t), b k (t), and d k (t) and their time derivatives grow at most polynomially with the system size. This work is organized as follows. In section I, we provide a brief introduction to free-fermion sampling, recalling established results in the literature and connecting them to the most general form of quadratic-linear Hamiltonians. In Section II, we derive three new algorithms allowing us to solve distinct classes of fermionic problems involving quadratic Lindblad jump operators and prove that these algorithms run in time that is polynomial in both L and the inverse of the distance from the exact distribution. In Section III, we establish generic class of systems that belong to the hard class and show their robustness to the presence of minor imperfections. In Section IV, we show how our complexity result applies to realistic experimental settings. We demonstrate that natural pair loss in cold atomic systems can be controlled and utilized to implement universal quantum computing. The dissipation-assisted gates provide an alternative to unitary gates with a potential advantage in the speed of two-qubit operations.

I. FREE-FERMION SAMPLING
In this Section, we discuss the noninteracting fermion problem in the absence of dissipation. We recap the work of Terhal and DiVincenzo [6], which shows that all output probabilities P (r|r ) in Eq. (3) and the marginal probabilities can be obtained using a classically tractable analytical formula. Before referring to this result, we need to incorporate the linear terms present in Eq. (2) into effective quadratic dynamics. In order to do so, we consider a slightly larger system containing an extra ancilla (L + 1)th mode [73], labeled as n = L. Next, we choose new effective dynamics such that the ancilla mode remains in the state |+ ≡ (|0 + |1 )/ √ 2 during the entire evolution, including the initial and final times, i.e.
To construct such dynamics, we consider a new Hamiltonian by replacing γ i → iγ i γ 2L , where γ 2L and γ 2L+1 are Majorana operators acting on the ancilla mode. It is straightforward to check that such a transformation results in a new purely quadratic Hamiltonian (without any linear terms) that keeps the state of the ancilla stationary and does not modify the dynamics of the original Hamiltonian. The new coefficients in Eq.
(2) are where we use by default β 2L = β 2L+1 = 0. Given that the modified initial and final conditions for the system and the ancilla are {r } → {r , s }, {r} → {r, s}, s, s ∈ {0, 1}, the probability P (r|r ) of obtaining outcome r for the original system can be computed from the probability P ({r, s}|{r , s }) for the system with the ancilla as follows: Summarizing, this method ensures that the dynamics of a linear-quadratic Hamiltonian can always be reduced to the dynamics of a quadratic one by expanding the system size by one mode. Therefore, we henceforth consider only quadratic Hamiltonians.
Let us derive the formula for the sampling probability. We start from a (backwards) time-evolved Majorana fermion op- Here T exp is the standard time-ordered exponential. Given the quadratic structure of the Hamiltonian, this evolution is ä is a unitary 2L × 2L matrix. One can use this expression to derive the time evolution of a fermion operator as where T nj ≡ R 2n,j + iR 2n+1,j are elements of a L × 2L transformation matrix T . Labeling the initially empty sites as l i and recalling that the initial fermion positions are r i and that the final positions are r i , the linearity allows to write the output probability in Eq. (3) at any time as This expression can be computed efficiently using Wick's theorem and written in a compact form. Let I be a subset of indices with increasing order and A[I] be the matrix whose elements satisfy Then the result can be written as [6] where Pf is the Pfaffian, and M is a 4L × 4L matrix where, in turn, the 2L × 2L matrix Λ is The expression in Eq. (10) can be efficiently evaluated on a classical computer using existing polynomial-time algorithms for computing Pfaffians [8]. Similarly, marginal probabilities can be efficiently computed conditioning on the output of a given fraction of sites, as in Ref. [6], which is enough to efficiently sample from the output probability distribution.

II. EASY CLASSES
Here we present three algorithms that allow simulating specific fermionic problems involving quadratic Hamiltonians and quadratic-linear Lindblad jump operators. All methods are based on stochastic unraveling, i.e., replacing dissipative dynamics by a stochastic free-fermion Hamiltonian without changing the outcome distribution (see also Ref. [19]). Since each stochastic realization can be simulated efficiently by a classical computer, as established in the previous section, a classical computer can serve as a black-box sampler that reproduces measured outcomes. In this section, we demonstrate that the classes of problems belonging to the aforementioned easy classes 1-3 are efficiently simulable. In particular, we show that these problems require computation resources C (number of operations) bounded as C ≤ poly L, t 2 / to sample from a distribution that is -close to the target distribution. Therefore, we establish efficient classical algorithms for approximate dissipative fermion sampling in the presence of certain classes of quadratic-linear Lindblad jump operators.

A. Efficient classical algorithms
Let us define Easy Class 1 (EC1) as problems that involve quadratic-linear self-adjoint Lindblad sets A(t) defined as follows. We assume that at any time one can divide the set as a union of two equal-size subsets, A = A 1 ∪ A 2 , where the Hermitian conjugate of every Lindblad operator in A 1 returns an operator in A 2 (and vice versa). Under this division, any Hermitian Lindblad operator must be included in both subsets A 1 and A 2 with normalization factor 1/ √ 2. The latter splitting can be seen as a transformation that keeps the Lindblad equation invariant, as defined earlier below Eq. (1). Examples from EC1 include several important physical models such as dephasing and classical fluctuations (see examples of sets in lines 1-4 in Table I).
In previous works, it was shown that such systems have two-point correlation functions that are classically simulable by solving a closed set of linear equations [67][68][69]. This is indeed a strong indication that the system can be simulable in the broader context of sampling complexity. However, as we noted previously, Wick's theorem is not applicable to non-Gaussian states (see Appendix D). This means that the scheme we utilized to obtain Eq. (10) does not work any more. We now show an alternative scheme using stochastic unraveling that leads us to the easiness result.
To efficiently simulate dynamics from EC1, we consider a stochastic Hamiltonian where θ k (t) is a complex random variable taking constant values θ k (t) = ξ nk / √ ∆τ during short time intervals t ∈ [n∆τ, (n + 1)∆τ ]. Later we also consider θ k (t) as operators. The discrete complex Gaussian variables ξ nk satisfy Eξ nk = 0, Eξ * nk ξ n k = δ nn δ kk δ ab , where E denotes the expectation value taken over the random variables. Then, given a stochastic Hamiltonian of the form in Eq. (13), the original dynamical map V(t 2 , t 1 ) generated by the Lindblad equation can be approximated as where δV(t 2 , t 1 )∆τ is the lowest-order correction (to be explicitly derived below) and V st is a stochastic unitary map In the above, ä encodes the time evolution due to H (t) in Eq. (13). The average E in Eq. (14) is taken over the stochastic fields θ k (t). The resulting output probabilities satisfy where P st (r|r ) is the output probability for the unitary dynamics in Eq. (15) obtained via the formula in Eq. (10). Therefore, a computer programmed to sample from the distribution for a randomly chosen set of unitary trajectories will produce outcomes with the same probabilities as the physical process following Lindbladian evolution, up to O(∆τ ) error.
The cost of suppression of this error in terms of computational resources will be discussed later in this section. Here we just specify that the correction to the dynamical map, which we treat as an error, can be expressed as where D(t) is a time-local superoperator Here, the operators D can be expressed as polynomials of degree four in terms of the Hamiltonian and Lindblad jump operators at time t. Therefore, D (i) α (t) can always be presented as a sum of terms, each being a product of no more than eight Majorana operators. The specific form of these operators and the derivation of Eq. (14) can be found in Appendix A.
Although the proposed unraveling scheme represents dynamics in terms of stochastic trajectories for Gaussian pure states, the resulting averaged mixed state is non-Gaussian, in contrast to previously studied problems [16,17]. Therefore, the overall dynamics of EC1 represents dissipative interactions of fermions, while the proposed method can be seen as a good choice of time-dependent density matrix decomposition.
Let us consider another class of problems, easy class 2 (EC2), that include unitary quadratic Lindblad jump oper- where Γ k (t) ≥ 0 are timedependent rates and Y k (t) = exp(−iG k (t)) are unitary operators generated by quadratic-linear Hamiltonians G k (t) of the form in Eq. (2). To classically simulate dynamics under EC2, we also consider discretized time evolution with sufficiently small timesteps ∆τ and set the unitary transformation U (t 1 , t 2 ) = n2 n=n1 U n , where the timestep unitaries U n are generated randomly according to the rule Here p k are probabilities that are used to generate the respective outcomes, U 0 , and t n ∈ [n∆τ, (n + 1)∆τ ] are random times generated from the uniform distribution.
Notwithstanding the slightly different stochastic unraveling, the procedure for approximating EC2 is the same as for EC1. In particular, the system dynamics is described by the expression in Eq. (14) leading to the distribution in Eq. (16), with the average taken over stochastic unitary realizations. The correction term has the form in Eq. (18), but the operators D Finally, let us consider easy class 3 described by generic linear Lindblad jump operators A k (t) = i b i k γ i + d k , which can be obtained by setting a k = 0 in Eq. (4) without assuming any additional restrictions on the set A(t). Previous works had already shown that linear jump operators can be simulated classically [16,17]. However, this proof applies only to Gaussian states and cannot be extended to, for example, Lindblad equations that also contain easy quadratic jump operators. Below we propose a different way of simulating linear jump operators similar to the one we used for EC1. This technique would allow us to combine EC1, EC2, and EC3 into a single easy class of Lindblad equations, including both quadratic and linear processes. Now let us show that simulation of linear jump operators is equivalent to simulating a unitary system extended by a number of ancilla modes. In particular, we require L a = t/∆τ ancilla fermion modes equal to the number of time steps after discretization. Let us enumerate the ancilla modes described by Majorana fermion operators γ 2n and γ 2n+1 using indices n = L, . . . , L + L a − 1. We also assume that the ancilla modes are initialized in the vacuum state and traced out after performing the evolution. The dynamics of both the system and the ancillas can be described as unitary evolution with the Hamiltonian in Eq. (13), with one important difference. Now, the quantities θ k (t) in the time interval t ∈ [n∆τ, (n + 1)∆τ ] are operators instead of numbers, and are given by The random variables ξ nk are the same as in EC1. The idea is that we pair a loss (gain) term on the system with a gain (loss) term on the ancilla to make the overall Hamiltonian quadratic. After discarding the ancilla modes, the evolution becomes equivalent to the target dissipative dynamics, up to a discretization error that originates from the approximation in Eq. (14) and Eq. (18), with D We note that the stochastic simulation method takes advantage of the system state being a convex mixture of Gaussian density operators [74,75]. This is a particular case of the more general property that any state of a fermionic system with well-defined parity has a representation as a convex distribution over generic (not necessarily Hermitian) Gaussian operators [76].

B. Performance of the classical algorithms
Let us quantify the error of the method of quantum trajectories used for easy classes 1-3, and then show that the sampled distribution can be made arbitrarily close to the exact one with an appropriate choice of the timestep ∆τ . In order to characterize the approximation error associated with sampling from a distributionP (r|r ) different from the ideal distribution P (r|r ), we utilize the total variation distance where the maximization is over all possible initial productstate configurations r . Using convexity of the absolute value and the expression for the correction in Eqs. (17)-(18), the error can be bounded as (see Appendix E), where and ρ r (t) = V st (t, 0)ρ r are operators transformed according to unitary evolution for a single stochastic trajectory, and the average E is taken over all trajectories. We now use the following lemma to further bound this expression.
Lemma. Consider two sparse operators O 1 and O 2 whose matrix elements satisfy (24) where d H (r, r ) is the Hamming distance between states with the positions of fermions r and r . Let ρ be a normalized positive semidefinite operator, ρ ≥ 0, Tr ρ = 1, then where O α max is the max-norm.
The proof of the Lemma can be found in Appendix E. The result of the Lemma allows us to simplify Eq. (23) as where k iα is the locality of the operator D where k m = 8 for EC1/EC3, and k m = 4 for EC2. We can also bound the max-norm by the (spectral) operator norm As a result, the error bound is given by Since the matrices D that keeps the error in Eq. (28) arbitrarily small, suppressed at least polynomially with the number of modes L.
Let us now estimate the amount of computational resources required to perform the above sampling procedure. For each sample, the algorithm randomly chooses a unitary trajectory according to the given prescription for each class EC1-EC3 and, according to the Terhal-DiVincenzo algorithm, produces outputs from the free-fermion distribution in Eq. (10). In particular, it samples the output at site i conditioned upon the outcomes sampled at sites j < i, for which the marginal probabilities should also be computed. Consider cases of EC1 and EC2 that do not require ancillas. Once the matrix T is obtained in Eq.
For EC3, the derivation is the same up to adjusting the system size to include the ancilla modes, L → L + t/∆τ . This case also has a similar polynomial upper bound on the classical runtime in the form of Eq. (30) as long as the evolution time t is polynomial. Finally, let us analyze the case when the conditions of classes 1-3 are violated. Strictly speaking, then the stochastic method fails as it generally maps the problem to a nonquadratic fermionic evolution, which is not believed to be simulable for arbitrarily long time using only polynomial overhead in the number of fermion modes and computational time. . However, we can still efficiently simulate the system after this mapping if the product of evolution time t and the correction rate δΓ (of processes violating easiness conditions) remains small. In particular, if the product is bounded as δΓt < c/L 2 for some constant c, the dynamics remains classically easy. To obtain this result, we consider a more general form of stochastic unraveling in Eq. (15) with nonunitary unraveling. This formula can be Taylor expanded as V st → V st = V st + δΓτ K 1 + (δΓτ ) 2 K 2 + . . . , where K n are local correction superoperators and τ is the evolution time. Therefore, we can update the bound for the operator norms in Eq. (27) α involves at most k m fermion Fock operators and the action of K n involves at most four fermion operators, we see that K n D

III. HARD CLASS
We have so far demonstrated cases when the probability distribution generated by the Lindblad equation is efficiently simulable on a classical computer. Can we extend these proofs to the most general case of quadratic A k 's? Since quadratic operators A k correspond to single-fermion jumps in many cases, one may expect that the problem can be solved in the single-particle sector, similar to unitary free-fermion dynamics. However, such an intuition is incomplete. A simple explanation can be obtained using the Fermi exclusion principle that requires the transition between two modes to depend on the occupation of the target mode; thus a quadratic Lindbladian jump operator can induce many-body correlations in the system that quickly become classically intractable.

A. Reduction to a generic quantum circuit
We now provide a rigorous argument for worst-case hardness based on the equivalence of dynamics under classes H(t) and A k (t) in Eq. (1) on the one hand and universal quantum computing on the other. Let us start with the simplest map utilizing quadratic Hamiltonians. We distribute all modes into L/2 pairs, each pair corresponding to a logical qubit in the state |0 L = |01 or |1 L = |10 . Then, utilizing only quadratic Hamiltonians and Lindblad jump operators, we can implement any quantum circuit with arbitrary precision. Thus, by showing the equivalence of the dynamics to universal quantum computation, we obtain hardness results for both estimating time-evolved local observables O(t) and sampling from the time-evolved state in any local basis. The obtained hardness result is therefore on par with the best complexity-theoretic evidence that simulating quantum circuits (in both senses) is hard.
First, using single-fermion hopping between the two sites of a qubit, we can reproduce arbitrary single-qubit operations [77]. Second, to approximately generate a desired two-qubit gate, we can use hopping combined with a quadratic Lindblad operator. In particular, assigning the two-qubit logical states |00 L = |0101 , |01 L = |0110 , |10 L = |1001 , and |11 L = |1010 , the control-Z gate can be implemented by simultaneously applying the hopping Hamiltonian H = J(c † 2 c 3 + c † 3 c 2 ) and pair-loss operator A = Γc 3 c 4 for time t = π/J, in the limit Γ J. This type of dynamics can be analyzed as follows. The logical states |01 L and |10 L remain invariant in the course of the evolution. At the same time, in the limit γ ≡ Γ/J → ∞, due to the quantum Zeno effect, the Lindblad operator's action disallows any coherent transition involving states where qubits 3 and 4 are both occupied (i.e. | · ·11 ). As a result, the logical state |00 L is unaffected by the evolution. Therefore, the only evolving logical state is |11 L , which acquires a phase factor exp(iπ) = −1 after time t = π/J. As a result, the effective transformation on the two logical qubits is the control-Z gate Together with arbitrary single-qubit operations, the control-Z gate is enough to obtain dynamics universal for quantum computing and hence hard to approximately sample from, assuming standard conjectures in complexity theory [78,79]. Importantly, the performance of the dissipative gate relies on the Zeno-effect blockade effective for γ → ∞. In the limit of large but finite γ, the two-qubit system has the probability = 2π/γ + O(γ −2 ) of ending up in states |0011 or |0000 , which could result in an error in the gate (see Appendix F). To avoid computational error, we can choose the ratio γ to be arbitrarily large by taking vanishing J → poly −1 (L) for any given Γ > 0. Therefore, we can keep the error below any given threshold at the cost of increased overall computation time, which remains polynomial in system size.
The proposed architecture is not unique and allows for modified and generalized realizations of logical qubits and gates. For example, if the pair decay is always present on any two neighboring modes, one may introduce an empty ancilla mode between two logical qubits in order to ensure that logical states do not decay. As another example, if the control Hamiltonian is linear in terms of Majorana operators, a logical qubit can be encoded using just a single mode. Moreover, for a reader focused on applications, we discuss below a practical modification of qubit encoding implementable in cold atoms. Now let us show that pair loss is not the unique dissipation present in the hard class. In fact, this class also includes any quadratic dissipation connected to pair loss by a time-dependent linear Bogoliubov transformation, A (t) = ΓY (t)c 1 c 2 Y † (t), where Y (t) = exp(−iG(t)) is a freefermion unitary transformation and G(t) is a Hermitian operator from the quadratic-linear class in Eq. (2). To demonstrate Complexity phase diagram for a fermionic system with simultaneous pair losses and gain. The plot illustrates the connection between complexity of simulation and the hypothetical dissipative control-Z gate error + (both axes use log scale). When the error is smaller than the best-known two-qubit error-correction threshold p0, the worst-case system dynamics is equivalent to that of a fault-tolerant quantum computer (blue shaded regions) and, according to existing complexity conjectures, is classically computationally hard to simulate. In contrast, when the rates of gain and loss are exactly equal, the problem belongs to EC1 (vertical red line) with the effective classical algorithm provided in the text. The result for the unshaded region remains inconclusive. The dashed line represents qualitative extrapolation. this equivalence, we consider the pair-loss scheme described above but simultaneously replace all pair losses A with A (t), the Hamiltonian H(t) with H (t) = Y † (t)H(t)Y (t), and instead of the initial and final states, choose states transformed by Y (0) † and Y (t), respectively. The resulting process has the same probability distribution; thus its complexity would be the same. As a result, Lindbladians such as incoherent transitions A = Γc † 1 c 2 or pair gains A = Γc † 1 c † 2 are also classically hard in combination with free-fermion dynamics (see Table I).

B. Robustness of the hardness result
The error associated with imperfect Zeno blockade cannot be arbitrarily suppressed by slowing down the computation if there are small generic corrections to the dissipative dynamics. These corrections can be viewed as the presence of additional Lindblad jump operators with total rate Γ . Such terms generate additional transitions with the probability ∼ πΓ /J, where Γ is the combined rate of added operators A and/or other errors. In contrast to the imperfect-Zenoblockade error, this type of error diverges for small J. Therefore, there is an optimal value J ∼ Γ /Γ that minimizes the overall gate error to + ∼ O(1), including, besides standard errors, leakage into states outside of the logical Hilbert space. For fixed Γ, there always exists a choice of Γ ∼ O(1) that keeps the error below any provided threshold, + < p 0 , where p 0 > 0. According to the leakage threshold theorem in Ref. [80], which is a generalization of earlier standard threshold results [81][82][83], a universal set of such gates can be used to implement fault-tolerant quantum computing. Therefore, there are instances of Lindblad evolutions that remain hard to simulate for arbitrarily long times.
One particular example of a dissipative correction to ideal dynamics is the presence of pair gain A ij = Γ c † i c † j that acts on exactly the same sites as pair loss A ij . In this case, the minimum error is + = 8π 2 Γ /Γ and the problem remains hard for a classical computer if Γ ≤ p 2 0 /8π 2 Γ. Since the entangling gate is also implementable using pair gain instead of loss, this inequality also works after replacing Γ by Γ . Thus, the problem of simulating the evolution in the regions Γ /Γ ≤ p 2 0 /8π 2 and Γ /Γ ≥ 8π 2 /p 2 0 is classically hard. The complexity for the rest of parameter space remains an open problem. Notably, there exists at least one point in this range, Γ = Γ , that is easily simulable by a classical computer since it is in EC1. Therefore, by changing the ratio Γ /Γ, we can potentially induce a complexity phase transition. Figure  2 illustrates the connection between gate error and sampling complexity.
Summarizing, we established quantum computational universality of quadratic dissipation combined with free-fermion dynamics, where dissipation replaces the unitary interactions between fermions. This result opens a possibility of using simple dissipation processes as a resource for quantum computing. In the following section, we illustrate the feasibility of this proposal by considering a system of cold atoms.

IV. APPLICATION TO COLD ATOMS
In this section, we discuss an experimentally relevant system where naturally occurring dissipation can be a source of computational hardness. In particular, we study trapped cold fermionic atoms and consider nonunitary pair collisions as a viable dissipation mechanism. Collision dynamics can be controlled by a magnetic Feshbach resonance; tuning into the resonance can suppress unitary interactions and amplify the loss rate, thus physically implementing the Zeno regime we studied in the previous section. Therefore, we show that cold atoms harbor a natural way of implementing quantum computing using dissipation-assisted operations. While one may combine dissipation with elastic interactions to increase the fidelity of the gates, our complexity analysis shows that dissipation alone is sufficient. Below we provide details on how to implement collision control and to construct the gates. We also estimate the resulting gate error for a particular system.

A. Feshbach resonance
The Feshbach resonance provides a perfect tool to manipulate interactions between trapped atoms. Several mechanisms are available for practical use including magnetic, optical, and orbital Feshbach resonance [34,37,38]. For concreteness, we study only magnetic resonance here. The other two mechanisms have a qualitatively similar effect on atomic interactions. We study magnetic Feshbach resonance since it does not involve laser transitions and potentially has smaller scope for error.
We also require a fermionic atom that can be cooled, trapped and prepared in specific spin states with the requisite interaction properties. A promising example we illustrate here is the 40 K atom in its 2 S atomic ground state, which has an electron spin S = 1/2 and nuclear spin I = 4, giving rise to total spin f = 9/2 or 7/2. The Zeeman substructure of the ground-state hyperfine manifold, shown in Fig. 3(a), gives rise to magnetically tunable Feshbach resonances in the interaction of two atoms for controlling elastic and dissipative collisions [37,38]. This example serves as a proper illustration of how dissipation and interaction rates can be tuned, and alternative schemes for both alkali and alkaline-earth (like) atoms can be proposed using other types of Feshbach resonance.
It is straightforward to set up and numerically solve for the scattering and bound states of two 40 K atoms, including the atomic electrons and nuclear spins, their mutual interactions, and the mass-scaled adiabatic Born-Oppenheimer molecular potentials for the 1 Σ + g and 3 Σ + u states [85]. We use the standard coupled channels method [37,86] to set up the full spin Hamiltonian and solve the matrix Schrödinger equation for the scattering states. Such models, when calibrated against bound state and scattering data, provide highly accurate predictions of the properties associated with magnetically tunable Feshbach resonance states used to tune the scattering properties of two ultracold K atoms [85,[87][88][89][90][91][92][93].
The collision of two 40 K atoms is characterized by the quantum numbers of the two separated atoms with resultant total spin projection m F = m f 1 + m f 2 and relative angular momentum, or "partial wave," and m . States with the same total angular momentum projection M tot = m F + m are coupled in the molecular Hamiltonian for two atoms. Two like fermions collide with odd partial waves, e.g., the p wave with = 1, whereas two unlike fermions can collide via even or odd partial waves, including the s wave with = 0. A collision "channel" is defined by the partial wave and the spin quantum numbers of the two atoms. Strong pair loss via spin-exchange interactions is only possible if there is a product channel available with the same M tot and as the entrance channel; otherwise weak spin-dipolar spin relaxation is possible where m F and m change by two units, conserving M tot . Strong s-wave spin-exchange relaxation is possible for the spin channels 1 + 4 or 1 + 5, but not for 1 + 2 or 1 + 3 channels; furthermore, only weak s-wave spin-dipolar spin relaxation is possible in the 1 + 3 channel. We also do not consider the much weaker p-wave spin relaxation in these channels at ultracold temperatures (see below). Consequently, we now concentrate on the 1 + 4 and 1 + 5 s-wave collisions for engineering dissipative collisions with weak on-site unitary interaction.
Very-low-energy s-wave elastic and dissipative collisions in the threshold regime are adequately described by a complex scattering lengthã 0 [94,95], defined as the k → 0 limit of the energy-dependent complex scattering length [96][97][98][99], Here k is the relative collisional momentum for a collision of two atoms with reduced mass µ at energy E = 2 k 2 /(2µ), a note also that we use G (Gauss) as the magnetic field unit in this paper because of its near-universal usage among groups working in this field, 1 G=10 −4 T. and S(k) is the diagonal element of the unitary S-matrix for the collision channel in question. In this subsection, we keep track of for added clarity. The coupling constant for the lowenergy zero-range regularized pseudopotential approximation for atomic interactions is g = 2π 2 Re(ã 0 )/µ [37,38,100]. The dissipative loss rateṅ 1 =ṅ 2 = −K 2 n 1 n 2 from colliding atoms in a gas with densities n 1 and n 2 is given by the rate constant K 2 = −4π Im(ã 0 )/µ [96,97] (since Im(ã 0 ) is zero or negative, K 2 is positive definite; g can be positive or negative). Using counterpropagating laser beams, it is possible to construct an array of trapping cells in an optical lattice structure [101]. Each cell is approximately harmonic and, in its ground state, may hold exactly zero, one, or two atoms. The scattering length formulation can readily be adapted to two atoms in an optical lattice cell to calculate the interaction energy or dissipative loss rate. For a harmonic trap with frequency ν = ω/(2π), the analytic interaction energy for the lattice ground state from the zero-range pseudopotential is (3/2 + (2/ √ π) Re(ã 0 )/d) ω, where the harmonic length d = /(µω) [102]. If the lattice zero-point energy 3 ω/2 is large enough, Re(ã k ) may need to be evaluated at the lattice eigenenergy instead of taking Re(ã 0 ) in the k → 0 limit [103,104]. The decay rate Γ of an atom from the cell is given by K 2n , wheren = dr|Ψ 0 (r)| 4 = 1/(π 3/2 d 3 ) can be interpreted as a mean local density in the ground state of the lattice cell with wave function Ψ 0 (r) [105].
The figure of merit for our dissipative quantum gate, the opposite requirement from that of Ref. [105], is that |Im(ã 0 )/Re(ã 0 )| 1. This is possible to achieve using two 40 K atoms in states 1 and 4 or states 1 and 5, as we now show from our coupled channel calculations. Using the mass scaled potentials of Ref. [85] and including s and d waves ( = 0 and 2) in the coupled channels expansion for unlike spin species gives the scattering lengths shown in Fig. 3(b).
There are two regimes where the interaction energy proportional to Re(ã 0 ) is small and the dissipation rate proportional to Im(ã 0 ) is large. These are in the "core" of the resonance, rounded into a dispersive shape by the decay, and near the zero crossing where Re(ã 0 ) = 0.
In order to get a sense of time scales, we can assume a harmonic length on the order of 100 nm, for whichn ≈ 2 × 10 14 cm −3 . If we take the van der Waals length R vdW of two 40 K atoms, 3.4nm [37], as a "typical" size for Im(ã 0 ), then K 2 ≈ 1.3 × 10 −10 cm 3 /s, giving a decay time of Γ −1 = 40 µs. The next subsection discusses how such magnitudes could enable the realization of dissipation-assisted quantum computing.
We note that there are other spin channels for 40 K and in other species where Feshbach tuning of a favorable ratio Im(ã k /d)/Re(ã k /d) 1 could be feasible. This may be possible for like fermions, where only p-wave channels are available. However, p-wave interactions, treated by Eq. (32) with a p-wave S-matrix element, are typically suppressed by a factor on the order of k 2 R 2 vdW relative to the range of s-wave processes, due to the threshold law for p waves [97,[106][107][108]. This suppression factor is on the order of 0.001 for 40 K atoms with an energy on the order of 1 µK, so it would be harder to find ranges suitable for experimental control.

B. Dissipation-assisted quantum computing
We now modify the idea from Sec. III to make it feasible for cold atomic systems. We limit our attention to two distinct states of trapped atoms, denoted here as µ 1 and µ 2 . A single trap is described using the following four basis states: the empty state |0 i , single-occupied states |µ α i = c † iµα |0 i , α = 1, 2, and the double-occupied state |µ 1 µ 2 i = −|µ 2 µ 1 i = c † iµ1 c † iµ2 |0 i . We consider the lattice Hamilto- The quantities J µ ij (t) are real tunneling amplitudes, and ∆ µ i (t) are on-site potentials, both of which can depend on the atomic state µ. Two distinct atoms located in the same trap are subject to elastic interactions V = E i n iµ1 n iµ2 , where E is the interaction energy and n iµ = c † iµ c iµ is the µ-occupation number of the ith trap. As shown in Sec. IV A, the interaction E may be made to vanish for a specific pair of states µ 1 and µ 2 by, for example, manipulating the magnetic field as shown in Fig. 3(c). Also, atoms in the same trap undergo pair loss with rate Γ, described by Lindblad jump operators A i = Γc iµ1 c iµ2 . For E = 0, the entire dynamics is described by the master equation Eq. (1).
The computational scheme utilizes pairs of sites to encode individual logical qubits. The logical qubit states are |0 L = |0 |µ j and |1 L = |µ j |0 , irrespective of the atom's type µ j . Single-qubit gates can be performed using the local potentials and coherent hopping between logical qubit sites.
Previously, we discussed how to use Feschbach resonance to induce and control the dissipation of alkali atoms using the example of 40 K. We now step aside to describe how a generic entangling gate works for a larger group of atoms, assuming that the same level of control can be applied to them as well. We consider two distinct ways of constructing entangling gates, depending on the atomic electronic structure. The first method we consider is designed for alkalineearth(like) atoms such as 87 Sr [33] and 173 Yb [31,32]. We can use nuclear-spin polarized metastable states 1 S 0 and 3 P 0 as the two species µ 1 and µ 2 [24] in order to apply a speciesdependent hopping term J µ ij . The scattering length between µ 1 and µ 2 can then be potentially tuned by optical or orbital Feshbach resonances to be purely imaginary. Alternatively, one could use 3 P 2 instead of 3 P 0 , in which case a magnetic Feshbach resonance is also an option. However, this method has limited applicability to alkali atoms, for example 87 Rb [26] or 40 K (described above), where the states µ i are encoded into different angular momentum projections m f . This is because a state-dependent lattice in alkali atoms [27] can exhibit significant single-atom dissipation rates. For alkali atoms, therefore, we propose an additional scheme that does not rely on an internal-state-dependent lattice and uses the same lattice potential for both states. In order to implement entangling gates, we make use of single-qubit rotations with single-site resolution. This can be achieved using two-photon Raman transitions induced by focused laser beams or other similar techniques [28,109,110]. Both schemes can be used interchangeably.
The performance of the gate can be disrupted by errors, including imperfect Zeno-effect error, the single-qubit phase gate error (for the second scheme), and background dissipation error. The background dissipation error can be bounded above by Γ t, where Γ is the background dissipation rate, and t = π/J is the gate performance time (neglecting the time taken for the single-qubit phase gate). The single-qubit phase gate error 0 is fundamentally limited by light scattering loss during the Raman transition, which depends on the characteristic linewidth γ of the excited levels and the detuning limited by fine structure splitting ∆. This error is estimated to be 0 ∼ γ/∆. Deviations from perfect single-site addressability during the single-qubit phase gate can also give rise to errors, which can nevertheless be greatly reduced by subwavelength addressability techniques [27,28,[110][111][112][113][114]. Finally, the error caused by an imperfect Zeno effect can be approximated by the first term in the Taylor expansion of the infidelity in (Γt) −1 (see Appendix F), leading to the total error where ζ = E/Γ is the loss-to-interaction ratio. The error can be minimized by making the choice t = 2π/ ΓΓ (1 + 4ζ 2 ), leading to the expression .
The dependence of the second term on the magnetic field is illustrated in Fig. 4(b) under the same choice of parameters as in Fig. 3. For a suitable choice of magnetic field strength, the theoretical upper bound for the error can be as low as ∼ 5 × 10 −4 for the background dissipation rate Γ = 10 −2 s −1 . For 40 K atoms, the optical transition error 0 can be estimated using the values γ 2π × 6.0 MHz and ∆ 2π × 1.7 THz [115], leading to the upper bound 0 ∼ 10 −6 , which is an insignificant contribution to the overall error. As a result, the theoretical bound for the gate error approaches the characteristic thresholds given by many errorcorrecting schemes.
From the comparison of elastic and inelastic rates in Fig. 3(c), one can see the advantage of dissipative schemes over unitary ones. While the gate time for both elastic and inelastic schemes is proportional to the corresponding inverse collision rate, the resonant dissipation rate can be significantly larger than the accessible elastic rates. Therefore, dissipationassisted gates can be faster and can thus experience smaller level of errors due to the background noise.
Summarizing, we have established that naturally occurring pair-loss processes in cold atoms can be enhanced and used as a resource for quantum computing. This conclusion opens a possibility to improve quantum computing with cold atoms [24] or, in absence of sufficient control, use them as a platform for quantum supremacy experiments.

V. DISCUSSION
In this work, we have demonstrated how simple forms of dissipation affect the complexity of simulation of noninteracting fermions. In particular, focusing on linear-quadratic Lindblad jump operators, we have shown the existence of two complementary complexity classes of Lindblad jump operators, easy and hard for simulation on a classical computer. Using the error-correction formalism, we showed that the hard class has a finite volume in the parameter space and tolerates the presence of small arbitrary corrections. At the same time, the easy classes may have small measure and could become hard even as a result of arbitrarily small corrections to the master equation.
We have expanded the region of classical simulability of free-fermions in the presence of Markovian errors from single-qubit loss/gain to more general quadratic-linear Lindblad jump operators. The algorithms we devise for EC1-EC3 based on the stochastic unraveling approach provably work in polynomial time. This shows that a large class of dissipation processes such as dephasing or single-fermion decay can be treated with the help of efficient classical algorithms.
At the same time, more complex processes are BQPcomplete, where BQP stands for the class Bounded-error Quantum Polynomial time. We show this fact by explic-itly constructing an entangling gate and showing the equivalence of the problem with universal quantum computation. We thus place limitations on the extent to which the simulability result may be extended, since we believe quantum computation is strictly more powerful than classical computation. Our detailed analysis shows that it is within the range of experimental feasibility to implement with cold atoms a quantum computer with purely dissipative atom-atom interactions, an exciting possibility for experiments in quantum computing. For example, dissipative quantum systems such as alkaline-earth atoms may serve in the next generation of quantum supremacy experiments. Also, our result suggests that simulating fermion dynamics may be hard for quantum particles experiencing dissipation, for example, quasiparticles in solid-state systems. Future work can explore the hardness of simulation of electronic systems with quasiparticle dynamics approximated with quadratic-linear Lindblad jump operators that include the effects of electron-electron, electronphonon, and electron-impurity scattering processes. Alternatively, physical systems following such dynamics with high accuracy may be a future platform for quantum-computing experiments.
It may be interesting to explore the connection of our results with the theory of matchgate (free-fermion) computations and the role played by non-Gaussianity. Quadratic fermionic Hamiltonians and single-fermion loss give rise to Gaussian operations and are hence easily simulable [17]. It is known [116] that any non-Gaussian fermionic state is a resource for fermionic computation, boosting the computational power of free fermions from being classically simulable to being universal for quantum computation. Our results suggest that quadratic-linear Lindblad jump operators are non-Gaussian in general. Therefore, it would be interesting to quantify the amount of non-Gaussianity (or "magic") for the Lindblad operations we study here.
Along the same lines, one can quantify a different resource for nonclassicality, such as a suitable measure of entanglement for open fermionic systems. Efficient sampling from the full output distribution in arbitrary bases can allow for efficient computation of certain measures of entanglement such as Rényi entropies [117]. Relatedly, Ref. [118] has studied the logarithmic negativity for free fermions with gain/loss Lindblad terms.
Further, one may also consider how the complexity of simulating dynamics under quadratic Lindbladians changes with time. Since the system starts off in a Fock state that is easy to sample from, and dynamics under quadratic Hamiltonians with quadratic Lindblad jump operators can generate states that are hard to sample from, one can see a dynamical transition in sampling complexity [13,119]. It is worthwhile to investigate whether these transitions are sharp or coarse (as defined in Ref. [119]) since this can identify what "universality class" free fermions with noise belong to. Techniques such as Lieb-Robinson-like bounds for the evolution of free particles with dissipation [120,121] would be relevant here.
Another exciting direction is the study of worst-to-averagecase equivalence in complexity, which seeks to understand the complexity of typical instances as opposed to worst-case in-stances [79,122,123]. It would be interesting to see if the Cayley path technique of Ref. [122] can be adapted to argue for average-case hardness of dissipative fermionic dynamics.
Lastly, let us address the case of bosonic particles. First, in the limit of strong dissipation, our result on computational universality applies readily to bosons. This conclusion combines with the fact that the proposed Feshbach realization scheme is generally similar for bosons and fermions. The major difference arises for identical atoms, since identical bosons also feature strong and broad resonances in s-wave scattering, whereas only p-wave resonances exist for fermions due to the Pauli principle (see, for example, Ref. [124]). Thus, the strong dissipation regime is easier to achieve for identical bosons in general. However, in the case of distinct fermions, e.g. in different spin states, there is essentially no difference between bosons and fermions as a class [125].
Next, in contrast to hardness results, classical easiness results for bosons are rather distinct from fermions. For Fockbasis measurements, free bosons are already believed to be hard to simulate in the sense of sampling from the output distribution [78]. This hardness equally holds for some problem variations such as using Gaussian initial states with Fock measurements [126,127] or homodyne measurements for non-Gaussian initial states [128]. However, since free boson evolution is not believed to be computationally universal, our results imply that quadratic dissipation boosts the computational power of free bosons in these settings to that of computational universality. Thus, it leads to another type of complexity transition between two classically hard classes. This transition is an interesting object for future studies.
Notably, there exist settings where the free-boson evolution is classically easy to simulate; an example is the situation when the system is accessed by homodyne measurements in combination with Gaussian initial states. Moreover, in this example, adding either single-boson or even many-body decay processes preserve classical easiness of free-boson evolution, unlike the fermionic cases considered in our work. However, it is also an open question whether or not generic dissipative processes involving boson creation can induce a complexity transition. The analysis of these processes can be done using previous works addressing the complexity of non-Gaussian bosonic states. For example, Ref. [129] suggests specific non-Gaussian states that can be useful as computational resources when making homodyne measurements, while Ref. [130] analyzes photon loss from the point of view of computational complexity (with Fock-basis measurements). We keep studying this problem as well as analyzing the robustness of the resulting phases for future research. In this Appendix, we analyze the convergence of the average unitary stochastic evolution to the exact Lindblad dynamics in the case of easy class 1 (EC1). First, we set the initial time to be zero and consider the final time t being an integer multiple of timestep ∆τ . This assumption holds without loss of generality since ∆τ may be adjusted appropriately to capture any particular final time. Then the overall evolution of unitary can be written as a product where the timestep unitary U n is expressed in terms of a timeordered exponential generated by the stochastic Hamiltonian H (t) in Eq. (13), Here H(t) is the original time-dependent Hamiltonian, and R(t) = k ξ nk A k (t)+ξ * nk A † k (t) is the normalized stochastic part, where ξ nk are independent complex Gaussian variables defined for times n∆τ ≤ t ≤ (n + 1)∆τ .
Let us consider the ordered exponential expansion of the timestep unitary in Eq. (A2): where we denote the discretized value of an operator O n and permutation sum, respectively, as The average over the stochastic field can be taken for each timestep independently. Therefore, the effect of the timestep unitary in Eq. (A4) is EU n ρU † n = I + L n ∆τ + 1 2 L 2 n ∆τ 2 ρ + D n ρ∆τ 2 + O(∆τ 3 ).

(A6)
In the equation above, L n is the generator of the original Lindblad equation, lim ∆τ →0 L n = L(n∆τ ), expressed as and D n represents the lowest-order correction occurring due to the timestep being nonzero: Here we use the notation The overall expression in Eq. (A8) can be written in a compact form, where D (i) αn = poly(H n , A kn ) are polynomials of degree less than four.
The averaged stochastic map in Eq. (A6) can be rewritten as a continuous evolution and then decomposed using Dyson series for the small parameter ∆τ , where the generators L(t) and D(t) are continuous versions of the operators in Eq. (A7) and Eq. (A8), in which the ∆τaveraged operators A kn and H n are replaced by the corresponding instantaneous values at time t, i.e. A(t) and H(t), respectively. To obtain the expression in Eq. (17), we recursively replace V(t 2 , t ) and V(t , t 1 ) on the right-hand side by their stochastic average and collect all O(∆τ 2 ) terms. In this Appendix, we analyze the convergence of the average stochastic unitary evolution to the Lindblad dynamics in the case of easy class 2 (EC2). The single timestep evolution averaged over stochastic unitaries in Eq. (19) is equivalent to the map where the target Liouville operator is The correction now takes the form denoting C kn = Γ n A † kn + i 2 [A † kn , H n ] and Γ n = ∆τ −1 n dt k Γ k (t). This expression has the form of Eq. (A10) with operators D (i) αn being a sum of products of at most four Majorana fermion operators.
Appendix C: Easy Class 3 In this Appendix, we analyze easy class 3 (EC3) and show the convergence of the system-ancilla stochastic evolution under the Hamiltonian in Eq. (13) using the stochastic operators in Eq. (20) to the dissipative dynamics with linear Lindblad jump operators. Let us start from a many-body pure state of the fermions occupying L modes of the system and L a ancilla modes at time t = n∆τ , denoting it as |Ψ n . At the nth timestep, the evolution acts on the system and the nth ancilla mode only. Thus, the state at time t = n∆τ is a product state of subsystem states: (1) correlated state of L system modes together with the first n ancilla modes and (2) the product states of the remaining L a − n ancilla modes, i.e. |Ψ n = |φ n L+n ⊗ |0 La−n . (C1) The evolution is governed by the Hamiltonian where the stochastic terms are K(t) = k f nk A k (t)(γ 2(L+n) + iγ 2(L+n) ) (C3) at times n∆τ ≤ t ≤ (n + 1)∆τ , and f nk are independent real Gaussian variables. Then, at the (n + 1)th step, the systemancilla state |Ψ n = U n |Ψ n is |Ψ n+1 = |Ψ n − i∆τ 1/2 K n |φ n |1 |0 La−n−1 − ∆τ iH n + 1 2 K † n K n |φ n |0 La−n − ∆τ 3/2 1 2 {H n , K n } − i 6 K n K † n K n |φ n |1 |0 La−n−1 where we used the discrete-time operator values H n and K n obtained as in Eq. (A4). The interpolated continuous-time evolution for the density matrix of the system can be presented in the form where x is the floor function. The target Liouville operator is and the correction is where As is the case for EC1 and EC2, the correction is described by Eq. (A10) with operators D where L is the system size and J is a real hopping coefficient. For simplicity, we set an equilibrium thermal state ρ(0) = Z −1 exp(−βH) as the (Gaussian) initial state, where Z is the partition function and β is an inverse temperature. As a model of quadratic jump dissipation from EC1, we consider 2L − 2 independent incoherent jump processes between adjacent sites, where Γ q is the jump rate. As required by EC1, we make these rates equal for the processes going in opposite directions (hopping to the left and to the right). As an example from EC2, we consider L − 1 unitary jump operators where J u is a real hopping coefficient and Γ u is the corresponding dissipation rate. Finally, as an example from EC3, we use the 2L linear dissipation operators with F s 2n = Γ s c n , F s 2n+1 = Γ s c † n , where Γ s = Γ s are fermion decay/gain rates.
For all three cases, we analyze the deviation of the 4-point correlation function from the prediction given by Wick's theorem, ε(t) = γ 1 γ 2 γ 3 γ 4 t − γ 1 γ 2 t γ 3 γ 4 t + γ 1 γ 3 t γ 2 γ 4 t − γ 1 γ 4 t γ 2 γ 3 t , where we define · t := Tr(·ρ(t)). The result is shown in Fig. 5. As we can see, Wick's theorem is violated for EC1 and EC2 at times t > 0. The process from EC3 alone remains Gaussian, as discussed previously [16,17]. The inapplicability of Wick's theorem for EC1 and EC2 implies that the simulability results of Refs. [67][68][69] do not lead to classical easiness of sampling from the full output distribution when measuring the time-evolved state in the fermion-number basis.
The restriction of the effective Hamiltonian to the subspace spanned by the basis {|01 , |P 1 , |P 2 , |EX } is where E is the interaction energy of a fermion pair on the same site. The non-Hermitian nature of the Hamiltonian reflects additional leakage to the fully empty state |V C due to pair loss. In the strong dissipation limit J Γ, the leakage error for the gate can be defined as where S = exp(−iH S t), t = π/J is the characteristic time of hopping between lattice sites, and ζ = E/Γ is the ratio between the interaction energy and the pair loss rate.