Quantum Simulation of the Sachdev-Ye-Kitaev Model by Asymmetric Qubitization

We show that one can quantum simulate the dynamics of a Sachdev-Ye-Kitaev model with $N$ Majorana modes for time $t$ to precision $\epsilon$ with gate complexity $O(N^{7/2} t + N^{5/2} t \,{\rm polylog}(N/ \epsilon))$. In addition to scaling sublinearly in the number of Hamiltonian terms, this gate complexity represents an exponential improvement in $1/\epsilon$ and large polynomial improvement in $N$ and $t$ over prior state-of-the-art algorithms which scale as $O(N^{10} t^2 / \epsilon)$. Our approach involves a variant of the qubitization technique in which we encode the Hamiltonian $H$ as an asymmetric projection of a signal oracle $U$ onto two different signal states prepared by state oracles, $A\left\vert{0}\right\rangle \mapsto \left\vert{A}\right\rangle$ and $B \left\vert{0}\right\rangle \mapsto \left\vert{B}\right\rangle$, such that $H = \left\langle{B}\right\vert U\left\vert{A}\right\rangle$. Our strategy for applying this method to the Sachdev-Ye-Kitaev model involves realizing $B$ using only Hadamard gates and realizing $A$ as a random quantum circuit.


I. INTRODUCTION
The AdS/CFT correspondence is a conjectured relationship between the quantum physics of correlated many-body systems and the classical physics of gravity in one higher dimension [1]. Holographic dualities such as AdS/CFT have become increasingly important tools for studying quantum gravity; however, it has generally been difficult to find simple models that capture exotic features such as black holes [2]. Introduced by Kitaev [3] based on a variant of an earlier model by Sachdev and Ye [4], the Sachdev-Ye-Kitaev (SYK) model has been widely studied in recent years as an example of a simple quantum many-body system which may have an interesting holographic dual. The SYK model can be expressed as H = 1 4 · 4! N −1 p,q,r,s=0 J pqrs γ p γ q γ r γ s (1) where the J pqrs are real-valued scalars drawn randomly from a normal distribution with variance σ 2 = 3!J 2 /N 3 and the γ p are Majorana fermion mode operators. A feature of the SYK model is that for large N and strong coupling J it is possible to sum all of the Feynman diagrams to obtain (among other properties) outof-order-time correlation functions [3]. Such methods reveal that the SYK model is maximally chaotic (having the same Lyapunov exponent as black holes in Einstein gravity) [3]. This property, in conjunction with the emergence of an approximate conformal symmetry at low temperatures, suggests that the holographic dual of the SYK model is a theory of Einstein gravity [5].
However, questions remain about the SYK model which have not yet been solved by analytical methods. Since quantum computers are widely conjectured to be capable of modeling otherwise intractable quantum systems [6], they may also prove useful tools for studying the * Corresponding author: babbush@google.com quantum many-body side of such holographic duals. For example, it has been analytically challenging to obtain the density of states associated with interesting states (e.g. the thermal state) of the SYK model. A quantum computer could efficiently sample such distributions by performing quantum phase estimation [7].
Earlier work [8] has proposed a scheme for quantum simulating SYK model dynamics using a Lie-Trotter based algorithm [9] with gate complexity O(N 10 t 2 /ǫ), where t is time and ǫ is target precision. Yet such scaling suggests that interesting (e.g. N > 100) SYK model quantum simulations would remain intractable for even a fault-tolerant quantum computer. In this paper we describe an algorithm to quantum simulate the SYK model with gate complexity O(N 7/2 t+N 5/2 t polylog(N/ǫ)). We compile the bottleneck components of our approach to Clifford + T gates and find that interesting simulations are possible with fewer than ten million Jt T gates. This relatively low T complexity makes the application attractive because T gates consume many logical qubits and take much longer to apply than any other operation within all known two-dimensional error-correcting codes with low thresholds (e.g. the surface code) [10].

II. ASYMMETRIC QUBITIZATION
Our method is to use a linear combinations of unitaries (LCU) approach combined with an "asymmetric" extension of the qubitization simulation framework [11]. The original formulation of qubitization (we refer to as "symmetric" qubitization) requires a state preparation with amplitudes proportional to the square roots of the coefficients J pqrs . Instead, our asymmetric form of qubitization permits a state preparation with amplitudes linearly proportional to the normally distributed J pqrs . Such a state can be prepared with low gate complexity by using a random quantum circuit.
The standard LCU query model [12] for Hamiltonian simulation uses two unspecified unitary operators, G and U . The operator G initializes a state, G|0 → |G in an ancilla register. Then U performs one of a set of unitaries on the target system controlled on the ancilla register. Finally, measurement of this ancilla register in the state |G gives the linear combination of unitaries applied to the target system. Conventional qubitization [11] requires that which then enables Hamiltonian simulation to be achieved via quantum signal processing [13]. In contrast to the LCU model, no postselection or amplitude amplification is required for quantum signal processing. In Eq. (2), U acts on the system and ancilla registers whereas G is a state only on the ancilla register. Therefore G|U |G gives an operator acting on the target system. We have included a scaling factor λ, since usually an operator proportional to the desired Hamiltonian is given. For our application, U is a controlled unitary operation between the ancilla and target systems, but the qubitization framework allows U to be more general. For our SYK model simulation we would like to encode the Hamiltonian using two different state preparation oracles; i.e., in contrast to Eq. (2), we would like to have state preparations A|0 → |A and B|0 → |B together with a unitary V such that B|V |A = H/λ. This is what we mean by "asymmetric" qubitization. In order to use qubitization with these asymmetric states, we can bundle the operations A and B together with V , so 0|B † V A|0 = H/λ. The difficulty now is that B † V A will not be self-inverse, but qubitization works best with self-inverse operators.
We will take V to be self-inverse. That means we can add an ancilla qubit and construct a self-inverse U as (in a block matrix representation) That is, the ancilla controls which state preparation is performed, then we perform V and a not on the ancilla qubit, then there is an inverse controlled state preparation. The operation V is the most costly and is only performed once. The complexity of A and B is logarithmic compared to the complexity of V , so the overall complexity is not significantly increased. Reference [11] proposes a method to construct self-inverse operators, but it doubles the complexity. An alternative realization of asymmetric qubitization is discussed in the Appendix A.
In our case, we take the Hamiltonian to be a sum of self-inverse terms and V will be a controlled H ℓ , The state preparation A will be a random orthogonal operation producing Gaussian distributed amplitudes, and B will just give an equal superposition.
To perform Hamiltonian simulation with asymmetric qubitization, we can use the same method as for symmetric qubitization, with the new operation U constructed above, and G just being a Hadamard on the new ancilla qubit. For quantum signal processing [13], one combines the operation U with a reflection R = 2 |G G| − 1 1 as W = RU . Applications of W and W † are controlled by a qubit, and these controlled-W operations are interspersed with rotations of the control qubit. Using M controlled operations can yield an overall polynomial in W of order M [13]. This polynomial is chosen to give the Hamiltonian evolution e −iHt . The order of the polynomial required scales as λ, which therefore governs the complexity. In the following we will examine the additional cost in λ due to the asymmetric qubitization, explain how to perform the state preparations and controlled operation V , then determine the overall cost.
Suppose that operators A and B prepare states |A and |B as Then, taking the register |ψ to be the target subsystem on which the operators H ℓ act, Thus, using Eq. (4) we obtain, so that for all ℓ, we require λα ℓ β * ℓ = w ℓ . Therefore, taking the absolute value and summing gives which comes from the Cauchy-Schwarz inequality applied to normalized states |A and |B . Hence, λ takes its minimum value of ℓ |w ℓ | in the case that |α ℓ | = |β ℓ |. This means that the most efficient case is that where the preparation and inverse preparation are symmetric. For some classes of Hamiltonians, the preparation can be made heavily asymmetric without much cost. In particular, consider the asymmetric case where Then using Eq. (10) gives Using · to indicate the mean over ℓ, this means that the overhead from using the asymmetric |A and |B over the symmetric case is |w ℓ | 2 / |w ℓ | . In the case where the w ℓ are drawn from a normal distribution (e.g. the SYK model), the additional complexity is only π/2 ≈ 1.25.

A. The State Preparation Circuits
Our strategy for implementing A and B oracles to simulate the SYK model is straightforward; B will consist of Hadamard gates which initialize the symmetric superposition state and A will be a random quantum circuit with orthogonal rotations. Let us assume for simplicity that the number of terms in our SYK model Hamiltonian (L = N 4 ) is a binary power. Then, one can initialize the state |B such that β ℓ = 1/ √ L by implementing the circuit B as a sequence of log L Hadamard gates.
Consider the state |A output by a random quantum circuit A with orthogonal rotations. We will use orthogonal rotations in order to ensure that the amplitudes remain real. Using the definition in Eq. (6), known properties of orthogonal random quantum evolutions [14] hold that the α ℓ are Gaussian distributed with zero mean and variance equal to the Hilbert space dimension. This approach to simulating the SYK model reveals a possibly interesting connection between chaos in random quantum circuits [15] and chaos in AdS 2 holography [16].
The asymmetric state preparation requires that the α ℓ are proportional to the desired weightings for the terms in the Hamiltonian. Because they have a normal distribution, they correctly generate the weights w ℓ . There is a difference in the variance, but the variance only affects an overall scaling of the values of the α ℓ , which is what we expect because the α ℓ correspond to the amplitudes of a normalized state. The scaling factor is taken into account in λ, giving the complexity of the calculation.
The value of λ for this asymmetric state preparation is given by Eq. (13), which only depends on w ℓ and the dimension. It does not explicitly depend on α ℓ , because that proportionality is already taken into account by the fact that the α ℓ are normalized. For the SYK model, where we have changed the index to pqrs to match the notation used for J. Then, the mean square value is There is an approximate equality here, because the mean of w pqrs indicates the mean summing over pqrs, not the expectation value of the probability distribution according to which J pqrs are chosen. Using Eq. (13) then gives An outstanding question is how large the orthogonal random quantum circuits should be in order to achieve sufficient convergence to the Gaussian distribution in the coefficients. There have been many theoretical results on related questions such as the convergence of random quantum circuits to t-designs of the Haar measure [17][18][19][20]. For a one dimensional random quantum circuit on log L qubits, the circuit approaches an ǫ-approximate 2design in gate complexity O(log 2 L+log(L/ǫ)) [17]. However, for circuits in higher dimension the gate complexity to achieve similar states is closer to O(log(L/ǫ)) [20].
While these works typically do not focus on the convergence of amplitudes to a Gaussian distribution, this topic was recently studied numerically in [15]. There, authors found rapid convergence of the probabilities (squared amplitudes) to the Porter-Thomas distribution [14], which corresponds to convergence of real and imaginary components of the amplitudes to a Gaussian distribution. Thus, to avoid an in depth discussion of the requisite circuit size, we will conservatively assume that to achieve amplitudes that are within ǫ distance of Gaussian distributed amplitudes, it suffices to use circuits of size O(polylog(L/ǫ)) = O(polylog(N/ǫ)).

B. The Hamiltonian Application Circuit
The circuit referred to in Eq. (5) as U should act as U |p |q |r |s |ψ → |p |q |r |s γ p γ q γ r γ s |ψ .
To implement this on a quantum computer, the Majorana operators are represented by strings of Pauli operators according to the Jordan-Wigner representation. We therefore need to apply the following transformation Circuits for exactly this transformation were introduced in Section 3B of [21], reproduced here in Figure 1. The T complexity of this circuit implementing the transformation of Eq. (18) is exactly 4N − 4 and the Clifford complexity is also O(N ). Clearly, four applications of this primitive are sufficient to implement U ; thus, our total implementation has gate complexity 16N − 16 and uses only log N additional ancillae. where N is the size of the target register. The circuit acts on N + 2 log N + 1 qubits, including log N ancillae. Note that this diagram uses notation for setting an ancilla to the logical and of two other qubits, described in the figure to the right. Our implementation of U consists of four of these circuits, which would require a total of 16N − 16 T gates and log N ancillae. This circuit is explained in detail in Section 3B of [21]. (Right): circuit for computing and uncomputing and operations [22], defined in terms of Toffoli gates and Clifford+T gates.

IV. HAMILTONIAN EVOLUTION
Now we give more details on how to perform the Hamiltonian evolution and determine its cost. A self-inverse operation U is equivalent to a reflection. Then R provides a second reflection, and a product of reflections has a spectrum given by Theorem 1 of [23] (a product of reflections gives a rotation in the same way as for Grover's algorithm [24]). The eigenvalues then correspond to the angle of rotation for each eigenstate. Applying that theorem shows that W = RU has eigenvalues of e ±i arccos(h/λ) for eigenvalue h of H [11]. One can show that (see the proof of Lemma 6 of [25] or Lemma 16 of [26]) where T n (·) is the n th Chebyshev polynomial of the first kind. This result originates from the relation T n (cos θ) = cos nθ. Here θ is equivalent to h/λ, and applying W n gives a rotation by nθ. Using the Jacobi-Anger expansion where J n (·) is the n th Bessel function of the first kind, and we have used cos(±n arccos(H/λ)) = T n (H/λ). Therefore a polynomial in W can be used to approximate e −iHt , and this polynomial can be generated via quantum signal processing [13]. The polynomial order needed to approximate e −iHt to within error ǫ is For large λt, it can be shown that the order required is approximately (see Appendix B) The complexity of simulating e −iHt via quantum signal processing is twice this, in terms of the number of applications of R and U . That can be seen from Theorem 1 of [11], where the order is N/2, and the total number of controlled operations is N .

V. CONCLUSION
Using quantum signal processing, our simulation requires a number of applications of A, B, and U given by Eq. (21). The total cost of simulation is thus, where C A , C B and C U are the cost of implementing A, B and U , respectively. In terms of gate complexities of the explicit implementations advocated for in this work, C A = O(polylog(N/ǫ)), C B = O(log N ), and C U = O(N ). We have found that λ = O(N 5/2 ). Thus, the total asymptotic complexity becomes O N 7/2 t + N 5/2 t polylog (N/ǫ) .
Since we have compiled all bottleneck components down to Clifford + T gates, we are also able to report the leading order scaling of the T count of the algorithm as 2 √ 6 N 7/2 Jt. Assuming reasonable precision goals, for N = 100 the leading order T count is less than ten million Jt and for N = 200, the leading order T count is less than one hundred million Jt. This should be compared to the roughly 10 12 T gates required to simulate the active site of FeMoco with 108 qubits (a molecule relevant to Nitrogen fixation) [27][28][29], 10 9 T gates required to simulate 100-200 qubit interesting problems in solid-state electronic structure [21,30] or the roughly 10 9 T gates required to simulate a 100 qubit onedimensional Heisenberg model for classically intractable durations [31]. While the particular value of Jt would depend on the application, this analysis reveals that simulation of the SYK model is among the most viable applications of the first surface code quantum computers.
holds for all steps in oblivious amplitude amplification, regardless of symmetry. Define the step of oblivious amplitude amplification where R 0 ≡ 2 |0 0| − 1 1 is a reflection on the control register. This operation can very easily be made controlled, just by making the reflection R 0 controlled. Then, Therefore, it is possible to generate both the even and odd terms required in Eq. (B1) by using steps of U or steps of U followed by B † U A. In order to generate the polynomial required by quantum signal processing, an ancilla qubit can be used, which controls whether the odd or even terms in the polynomial in Eq. (B1) are generated. This qubit would primarily control the qubit rotations used in quantum signal processing. Ultimately this qubit would also control the operations B † U A, applying them to produce the odd terms.
In this approach, combining the odd and even terms via a linear combination of unitaries gives a success amplitude that is not unity. Amplitude amplification requires repeating the sequence of operations at least three times, giving an additional multiplicative overhead. In the scheme proposed in the main text, U is applied only once. There is a factor of two for the preparations, but those have complexity that is logarithmic compared to the complexity of U . Thus, the scheme in the main text is the most gate-efficient of the alternatives considered.
For completeness, we derive the equations in Eq. (A6) and Eq. (A5). Let P = |G G|, and R = 2P − 1 1, and we require that G|U |G = H. We claim that Let us consider Eq. (A5) first. For the case m = 0 the right hand side is Therefore this expression is clearly true for m = 0. Then, if it is true for m − 1, we get Therefore, We will put τ = λt for simplicity. We are interested in the regime where the cutoff K is large, and slightly larger than τ . This is known as the transition region, and the asymptotic form is given in Eq.