Limits on Classical Simulation of Free Fermions with Dissipation

Complexity of Fermionic Dissipative Interactions and Applications to Quantum Computing Oles Shtanko,1, 2 Abhinav Deshpande,1, 2 Paul S. Julienne,2 and Alexey V. Gorshkov1, 2 arXiv:2005.10840v2 [quant-ph] 17 Sep 2021 1 Joint Center for Quantum Information and Computer Science, NIST/University of Maryland, College Park, Maryland 20742, USA 2 Joint Quantum Institute, NIST/University of Maryland, College Park, Maryland 20742, USA 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–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–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 H(t), Ak(t) ψr' 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 Ak (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–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]. One motivation behind this work is the existence of a variety of accessible fermionic physical systems involving inelastic processes. Examples of dissipative processes described by quadratic Lindblad jump operators include two-body loss in 2 trapped alkali atoms [26–28], alkaline-earth atoms [29–34], and cold molecules [35, 36]. As we show, Feshbach resonances [37, 38] can be used to significantly suppress coherent interactions between cold atoms, simultaneously increasing the rate of atom-pair trap losses. More general types of dissipation can be created by adding a source of atoms [39–41] or inelastic photon scattering [42–44]. In solid-state physics, examples of processes described by quadratic jump operators in Markovian approximation are Cooper-pair loss [45, 46] and phonon-induced dephasing [47]. Recent progress in the control of dissipative electronic systems has brought them into focus in condensed-matter physics. Some of the novel effects in noninteracting and mean-field fermionic systems include dissipation-induced magnetism [48–50], dissipative superfluids and superconductors [51, 52], dissipative Kondo effect [53, 54], non-Hermitian topological phases [55–62], and non-Hermitian localization [63–65]. 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– 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] k A X dρ 1 = −i[H(t), ρ] + Ak (t)ρA†k (t) − {A†k (t)Ak (t), ρ}, dt 2 k=1 (1) where {X, Y } ≡ XY +Y X is the anticommutator, ρ(t) is the density matrix of the system, H(t) is a noninteracting Hamiltonian, and Ak (t) ∈ A(t) form a set of kA 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 ρ(t2 ) = V(t2 , t1 )ρ(t1 ) between arbitrary times t1 and t2 ≥ t1 satisfies V(t2 , t1 ) = V(t2 , τ )V(τ, t1 ) for any t2 ≥ τ ≥ t1 . 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 Ak ⇄ { pAk , 1 − pAk }, 0 ≤ p ≤ 1. 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 3 Type Examples of Ak Complexity Dephasing Particle shuffle Classical fluctuations Classical pair fluctuations c†1 c1 c†1 c2 & c†2 c1 c†1 & c1 c†1 c†2 & c1 c2 Easy (EC1) Mixing unitaries 2c†1 c1 − 1 + i(c†2 + c2 ) Easy (EC2) Single-particle loss/gain c1 OR c†1 Incoherent hopping Pair loss/gain c†1 c2 Easy (EC3) Hard c1 c1 OR c†1 c†1 TABLE I. Comparison between different types of noninteracting fermion dynamics with additional dissipation. For simplicity, we provide examples for two modes out of L, denoted by numbers 1 and 2. The symbol & means that both operators are present in the set A(t) with factors equal in absolute value. Abbreviations EC1, EC2, EC3 stand for easy class 1, 2, and 3 described in the text. Note: Different classes can be combined by summing the right-hand sides of the corresponding master equations. However, the jump operators cannot be combined without precaution: the sum of two easy jump operators from different classes does not necessarily produce a simulable jump operator. any operator, including the Hamiltonian and Lindblad jump operators, in terms of fermionic Fock operators c†n and cn , n ∈ {0, 1, . . . L − 1}. Fock operators create and annihilate a single fermion in a particular mode and satisfy the canonical commutation relations {cn , cm } = 0, {c†n , cm } = δ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 = cn + c†n and γ2n+1 = −i(cn − 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] H(t) = i 2 2L−1 X αij (t)γi γj + i,j=0 2L−1 X Ak (t) = βi (t)γi , (2) i=0 where ak (t) are antisymmetric 2L × 2L matrices, bk (t) are 2L vectors, and dk (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, P (r|r′ ) = hψr |ρ(t)|ψr i, 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 Ak (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 ρ(0) = |ψr′ ihψr′ |, (3) where the vectors r′ and r denote the positions of occupied modes in the initial and final (measured) product-state configurations, respectively, and |ψr i is a product state defined as |ψr i = c†r1 . . . c†rN |0i = γ2r1 . . . γ2rN |0i. Here |0i is the vacuum state defined as the state satisfying cn |0i = 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 6= Ñ . 2L−1 2L−1 X i X ij bik (t)γi + dk (t), ak (t)γi γj + 2 i,j=0 i=0 (4) where ak (t) and bk (t) are arbitrary complex-valued 2L × 2L matrices and 2L-vectors respectively, and dk (t) is a number. Notably, the same master equation allows representation with more than one set of jump operators A = {Ak }. A set can be reduced to a smaller one if its jump operators are linearly dependent [15]. Therefore, the number kA 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 ak (t), bk (t), and dk (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 4 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 |+i ≡ (|0i + |1i)/ 2 during the entire evolution, including the initial and final times, i.e. |ψr′ i → |ψr′ i ⊗ |+i, |ψr i → |ψr i ⊗ |+i. (5) 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 αij → α̃ij = αij + δi2L βj − δj2L βi , (6) 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: P (r|r′ ) = 1 2 X s,s′ ∈{0,1} P ({r, s}|{r′ , s′ }). (7) 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 Ä R tfermion opä † erator γi (t) = Ut γi Ut , where Ut = T exp −i 0 H(t′ )dt′ . Here T exp is the standard time-ordered exponential. Given the quadratic structure of the Hamiltonian, this evolution is P a linearÄ transformationä γi (t) = i Rij (t)γj , where R = Rt ′ ′ T exp −2 0 α(t )dt is a unitary 2L × 2L matrix. One can use this expression to derive the time evolution of a fermion operator as Ut cn Ut† = X 1 Tnj γj , Ut (γ2n + iγ2n+1 )Ut† = 2 j (8) where Tnj ≡ R2n,j + iR2n+1,j are elements of a L × 2L transformation matrix T . Labeling the initially empty sites as li′ and recalling that the initial fermion positions are ri′ and that the final positions are ri , the linearity allows to write the output probability in Eq. (3) at any time as P (r|r′ ) = hψr | Ut |ψr′ ihψr′ | Ut† |ψr i = hψr |Ut c†r′ cr1′ . . . cl′ c†l′ Ut† |ψr i L−N 1 L−N X ∗ ∗ ′ = Tr1′ n1 Tr1′ m1 . . . TlL−N m L Tl ′ L−N nL n1 ,...nL ;m1 ,...mL−N × h0|γ2rN . . . γ2r1 γn1 γm1 . . . γmL γnL γ2r1 . . . γ2rN |0i. (9) 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 A[I]ij ≡ AIi ,Ij . Consider the set I = {ri′ , L + lj′ , 2L + 2rk }, where i ∈ {1, 2, . . . N }, j ∈ {1, 2, . . . L − N }, and k ∈ {1, 2, . . . Ñ } take all possible values. Then the result can be written as [6] P (r|r′ ) = Pf M [I], where Pf is the Pfaffian, and M is a 4L × 4L matrix é Ñ T ΛT T T ΛT † T Λ M = T ∗ ΛT T T ∗ ΛT † T ∗ Λ , ΛT T ΛT † I where, in turn, the 2L × 2L matrix Λ is Å ã 1 i . Λ = IL×L ⊗ −i 1 (10) (11) (12) 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 5 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, t2 /ǫ 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 = A1 ∪ A2 , where the Hermitian conjugate of every Lindblad operator in A1 returns an operator in A2 (and vice versa). Under this division, any Hermitian Lindblad operator must be included in both subsets √ A1 and A2 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–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 nonGaussian 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 X θk (t)Ak (t) + θk† (t)A†k (t), (13) H ′ (t) = H(t) + Ak ∈A1 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(t2 , t1 ) generated by the Lindblad equation can be approximated as V(t2 , t1 ) = EVst (t2 , t1 ) + δV(t2 , t1 )∆τ + O(∆τ 2 ), (14) where δV(t2 , t1 )∆τ is the lowest-order correction (to be explicitly derived below) and Vst is a stochastic unitary map Vst (t2 , t1 )ρ = U (t2 , t1 )ρU † (t2 , t1 ). (15) ä Ä Rt In the above, U (t2 , t1 ) = T exp −i t12 dt′ H ′ (t′ ) 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 P (r|r′ ) = E Pst (r|r′ ) + O(∆τ ), (16) where Pst (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 Z t2 (17) δV(t2 , t1 ) = E dt′ Vst (t2 , t′ )D(t′ )Vst (t′ , t1 ), t1 where D(t) is a time-local superoperator X D(t)ρ = Dα(1) (t)ρDα(2) (t). (18) α (i) Here, the operators Dα (t) = poly4 (H(t), Ak (t)) can be expressed as polynomials of degree four in terms of the Hamiltonian and Lindblad jump operators at time t. Therefore, (i) Dα (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 operp ators Ak = Γk (t)Yk (t), where Γk (t) ≥ 0 are timedependent rates and Yk (t) = exp(−iGk (t)) are unitary operators generated by quadratic-linear Hamiltonians Gk (t) of the form in Eq. (2). To classically simulate dynamics under EC2, we also consider discretized time evolution with sufficiently smallQ timesteps ∆τ and set the unitary transformation n2 Un , where the timestep unitaries Un are U (t1 , t2 ) = n=n 1 generated randomly according to the rule ® Yk (tn ), pk = Γk (tn )∆τ, 0 P Un = Un × (19) I, p0 = 1 − k Γk (tn )∆τ. 6 Here pk are probabilities that are to generate the Ä used ä reR (n+1)∆τ spective outcomes, Un0 = T exp −i n∆τ H(t)dt , and tn ∈ [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 oper(i) ators Dα (t) here are degree-two polynomials in the Hamiltonian and Lindblad jump operators. The detailed form of the operators along with the derivation can be found in Appendix B. Finally, let us consider easy class 3 described by generic P linear Lindblad jump operators Ak (t) = i bik γi + dk , which can be obtained by setting ak = 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 La = 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 + La − 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 θk (t) = ξnk ∆τ −1/2 (γ2(L+n) + iγ2(L+n) ). (20) 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 (i) Eq. (14) and Eq. (18), with Dα (t) expressed as degree-four polynomials in the Hamiltonian and Lindblad jump operators, as shown in Appendix C. 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 distribution P̃ (r|r′ ) different from the ideal distribution P (r|r′ ), we utilize the total variation distance ǫ= X 1 P̃ (r|r′ ) − P (r|r′ ) , max ′ 2 r r (21) 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), XZ t ∆τ ǫ≤ (22) dt′ Crα′ (t, t′ ) + O(∆τ 2 ), max 2 r′ α 0 where Crα′ (t, t′ ) = E X r hψr |Dα(1) (t, t′ )ρr′ (t)Dα(2) (t, t′ )|ψr i . (23) Here = Vst (t, t and ρr′ (t) = Vst (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. (i) Dα (t, t′ ) ′ (i) )Dα (t′ ) Lemma. Consider two sparse operators O1 and O2 whose matrix elements satisfy hψr |Oα |ψr′ i = 0 if dH (r, r′ ) ≥ kα , α ∈ {1, 2}, (24) where dH (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 X r |hψr |O1 ρO2 |ψr i| ≤ Lk1 +k2 kO1 kmax kO2 kmax , (25) k1 !k2 ! where kOα kmax 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 Crα′ (t, t′ ) ≤ Lk1α +k2α EkDα(1) (t, t′ )kmax kDα(2) (t, t′ )kmax , k1α !k2α ! (26) (i) where kiα is the locality of the operator Dα (t, t′ ), i.e. the maximum number of Majorana operators in its decomposition. Because Vst is a map describing free-fermion evolution, 7 (i) the locality kiα of the operator Dα (t, t′ ) is equal to the local(i) ity of Dα (t′ ). At the same time, as analyzed in the previous (i) section, the localities of operators Dα (t′ ) satisfy kµα ≤ km , where km = 8 for EC1/EC3, and km = 4 for EC2. We can also bound the max-norm by the (spectral) operator norm kDα(i) (t, t′ )kmax ≤ kDα(i) (t, t′ )k = kDα(i) (t′ )k. As a result, the error bound is given by Z ∆τ L2km X t ′ (1) ′ E dt kDα (t )kkDα(2) (t′ )k. ǫ≤ 2 (km !)2 α 0 (27) (28) (i) Since the matrices Dα = poly(H, Ak ) are generated by a quadratic-linear Hamiltonian H and set of Lindbladians Ak , we can always find a polynomially large bound for the norm (i) kDα (t′ )k ≤ poly(L). Thus, there always exists a discretization step ∆τ ≤ ǫ t poly(L) (29) 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. (8), the number of steps to compute the distribution is equal to C ′ = L × O(L3 ) = O(L4 ), where the factor O(L3 ) is the upper bound on the time it takes to compute a Pfaffian of an O(L) × O(L) matrix. Further, the runtime for obtaining the matrix T is proportional to t/∆τ × M (2L), where M (n) . O(n3 ) is the time for n × n matrix multiplication. In sum, the total bound on the runtime for
each trajectory is bounded as C ∼ O(L4 ) + O L3 t/∆τ . Choosing ∆τ = ǫ/(t × poly(L)), the runtime is Å ã t2 C ≤ poly L, . (30) ǫ 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/L2 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 ′ Vst → Vst = Vst + δΓτ K1 + (δΓτ )2 K2 + . . . , where Kn are local correction superoperators and τ is the evolution time. Therefore, we can update the bound for the op(i) (i) erator norms in Eq. (27) as kDα (t, t′ )k ≤ kDα (t′ )k + (i) ′ (i) δΓτ kK1 Dα (t )k + (δΓτ )2 kK2 Dα (t′ )k + . . . . Since the (i) operator Dα involves at most km fermion Fock operators and the action of Kn involves at most four fermion opera(i) tors, we see that kKn Dα (t)k ≤ O(Lkm +2n ). Therefore, the (i) norm kDα (t, t′ )k is always bounded by a poly(L) value if δΓ < c/(L2 τ ), where τ = t − t′ ≤ t. This result leads to Eq. (30). As a result, if the dissipation is close enough to the symmetric point, the evolution remains classically easy. This result may be helpful for analyzing the precision needed for implementing this dynamics in intermediate-scale quantum devices. 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 Ak ’s? Since quadratic operators Ak 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 Ak (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 |0iL = |01i or |1iL = |10i. 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 hO(t)i 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 |00iL = |0101i, |01iL = |0110i, |10iL = |1001i, and |11iL = |1010i, the control-Z gate can be implemented by simultaneously applying the hopping Hamiltonian H = J(c†2 c3 + c†3 c2 ) and pair-loss operator A = Γc3 c4 for time t = π/J, in the limit Γ ≫ J. This type of dynamics can be analyzed as follows. The logical states |01iL and |10iL 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. | · ·11i). As a result, the logical state |00iL is unaffected by the evolution. Therefore, the only evolving logical state is |11iL , 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 è Ö 1 0 0 0 0 1 0 0 . (31) |ψi → Uπ |ψi, Uπ = 0 0 1 0 0 0 0 −1 Together with arbitrary single-qubit operations, the controlZ 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 |0011i or |0000i, 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)c1 c2 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 Error ϵ+ϵ' 8 Hard Easy Hard p0 p02/8π2 1 8π2/p02 Pair loss-to-gain ratio Γ/Γ FIG. 2. 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 c2 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. Therep fore, 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, ǫ+ǫ′ < p0 , where p0 > 0. According to the leakage threshold theorem in Ref. [80], which is a generalization of earlier standard threshold results [81–83], a universal set of such gates can be used to implement fault-tolerant quantum computing. Therefore, 9 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 p pair loss Aij . In this case, the minimum error is ǫ + ǫ′ = 8π 2 Γ′ /Γ and the problem remains hard for a classical computer if Γ′ ≤ p20 /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 Γ′ /Γ ≤ p20 /8π 2 and Γ′ /Γ ≥ 8π 2 /p20 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 3 + potentials for the 1 Σ+ g and Σ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–93]. The collision of two 40 K atoms is characterized by the quantum numbers of the two separated atoms with resultant total spin projection mF = mf 1 + mf 2 and relative angular momentum, or “partial wave,” ℓ and mℓ . States with the same total angular momentum projection Mtot = mF + 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 Mtot and ℓ as the entrance channel; otherwise weak spin-dipolar spin relaxation is possible where mF and mℓ change by two units, conserving Mtot . 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–99], ãk = 1 1 − S(k) . ik 1 + S(k) (32) Here ~k is the relative collisional momentum for a collision of two atoms with reduced mass µ at energy E = ~2 k 2 /(2µ), 10 (b) 2 E/h (GHz) 1.5 f=7/2 1 Scattering length (rB) mf=-7/2 mf=7/2 0.5 mf=9/2 f=9/2 0 mf=-9/2 -0.5 0 100 200 6 45 213 300 500 400 B (Gauss) Magnetic field B (Gauss) (c) 800 1+4 (Re) 1+4 (-Im) 1+5 (Re) 1+5 (-Im) 600 400 200 0 220 800 Re -Im 600 400 200 0 200 200 600 Scattering length (rB) (a) 240 260 280 300 Magnetic field B (Gauss) 320 240 250 260 Magnetic field B (Gauss) FIG. 3. Magnetic Feshbach resonance for 40 K atoms. (a) The hyperfine-Zeeman energy levels E/h (GHz) of the f = 9/2 and 7/2 manifolds versus magnetic field B (Gauss) a , labeled 1, 2, . . . in order of increasing energy. The levels labeled 1, 2, . . . have spin projections mf = −9/2, −7/2, . . . , respectively. The spin and Zeeman coupling parameters are taken from Ref. [84]. (b) Scattering lengths of two 40 K atoms. The solid and dashed lines represent − Im(ã0 ) ∝ K2 and Re(ã0 ), respectively; blue and red represent the 1 + 4 and 1 + 5 channels, respectively. (c) Magnetic Feshbach resonance for the 1 + 4 channel. The shaded regions depict magnetic fields where the elastic interaction between two atoms is smaller than the pair-loss rate, marking the regime of dissipative fermionic dynamics. All lengths are provided in Bohr-radius units rB = 5.29 × 10−9 cm, and the collision energy is E/kB = 1µK, where kB is the Boltzmann constant. 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 = −K2 n1 n2 from colliding atoms in a gas with densities n1 and n2 is given by the rate constant K2 = −4π~ Im(ã0 )/µ [96, 97] (since Im(ã0 ) is zero or negative, K2 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 p 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 R rate Γ of an atom from the cell is given by K2 n̄, where n̄ = dr|Ψ0 (r)|4 = 1/(π 3/2 d3 ) 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 which n̄ ≈ 2 × 1014 cm−3 . If we take the van der Waals length RvdW of two 40 K atoms, 3.4nm [37], as a “typical” size for Im(ã0 ), then K2 ≈ 1.3 × 10−10 cm3 /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 2 factor on the order of k 2 RvdW relative to the range of s-wave processes, due to the threshold law for p waves [97, 106–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 |0ii , single-occupied states |µα ii = c†iµα |0ii , α = 1, 2, and the double-occupied state |µ1 µ2 ii = −|µ2 µ1 ii = c†iµ1 c†iµ2 |0ii . We consider the lattice Hamilto- Minumum error 10 01 11 μ1 μ2 40 K 10 00 3 (b) 10 (a) 2 11 f = 9/2 mf = −9/2 f = 9/2 mf = −3/2 240 250 260 B (Gauss) FIG. 4. Dissipative gate utilizing cold atoms. (a) Two-qubit logical states encoded using a pair of atoms in distinct states occupying four neighboring sites. (b) The upper bound for the total dissipative gate error as a function of magnetic field for 40 K, given the background error rate Γ′ = 10−2 s−1 . The shaded areas describe the weak interactions regime shown in Fig. 3(c). P µ (t)(c†iµ cjµ + nian H = Hµ1 +Hµ2 +V , where Hµ = hiji Jij P µ µ (t) are real tunh.c.) + i ∆i (t)c†iµ ciµ . The quantities Jij µ neling amplitudes, and ∆i (t) are on-site potentials, both of which can depend on the atomic state µ. Two distinct atoms located in P the same trap are subject to elastic interactions V = E i niµ1 niµ2 , where E is the interaction energy and niµ = c†iµ ciµ 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 Ai = Γciµ1 ciµ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 |0iL = |0i|µj i and |1iL = |µj i|0i, 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 S0 and 3 P0 as the two species µ1 and µ2 [24] in order to apply a speciesµ dependent hopping term Jij . 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 P2 instead of 3 P0 , 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 mf . 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. We now give details of the two schemes. Consider the four two-qubit logical states |00iL = |0, µ1 , 0, µ2 i, |01iL = |0, µ1 , µ2 , 0i, |10iL = |µ1 , 0, 0, µ2 i, and |11iL = |µ1 , 0, µ2 , 0i [see Fig. 4(a)], where the comma separates states of individual traps. For the first scheme, an entangling control-Z gate is performed in a single step by applying the hopping H = J(c†2µ2 c3µ2 + h.c.) for state µ2 between traps 2 and 3 for time t = π/J. As a result, the states |00iL and |01iL remain invariant under the evolution, while any transitions involving the state |01iL are blocked by the quantum Zeno effect. In the limit Γ/J → ∞, the overall unitary operation in the logical Hilbert space is described by the control-Z gate in Eq. (31). For the second scheme, the control-Z gate can be applied in three steps: (1) apply the state-independent P hopping H = J i (c†2µi c3µi + h.c.) between traps 2 and 3 for time t = π/(2J); (2) apply a single-qubit phase gate |µ1 i3 → eiπ |µ1 i3 , |µ2 i3 → |µ2 i3 on site 3; (3) repeat the first step. As a result, states |01iL and |10iL remain stationary, state |00iL acquires a total phase 2π (from the phase gate and hopping), and state |11iL acquires phase π. Thus, the second scheme also implements a control-Z gate. 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–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 ǫ = ǫ 0 + Γ′ t +  1  8π 2 1 +O 2 2 , 2 Γt 1 + 4ζ Γ t (33) where ζ = E/Γ is the loss-to-interaction ratio. p The error can be minimized by making the choice t = 2π/ ΓΓ′ (1 + 4ζ 2 ), 12 leading to the expression ǫ = ǫ0 + 4π Γ′ . Γ(1 + 4ζ 2 ) (34) 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 nonGaussian 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- 13 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 nonGaussian 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 nonGaussian 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. ACKNOWLEDGMENTS We thank Bill Fefferman, Mohammad Hafezi, Trey Porto, and Daniel Gottesman for fruitful discussions. This work was supported by DoE ASCR Accelerated Research in Quantum Computing program (award No. DE-SC0020312), U.S. Department of Energy Award No. DE-SC0019449, NSF PFCQC program, DoE ASCR Quantum Testbed Pathfinder program (award No. DE-SC0019040), AFOSR, ARO MURI, ARL CDQI, and NSF PFC at JQI. Appendix A: Easy Class 1 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 t/∆τ U (t) = Y Un , (A1) n=0 where the timestep unitary Un is expressed in terms of a timeordered exponential  Z (n+1)∆τ  Un = T exp −i dtH ′ (t) (A2) n∆τ generated by the stochastic Hamiltonian H ′ (t) in Eq. (13), R(t) H ′ (t) = H(t) + √ . ∆τ (A3) Here H(t) is the original time-dependent Hamiltonian, and P ∗ A†k (t) is the normalized stochastic R(t) = k ξnk Ak (t)+ξnk 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):  1  Un = I − i∆τ 1/2 Rn − ∆τ iHn + Rn2 2  i 3 3/2 − ∆τ PHn Rn − Rn 6 (A4)  i 1 4 1 2 2 − ∆τ Hn − PHn Rn Rn − Rn 2 3 12 + O(∆τ 5/2 ), where we denote the discretized value of an operator On and permutation sum, respectively, as Z Z (n+1)∆τ 1 1 On = dtO(t), dtO(t) ≡ ∆τ n ∆τ n∆τ (A5) X Oσ(1) . . . Oσ(m) . PO1 . . . Om = σ∈Sm The average over the stochastic field can be taken for each timestep independently. Therefore, the effect of the timestep unitary in Eq. (A4) is   1 EUn ρUn† = I + Ln ∆τ + L2n ∆τ 2 ρ 2 (A6) + Dn ρ∆τ 2 + O(∆τ 3 ). 14 In the equation above, Ln is the generator of the original Lindblad equation, lim∆τ →0 Ln = L(n∆τ ), expressed as  X 1 Ln ρ = −i[Hn , ρ] + Akn ρA†kn − {A†kn Akn , ρ} 2 k  X † 1 Akn ρAkn − {Akn A†kn , ρ} , + 2 k (A7) and Dn represents the lowest-order correction occurring due to the timestep being nonzero: 1 X † Dn ρ = Ak′ n A†kn ρ[Ak′ n , Ak ] + A†k′ n Akn ρ[Ak′ n , A†kn ] 4 ′ Appendix B: Easy Class 2 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 Z    X † 0 EUn ρUn = Un ρ + dt Γk (t) Yk (t)ρYk† (t) − ρ Un0† n k 1  = I + Ln + L2n ρ + Dn ρ∆τ 2 + O(∆τ 3 ), 2  (B1) kk  where the target Liouville operator is + Ak′ n A†kn ρ[A†k′ n , Ak ] + Ak′ n Akn ρ[A†k′ n , A†kn ] X  X Ln ρ = −i[Hn , ρ] + Akn ρA†kn − Γkn ρ. † † † † + Akn ρVkn + Vkn ρAkn + Vkn ρAkn + Akn ρVkn k (B2) k + Wn ρ + ρWn† . (A8) Here we use the notation X1 † 1 {Akn , {A†k′ n Ak′ n }} − PA†kn A†k′ n Ak′ n , Vkn = 4 6 k′ i  h iX [Hn , Akn ], A†kn +{Hn , A†kn Akn } Wn = − 6 k 2 1 X  1 X † {Akn , Akn } + PA†kn Akn A†k′ n Ak′ n . − 8 48 ′ k kk The correction now takes the form  X † Akn ρCkn + Ckn ρA†kn Dn ρ = k − 1X 1 Akn Ak′ n ρA†k′ n A†kn − Γ2n ρ, 2 ′ 2 (B3) kk denoting Ckn = Γn A†kn + 2i [A†kn , Hn ] and Γn = R P ∆τ −1 n dt k Γk (t). This expression has the form of (i) Eq. (A10) with operators Dαn being a sum of products of at most four Majorana fermion operators. (A9) The overall expression in Eq. (A8) can be written in a compact form, X (2) (1) Dn ρ = Dαn ρDαn , (A10) α (i) where Dαn = poly(Hn , Akn ) 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 ∆τ , Z t 2  EVst (t2 , t1 ) = T exp dt′ (L(t′ ) + D(t′ )∆τ ) +O(∆τ 2 ) t1 = V(t2 , t1 ) + Z t2 t1 dt′ V(t2 , t′ )D(t′ )V(t′ , t1 )∆τ + O(∆τ 2 ), (A11) 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 Akn and Hn 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(t2 , t′ ) and V(t′ , t1 ) on the right-hand side by their stochastic average and collect all O(∆τ 2 ) terms. 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 La ancilla modes at time t = n∆τ , denoting it as |Ψn i. 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 La − n ancilla modes, i.e. |Ψn i = |φn iL+n ⊗ |0iLa −n . (C1) The evolution is governed by the Hamiltonian H ′ (t) = H(t) ⊗ IA + √  1  K(t) + K † (t) , ∆τ where the stochastic terms are X K(t) = fnk Ak (t)(γ2(L+n) + iγ2(L+n) ) k (C2) (C3) 15 |Ψn+1 i = |Ψn i − i∆τ 1/2 Kn |φn i|1i|0iLa −n−1   1 − ∆τ iHn + Kn† Kn |φn i|0iLa −n 2   i 3/2 1 {Hn , Kn } − Kn Kn† Kn |φn i|1i|0iLa −n−1 − ∆τ 2 6  1 i 2 2 − ∆τ Hn − {Hn , Kn† Kn } 2 3  i 1 − Rn† Hn Kn − (Kn† Kn )2 |φn i|0iLa −n 3 12 + O(∆τ 5/2 ), (C4) where we used the discrete-time operator values Hn and Kn obtained as in Eq. (A4). The interpolated continuous-time evolution for the density matrix of the system can be presented in the form   d 1 ρ= E TrA |Ψn+1 ihΨn+1 | − |Ψn ihΨn | dt ∆τ n=⌊t/∆τ ⌋  1 2 = I + Ln + Ln ρ + Dn ρ∆τ 2 + O(∆τ 3 ), 2 (C5) where ⌊x⌋ is the floor function. The target Liouville operator is X 1 Ln ρ = −i[Hn , ρ] + Akn ρA†kn − {A†kn Akn , ρ} (C6) 2 k and the correction is  X 1 † 1 Dn ρ = Ak Ak′ ρA†k′ Ak − Ak Ak′ ρA†k′ A†k 4 2 kk′  X + Ak ρQk + Q†k ρA†k +M ρ + ρM † , 0.08 EC1: quadratic jumps EC2: unitary jumps EC3: linear jumps 0.06 0.04 0.02 0.00 0.0 0.2 0.4 0.6 0.8 1.0 Time t FIG. 5. Performance of Wick’s theorem. The plot shows the dependence of error in Eq. (D5) as function of time for three processes, involving quadratic jump operators from EC1 (blue, circles), unitary jump operators from EC2 (orange, squares), and linear jump operators from EC3 (green, triangles) as discussed in the text of Appendix D. The systems size is L = 4. The parameters are chosen as βJ = βΓq = βΓs = βJu = βΓu = 1, βΓ′s = 2. The plot demonstrates that jump operators from EC1 and EC2 violate Wick’s theorem, in contrast to the linear jump operators. processes violating Wick’s theorem. First, consider a 1D nearest-neighbor hopping Hamiltonian to be the same for all examples below, H=J L−1 X cn+1 cn + h.c., (D1) n=1 (C7) k where 1 {A† , A† ′ Ak′ } 12 k k  i X † 1 M= Ak HAk − {H, A†k Ak } 6 2 k   X 1 1 − A†k′ Ak′ A†k Ak − A†k′ Ak A†k Ak′ . 12 ′ 2 0.10 Wick\'s theorem error (t) at times n∆τ ≤ t ≤ (n + 1)∆τ , and fnk are independent real Gaussian variables. Then, at the (n + 1)th step, the systemancilla state |Ψn i = Un |Ψn i is Qk = (C8) kk As is the case for EC1 and EC2, the correction is described by (i) Eq. (A10) with operators Dαn being a sum of products of at most eight Majorana fermion operators. Appendix D: Non-Gaussianity of EC1 and EC2 In this Appendix, we show a basic example where jump operators from easy classes (EC) 1 and 2 may be non-Gaussian 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, p p q q F2n = Γq c†n+1 cn F2n+1 = Γq c†n cn+1 , (D2) 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 p
Fn = Γu exp −iJu (c†n cn+1 + h.c.) , (D3) where Ju 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 p p s s F2n = Γs c n , F2n+1 = Γ′s c†n , (D4) where Γs 6= Γ′s are fermion decay/gain rates. 16 For all three cases, we analyze the deviation of the 4-point correlation function from the prediction given by Wick’s theorem, ε(t) = hγ1 γ2 γ3 γ4 it − hγ1 γ2 it hγ3 γ4 it + hγ1 γ3 it hγ2 γ4 it − hγ1 γ4 it hγ2 γ3 it , r (D5) where we define h·it := 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–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. Appendix E: Error analysis In this Appendix, we first derive Eq. (22) and then provide the proof of the lemma in the main text. The error can be formally expressed in terms of evolution superoperators as ǫ= X 1 max |hψr |EVst (t, 0)ρr′ (0) − V(t, 0)ρr′ (0)|ψr i , ′ 2 r r (E1) where V(t2 , t1 ) is the Markovian map generated by Eq. (1) in the main text, Vst (ξ, t2 , t1 ) is a unitary trajectory map depending on either a realization of the discrete stochastic field ξkn (EC1 and EC3) or a random choice of unitaries (EC2). We use the Dyson-like expansion in Eq. (17) and the convexity of the absolute value to upper bound the error as Z t X ∆τ ′ dt ǫ≤ E max hψr |Vst (t, t′ )D(t′ )Vst (t′ , 0)ρr′ |ψr i r′ 2 0 r + O(∆τ 2 ). (E2) Using the fact that Vst is a unitary map, we can rewrite −1 ′ Vst (t′ , 0) = Vst (t , t)Vst (t, 0), (E3) where the inverse of a unitary map is well defined through the inverse unitary transformations. This expression leads directly to Eq. (22), taking into account that X −1 ′ −1 ′ Vst (t , t)D(t)ρVst (t , t) = Dα(1) (t, t′ )ρDα(2) (t, t′ ), α (E4) (i) ity as X X X |hψr |O1 ρO2 |ψr i| = pµ φµr1 φµ∗ r2 hψr |O|ψr1 ihψr2 |O|ψr i (i) where Dα (t, t′ ) = Vst (t, t′ )Dα (t′ ). 1. Proof of the lemma Let us rewrite the left-hand P side of Eq. (25) using the spectral decomposition ρ = µ pµ |φµ ihφµ | and triangle inequal- r ≤ X µ pµ XX r r1 r2 µ,r1 r2 |φµr1 ||φµr2 ||hψr |O|ψr1 i||hψr2 |O|ψr i| ≤ kO1 kmax kO2 kmax X X X µ,r r1 ∈D(k1 ,r) r2 ∈D(k2 ,r) pµ |φµr1 ||φµr2 |, (E5) where we denote kOkmax = maxij |Oij | to be the max-norm of the matrix O, and Dk (r) is a sphere with radius k with respect to Hamming distance. Using the inequality  1 µ 2 |φµr1 ||φµr2 | ≤ |φr1 | + |φµr2 |2 (E6) 2
and the property that the sphere D(k, r) contains Lk ≤ Lk /k! states, we obtain X 1 kO1 kmax kO2 kmax Lk1 +k2 , |hψr |O1 ρO2 |ψr i| ≤ k !k ! 1 2 r (E7) where we use the fact that the density matrix is properly normalized, Tr ρ = 1. Appendix F: Dissipative gates errors In this Appendix, we derive the error of dissipative gates analyzed in Sections III and IV. In the case of imperfect Zeno blockade, the major source of error is associated with leakage to the out-of-logic states. In the scheme proposed in Section III, there are two relevant out-of-logic states into which leakage occurs from the state |00iL = |0101i, namely |0011i and |0000i. The first of these states (|0011i) is accessed via a unitary channel, while the second of these states (|0000i) is accessed via a dissipative channel. The simplest way to describe leakage is to consider unitary evolution of basis states {|0101i, |0011i} and including the second-channel leakage using a non-Hermitian term. The resulting non-Hermitian Hamiltonian is Å ã 0 J HS = . (F1) i J −2Γ The leakage error can be computed as |h00|L S|00iL |2 ≡ 1 − 2ǫ = 1 − 1 4π +O 2 , γ γ (F2) where S = exp(−iπHS /J) and γ = Γ/J. For the scheme involving cold atoms in Section IV, the relevant out-of-logic states are |P 1i = |0, µ2 µ1 , 0, 0i, |P 2i = |0, 0, µ1 µ2 , 0i, |EXi = |0, µ2 , µ1 , 0i, |V Ci = |0, 0, 0, 0i. (F3) 17 The restriction of the effective Hamiltonian to the subspace spanned by the basis {|01i, |P 1i, |P 2i, |EXi} is Ü ê 0 J J 0 J E − 2i Γ 0 J HS′ = , (F4) J 0 E − 2i Γ J 0 J J 0 pair loss. In the strong dissipation limit J ≪ Γ, the leakage error ǫ for the gate can be defined as 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 Ci due to where S ′ = exp(−iHS′ 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. [1] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information: 10th Anniversary Edition (Cambridge University Press, 2010). [2] R. P. Feynman, Simulating physics with computers, Int. J. Theor. Phys 21, 467 (1982). [3] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996). [4] E. Knill, Fermionic linear optics and matchgates, arXiv:quantph/0108033 (2001). [5] L. G. Valiant, Quantum circuits that can be simulated classically in polynomial time, SIAM J. Comput. 31, 1229 (2002). [6] B. M. Terhal and D. P. DiVincenzo, Classical simulation of noninteracting-fermion quantum circuits, Phys. Rev. A 65, 032325 (2002). [7] S. Aaronson and A. Arkhipov, Bosonsampling is far from uniform, arXiv:1309.7460 (2013). [8] M. Wimmer, Algorithm 923: Efficient numerical computation of the pfaffian for dense and banded skew-symmetric matrices, ACM T. Math. Software 38, 30 (2012). [9] N. Schuch and F. Verstraete, Computational complexity of interacting electrons and fundamental limitations of density functional theory, Nat. Phys. 5, 732 (2009). [10] R. Jozsa and M. V. d. Nest, Classical simulation complexity of extended clifford circuits, arXiv:1305.6190 (2013). [11] B. Fefferman and C. Umans, The power of quantum fourier sampling, arXiv:1507.05592 (2015). [12] M. J. Bremner, A. Montanaro, and D. J. Shepherd, Averagecase complexity versus approximate simulation of commuting quantum computations, Phys. Rev. Lett. 117, 080501 (2016). [13] A. Deshpande, B. Fefferman, M. C. Tran, M. Foss-Feig, and A. V. Gorshkov, Dynamical phase transitions in sampling complexity, Phys. Rev. Lett. 121, 030501 (2018). [14] G. Muraleedharan, A. Miyake, and I. H. Deutsch, Quantum computational supremacy in the sampling of bosonic random walkers on a one-dimensional lattice, New J. Phys. 21, 055003 (2018). [15] H.-P. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, 2007). [16] T. Prosen, Spectral theorem for the lindblad equation for quadratic open fermionic systems, J. Stat. Mech.-Theory. E. 2010, P07020 (2010). [17] S. Bravyi and R. König, Classical simulation of dissipative fermionic linear optics, Quantum Information & Computation 12, 925 (2012). [18] C. Gardiner, P. Zoller, and P. Zoller, Quantum noise, Vol. 56 (Springer, 2004). [19] A. Chenu, M. Beau, J. Cao, and A. del Campo, Quantum simulation of generic many-body open system dynamics using classical noise, Phys. Rev. Lett. 118, 140403 (2017). [20] B. Misra and E. G. Sudarshan, The zeno’s paradox in quantum theory, J. Math. Phys. 18, 756 (1977). [21] A. Degasperis, L. Fonda, and G. C. Ghirardi, Does the lifetime of an unstable system depend on the measuring apparatus?, Il Nuovo Cimento A (1965-1970) 21, 471 (1974). [22] J. D. Franson, B. C. Jacobs, and T. B. Pittman, Quantum computing using single photons and the zeno effect, Phys. Rev. A 70, 062302 (2004). [23] Y.-Z. Sun, Y.-P. Huang, and P. Kumar, Photonic nonlinearities via quantum zeno blockade, Phys. Rev. Lett. 110, 223901 (2013). [24] A. J. Daley, M. M. Boyd, J. Ye, and P. Zoller, Quantum computing with alkaline-earth-metal atoms, Phys. Rev. Lett. 101, 170504 (2008). [25] C. W. J. Beenakker, D. P. DiVincenzo, C. Emary, and M. Kindermann, Charge detection enables free-electron quantum computation, Phys. Rev. Lett. 93, 020501 (2004). [26] M. Anderlini, P. J. Lee, B. L. Brown, J. Sebby-Strabley, W. D. Phillips, and J. V. Porto, Controlled exchange interaction between pairs of neutral atoms in an optical lattice, Nature 448, 452 (2007). [27] P. J. Lee, M. Anderlini, B. L. Brown, J. Sebby-Strabley, W. D. Phillips, and J. V. Porto, Sublattice addressing and spindependent motion of atoms in a double-well lattice, Phys. Rev. Lett. 99, 020402 (2007). [28] C. Zhang, S. L. Rolston, and S. Das Sarma, Manipulation of single neutral atoms in optical lattices, Phys. Rev. A 74, 042316 (2006). [29] A. D. Ludlow, N. D. Lemke, J. A. Sherman, C. W. Oates, G. Quéméner, J. von Stecher, and A. M. Rey, Cold-collisionshift cancellation and inelastic scattering in a yb optical lattice clock, Phys. Rev. A 84, 052724 (2011). [30] R. Zhang, Y. Cheng, H. Zhai, and P. Zhang, Orbital feshbach resonance in alkali-earth atoms, Phys. Rev. Lett. 115, 135301 (2015). [31] L. Riegger, N. Darkwah Oppong, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Localized magnetic moments with tunable spin exchange in a gas of ultracold fermions, Phys. Rev. Lett. 120, 143601 (2018). [32] F. Scazza, C. Hofrichter, M. Höfer, P. C. De Groot, I. Bloch, and S. Fölling, Observation of two-orbital spin-exchange interactions with ultracold su(n)-symmetric fermions, Nat. Phys. 10, 779 (2014). [33] R. Le Targat, X. Baillard, M. Fouché, A. Brusch, O. Tcherbakoff, G. D. Rovera, and P. Lemonde, Accurate optical lattice clock with 87 Sr atoms, Phys. Rev. Lett. 97, 130801 (2006). 1 − 2ǫ ≡ |h01|L S ′ |01iL |2 =1−  1  8π 2 1 , + O Γt 1 + 4ζ 2 Γ2 t 2 (F5) 18 [34] G. Cappellini, L. F. Livi, L. Franchi, D. Tusi, D. Benedicto Orenes, M. Inguscio, J. Catani, and L. Fallani, Coherent manipulation of orbital feshbach molecules of two-electron atoms, Phys. Rev. X 9, 011028 (2019). [35] N. Syassen, D. M. Bauer, M. Lettner, T. Volz, D. Dietze, J. J. Garcı́a-Ripoll, J. I. Cirac, G. Rempe, and S. Dürr, Strong dissipation inhibits losses and induces correlations in cold molecular gases, Science 320, 1329 (2008). [36] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Observation of dipolar spin-exchange interactions with lattice-confined polar molecules, Nature 501, 521 EP (2013). [37] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Feshbach resonances in ultracold gases, Rev. Mod. Phys. 82, 1225 (2010). [38] S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of ultracold atomic fermi gases, Rev. Mod. Phys. 80, 1215 (2008). [39] G. Barontini, R. Labouvie, F. Stubenrauch, A. Vogler, V. Guarrera, and H. Ott, Controlling the dynamics of an open manybody quantum system with localized dissipation, Phys. Rev. Lett. 110, 035302 (2013). [40] R. Labouvie, B. Santra, S. Heun, S. Wimberger, and H. Ott, Negative differential conductivity in an interacting quantum gas, Phys. Rev. Lett. 115, 050601 (2015). [41] R. Labouvie, B. Santra, S. Heun, and H. Ott, Bistability in a driven-dissipative superfluid, Phys. Rev. Lett. 116, 235302 (2016). [42] H. P. Lüschen, P. Bordia, S. S. Hodgman, M. Schreiber, S. Sarkar, A. J. Daley, M. H. Fischer, E. Altman, I. Bloch, and U. Schneider, Signatures of many-body localization in a controlled open quantum system, Phys. Rev. X 7, 011034 (2017). [43] Y. S. Patil, S. Chakram, and M. Vengalattore, Measurementinduced localization of an ultracold lattice gas, Phys. Rev. Lett. 115, 140402 (2015). [44] J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Observation of parity-time symmetry breaking transitions in a dissipative floquet system of ultracold atoms, Nat. Commun. 10, 855 (2019). [45] R. Bonifacio, F. Casagrande, and M. Milani, Superradiance and superfluorescence in josephson junction arrays, Phys. Lett. A 101, 427 (1984). [46] M. C. Cassidy, A. Bruno, S. Rubbert, M. Irfan, J. Kammhuber, R. N. Schouten, A. R. Akhmerov, and L. P. Kouwenhoven, Demonstration of an ac josephson junction laser, Science 355, 939 (2017). [47] P. E. Dolgirev, J. Marino, D. Sels, and E. Demler, Non-gaussian correlations imprinted by local dephasing in fermionic wires, Phys. Rev. B 102, 100301 (2020). [48] T. E. Lee and C.-K. Chan, Heralded magnetism in nonhermitian atomic systems, Phys. Rev. X 4, 041001 (2014). [49] T. E. Lee, S. Gopalakrishnan, and M. D. Lukin, Unconventional magnetism via optical pumping of interacting spin systems, Phys. Rev. Lett. 110, 257204 (2013). [50] C. Joshi, F. Nissen, and J. Keeling, Quantum correlations in the one-dimensional driven dissipative xy model, Phys. Rev. A 88, 063835 (2013). [51] S. Diehl, W. Yi, A. J. Daley, and P. Zoller, Dissipation-induced d-wave pairing of fermionic atoms in an optical lattice, Phys. Rev. Lett. 105, 227001 (2010). [52] K. Yamamoto, M. Nakagawa, K. Adachi, K. Takasan, M. Ueda, and N. Kawakami, Theory of non-hermitian fermionic superfluidity with a complex-valued interaction, Phys. Rev. Lett. 123, 123601 (2019). [53] J. A. S. Lourenço, R. L. Eneias, and R. G. Pereira, Kondo effect in a PT -symmetric non-hermitian hamiltonian, Phys. Rev. B 98, 085126 (2018). [54] M. Nakagawa, N. Kawakami, and M. Ueda, Non-hermitian kondo effect in ultracold alkaline-earth atoms, Phys. Rev. Lett. 121, 203001 (2018). [55] M. S. Rudner and L. S. Levitov, Topological transition in a non-hermitian quantum walk, Phys. Rev. Lett. 102, 065703 (2009). [56] T. E. Lee, Anomalous edge state in a non-hermitian lattice, Phys. Rev. Lett. 116, 133903 (2016). [57] D. Leykam, K. Y. Bliokh, C. Huang, Y. D. Chong, and F. Nori, Edge modes, degeneracies, and topological numbers in nonhermitian systems, Phys. Rev. Lett. 118, 040401 (2017). [58] S. Yao and Z. Wang, Edge states and topological invariants of non-hermitian systems, Phys. Rev. Lett. 121, 086803 (2018). [59] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Topological phases of non-hermitian systems, Phys. Rev. X 8, 031079 (2018). [60] H. Shen, B. Zhen, and L. Fu, Topological band theory for nonhermitian hamiltonians, Phys. Rev. Lett. 120, 146402 (2018). [61] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-hermitian physics, Phys. Rev. X 9, 041015 (2019). [62] T. Yoshida, K. Kudo, and Y. Hatsugai, Non-hermitian fractional quantum hall states, Scientific Reports 9, 16895 (2019). [63] K. B. Efetov, Directed quantum chaos, Phys. Rev. Lett. 79, 491 (1997). [64] I. Y. Goldsheid and B. A. Khoruzhenko, Distribution of eigenvalues in non-hermitian anderson models, Phys. Rev. Lett. 80, 2897 (1998). [65] S. Longhi, D. Gatti, and G. Della Valle, Non-hermitian transparency and one-way transport in low-dimensional lattices by an imaginary gauge field, Phys. Rev. B 92, 094204 (2015). [66] To be precise, for the first task (estimating few-body observables), the error is measured by the maximum difference in the estimated and ideal expectation values of a unit-norm observable O. For the second task (sampling), the error is measured through the maximum variation distance between the ideal and the sampled probability distributions. [67] B. Horstmann, J. I. Cirac, and G. Giedke, Noise-driven dynamics and phase transitions in fermionic systems, Phys. Rev. A 87, 012108 (2013). [68] M. Žnidarič, Exact solution for a diffusive nonequilibrium steady state of an open quantum chain, J. Stat. Mech.–Theory E. 2010, L05002 (2010). [69] V. Eisler, Crossover between ballistic and diffusive transport: the quantum exclusion process, J. Stat. Mech.–Theory E. 2011, P06007 (2011). [70] S. Bravyi, B. M. Terhal, and B. Leemhuis, Majorana fermion codes, New J. Phys. 12, 083039 (2010). [71] G. Lindblad, On the generators of quantum dynamical semigroups, Comm. Math. Phys. 48, 119 (1976). [72] E.-M. Laine, J. Piilo, and H.-P. Breuer, Measure for the nonmarkovianity of quantum processes, Phys. Rev. A 81, 062115 (2010). [73] J. Colpa, Diagonalisation of the quadratic fermion hamiltonian with a linear part, J. Phys. A-Math. Gen. 12, 469 (1979). [74] M. Oszmaniec, J. Gutt, and M. Kuś, Classical simulation of fermionic linear optics augmented with noisy ancillas, Phys. Rev. A 90, 020302 (2014). [75] F. de Melo, P. Ćwikliński, and B. M. Terhal, The power of noisy fermionic quantum computation, New J. of Phys. 15, 013015 (2013). [76] J. F. Corney and P. D. Drummond, Gaussian phase-space representations for fermions, Phys. Rev. B 73, 125112 (2006). 19 [77] M. S. Underwood and D. L. Feder, Bose-Hubbard model for universal quantum walk-based computation, Phys. Rev. A 85, 052314 (2012). [78] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, in Proceedings of the forty-third annual ACM symposium on Theory of computing (ACM, 2011) pp. 333– 342. [79] A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani, On the complexity and verification of quantum random circuit sampling, Nat. Phys. 15, 159 (2019). [80] P. Aliferis and B. M. Terhal, Fault-tolerant quantum computation for local leakage faults, Quantum Information & Computation 7, 139 (2007). [81] D. Aharonov and M. Ben-Or, Fault-tolerant quantum computation with constant error, in Proceedings of the Twenty-Ninth Annual ACM Symposium on Theory of Computing - STOC ’97 (ACM Press, El Paso, Texas, United States, 1997) pp. 176– 188. [82] E. Knill, R. Laflamme, and W. H. Zurek, Resilient quantum computation, Science 279, 342 (1998). [83] B. M. Terhal, Quantum error correction for quantum memories, Rev. Mod. Phys. 87, 307 (2015). [84] E. Arimondo, M. Ingucio, and P. Violino, Experimental determination of the hyperfine structure in the alkali atoms, Rev. Mod. Phys. 49, 31 (1977). [85] S. Falke, H. Knöckel, J. Friebe, M. Riedmann, E. Tiemann, and C. Lisdat, Potassium ground-state scattering parameters and born-oppenheimer potentials from molecular spectroscopy, Phys. Rev. A 78, 012503 (2008). [86] H. T. C. Stoof, J. M. V. A. Koelman, and B. J. Verhaar, Spin-exchange and dipole relaxation rates in atomic hydrogen: rigorous and simplified calculations, Phys. Rev. B 38, 4688 (1988). [87] T. Loftus, C. A. Regal, C. Ticknor, J. L. Bohn, and D. S. Jin, Resonant control of elastic collisions in an optically trapped fermi gas of atoms, Phys. Rev. Lett. 88, 173201 (2002). [88] C. A. Regal and D. S. Jin, Measurement of positive and negative scattering lengths in a fermi gas of atoms, Phys. Rev. Lett. 90, 230404 (2003). [89] C. D’Errico, M. Zaccanti, M. Fattori, G. Roati, M. Inguscio, G. Modugno, and A. Simoni, Feshbach resonances in ultracold 39 k, New J. Phys. 9, 223 (2007). [90] A. Ludewig, Feshbach Resonances in 40 K, Ph.D. thesis, University of Amsterdam (2012). [91] J. S. Krauser, J. Heinze, S. Götze, M. Langbecker, N. Fläschner, L. Cook, T. M. Hanna, E. Tiesinga, K. Sengstock, and C. Becker, Investigation of feshbach resonances in ultracold 40 K spin mixtures, Phys. Rev. A 95, 042701 (2017). [92] L. Tanzi, C. R. Cabrera, J. Sanz, P. Cheiney, M. Tomza, and L. Tarruell, Feshbach resonances in potassium bose-bose mixtures, Phys. Rev. A 98, 062712 (2018). [93] R. Chapurin, X. Xie, M. J. Van de Graaff, J. S. Popowski, J. P. D’Incao, P. S. Julienne, J. Ye, and E. A. Cornell, Precision test of the limits to universality in few-body physics, Phys. Rev. Lett. 123, 233402 (2019). [94] N. Balakrishnan, V. Kharchenko, R. Forrey, and A. Dalgarno, Complex scattering lengths in multi-channel atomñmolecule collisions, Chem. Phys. Lett. 280, 5 (1997). [95] J. L. Bohn and P. S. Julienne, Prospects for influencing scattering lengths with far-off-resonant light, Phys. Rev. A 56, 1486 (1997). [96] J. M. Hutson, Feshbach resonances in ultracold atomic and molecular collisions: threshold behaviour and suppression of poles in scattering lengths, New J. Phys. 9, 152 (2007). [97] Z. Idziaszek and P. S. Julienne, Universal Rate Constants for Reactive Collisions of Ultracold Molecules, Phys. Rev. Lett. 104, 113202 (2010). [98] M. D. Frye, P. S. Julienne, and J. M. Hutson, Cold atomic and molecular collisions: approaching the universal loss regime, New J. Phys. 17, 045019 (2015). [99] T. L. Nicholson, S. Blatt, B. J. Bloom, J. R. Williams, J. W. Thomsen, J. Ye, and P. S. Julienne, Optical feshbach resonances: Field-dressed theory and comparison with experiments, Phys. Rev. A 92, 022709 (2015). [100] K. Huang and T. D. Lee, Quantum mechanical many-body problem with hard sphere interaction, Phys. Rev. 105, 767 (1957). [101] I. Bloch, Ultracold quantum gases in optical lattices, Nat. Phys. 1, 23 (2005). [102] T. Busch, B.-G. Englert, K. Rzažewski, and M. Wilkens, Two cold atoms in a harmonic trap, Found. of Phys. 28, 549 (1998). [103] D. Blume and C. H. Greene, Fermi pseudopotential approximation: Two particles under external confinement, Phys. Rev. A 65, 043613 (2002). [104] E. L. Bolda, E. Tiesinga, and P. S. Julienne, Effectivescattering-length model of ultracold atomic collisions and feshbach resonances in tight harmonic traps, Phys. Rev. A 66, 013403 (2002). [105] E. Tiesinga, C. J. Williams, F. H. Mies, and P. S. Julienne, Interacting atoms under strong quantum confinement, Phys. Rev. A 61, 063416 (2000). [106] B. DeMarco, J. L. Bohn, J. P. Burke, Jr., M. Holland, and D. S. Jin, Measurement of p-wave threshold law using evaporatively cooled fermionic atoms, Phys. Rev. Lett. 82, 4208 (1999). [107] C. A. Regal, C. Ticknor, J. L. Bohn, and D. S. Jin, Tuning p-wave interactions in an ultracold fermi gas of atoms, Phys. Rev. Lett. 90, 053201 (2003). [108] Z. Idziaszek and T. Calarco, Pseudopotential method for higher partial wave scattering, Phys. Rev. Lett. 96, 013201 (2006). [109] D. D. Yavuz, P. B. Kulatunga, E. Urban, T. A. Johnson, N. Proite, T. Henage, T. G. Walker, and M. Saffman, Fast ground state manipulation of neutral atoms in microscopic optical traps, Phys. Rev. Lett. 96, 063001 (2006). [110] A. V. Gorshkov, L. Jiang, M. Greiner, P. Zoller, and M. D. Lukin, Coherent quantum optical control with subwavelength resolution, Phys. Rev. Lett. 100, 093005 (2008). [111] K. D. Stokes, C. Schnurr, J. R. Gardner, M. Marable, G. R. Welch, and J. E. Thomas, Precision position measurement of moving atoms using optical fields, Phys. Rev. Lett. 67, 1997 (1991). [112] D. Schrader, I. Dotsenko, M. Khudaverdyan, Y. Miroshnychenko, A. Rauschenbeutel, and D. Meschede, Neutral atom quantum register, Phys. Rev. Lett. 93, 150501 (2004). [113] J. R. Gardner, M. L. Marable, G. R. Welch, and J. E. Thomas, Suboptical wavelength position measurement of moving atoms using optical fields, Phys. Rev. Lett. 70, 3404 (1993). [114] D. D. Yavuz and N. A. Proite, Nanoscale resolution fluorescence microscopy using electromagnetically induced transparency, Phys. Rev. A 76, 041802 (2007). [115] T. Tiecke, Feshbach resonances in ultracold mixtures of the fermionic quantum gases, Ph.D. thesis (2009). [116] M. Hebenstreit, R. Jozsa, B. Kraus, S. Strelchuk, and M. Yoganathan, All pure fermionic non-Gaussian states are magic states for matchgate computations, Phys. Rev. Lett. 123, 080503 (2019). [117] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, 20 [118] [119] [120] [121] [122] [123] [124] B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing Rényi entanglement entropy via randomized measurements, Science 364, 260 (2019). V. Alba and F. Carollo, Spreading of correlations in Markovian open quantum systems, arXiv:2002.09527 (2020). N. Maskara, A. Deshpande, M. C. Tran, A. Ehrenberg, B. Fefferman, and A. V. Gorshkov, Complexity phase diagram for interacting and long-range bosonic Hamiltonians, arXiv:1906.04178 (2019). D. Poulin, Lieb-Robinson Bound and Locality for General Markovian Quantum Dynamics, Phys. Rev. Lett. 104, 190401 (2010). T. Barthel and M. Kliesch, Quasilocality and Efficient Simulation of Markovian Quantum Dynamics, Phys. Rev. Lett. 108, 230504 (2012). R. Movassagh, Cayley path and quantum computational supremacy: A proof of average-case #P -hardness of Random Circuit Sampling with quantified robustness, arXiv:1909.06210 (2019). A. Bouland, B. Fefferman, Z. Landau, and Y. Liu, Noise and the frontier of quantum supremacy, arXiv:2102.01738 (2021). D. J. Ahmed-Braun, K. G. Jackson, S. Smale, C. J. Dale, B. A. Olsen, S. J. Kokkelmans, P. S. Julienne, and J. H. Thywissen, Probing open-and closed-channel p-wave resonances, arXiv:2101.02700 (2021). [125] We note that individual species do show a wide range of variation due to their specific molecular physics. An example of such variation is provided by 87 Rb and 85 Rb. The former has no broad s-wave Feshbach resonances whereas the latter has a very prominent one. The difference is due to the nature of the singlet and triplet scattering lengths for the two different mass isotopes: nearly identical values suppress the Feshbach resonance for 87 Rb, but quite different values permit a strong Feshbach resonance for 85 Rb. [126] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Gaussian Boson Sampling, Phys. Rev. Lett. 119, 170501 (2017). [127] R. Kruse, C. S. Hamilton, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, A detailed study of Gaussian Boson Sampling, Phys. Rev. A 100, 032326 (2019). [128] L. Chakhmakhchyan and N. J. Cerf, Boson sampling with gaussian measurements, Phys. Rev. A 96, 032326 (2017). [129] U. Chabaud, T. Douce, D. Markham, P. van Loock, E. Kashefi, and G. Ferrini, Continuous-variable sampling from photonadded or photon-subtracted squeezed states, Physical Review A 96, 062307 (2017). [130] A. Deshpande, A. Mehta, T. Vincent, N. Quesada, M. Hinsche, M. Ioannou, L. Madsen, J. Lavoie, H. Qi, J. Eisert, D. Hangleiter, B. Fefferman, and I. Dhand, Quantum Computational Supremacy via High-Dimensional Gaussian Boson Sampling, arXiv:2102.12474 (2021).