Significant progress has been made in the last few years toward building noisy intermediate scale quantum (NISQ) computers [1]. These quantum devices have fewer qubits ranging from dozens to a few hundred and despite their limitations, they have been used to tackle certain classically challenging problems and emulate numerous many-body phases of matter [2,3]. However, more conventional quantum algorithms designed for ideal quantum computers do not in general perform meaningful computations when run on a NISQ device, suggesting that general usage will require quantum error correction.
Until such a time as quantum error correction is achieved, it remains important to understand the role of noise in creating issues for known quantum algorithms. One such algorithm is Grover's algorithm for unstructured database search [4]. It has been the subject of intense study and has, in fact, been implemented on various NISQ devices [5][6][7]. The effects of noise on Grover's search have been studied using different models and efforts have been made towards putting limits on the performance of the algorithm. For example, Long et al. [8] determined that the size of the database is constrained by O(1/δ 2 ) where δ is an error parameter representing gate imperfection. In the presence of random phase error in the oracle gate, Shenvi et al. [9] have shown that in order to maintain the probability of success, the error must decrease as k 1/4 when the database size is increased by a factor of k. In the case of depolarizing noise channel, Salas [10] found that the noise threshold behaves as N -1.1 , where N is the size of the database. Norman and Ataba [11] determined that the allowed noise rate scales as N -2/3 if Gaussian noise is added at each step of the algorithm. Finally, Shapira et al. [12] showed that for Grover's algorithm with L qubits, with unitary noise characterized by standard deviation , must satisfy < O(L -1/2 N -1/4 ) where N = 2 L , in order for the algorithm to maintain a significant efficiency. It was also shown that any NISQ algorithm will fail to achieve the quadratic speedup over classical algorithms in the unstructured database search problem [13].
These cases primarily studied unphysical noise models on the long-range many-qubit Grover gates and consistently found that exponentially small error rates prevent the algorithm from succeeding. These error models make generalizing the results rather challenging, leaving open the question of whether other physically realistic error models might have improved performance bounds. Furthermore, the standard Grover algorithm involves a call to an oracle which is a nonk-local operation that is unphysical in most realizations. In practice, such non-k-local gates are implemented by decomposing it in terms of k-local gates (k = 1, 2 in most cases that act on one qubit, or two qubits that are not necessarily adjacent to each other). We will consider the impact of errors on a circuit in which the oracle is unraveled into single-qubit gates and 2-local unitary gates (cf. Fig. 2), which comes at the cost of larger depth.
In order to get a better analytical handle on the error tolerance of Grover's algorithm, we consider here the effects of systematic gate noise, i.e., undesired noise in the quantum gates that is repeated systematically throughout the course of the algorithm. Systematic noise naturally occurs in most experiments due to calibration issues. For our purposes, it imbues the noisy algorithm with a time-periodic (Floquet) structure, such that we can find the exact scaling of target FIG. 1. Circuit for the four-qubit Grover operator G 0 . X, Z, and H are the Pauli X , Pauli Z, and the Hadamard gate respectively. state dynamics, noise thresholds, and other physically relevant quantities by employing both analytical and numerical approaches. We show that these signatures can be captured by a simple random-matrix Hamiltonian whose perturbative effect on the Floquet spectrum describes gap-closing transitions that determine the response to noise. We find that for a database of size N = 2 L , the algorithm loses its computational power at a noise strength that scales as O(N -1/2 ). Additionally, an ergodicity transition occurs at a noise strength that scales as O(L -1 ) with the system size. We show that these analytical predictions align perfectly with the exact numerical solution and comment on the relevance to using Grover's algorithmas well as related amplitude amplification techniques-in the NISQ era.
We start by reviewing the basic formulation of Grover's algorithm. Suppose we want to search for an item w in an unstructured database with N = 2 L items. For simplicity, let us assume that the target item w is unique. The search problem can be expressed by using an indicator function f (u) that takes integer inputs u, with 0 u N -1, such that
1 if u is the solution, i.e., u = w 0 otherwise.
In the case of Grover's algorithm, this indicator function is implemented using a unitary operator U w , known as the oracle, whose action on the computational basis states is
The complexity of an algorithm can be measured by the number of times the function f (u) is evaluated. Classical algorithms require ∼N/2 "time steps"-calls to a circuit that evaluates the function f (u)-on average to find the target state. Grover's algorithm requires O( √ N ) evaluations of the oracle U w to obtain the target state, and thus provides FIG. 2. Example decomposition of the four-qubit-controlled X gate in Fig. 1 in linear depth using the algorithm described in [21]. The gates with the same control qubit are shown together. quadratic speed up of query complexity over the classical counterparts.
Instead of the elements of the database directly, Grover's algorithm focuses on the indices of the elements, such that the database can be represented using L qubits as {|0 , |1 , . . . , |2 L -1 }. The algorithm starts with the uniform superposition state |x = 1 √ N N-1 j=0 | j , and then amplifies the amplitude of the target state |w using the action of the two unitary operators. First, it applies the oracle, which flips the sign of all the states except for the state |w ,
Second is the "diffusion" operator U x , which flips the amplitude of all the states about their mean,
The combined action of these two reflections is written in terms of the Grover operator as G w = U x U w , which gradually rotates the initial state |x towards the target state |w . The action of G w on the state of the system can be interpreted geometrically by writing the state of the system at t = 0 as
where
and sin θ = 1/ √ N. Then the state of the system at time t can be written as
In other words, G w rotates the wave function in the twodimensional subspace spanned by |x w and |w by angle 2θ. From Eq. ( 6), we get the probability of measuring the target |w after t steps as
Therefore, to obtain the target state with maximum probability, one performs a measurement on the system at t = π 4 √ N steps [14]. Not only does this scale faster than a classical computer; it has in fact been shown that Grover's algorithm is asymptotically optimal on an ideal quantum computer [15].
Since Grover's algorithm involves the repeated application of an identical unitary G w , it can be achieved by an appropriately chosen time-periodic Hamiltonian H (t + T ) = H (t ). In other words, Grover's algorithm is an example of a Floquet problem [16]. This allows us to define the time-independent Floquet Hamiltonian H F as [17,18]
where we choose units with h = 1 throughout. It follows that eigenstates |u α of G w satisfy
The quasienergies φ α are defined modulo 2π/T , where period T ≡ 1 represents the time required to apply G w on the qubits. The actual value of time T is irrelevant to this analysis because we can express everything using dimensionless phases rather than working with energies. Before we write the exact eigenstates and quasienergies for the Grover operator, we state a useful theorem that simplifies our analysis: the spectrum of the Grover operator is independent of the target state w. To prove this, consider conjugation by the unitary matrix
where the notation j ∈ w indicates to take a product over Pauli X operators acting only on sites where w is in state 1. Then
Since G 0 and G w are unitarily equivalent, their eigenvalues are identical. We will therefore focus on target state |0 ⊗L ≡ |0 for the remainder of the work. The Grover operator G 0 has two eigenstates with nonzero quasienergies, which we will refer to as the special states,
where
and
The rest of the eigenstates, which we will refer to as the bulk states, are all degenerate at zero quasienergy,
Noting that |x is a "zero-momentum" superposition of the nontarget states, we choose the following basis for the degenerate bulk subspace:
where k n = (2π n)/(N -1), and n = 1, 2, . . . , N -2. In the absence of noise, these states are gapped out, i.e., there is an energy difference (gap) between these bulk states and the states involved in the dynamics, and the evolution of the system takes place in the subspace that does not involve these states.
For any practical realization of Grover's algorithm, the Grover operator G 0 must be decomposed into one-and twoqubit gates. This can be done for any unitary operator with an appropriate universal set of quantum gates [19], but the decomposition is not unique. Therefore, choosing an efficient gate decomposition will be crucial to the performance of algorithm upon adding noise. In Fig. 1 a quantum circuit for Grover's operator G 0 for four qubits is shown. The oracle and the Grover's diffusion operators are written in terms of a multicontrolled X , a.k.a. L-qubit Toffoli, gate for ease of decomposition, as such gates are more commonly studied in the literature than the desired multicontrolled Z gate. Decomposing the L-qubit Toffoli can be done in many different ways, such as using automated transpilation methods (e.g., "transpile" in Qiskit) or via the general algorithm described in the classic paper by Barenco et al. [20]. However, these decompositions result in circuits whose depth scales nonoptimally with system size; it appears to scale exponentially with L using Qiskit transpilation and is known to scale as L 2 using the Barenco algorithm. In anticipation of improving performance on the noisy algorithm, we instead use the linear depth decomposition algorithm described in [21]. For a multicontrolled X gate with L qubits, this algorithm provides a decomposition with circuit depth 8L -20 and with 2L 2 -6L + 5 controlled x-rotation gates [22], hereafter denoted by CR x (θ ). An example for the case of L = 4 is shown in Fig. 2. Note that, while the depth is linear, it relies on the ability to perform multiple finite-range [23] control gates in parallel. Note also that here a single Floquet unitary already has extensive depth ∼L, as opposed to more conventional Floquet operators whose depth is usually L-independent.
We consider a coherent systematic noise model with random under or over-rotation of each single or two-qubit quantum gate,
where the index k refers to the kth gate in the decomposition of the Grover operator, k is drawn from the uniform distribution [-1, 1], and δ is the noise strength. Importantly, each Floquet cycle has an identical noise configuration, so that the operator G w is identical across time steps. This ensures that each iteration of the Grover algorithm is a (random) Floquet problem.
To implement this systematic noise model for our onequbit gates, we start from unitaries of the form U k (0) ≡ σ k ∈ {X, Z, H} which represent single-qubit π rotations around the x, ẑ, and ( x + ẑ)/ √ 2 axes, respectively. We can write this as
represents a rotation by angle θ about the σ k axis. Then we define the noisy gates as a random over-or under-rotation, Similarly, we write the noisy two-qubit controlled rotations as
where the convention is chosen such that
Note that an identical noise amplitude is added to each control gate independent of the magnitude of the desired rotation θ k . For small θ k δ, this means that the gate is predominately noise, suggesting that further improvements can be made by simply neglecting gates that fall below this threshold. Doing so would complicate the analysis below, but will be worth investigating in future research.
To begin exploring the effect of system noise on the performance of the algorithm, we numerically calculate the Floquet quasienergies and eigenstates of the Grover operator in the presence of systematic noise. To probe the thermalization of the bulk states, we calculate the bipartite von Neumann entanglement entropy S A of each eigenstate. We choose a subsystem A of L/2 consecutive qubits (with periodic boundary conditions) and calculate the entropy S A = -Tr(ρ A ln ρ A ), where ρ A = Tr Ā (|ψ ψ|) is the reduced density matrix in eigenstate |ψ . For each eigenstate, the entanglement entropy is averaged over the L/2 independent choices of subsystem A.
An example of the quasienergies and their entanglement entropies are shown for a single noise realization in Fig. 3. The bulk quasienergies have their degeneracy broken by the noise. Moreover, the quasienergies of these bulk states show a clear linear dependence on the noise strength δ, which comes from the exponential degeneracy of the unperturbed spectrum. This suggests that the bulk properties are well described by first-order perturbation theory, even at relatively large δ, which we use in the next section to derive an effective Hamiltonian for the quasienergy spectrum.
In addition to this linear δ dependence of quasienergy, the bulk eigenstates show entanglement characteristic of a thermal quantum system, where they should match the thermal entropy S(E ). The entanglement peaks in the middle of the bulk spectrum at approximately its maximum value [24]. As expected, the entropy density approaches zero at either edge of the bulk spectrum. Further evidence for thermalizing (random matrix) behavior will be given in Sec. III B 1.
Finally, below some critical value of noise strength δ < δ c,gap , the special states remain gapped from the bulk. Since Grover's algorithm involves coherent oscillations within the special state manifold, this suggests that there is a "computing" regime at small δ where finite probability of finding the target state will survive the systematic noise. We will indeed show that this is true at small enough δ though, as we will see in our analysis of the special states (Sec. III D), the computing transition happens at a much smaller critical value δ c,comp δ c,gap . Anticipating the results that we are about to show, δ c,gap will in fact correspond to an ergodicity transition;
for δ < δ c,gap , the initial states within the special states subspace do not reach their thermal equilibrium expectation values, i.e., retain memory of their initial conditions even after a long time, while for δ > δ c,gap , all states obtain their thermal expectation values.
We start by using our observation of clear linear δdependence of the bulk spectrum to derive an effective Hamiltonian from first-order degenerate perturbation theory. To achieve this goal, let us begin by writing the Grover operator with noise strength δ as
where we have written each gate in the decomposition of the Grover operator as generalized rotations in Hilbert space [see Eq. (19)]. Here M = 4L 2 -6L + 10 is the total number of gates in the decomposition of the L-qubit Grover operator. Individual commuting gates are divided into separate terms in the product. Expanding the exponentials and keeping terms only up to first order in δ, we obtain
FIG. 4. Plot of noise-averaged difference between the exact quasienergy of G(δ) and eigenvalues of δH eff for L = 8 using Eq. ( 25). Error bands represent the standard error of mean for five noise realizations.
where we have defined the effective Hamiltonian as
We identify H eff as the generator of gate perturbations in the Heisenberg picture [25]. To demonstrate that H eff correctly describes the exact bulk quasienergies, let us consider the action of G(δ) on these eigenstates. Using Eq. ( 15), and keeping terms within the first order in δ, we get
Thus, after discarding the two special states from H eff , its eigenvalues should match with the exact bulk quasienergies of G(δ). In Fig. 4 we plot the noise-averaged difference between the exact and the effective energies
as a function of noise strength δ. For small δ, the kth quasienergy can be written as
. From this, we can write φ exact,k /δ -E eff,k = O(δ 2 )/δ. Then plotting φ exact,k /δ -E eff,k as a function of δ will asymptote to zero as δ → 0 if the higher order terms δ 3 , δ 4 are small. In Fig. 4 we plot the disorder averaged sum of all these k-quasienergies.
We see that deviations vanish linearly as δ → 0, confirming the correctness of first-order perturbation theory. Indeed, we numerically see that first-order perturbation accurately captures the spectrum nearly up to the value of δ when the bulk spans the entire Floquet Brillouin zone from -π to π . It is worth noting that our effective Hamiltonian is defined to be independent of the scaling parameter δ, which allows us to analyze the statistical properties of H eff for different noise strengths and therefore determine many important properties such as thermalization of the bulk states and scaling of δ c,gap .
Having shown the accuracy of H eff in modeling our spectrum, we now turn to analyzing its properties with a focus on the bulk states. H eff consists of random terms which start as 2-local terms before Heisenberg evolution causes operator spreading. Some remnant of the locality remains, as we argue numerically in Appendix A. However, lacking any other structure, we would expect this Hamiltonian to exhibit generic, thermalizing behavior. An exact analytical treatment of this Hamiltonian is not feasible, therefore we resort to analyzing its properties statistically using random matrix theory, which is a natural and powerful framework to study disordered many-body systems. One way to test if the Hamiltonian belongs to a random matrix ensemble is to calculate its energy level statistics and the Kullback-Leibler divergence that captures the level repulsion, which is a characteristic feature of a thermalized phase. On the other hand, if the energy level statistics is close to the Poisson distribution, the Hamiltonian will display properties of a nonergodic system. We now show that H eff thermalizes by numerically demonstrating its consistency with random matrix theory.
In the thermal phase, the energy levels of a many-body Hamiltonian show level repulsion, and their statistics follow that of a random matrix, the underlying ensemble of which depends on the symmetries of the Hamiltonian in the problem. Here we examine two properties of H eff to demonstrate signatures of random matrix theory (RMT). First, to indicate the existence and amount of level repulsion (which is characteristic of the RM ensemble), we examine the level spacing ratio r n = min( E n , E n+1 )/max( E n , E n+1 ), where E n = E n -E n-1 is difference between the ordered eigenenergies of H eff [26]. Second, to study the ergodic structure of the eigenstates, we calculate the Kullback-Leibler (KL) divergence (also called the relative entropy) in the computational basis. It is defined as KLd n = dimH i=1 p i ln(p i /q i ), where p i = | i|n | 2 and q i = | i|n + 1 | 2 are the occupation probabilities of two energetically neighboring eigenstates |n and |n + 1 in the Z (computational) basis {|i } [27,28]. To calculate these two metrics, we sample H eff for different noise realizations, and calculate the level spacing ratio and KL divergence for each, and then average over noise realizations. The results, shown in Fig. 5, indicate that the noise-averaged level spacing ratio r n and KL divergence KLd n are consistent with that of Gaussian Unitary Ensemble (GUE): r GUE = 0.599 [29], KLd GUE = 1 [30]. This confirms that the system will thermalize within the bulk state manifold that is well described by H eff ; note, however, that the data shown actually include all states, not just the bulk. Furthermore, GUE indicates that thermalization will occur with no conserved quantities besides energy, as is perhaps expected from the unstructured form of H eff . Note, however, that both the level spacing ratio and K-L divergence probe only properties of states that are nearby one another in energy. Therefore, other RMT properties that rely on the complete spectrum and full random matrix form-such as the Wigner semicircle law-are not in general expected to hold. In fact, we will see that the locality structure that exists in H eff prevents the formation of a Wigner semicircle, as is known to happen in other local ergodic Hamiltonians [31].
Having shown that H eff obeys certain properties of RMT, we now use that to solve for other properties of the noisy algorithm, starting with the critical noise δ c,gap and its scaling with system size. We start by noting that the algorithm stops giving meaningful answers once the two special states enter into the bulk, i.e., once the gap closes. As evident from Fig. 3, this happens near the point where the bulk energies span from -π to π due to weak dependence of the special state energies on δ. We will study δ dependence of the special state energies in Sec. III D. We define the noise value at this gap-closing point as the critical noise δ c,gap . As indicated by our numerical findings, the quasienergy of the noisy Grover operator is well captured by the effective Hamiltonian up to the point of merger of the special states, therefore we will study the eigenvalues of the latter to determine the critical noise and its scaling with system size.
The relevant property of H eff that determines δ c,gap is the second moment of its energy distribution. While each realization of H eff will have nonzero trace (first moment), the fact that H eff is generated by Paulis implies that this term must be proportional to the identity matrix and therefore creates a uniform energy shift that does not affect the gap-closing transition [32]. To determine the second moment, it is convenient to first shift the spectrum to make H eff traceless: which can be done by individually removing the traces from the generators, Hj = H j -Tr(H j )/N. Then the variance is given by
where we used that E ( j k ) = 1 3 δ jk . Notice that E 0 ∼ L comes from the overall gate count, M ∼ L 2 ; by reducing scaling of the gate count, we would modify this system-size dependence. This scaling is corroborated numerically in Fig. 6.
Having obtained an accurate description of the bulk spectrum, we can derive the scaling of the critical noise δ c,gap . Since the Gaussian distribution does not have a sharp cutoff, let us instead consider when the eigenvalues within, say, 3σ FIG. 7. Critical noise strength δ c,gap as a function of system size shown in a log-log plot. The exact critical noise is calculated at the point where the gap between the two special states and the bulk closes and is averaged over 20 noise realizations. The dashed line on the plot represents the best-fit line for these exact noises. Predicted critical noises are averaged over 512 noise realizations and are calculated when the two special states merge into 67%, (-2E 0 , 2E 0 ), and 99.7%, (-3E 0 , 3E 0 ) of the bulk eigenvalues spanning the Floquet Brillouin zone.
span the entire Floquet Brillouin zone. Then we can write
which gives
It is easy to see that irrespective of the fraction of eigenvalues we consider that span the entire Brillouin zone, the scaling of the critical noise does not change. In Fig. 7, we illustrate this scaling behavior by considering the bulk eigenvalues within 2σ and 3σ intervals. The parameter δ c,gap does not represent a conventional phase transition since it scales as a power law with system size L. However, one can define the observable: the late-time average of occupation in the state |x , which vanishes at δ = δ c,gap in a manner characteristic of conventional second-order phase transitions. An alternative perspective is to consider δ eff = Lδ, which represents the aggregate error rate across all qubits. In terms of this rescaled noise parameter, the transition occurs at δ eff = O(1). This type of parameter rescaling is common in both mean-field theory models and nonstandard thermodynamic limits [33].
Our analysis gives the correct system-size scaling of δ c,gap , but also shows that a full random matrix treatment is incorrect. Specifically, the drastically different variances of the diagonal and the off-diagonal elements of H eff in Fig. 6 deviate from the prediction of a Gaussian Unitary Ensemble (GUE):
for 1 m, n N, and where E GUE [•] denotes the average over different realizations of the GUE and λ sets the units of energy [34]. As a result of this discrepancy in the scaling of elements, the distribution of the eigenvalues of H eff deviates from Wigner's semicircle law, instead following a normal distribution
Similar deviations from the semicircle law in the energy spectra of many-body systems have been noted in other interacting many-body systems [35][36][37], leading to the conclusion that completely random matrices do not accurately represent realistic many body systems. The underlying assumption of these matrices-that each particle interacts simultaneously with all other particles-contradicts the nature of most systems where interactions are predominantly few-body. The presence of k locality here is not immediately obvious due to operator spreading of the 2-local perturbations when taken in the Heisenberg picture (see Eq. 23). We find strong evidence for approximate k locality from a detailed analysis of matrix elements in H eff , which is done in Appendix A.
It is worth noting that the exact value of the gap closing between bulk and special states will not happen precisely at δ c,gap derived above, but rather a slightly smaller value. This is because the edges of the spectrum of H eff do not scale identically as its variance. For example, conventional local models have an extensive energy range but have energy width scaling as √ L. However, the energy width E 0 ∼ L is the most important parameter for determining ergodicity, as this is the energy scale on which a finite density of states will come in contact with the special states. The precise scaling of the range of bulk energies, which determines initial gap closing but not ergodicity breaking, will be left for future work.
Finally, we turn our attention to the two special states. In the absence of noise, the actions of the Grover operator on the states |0 and |x are
Then we can write the noiseless Grover operator in the {|0 , |x } block as
Now, adding back in the systematic noise, we can repeat our analysis using the effective Hamiltonian within first-order perturbation theory and apply it to these special states. Restricting Eq. ( 24) to the special states manifold by writing it in |0 and |x basis,
we obtain
FIG. 8. Plot of difference between the exact special states quasienergy of G spec (δ) and eigenvalues of -θτ y + δH spec eff for L = 8. Error bands represent the standard error of mean for five noise realizations.
The eigenvalues of -2θτ y + δH spec eff -π capture the quasienergies of two special states, as shown in Fig. 8, once again confirming the effectiveness of first-order perturbation theory.
Next, let us analyze the statistical properties of H spec eff to understand the behavior of the two special states in the presence of noise. Taking the logarithm on both sides of Eq. ( 40):
we can expand in the Pauli basis to obtain
where b j = tr(H spec eff τ j )/2 with τ x,y,z,0 representing the Pauli matrices (τ 0 = 1). Then the square of the energy gap between the two special states is given by
which sets the oscillation frequency between the special states. Averaging over noise configurations and noting that E[b x,y,z,0 ] = 0, we see that the gap is set by their variance:
If H eff were an N × N GUE random matrix, then H spec eff would be a 2 × 2 GUE random matrix, for which E (b 2
x,y,z ) ∼ 1/N. However, as we have seen that non-GUE scaling of the diagonals exists in H eff , we need to more carefully examine whether GUE scaling holds within the special state manifold. As shown in Fig. 9, indeed we do not have exponential scaling of these Pauli coefficients nor the gap. While the off-diagonal terms do decrease exponentially fast θ 2 , b 2
x , b 2 y ∼ 1/N, the diagonal term instead increases as b 2 z ∼ L. Putting these scalings into Eq. ( 42), we see that E[ 2 ] ∼ A/N + Bδ 2 /N + CLδ 2 for some order-1 coefficients A, B, and C. Then, for δ δ c,comp ≡ √ A/(CLN ) ∼ 1/ √ L2 L , the linear b z term dominates over the θ 2 term giving a power-law-large gap. While this increases the oscillation frequency within the special states, it comes at the cost of exponentially decreasing the oscillation amplitude. In the limit of δ δ c,comp , the states |0 and |x become eigenstates of G spec (δ) and therefore their probabilities remain at at all times T , up to exponentially small oscillations. Since there is no longer an order-1 probability of finding the target state at any times, this implies that the model has lost computing power. Therefore, our special state analysis unfortunately implies that at this much smaller critical noise δ c,comp = O(2 -L/2 ) we have a "computing transition," because the algorithm loses computing power for δ δ c,comp . By contrast, the gap-closing critical noise δ c,gap = O(L -1 ) is much larger, but occurs well after computational power has been lost.
Given that the large coefficient b z in H spec eff implies a special basis for the special state manifold, one might hope to recover an exponentially small gap by randomly rotating the basis, e.g., by choosing a different noise model. However, we now argue that this is not possible for any local or k-local noise model. The key point is that the relatively larger magnitude of the x,y for generic k-local error models.
We analyzed Grover's algorithm with systematic noise by modeling the Grover operator as a Floquet unitary. We showed that the bulk energy spectrum is well described by an effective Hamiltonian whose spectral statistics match the Gaussian unitary ensemble but whose Gaussian eigenvalue distribution reflects the underlying k-local nature of the noise model. We use this structure to derive sharp crossovers between behaviors, finding a gap-closing ergodicity transition at relatively large noise strength δ c,gap = O(L -1 ) and a separate transition where computational power is lost above an exponentially small noise strength δ c,comp = O(2 -L/2 ). The random matrix theory model enables analytical predictions for many quantities, such as these noise thresholds and dynamics of the initial state.
The machinery we have developed for this simplest case of Grover's algorithm with a single solution should provide insights into more generally applying many-body physics approaches to noisy quantum algorithms. For example, we can readily see that Grover's algorithm with multiple solutions exhibits the same phenomenology since the core aspects of a degenerate gapped bulk and cat-like special eigenstates hold there as well. While we specifically focus on unraveling both the oracle and Grover diffusion operator, the same effect will occur if we just unravel the Grover's diffusion operator (reflection about |x ) since it is identical to the oracle for a single target state. Introducing noise to only one of these operators, rather than both, would merely adjust the critical disorder thresholds by approximately a factor of 2, affecting only the scaling prefactor rather than the fundamental behavior. A more interesting variant is to instead consider decompositions of the multiqubit Toffoli gate using ancilla qubits. If the ancillas are erased before reuse, they will behave similar to a conventional noise channel, which previous results suggest will yield an exponentially small error threshold. On the other hand, if the ancillas are uncomputed and the measured and postselected, the resulting dynamics is noisy but nonunitary, similar to dynamics generated by non-Hermitian Hamiltonians. Non-Hermitian dynamics are known to have qualitatively different phase structures than their Hermitian counterparts. How this affects the phase structure and noise tolerance of Grover's algorithm remains a fascinating open question.
Another important avenue to investigate for Grover's algorithm is the role of non-k-local gating. While most controllable quantum computers rely on k-local gates-even if swap operations or atomic repositioning relax the assumption of strict spatial locality-there are potential routes to remove this restriction by using a nonlocal cavity mode or the longrange tails of Rydberg blockade [38]. This not only allows more efficient decomposition of the multicontrolled gates, but also implies that systematic gate errors will be nonlocal as well. Nonlocal errors in H eff will cause it to behave more like a proper GUE random matrix, for which a power-law noise tolerance δ ∼ 1/L α can in principle be admitted for both the ergodicity transition and the computational transition. In the NISQ era of finite noise, it is crucial to know these noise scalings in order to determine near-term experimental feasibility. We therefore suggest that nonlocal gating may be a pathway to drastically improved noise tolerance for Grover's algorithm, based on our analytical understanding through random matrix theory.
In the context of many-body Floquet systems, decorrelating unitary errors between Floquet cycles often has the effect of destabilizing the phase of matter rather than slowing its decay. For example, in the topological Floquet-Anderson insulator, one of us showed in previous work [39] that it is stable to weak systematic noise, whereas uncorrelated noise between Floquet cycles leads to eventual decay. A similar situation occurs in many other Floquet phases of matter, such as time crystals, which are absolutely stable to weak systematic perturbations as one of us showed in [40] but unstable to nonsystematic noise [41]. However, our analysis suggest that the application of twirling operations to decorrelate the noise could significantly mitigate the impact of errors, improving the scaling from 2 -L/2 to 2 -L/4 . This observation is consistent with the results of Shapira et al. [12], although their analysis did not fully consider the effect of unitary noise in the oracle or the Grover diffusion operator.
Finally, while Grover's algorithm provides a concrete example of Floquet analysis in the context of useful quantum algorithms, the eventual goal is to apply similar tools to a wider class of algorithms. Rather than approaching such problem on an individual basis, we would hope to find inspiration for doing so in "grand unified" theories of quantum computation [42]. While no Floquet structure is immediately apparent in such prescriptions, we might hope to find a constituent block on which similar analysis can be applied and use that insight to understand systematic noise tolerance of a wider class of quantum algorithms.
In addition to gap-closing and computational transitions, the random matrix treatment gives analytical predictions for nearly all properties of the algorithm's behavior. In this Appendix, we describe the behavior of special state dynamics, which is the most relevant dynamical feature that allows the unperturbed algorithm G 0 to be successful. While we can make predictions for how systematic noise modifies these dynamics, their computational relevance is limited to the very small regime of δ δ c,comp . However, we can also make concrete predictions for ergodicity breaking as long as δ < δ c,gap .
Before directly examining dynamics, it is worth briefly thinking about the behavior of perturbation theory in this degenerate many-body setting. In particular, given that we start within the special state manifold, does a finite probability of staying in this manifold survive at late time?
In order to address this question, we define the effective Hamiltonian within the special state manifold as
where θ ≡ θ (δ). Then we write the eigenstates of H (δ) as
with energies (0) 1,2 near π. Then we can write
where H eff is divided into H spec eff , a block containing the two special states, and the remainder H eff = H eff -H spec eff such that G|ψ expanding both sides up to the first order in δ, we get
Therefore, the net bulk occupation in the dressed eigenstate |ψ j is approximately
While the exact expression for this off-diagonal matrix element, | xk |H eff |ψ (0) j | 2 , depends on the particular choice of k, the sum over these matrix elements is constrained to be no larger than E 2 0 (at least upon averaging over j) due to the fact that Tr[H 2 eff ] = E 2 0 N can be broken up into diagonal and off-diagonal pieces and the diagonals already given the E 2 0 N scaling. Therefore, the off-diagonal sum must give E 2 0 N or smaller, and the lack of summing over j in Eq. (B6) removes the multiplicative factor N. Assuming that the sum saturates this scaling, we conclude that
Therefore, the bulk occupation remains finite below δ c,gap , meaning that the special states will also remain finitely occupied at late times and thus break ergodicity as we have claimed.
Next, let us analyze the behavior of the probabilities P 0 (t ) and P x (t ) in the presence of noise. Consider the following ansatz for the two special states within first order in noise:
where we will see that P b is a related parameter that describes late-time occupation of the bulk given our actual starting state. This is an unusual ansatz, as it dresses the zeroth-order result from degenerate perturbation theory by a nonperturbative overall rescaling of the amplitude, (1 -P b ) 1/4 . This comes from the finite "leakage" into the bulk at late times but, as we'll show, seems to give a self-consistent picture of the full dynamics.
With this ansatz, we can write the initial state of the system as
where k labels the N -2 bulk energy eigenstates. But |x is normalized, so
Having decomposed the initial state into eigenstates [Eq. (B8)], we can now easily calculate its time evolution as
where t is the number of time steps. Then we simply need to calculate observables. One important observable is the probability to find ourselves in |x , which is approximately the Loschmidt echo: (B9)
Let us analyze Eq. (B9) term by term.
The existence of a non-exponentially small occupation in state |x is a clear indication of nonergodicity. It will differ between noise realizations, but we can estimate its noise average (which is equal to the participation ratio in the eigenstate Before continuing, it's worth doing a similar calculation for the target state:
≈ (1 -P b ) 2 cos 2 θ 2 sin 2 θ 2 ⇒ P x + P 0 = 1 -P b ;
as expected, P b is the late-time probability to leave the special state manifold.
These are the oscillations that occur at frequency ∼ 1/ √ N in the absence of noise and with notable modifications in both amplitude and frequency in the range δ c,comp < δ < δ c,gap . Specifically, it is clear that, when P b > 0 and/or θ = π/2 (as happens for finite δ), oscillation amplitude is less than unity.
Due to the π gap between bulk and special states, there are fast (period 2) oscillations in the occupation. However, since the bulk has finite spectral width, these oscillations are damped. The gap between E 1 and E 2 will not be resolvable on the decay timescale, so for all intents and purposes we can just set both of these to π . Using the Gaussian density of states obtained in the main text,
e -E 2 /(2δ 2 E 2 0 ) , we get
which is simply Gaussian relaxation on the timescale τ ∼ 1/(δE 0 ) ∼ δ c,gap /δ that we expect to see. Notice that this is a fast O(1) timescale as long we keep δ/δ c,gap fixed as we take L large.
Finally, the last term tells us how Loschmidt echo of a given state-in this case |x -decays when propagated by the random bulk matrix H eff . Again, using the Gaussian FIG. 11. Plot of average probability P 0 (t ) and P x (t ) for L = 10 as a function of time steps for different noise strengths: 0.01 < δ c,comp ≈ 0.02 < 0.06 < 0.08 < 0.1 < δ c,gap ≈ 0.15 < 0.17. When δ < δ c,gap , the probabilities oscillate with a larger frequency but smaller amplitude as compared to the noiseless case. The noises δ = 0.06, 0.08, 0.1 lie between δ c,gap and δ c,comp , for which there is no significant change in the probabilities as anticipated from the above analytical calculation. Noise δ = 0.17 > δ c,comp represents the noise-dominated dynamics. Error bands represents the standard error of mean for 256 noise realizations. approximation for the density of states, we find
= (1 -1 -P b ) 2 e -δ 2 E 2 0 t 2 . This is similar to the decay rate of the 2T -periodic oscillations (term C) with a rate that is √ 2 times faster. Finally, as a sanity check, let us confirm that P x (0) = 1: The predicted behaviors above are seen in Fig. 11. Note that the period-2 oscillations do appear on close inspection of the short-time behavior (cf. the sharp jump in the curve for δ = 0.17 in orange), but are not clear on the longer timescale plotted here.