title: Concentration for random product formulas
authors: Chi-Fang Chen Hsin-Yuan Huang Richard Kueng Joel A. Tropp
text: 
  
  
  
  
  
  These authors contributed equally.
  Department of Physics, Caltech, Pasadena, CA, USA
  
  These authors contributed equally.
  Institute for Quantum Information and Matter, Caltech, Pasadena, CA, USA
  Department of Computing and Mathematical Sciences, Caltech, Pasadena, CA, USA
  
  Institute for Quantum Information and Matter, Caltech, Pasadena, CA, USA
  Department of Computing and Mathematical Sciences, Caltech, Pasadena, CA, USA
  Institute for Integrated Circuits, Johannes Kepler University Linz, Austria
  
  Department of Computing and Mathematical Sciences, Caltech, Pasadena, CA, USA
  
  
   Quantum simulation has wide applications in quantum chemistry and physics.
      Recently, scientists have begun exploring the use of randomized methods for accelerating quantum simulation.
      Among them, a simple and powerful technique, called qDRIFT, is known to generate random product formulas
      for which the average quantum channel approximates the ideal evolution.
      qDRIFT achieves a gate count that does not explicitly depend on the number of term in the Hamiltonian, which contrasts with Suzuki formulas.
      This work aims to understand the origin of this speed-up by comprehensively analyzing a single realization of the random product formula produced by qDRIFT.
      The main results prove that a typical realization of the randomized product formula approximates the ideal unitary evolution up to a small diamond-norm error.
      The gate complexity is already independent of the number of terms in the Hamiltonian, but it depends on the system size and the sum of the interaction strengths in the Hamiltonian.
      Remarkably, the same random evolution starting from an arbitrary, but fixed, input state yields a much shorter circuit suitable for that input state. In contrast, in deterministic settings, such an improvement usually requires initial state knowledge. 
      
      The proofs depend on concentration inequalities for vector and matrix martingales, and the framework is applicable to other randomized product formulas. Our bounds are saturated by certain commuting Hamiltonians.
      
  
      
      
      
      
      
  
      
      
      
      
      
      
  Concentration for random product formulas
      Joel A. Tropp
      December 30, 2023
  =========================================
  
  
  
  
  
  § INTRODUCTION
  
  
  Simulating complex quantum systems is one of the most promising applications for quantum computers.  This task has many applications, such as developing new pharmaceuticals, catalysts, and materials <cit.>, as well as solving linear algebra problems <cit.>. The task of digital quantum (dynamics) simulation can be phrased as a compiling problem:
  approximate a given unitary, say a Hamiltonian evolution U = e^-  Ht, by a product of `simple' unitaries g_k:
  
      U=e^-  H t≈ V = g_1 ⋯ g_N.
  
  Such a decomposition into elementary gates should obey two conditions: (i) It should accurately approximate the target unitary.
  In this work, we require that the error in operator norm[For measuring time-evolved observables (ρ(t)O), the operator norm suffices; for quantum phase estimation of ground state energies things can get more complicated, though. There, the estimates may depend on the details and vary orders of magnitude <cit.> and the
  assessment of real costs is an on-going research direction that is beyond the scope of this work.]
  satisfies U-V ≤ϵ for some specified accuracy parameter ϵ.
  Moreover, (ii) the decomposition should cost as little as possible.  Several cost functions make sense in this context, but most of them can be approximated by the gate complexity, i.e., the number N of simple gates[In this work, we will count ^- h_jt as a single gate,  as in <cit.>. Strictly speaking, in a fault-tolerant quantum computer, we would further decompose each ^- h_jt into universal 2-qubit gates.
  According to the Solovay–Kiteav Theorem <cit.> this incurs at most a constant (multiplicative) overhead, but the exact cost depends on the type of hardware.
  
  ] on the right-hand side of Eq. (<ref>).
  
  
  
  
  
  
  
  
  
  
  The earliest compilation procedures for quantum simulation were based on product formulas <cit.>, also known as Trotterization, or splitting methods. They depend on the idea of approximating the matrix exponential of a sum by a product of matrix exponentials.
  
  
  
  
  
  
  
  We will review one such construction below in Eq. (<ref>).
  Subsequently, alternative principles have led to the development of other quantum simulation algorithms.
  These include linear combination of unitaries <cit.>, quantum signal processing <cit.> and qubitization <cit.>.
  By and large, these algorithms rely on more powerful quantum computing primitives to yield improved performance in accuracy and cost. We refer to <cit.> for a systematic review.
  
  
  
  
  
  
  Despite these advanced simulation techniques, product formulas
  
  have recently undergone a renaissance <cit.>.
  They are not only simple and (comparatively) easy to implement on near-term devices, but they also
  
  remain very competitive <cit.> provided that they incorporate information about the structure of the problem, such as
  
   initial state knowledge <cit.>, locality <cit.> or the commutator structure of Hamiltonian <cit.>.
  The purpose of this paper is to explore the randomized aspects of constructing product formulas.
  
  
  
  We begin by reviewing the first-order Lie–Trotter formula, which will be later contrasted with a randomized variant (qDRIFT). To that end, consider a quantum many-body Hamiltonian  H = ∑_j=1^L h_j composed of L simple terms h_j. The first-order formula cycles through each term in the Hamiltonian
  
      U ≈(exp( -i (tL / N) h_L ) ⋯exp( -i (tL / N) h_1 ))^N/L.
  
   To approximate the target unitary U up to accuracy ϵ in operator norm, a total gate count N = 𝒪(L λ^2 t^2 / ϵ) suffices[To finally compile into universal gates, run a standard gate synthesis algorithm (such as Solovay–Kiteav) for each simple term. This yields another multiplicative factor on the gate count.] <cit.>.  In this expression, t is the simulation time, L is the number of terms, and λ = ∑_j h_j summarizes the interaction strengths within H.
  
   The main idea is to cancel out the leading-order term in the Taylor expansion by including each of the L terms. Owing to this construction, the factor L remains in higher-order Suzuki formulas where the gate count N = 𝒪(L (λ t)^1+ o(1) / ϵ^o(1)), even though the time-dependence becomes nearly optimal. We refer to  Table <ref> for a sharper gate count incorporating commutators.
  
  
  
  
  
  
  
  
   
  
  
  
  
  
  
  
      
  
  
  
  
  
  
  
  
  Recently, researchers started using randomization to improve the performance of product formulas <cit.>.  Campbell <cit.> introduced the qDRIFT algorithm,
  
  
  which approximates the target evolution U=^- Ht by a quantum channel that results from averaging products V_N ⋯ V_1 of random unitaries. Each V_k corresponds to a short-time evolution based on a single term h_K from the Hamiltonian.
  The index K is selected randomly, according to an importance sampling distribution (p_1, …, p_L), constructed to match the leading order of time step
  
      𝔼[ V_k ]
      ≈exp( - i (t/N)𝔼[h_K/p_K ] ) = U^1/N.
  
  This approximation is achieved by averaging over a single (random) unitary.
  (In contrast, the first order Suzuki formula (<ref>) must cycle through all terms which incurs an extra L-factor.) Independence among the V_k's then ensures
  
  
      𝔼[ V_N ⋯ V_1 ] = 𝔼[ V_N ] ⋯𝔼[ V_1 ] ≈(U^1/N)^N = U.
  
  Campbell proved that the averaged Hamiltonian approximates the true Hamiltonian
  when the gate count satisfies
  
      N = 𝒪(λ^2 t^2 / ϵ).
  
  Observe that the gate count is independent of the number of terms L in the Hamiltonian.
  We refer to Figure <ref> for a summary of this procedure.
  
  
  
  
  
  
  
  
  
  
  
  Operationally, Campbell considers a black box that applies a new random product V_N ⋯ V_1 of unitaries every time it is invoked. The average of all possible product formulas is 𝒱^(N)(X) = [V_N ⋯ V_1 X V_1^†⋯ V_N^†] and forms a completely positve trace-preserving (CPTP) map. The ideal unitary also forms a CPTP map given by 𝒰(X) = U X U^†.
  Campbell proves that (<ref>) ensures that the two CPTP maps are ϵ-close in diamond distance.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  It is interesting to compare the gate count (<ref>) achieved by qDRIFT with very recent and powerful results for (deterministic) product formulas <cit.>, see Table <ref>. Broadly speaking, deterministic product formulas can achieve gate counts that are linear in time t and the number L of terms. The qDRIFT gate count, on the other hand, is independent of L, but quadratic in t. This implies that qDRIFT should be favored for short-time simulations rather than long-time simulations. For systems with many terms in the Hamiltonian (L ≫ 1), such as in quantum chemistry and the SYK model <cit.>[While qDRIFT provides performance guarantee on rather flexible choices of models, there is a specialized quantum algorithm for simulating the SYK models using much fewer gates <cit.>.], there are indeed physically relevant time-scales for which qDRIFT outperforms Suzuki formulas; see Table <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  To summarize, there are interesting quantum simulation problems where the qDRIFT gate count (<ref>) compares favorably with very recent and powerful results about high-order Suzuki formulas. This advantage is a consequence of randomization.
  But we may ask whether the benefit arises from the random choice of an individual product formula or whether it is due to the averaging of many product formulas together.  Which aspect produces the speed-up?
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  In this work, we analyze why random product formulas are efficient.
  
  
  To do so, we study the performance of a random instance of the product formula V_N ⋯ V_1.  By contrasting this instance with the average channel 𝒱^(N)(X), we deduce that random sampling alone allows us to avoid the dependency on the number L of terms in the Hamiltonian; it is not necessary to average many different product formulas.  Our results establish strong concentration bounds for two cases, depicted in Figure <ref>.
  
  First, let n denote the number of system constituents, e.g., qubits.
  If the gate count N obeys
  
      N ≥Ω( n t^2 λ^2 / ϵ^2),
  
  then, with high probability, a single realization of the random product formula
  approximates the ideal target unitary up to accuracy ϵ in operator norm.
  
  
  
  
  By the probabilistic method, this result also establishes, for the first first time, the existence of product formulas whose gate count N is independent of the number L of terms in the Hamiltonian but depends on the system size n instead.  On the other hand, we cannot easily verify the quality of any given instance.
  
  
  
  
  
  
  
  
  
  In practice,  we often wish to evolve a fixed input quantum state ρ, which may be arbitrary and unknown.
  
  This change in the problem statement has profound implications for randomized quantum simulation.
  With high probability, a random product formula with
  
      N ≥Ω( t^2 λ^2 / ϵ^2)
  
  terms suffices to achieve an ϵ-approximation V_N ⋯ V_1 ρ V_1^†⋯ V_N^† of the ideal time-evolved state U ρ U^† with respect to trace distance. Roughly speaking, this result implies that each input state has a set of product formulas that are n times shorter than a “general-purpose” product formula that works for all input states simultaneously.
  Although the set of effective product formulas depends on the choice of state and observable,
  the formulas are all produced by the same randomized procedure.  Remarkably, this procedure
  does not exploit any knowledge of the particular input state.
  
  Note that the gate count required for the original qDRIFT protocol (<ref>) and Rel. (<ref>) only differ in the scaling with accuracy: order 1/ϵ for the full qDRIFT channel vs. order 1/ϵ^2 for a single random product formula with fixed input. This discrepancy is reminiscent of the mixing Lemma by Hastings <cit.> and Campbell <cit.>: mixing unitaries yields a “quadratic" improvement in accuracy.   See Lemma <ref> below for a statement.
  
  
  
  In this work, we analyze several classes of randomized product formulas, including qDRIFT as a particular example.   The underlying theory is far more general because it relies on powerful concentration results for matrix and vector martingales.  The approach yields strong concentration results for any product of random unitary matrices, such as the ones that arise from randomized compiling <cit.>.   We are confident that these ideas are applicable to other problems that involve functions of random matrices, such as  <cit.>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    
  *Roadmap
  The rest of this paper is organized as follows.
  Section <ref> presents the main theoretical contributions to this work in detail. These are further supported and illustrated by numerical experiments presented in Section <ref>.
  Section <ref> contains two instructive examples, as well as a non-technical illustration of the underlying proof idea. We introduce related background for martingales in Appendix <ref>. Details and rigorous arguments are provided in Appendix <ref>).
  The final appendix section  (Section <ref>) establishes
  asymptotic tightness for time-evolving two simple (commuting) Hamiltonians.
  
  
  
  
  § MAIN RESULTS
   
  
  
  This section gives rigorous results for the error incurred by a randomly sampled product formula V_N ⋯ V_1,
  as compared with the ideal unitary evolution operator U = exp(-i H t).
  The results depend on the distance measure and the particular setup,
  which we discuss separately in the following subsections.
  
  
  
   §.§ Error bound in diamond distance: Worst-case error over all input states and observables
  
  
  The diamond distance is a standard distance measure for quantum channels. To compare two unitaries U_1 and U_2,
  the diamond distance is equivalent to
  
      dist_♢(U_1, U_2 )    = max_|ψ⟩:  stateU_1 |ψ⟩⟨ψ| U_1^† - U_2 |ψ⟩⟨ψ| U_2^†_1
         = max_|ψ⟩:  state   max_O: O≤ 1| ⟨ O ⟩_U_1 |ψ⟩ -  ⟨ O ⟩_U_2 |ψ⟩|,
  
  where ·_1 is the trace norm and ⟨ O ⟩_|ϕ⟩ = ⟨ϕ| O |ϕ⟩ is the expectation value of an observable O for the quantum state |ϕ⟩.
  
  Operationally, this means that the expectation value of O evaluated on the output state would differ at most by the diamond distance between U_1, U_2 for any input quantum state |ψ⟩ and any observable O with eigenvalues in [-1, 1].
  
  
  
  
  
  
  
  Our first main result bounds the gate complexity that suffices to guarantee that the randomly sampled product formula V_N ⋯ V_1 is close to the ideal evolution exp(-i t H) in this worst-case error metric.
  
  
  
  Consider an n-qubit Hamiltonian H = ∑_j h_j with λ = ∑_j h_j.
  Draw a randomized product formula V_N ⋯ V_1 from (<ref>) with gate count
   
      N ≥Ω((n+ log(1 / δ)) t^2 λ^2 / ϵ^2 ).
  
  With probability at least 1 - δ, the diamond distance error satisfies
   
      max_|ψ⟩:  state   max_O: O≤ 1| ⟨ O ⟩_V_N ⋯ V_1 |ψ⟩ -  ⟨ O ⟩_exp(-i t H) |ψ⟩| < ϵ.
  
  
  
  To keep notation as simple as possible, we have formulated this result for pure input states |ψ⟩.
  Convexity readily allows for extending the bound to mixed input states ρ = ∑_i p_i |ψ_i ⟩⟨ψ_i| as well.
  
  
  A complementary result bounds the expected approximation error in terms of gate count N.
  
  
  
  
  
  
  
  
  
  Consider an n-qubit Hamiltonian H = ∑_j h_j with λ = ∑_j h_j.
  A randomized product formula V_N ⋯ V_1, drawn from (<ref>),
  has expected diamond-distance error
  
  
      𝔼[ max_|ψ⟩:  state   max_O: O≤ 1| ⟨ O ⟩_V_N ⋯ V_1 |ψ⟩ -  ⟨ O ⟩_exp(-i t H) |ψ⟩| ] 
      ≲   √(n t^2 λ^2/N).
  
  
  
  
  The symbol ≲ applies in the large-N regime and suppresses constants.
  The proof sketch is presented in Section <ref>, and the detailed proof is given in Section <ref>.
  
  For comparison, the error bounds for the average quantum channel, established in <cit.>, imply that
  
      max_|ψ⟩:  state   max_O: O≤ 1| 𝔼_V_1, …, V_N [⟨ O ⟩_V_N ⋯ V_1 |ψ⟩] -  ⟨ O ⟩_exp(-i t H) |ψ⟩| 
      
      ≲   t^2 λ^2/N.
  
  
  As we can see, the error bound of the average over all product formulas is smaller than the error for an individual random product formula.
  
  To understand the discrepancy, it is valuable to think about a randomly sampled product formula as a random walk on the group of 2^n × 2^n unitary matrices.  Figure <ref> indicates why a single realization of the random walk V_N ⋯ V_1 has much greater error than the average of the random walks.
  
  In the error bound 𝒪(√(nt^2 λ^2 / N)), the square root reflects the statistical nature of the fluctuations in the random walk around its expectation.  The diamond norm requires us to control the behavior of the random product formula when applied to every 2^n-dimensional input state.  Remarkably, we only pay for the logarithm of the dimension, which coincides with the number n of qubits. We will discuss an intuitive example, where log(2^n) naturally arises from a union bound in Section <ref>. Formally, this pre-factor is a general feature of concentration inequality for matrix martingales (Sec. <ref>). Similar proof techniques could potentially be useful for controlling stochastic errors in other quantum computing applications.
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Faster simulation for a fixed input state
  
  
  
  
  
  
  
  
  
  
  
  
  
  In practice, it is common to perform the quantum simulation starting from a particular input state.
  The distinction with the previous setting is that the product formula only needs to work for one (arbitrary and possibly unknown) input state,
  not for all states simultaneously.  The next theorem asserts that much shorter product formulas suffice in the easier setting.
  
  
   
   Consider an n-qubit Hamiltonian H = ∑_j h_j with λ = ∑_j h_j and any input quantum state |ψ⟩.
  Draw a randomized product formula V_N ⋯ V_1 from (<ref>) with the number of gates
   
      N ≥Ω(t^2 λ^2 log(1 / δ) / ϵ^2).
  
   With probability at least 1 - δ, the output state V_N ⋯ V_1 |ψ⟩ satisfies
   
      max_O: O^† = O, O≤ 1| ⟨ O ⟩_V_N ⋯ V_1 |ψ⟩ -  ⟨ O ⟩_exp(- i t H) |ψ⟩| < ϵ,
  
   where ⟨ O ⟩_|ψ⟩ = ⟨ψ| O |ψ⟩.
   This is equivalent to the output state V_N ⋯ V_1 |ψ⟩ being ϵ-close to the ideal output state exp(-i t H) |ψ⟩ in trace distance.
  
  
  Again, convexity allows us to extend this bound to a (fixed) mixed state ρ = ∑_i p_i |ψ_i ⟩⟨ψ_i| as well. And, a complementary result bounds the expected approximation error in terms of gate count N.
  
  
  
  Consider an n-qubit Hamiltonian H = ∑_j h_j with λ = ∑_j h_j.
  A randomized product formula V_N ⋯ V_1, drawn from (<ref>),
  has expected trace distance error for any fixed input
  
  
      max_|ψ⟩:  state   𝔼[  max_O: O≤ 1| ⟨ O ⟩_V_N ⋯ V_1 |ψ⟩ -  ⟨ O ⟩_exp(-i t H) |ψ⟩| ] 
      ≲   √( t^2 λ^2/N).
  
  
  
  
  
  
  
  Theorem <ref> yields an n-fold improvement over the gate count from Theorem <ref>.
  
  So, for a 100-qubit system, focusing on a single input state leads to a 100× reduction
  in gate complexity over a simulation that works for all input states.
  The product formulas that work for a fixed input state may vary with the choice of state,
  but they are all generated using the same qDRIFT procedure.
  Taking a union bound, we can also construct short product formulas that work
  for a moderate number of input states by increasing the gate complexity slightly.
  
  
  
  The proof of Theorem <ref> is similar in spirit to the proof of Theorem <ref>. We construct a random walk from the (fixed) starting state |ψ⟩.
  We show that, with high probability, the output state is close to the ideal output state U |ψ⟩. However, what marks the difference from the unitary results (Theorem <ref>) is that
  vector concentration inequalities suffice.
  More precisely, we analyze the vector random walk using a geometric tool, called uniform smoothness <cit.>. The details appear in Section <ref>.
  The resulting n-fold improvement can also be understood as a result of switching the order of quantifiers, see <ref> for an explicit example.
  
  
  
  
  
  
  
  
  
  
  
   §.§ Even faster simulation for fixed input states and fixed observables
  
  
  
  We have seen that restricting attention to a fixed (but arbitrary) input state can yield considerable gate count improvements for random product formulas.
  It should not come as a surprise that additional savings are possible if one restricts attention to predicting a finite collection of outcome observables only (instead of demanding an accurate approximation for all observables).
  
  
  
  In this section, we consider the task of accurately estimating M expectation values ⟨ O_j ⟩_exp(-i t H) |ψ_j⟩ associated with the time evolutions of (potentially distinct) input states |ψ_j⟩.
  
  Estimating expectation values necessitates quantum measurements. In turn, the probabilistic nature of the measurement outcomes requires multiple samples
  (and, thus, repetitions of the time evolution procedure) for each pair of input and observable.
  
  This unavoidable bottleneck motivates the following procedure:
  
  We first sample R random product formulas: V^(1)_N ⋯ V^(1)_1, …, V^(R)_N ⋯ V^(R)_1. For the jth pair of input state |ψ_j⟩ and observable O_j, we perform R measurement[We perform single shot measurement for each repetition rather than estimate the average of the observable O_j. This is equivalent to measuring in the eigenbasis of O_j to obtain a single eigenvalue.]
  repetitions:
  
      repetition 1:|ψ_j⟩→    V^(1)_N ⋯ V^(1)_1 |ψ_j⟩
      →   Measure  O_j  to get ô_j^(1).
      ⋮
      repetition R:|ψ_j⟩→    V^(R)_N ⋯ V^(R)_1 |ψ_j⟩
      →   Measure  O_j  to get ô_j^(R).
  
  Subsequently, we approximate each target expectation value by the corresponding empirical average:
  
      ⟨ O_j ⟩_exp(-i t H) |ψ_j⟩≈1/R∑_r = 1^R ô_j^(r) =: ô_j.
  
  Naïvely applying Theorem <ref> would require gate count N ≥Ω(t^2 λ^2 / ϵ^2) and R ≥log (M)/ϵ^2 repetitions
  [The factor 1 / ϵ^2 in R comes from the statistical noise that accompanies Monte Carlo averaging within measurements. The factor log(M) in R comes from a union bound that ensures that all M estimates are accurate.]
  to ensure that all M estimators ô_1, …, ô_M have additive error
  less than ϵ.
  However, the fact that we are only measuring a particular set of observables O_1, …, O_M
  allows us to reduce the gate complexity N to
  
  Ω(t^2 λ^2 / ϵ).
   
   Consider a 
   Hamiltonian H = ∑_j h_j with λ = ∑_j h_j and M (arbitrary) input–observable pairs
   
   (|ψ_j⟩, O_j), 1 ≤ j ≤ M.
   Set
   
      N ≥Ω(t^2 λ^2 / ϵ)  R ≥Ω(log(M / δ) / ϵ^2).
  
   Draw R randomized product formulas of length N according to the qDRIFT procedure (<ref>).
   Then, with probability at least 1 - δ, the corresponding empirical averages ô_j defined in Eq. (<ref>) obey
   
      | ô_j -  ⟨ O_j ⟩_exp(-it H) |ψ_j⟩| < ϵ  for all 1 ≤ j ≤ M.
  
  
  The gate complexity achieved by this result exactly reproduces prior results <cit.> regarding the concentration of the average channel formed by random product formulas.
  
  This is not a coincidence, as Theorem <ref> may be viewed as a nontrivial consequence of the existing result. We refer to Section <ref> for a derivation.
  
  
  
  
  This estimation procedure is inspired by <cit.>, and it has a similar flavor.
  The paper <cit.> shows that order log (M) random measurements of a quantum system allow for the accurate estimation of M (arbitrary) observables.
  
  
  In stark contrast, a deterministic/structured product formula must make explicit use of the particular input state and measurement observable to guarantee its accuracy.
  We refer to Section <ref> for an example that illustrates this discrepancy.
  
  
  
  
   §.§ 
  Extension to general products of random unitaries
  
  
  
  So far, we have presented our results exclusively for qDRIFT. But the underlying concentration analysis readily extends to more general random walks on the unitary group.
  Let V_N ⋯ V_1 ∈ U(2^n) be a random product designed to approximate a target unitary U=U_N ⋯ U_1. We need two properties. (i) Causality: for 1 ≤ k ≤ N the random selection of V_k can only depend on previous choices for V_1,…,V_k-1:
  
      Pr[ V_k |V_N… V_k+1 V_k-1,… V_1 ]
      =   Pr[ V_k| V_k-1,… ,V_1 ]
  
  
  
  
  
  
  (ii) accurate approximation: Each realization of V_k (and their conditional expectation) must be close to the ideal unitary U_k. More precisely, let R,b_k >0 be constants such that 
      V_k-_k-1 V_k≤ R   and  𝔼_k-1 V_k - U_k ≤ b_k, 
      where    𝔼_k-1V_k:=𝔼[ V_k| V_k-1,…,V_1 ]
  
  almost surely for all 1 ≤ k≤ N.
  All concentration results from this work, as well as the main result from <cit.>, do readily extend to random products that do obey these two properties.
  
  
  
  Let V=V_N · V_1 ∈ U(2^n) be a random product that approximates a target product U=U_N ⋯ U_1 in a causal (Eq. (<ref>)) and small-step (Eq. (<ref>) fashion.
  Then, the associated unitary channels 𝒱(X) = VXV^† and 𝒰(X) = UXU^† obey
  
  
  
  
  
  
      -_♢   ≤ 2∑_k=1^N a_k    (worst case),
      𝔼-_♢   ≲√(n∑_k=1^N a_k^2)+ 2∑_k=1^N b_k    (typical case),
      
      
      𝔼(ρ) - (ρ) _1
         ≲√(∑_k=1^N a_k^2)+ 2∑_k=1^N b_k
         (fixed input),
      
      
      -𝔼_♢   ≤ 2 ∑_k=1^N b_k    (average case)
      ,
  
  where ≲ suppressed absolute constants. 
  
  
      max_ρ state U ρ U^† - 𝔼 V_N ⋯ V_1 ρ V_1^†⋯ V_N^†_1 ≤    2 N b+ Na^2
         (fixed input & fixed observable, see Rel. (<ref>)), 
      max_ρ state𝔼 U ρ U^† - V_N ⋯ V_1 ρ V_1^†⋯ V_N^†_1
      ≲   √(N)a+ 2Nb
         (fixed input & arbitrary observable), 
      𝔼max_ρ state U ρ U^† - V_N ⋯ V_1 ρ V_1^†⋯ V_N^†_1
      ≲   √(nN)a+ 2Nb    (arbitrary input & arbitrary observable). 
      max_ρ state U ρ U^† - V_N ⋯ V_1 ρ V_1^†⋯ V_N^†_1
      ≤    2Na    (deterministic bound).
  
  
  
  
  
  
  
    §.§.§ Concentration of randomly permuted Suzuki formulas
  
  So far, we have only introduced the first order Lie-Trotter formulas (<ref>). The 2nd
  order formula is typically called the Suzuki formula. For some short time τ >0, define
  
  
      S_2(τ) := ∏^L_j=1exp(-(τ/2) h_j ) ∏^1_j=Lexp(-(τ/2) h_j ).
  
  Higher order formulas, also named after Suzuki, are constructed recursively from S_2 (τ):
  
      S_2p(τ) := S_2p-2(q_pτ)^2 S_2p-2((1-4q_p)τ)S_2p-2(q_pτ)^2,
  
  where q_p:=1/(4-4^1/(2p-1)) <cit.>. We can combine r rounds of 2p-th order Suzuki formulas with time τ = t/r each to approximate a desired quantum evolution. This yields
  
      U = ^- t H =   ^- (t/r) H⋯^- (t/r) H
      ≈    S_2p(t/r) ⋯ S_2p(t/r) = S_2p(t/r)^r
  
  
  
  
  
  For ϵ fixed, we choose an appropriate number of rounds r and obtain a gate count N≈ rL that scales as
  
  
      N_det≳ tΛ L^2 (tΛ L/ϵ)^1/2p  with  Λ = max_j h_j,
  
  simple gates e^- (t/r) h_j to ensure worst-case approximation error 𝒰-𝒱_♢≤ϵ <cit.>.
  Note that the 2p-th order Suzuki formula S_2p(τ) does not specify a preferred order in which the terms in Hamiltonian h_j should appear.
  Childs, Ostrander and Su observed that randomly permuting the order of Hamiltonian terms
  within each block S_2p(t/r) can suppress the approximation error for low order Suzuki formulas <cit.> considerably. More precisely,
  we reshuffle the ordering in each Suzuki block by applying uniformly random permutations π_1,…,π_r to the term indices (h_j ↦ h_π_1 (j) for the first Suzuki block, etc.). Denote the resulting randomized Suzuki formula by V_2p^π_r⋯π_1= S_2p^π_r(t/r) ⋯ S_2p^π_1(t/r).
  Childs, Ostrander and Su proved that
  a total gate count
  
  
  
      N_avg≳     t Λ L^2 (tΛ/ϵ)^1/2p
  
  ensures that the average channel (averaged over all possible permutations) obeys
  -𝔼_π_r,…,π_1_2p^π_r ⋯π_1_♢≤ϵ.
  
  We note in passing that Eq. (<ref>) is tighter than the original bound presented in <cit.>. This slight improvement follows from replacing the original mixing Lemma <cit.>  in the proof of Childs et al. by a stronger statement derived in this work (Lemma <ref> below)
  
  
  
  
  Applying Theorem <ref> to the
  problem at hand supplies strong concentration around this expected behavior.
  
  
  
  
  
  Consider a n-qubit Hamiltonian H=∑_j h_j with Λ = max_j h_j, an associated time-t evolution U=^- t H and a Suzuki order 2p.
  Then, a total gate count of
  
      N_typ≳  tΛ L^2 (nΛ Lt/ϵ^2)^1/(4p+1)+t Λ L^2 (tΛ/ϵ)^1/2p,
  
  ensures that a randomly permuted Suzuki formula with r terms obeys 𝔼_π_1,…,π_r𝒰 - 𝒱_2p^π_r ⋯π_1⋯ S_2p^π_1(t/r) ≤ϵ (typical case). For a fixed (but arbitrary) input state, the gate count can be further reduced to
  
      N_fix≳  tΛ L^2(Λ Lt/ϵ^2)^1/(4p+1)+t Λ L^2 (tΛ/ϵ)^1/2p.
  
  
  
  
  
  
  
  
  
  
       
  
  For simplicity, we omitted the logarithmic dependence on failure probability δ.
  
  
  
  
  In contrast to qDRFIT, the parameters n, ϵ, L parameters now all raised to the 1/(4p+1)st power, and the different randomized settings yield very much the same gate count tΛ^2 L for large p. This is in accordance with observations provided  in <cit.>.
  
  For illustration, we have chosen to directly import older results (<ref>) to compare with the averaged channel bounds (<ref>). However,  (<ref>) can be sharped to scale with commutator <cit.>. Future work remains to carry out commutator analysis for the averaged channel(<ref>) and then swiftly upgrade Corollary <ref>.
  
  We refer to the appendix (Sec. <ref>) for detailed statements and proofs.
  
  
  
    §.§.§ Compiling
  
  Beyond Hamiltonian simulation, random product also help suppressing error in compiling. The task of compiling turns a perfect circuit diagram into a sequence of executable gates. However, since the gate set is discrete (or imperfect), the compiler can only return an approximate circuit
  
      U_N⋯ U_1 - V_N⋯ V_1≤ϵ,
  
  where each V_k will be further decomposed into universal gates (such as using the Solovay-Kiteav Theorem, see e.g. <cit.>) up to some individual accuracy U_k-V_k≤η. In the worst case, the local error U_k-V_k≤η accumulates linearly and we must conclude
  
      U-V≤ Nη,
  
  by means of the triangle inequality. This means that individual accuracy η =ϵ/N is required
  to ensure an overall approximation error of (at most) ϵ.
  This in turn, requires more gates for each approximation, because η-approximating V_k requires a gate count[For the Solovay-Kiteav Theorem, 3≤ c ≤ 4 and it is also known that c=1 is the best possible exponent.] proportional to log^c(η).
  
  Randomized compiling <cit.> addresses precisely this issue.
  It uses random gate synthesis to avoid that individual approximation errors add up “coherently” over the entire compilation. The resulting “incoherent” error accumulation can be much more favorable and the mixing Lemma <cit.> makes this intuition precise. In the following, we present a sharpened version of this statement that seems to be novel.
  
  
  Let U_k be a fixed unitary and V_k a random approximation thereof that obeys U_k-  V_k≤ b, for some b >0. Then, the associated channels obey
  
  
      𝒰_k-𝒱_k _♢≤ 2b.
  
  
  We refer to Appendix <ref>
   for a self-contained proof.
  
   Back to compiling, V_k-U_k can as small as (η^2) <cit.> by mixing appropriate synthesis protocols for V_k. This yields an improved bound for the average channel,
  
      𝒰-𝒱_♢≤(Nη^2),
  
  In words, the individual accuracy needed is suppressed quadratically to η = Ω(√(ϵ/N)).
  In <cit.>, Campbell pointed out that this square root improvement may turn into a multiplicative overhead reduction, depending on the gate synthesis efficiency log^c(√(η))=log^c(η)/2^c.
  Using Theorem <ref>, we can immediately bound the performance of randomized compiling without mixing.
  
  
  Suppose that we wish to approximate U=U_N ⋯ U_1 by a gate collection V=V_N ⋯ V_1 such that each V_k is random and obeys U_k-V_k≤η almost surely. Then,
  
      𝒰 - 𝒱≲√(nN)η.
  
  What is more, the √(n)-factor on the r.h.s. disappears if we restrict attention to a fixed input state.
  
  
  
  
  
  
  
  This translates to individual accuracy η = (ϵ/√(nN)), and η = (ϵ/√(N)) respectively. Both results interpolate between the worst case (<ref>) (“coherent” error accumulation) and best case (<ref>) (“incoherent” error accumulation).
  
  
  
  
  
  
  
   §.§ Numerical experiments
   
  
  
  
  In this section, we perform numerical experiments for simulating a simple Heisenberg model on a one-dimensional chain with a randomly sampled product formula.
  For n qubits, H = 1/n-1∑_i=1^n-1 X_i X_i+1 + Y_i Y_i+1 + Z_i Z_i+1 and we view this as a sum of 3(n-1) simple terms.
  The normalization keeps the interaction strength λ
  =3 as a constant and we consider a constant time evolution t=2.
  
  
  The numerical experiments for the error under various setups using different gate count N and different system sizes n are shown in Figure <ref>.
  We can see that the error ϵ when we consider all input states scales as √(n).
  In contrast, the error ϵ stays roughly the same when we only consider a single input state.
  This is in accordance with the theoretical predictions presented in Theorem <ref> and <ref>.
  
  
  
  
  
  
  
  
  
  
  
  § INSTRUCTIVE EXAMPLES AND PROOF IDEA
  
  
  
  
  
   §.§ Comparison between stochastic averages of product formulas and concrete instances
  
  
  
  This section considers an extremely simple Hamiltonian to pinpoint important differences between averaging random product formulas (that is, Campbell's black box) and concrete instances of product formulas.
  The example Hamiltonian is a 1-local non-interacting Hamiltonian with a Pauli-Z operator acting on each qubit:
  
      H =   1/n∑_k=1^n Z_k,
      
        where
      
      Z_k =   (k-1)-times𝕀⊗⋯⊗𝕀⊗ Z ⊗(n-k)-times𝕀⊗⋯⊗𝕀
  
  for 1 ≤ k ≤ n.
  
  The relevant parameters are L=n (number of terms), λ =1n∑_k=1^n Z_k   = 1 (interaction strength) and we fix the evolution time to t=π.
  
  Stochastic averages of random product formulas can accurately approximate the associated unitary evolution U=exp(- π H) after only a few iterations.
  The following observation is an immediate consequence of Campbell's main result <cit.>, see also Proposition <ref> below.
  
   
  Fix a target accuracy ϵ and set N = t^2 λ^2/ϵ≈ 10/ϵ.
  
  Then, N successive applications of the qDRIFT single-step average 𝒱(ρ)=1n∑_k exp(- iπN Z_k ) ⊗𝕀^(else)ρexp( iπN Z_k) ⊗𝕀^(else) (Campbell's black box)
  approximate the target unitary channel 𝒰(ρ)=Uρ U^† up to accuracy ϵ in diamond distance. In particular, 12𝒱^(N)(|ψ⟩⟨ψ|) - U |ψ⟩⟨ψ| U^†_1 ≤ϵ for all input states |ψ⟩⟨ψ|.
  
  
  This assertion seems remarkably strong. In particular, the sequence length N does not depend on the number of qubits n. Once n is sufficiently large it becomes impossible for concrete product formulas to achieve comparable results. The problem is that the sequence length N is too small to address all n qubits. This necessarily leads to substantial discrepancies between the simulated time evolution V_N ⋯ V_1 and the actual target U, see Figure <ref> for an illustration.
  
  
  
  Assume that n is an even number.
  It is impossible to accurately approximate the time evolution U defined in Eq. (<ref>) with  N < n/2 elementary gates of the form V_i=exp (-πNZ_k(i)) ⊗𝕀^(else).
  More precisely, for each product formula V=V_N ⋯ V_1, there exists an input state |ψ⟩⟨ψ| such that 12V |ψ⟩⟨ψ| V^† - U|ψ⟩⟨ψ| U^†_1 =1.
  
  
  
  
  
  All terms in the Hamiltonian (<ref>) commute. Hence, the associated target evolution factorizes nicely into tensor products: U = exp (-π H) = exp (-πnZ_1) ⊗⋯⊗exp (-πnZ_n). Up to a global phase, each single-qubit unitary affects the computational basis in the following fashion: exp( -πnZ) |0 ⟩ = |0 ⟩ and exp (-πn Z) |1 ⟩ = exp (2πn) |1 ⟩.
  These small phase shifts can add up for states that are in superposition. Consider the tensor product of a GHZ state on n/2 qubits with the all-zeroes state on the remaining half: |ψ̃⟩ = 1√(2)(|0 ⟩^⊗ n/2 + |1 ⟩^⊗ n/2) ⊗ |0 ⟩^⊗ n/2.
  Then,
  
      U |ψ̃⟩ =   exp(-2πn Z)^⊗ n|ψ̃⟩
      
      =   1√(2)( |0 ⟩ + (exp (2πn))^n/2 |1 ⟩) ⊗ |0 ⟩^⊗ n/2
      
      =   1√(2)(|0 ⟩^⊗ n/2- |1 ⟩^⊗ n/2) ⊗ |0 ⟩^⊗ n/2
  
  and we can easily check that input and output are orthogonal to each other: 12U |ψ̃⟩⟨ψ̃| U^† - |ψ̃⟩⟨ψ̃| _1 = 1.
  These features do not change if we permute the qubits in the input state |ψ̃⟩. Any combination of a GHZ state on one half of the qubits with computational |0⟩-states on the remaining ones obeys the same orthogonality relation. We can use this freedom to construct a worst-case input |ψ⟩ for a fixed product formula V=V_N ⋯ V_1 comprised of fewer than n/2 single-qubit gates. Simply initialize the (at most) n/2 qubits on which the product formula acts nontrivially in the computational 0-state and hide the GHZ component among the remaining qubits. By construction, the product formula V does not affect this input state at all.
  This is a worst case, because the target unitary U does rotate the hidden GHZ component into an orthogonal configuration:
  U |ψ⟩⟨ψ|U^† - V|ψ⟩⟨ψ| V^†_1 =  U |ψ⟩⟨ψ| U^† - |ψ⟩⟨ψ| _1 =1.
  
  
  This negative statement highlights that the gate count of (worst case) accurate product formulas must in general depend on the number of qubits and justifies the appearance of n in Theorem <ref>. Note, however, that Lemma <ref> is contingent on identifying a worst-case input state for a fixed (and known) product formula. If the input state is fixed, the situation can change dramatically. For instance, we could use explicit knowledge of the input to construct a product formula that accurately approximates its time evolution. Identifying an optimal product formula seems like a daunting task, but randomness can help.
  Theorem <ref> asserts that a collection of N ≳π^2/ϵ^2 randomly selected single-qubit gates approximate the time evolution (<ref>) of any fixed input state |ψ⟩ up to accuracy ϵ in trace distance. While this gate count is considerably larger than the one put forth in Corollary <ref>, it is still independent of the number of qubits. What is more, this assertion applies with high probability to any fixed input state.
  This capitalizes on another advantage of generating unstructured product formulas according to a randomized procedure: it is extremely difficult to fool a randomized compiling procedure with an already fixed 
   input.
  
  
  
  
  
  
  
   §.§ Instructive concentration argument for a simple Hamiltonian
  
  
  
   
  
  This section provides intuition for the concentration effects that ultimately imply Theorem <ref>
  by means of another example Hamiltonian that is composed of (commuting) Pauli-Z terms only:
  
  
      H =   1/2^n∑_p∈{0, 1}^nα_p Z_p  where
       Z_p=   Z_(p_1,…,p_n)=Z^p_1⊗⋯⊗ Z^p_n
  
  (with the convention that Z^0 =𝕀) and α_p∈{-1,1}.
  That is, the Hamiltonian is a sum of 2^n signed Pauli strings that are comprised of Z and 𝕀, as well as a global sign. A high-order Suzuki formula would require a gate complexity of 𝒪(L) = 𝒪(2^n) by putting down every term. In contrast, by random sampling, Theorem <ref> yields a gate complexity of 𝒪(n / ϵ^2) (Theorem <ref> yields (1/ϵ^2)). This is an exponential improvement in system size.
  
  
  
  The physical intuition is that all the terms in the Hamiltonian act on the same system with n qubits (a 2^n-dimensional Hilbert space), so their actions must “overlap” with one another, and can be efficiently estimated by random sampling.
  To see this effect more clearly, let us write down the unitary evolution exp(-i H) in the computational basis |b⟩ with multi-index b=(b_1,…,b_n) ∈{0, 1}^n.
  Note that all terms in the Hamiltonian (<ref>) are diagonal in the computational basis. This implies
  
      exp(-i H) |b⟩ =   exp(-i1/2^n∑_p∈{0, 1}^nα_p Z_p) |b⟩
      
      =   exp(-i1/2^n∑_p∈{0, 1}^n c_p(b)) |b⟩
      
      :=   e^-i S(b)|b⟩,
  
  where c_p(b) = α_p⟨b| Z_p|b⟩∈{-1, 1}. When we instead sub-sample an effective Hamiltonian Ĥ comprised of N randomly selected terms, the constructed product formula would be
  
      exp(-iĤ)   = exp(-i1/Nα_p_N Z_p_N) ⋯exp(-i1/Nα_p_1 Z_p_1) |b⟩
      
          =   exp(-i1/N∑_k c_p_k(b)) |b⟩ := e^-iŜ(b)|b⟩.
  
  By Hoeffding's inequality, Ŝ(b) = 1/N∑_k c_p_k(b) should concentrate around S(b) = 2^-n∑_p∈{0, 1}^n c_p(b) with standard deviation 1 / √(N) and an exponentially decaying tail (think Monte Carlo). An illustration and some helpful facts can be found in Figure <ref>.
  
  Let us now fix an (arbitrary) n-qubit basis state. Then, the probability of an ϵ-deviation (or larger), i.e.  |Ŝ(b) - S(b)| > ϵ,
  can be bounded by 1/e if we choose N = 1 / ϵ^2. However, if we wish the effective Hamiltonian to work for all basis states simultaneously, we would choose a larger N = n/ ϵ^2 to ensure that the deviation probability is further suppressed exponentially to 1 / e^n. By a union bound, |Ŝ(b) - S(b)| ≤ϵ for all 2^n computational basis states |b⟩ with probability at least 1 - 2^n / e^n. Albeit in the simplest example (commuting Hamiltonian), this demonstrates that a random product formula exp(-iĤ) can accurately simulate exp(-i H) up to error ϵ with only N = n / ϵ^2 gates, which further reduces to (1/ϵ^2) for an arbitrary basis state. The powerful tool of concentration for matrix (vector) martingales allows us to prove the same statement for any (non-commuting) many-body Hamiltonian.
  
  We will return to this example Hamiltonian in Section <ref> to show that this more general analysis yields an essentially optimal parameter dependence:
  dimension dependence that is tight: the scaling N ≥Ω( n t^2 λ^2/ϵ^2) in Theorem <ref> is unavoidable in general.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Proof idea for Theorem <ref> and <ref>
  
  
  
  This section sketches the main ideas and tools required to establish Theorem <ref> and Theorem <ref>. The other results follow from more elementary
  arguments. Detailed arguments and rigorous statements are provided in Appendix <ref> below.
  
  Consider an n-qubit Hamiltonian H = ∑_j=1^L h_j and an evolution time t.  The associated unitary evolution defines a (unitary) channel on n-qubit states:
  
      𝒰(| ψ⟩⟨ψ| ) =    U | ψ⟩⟨ψ| U^†  where
      
      U =   exp( -i t H )
      = exp( - i t ∑_j=1^L h_j ).
  
  Fix a number of steps N and set λ = ∑_j=1^L h_j. The task is to accurately approximate the target unitary U by a product formula, i.e., the composition of N simple unitary evolutions:
  
      𝒱^(N) (|ψ⟩⟨ψ|)  =   𝒱_N∘⋯∘𝒱_1 (| ψ⟩⟨ψ|) 
      
      =    V_N ⋯ V_1 | ψ⟩⟨ψ| V_1^†⋯ V_N^†.
  
  
  We quantify the difference between ^(N) and 𝒰 in diamond distance. That is, the worst case approximation error over all possible input states ρ in the presence of an unaffected quantum memory.
  
  Let ℰ,ℱ be two quantum channels, and let ℐ(|ψ⟩⟨ψ|) = |ψ⟩⟨ψ| denote the identity channel on an equally large ancilla system. The diamond distance between ℰ and ℱ is defined as
  
      12ℰ-ℱ_♢ = 12max_|ψ⟩⟨ψ|ℰ⊗ℐ(|ψ⟩⟨ψ|) - ℱ⊗ℐ (|ψ⟩⟨ψ|) _1,
  
  where the maximization ranges over all pure[Convexity ensures that the worst-case discrepancy is attained at a pure state |ψ⟩⟨ψ|. Hence, it is not necessary to consider mixed states ρ in this definition. We refer to <cit.> for details.]
  input states | ψ⟩⟨ψ| and
  ·_1 denotes the trace norm.
  First, we relate the diamond distance (<ref>) between the channels ^(N) and 𝒰, which are both unitary, to an operator norm distance of the associated matrices:
  
      12𝒱^(N)-𝒰_♢
      =   
      12max_|ψ⟩⟨ψ|𝒰(|ψ⟩⟨ψ|) - 𝒱^(N) (|ψ⟩⟨ψ|) _1
      
      
      ≤    V_N ⋯ V_1 - U ,
  
  
  see Lemma  <ref> below.
  This relation exploits the fact that stabilization
   (i.e., tensoring with the identity channel) is not
  necessary for computing the diamond distance of two unitary channels <cit.>.
  
  
  
  
  
  
  Now, we can deal with the i.i.d. random matrices V_N,…,V_1 in the more familiar operator norm. Add and subtract the expected product 𝔼[ V_N ⋯ V_1 ] = 𝔼[V_N ] ⋯𝔼[ V_1 ]= ( 𝔼V )^N to decompose the operator-norm difference into two qualitatively different contributions:
  
      V_N ⋯ V_1 - U 
      ≤   deterministic bias (𝔼 V )^N - U 
      +
      random fluctuation V_N ⋯ V_1 - 𝔼[ V_N ⋯ V_1 ] .
  
  
  These two contributions can be analyzed separately:
  
  
    i. Deterministic bias: Most product formulas arise from first decomposing the target unitary into a sequence of many small steps: U=(U^1/N)^N, where U^1/N = exp( - i (t/N) H ) is close to the identity matrix.
  This allows for approximating U^1/N by another process that is easier to implement.
  The random importance sampling model (<ref>) over individual Hamiltonian terms
  
  is designed to achieve this goal. The average approximation error scales inverse quadratically in the number of steps:
  (𝔼 V ) - U^1/N≤ t^2 λ^2 / N^2;
  see Lemma <ref> below.
  While small, this expected error does constitute a bias that is present in each of the N approximation steps. It can, and in general will, accumulate across different time steps:
  
  
      𝔼[V_N ⋯ V_1 ] - U   =    (𝔼 V )^N - (U^1/N)^N 
      ≤    N  (𝔼V)-U^1/N
      ≤   t^2 λ^2/N,
  
  see Lemma <ref> below.
  The first inequality is obtained from a telescoping sum. This upper bound diminishes as the number of steps N increases. For ϵ >0,
  
      N ≥2 t^2 λ^2/ϵ  ensures  
      (𝔼V)^N-U ≤ϵ/2.
  
  
  
    ii. Random fluctuation:
  We also need to control the deviation of a product of i.i.d. unitaries V_N ⋯ V_1 from its expectation 𝔼[V_N ⋯ V_1]= (𝔼 V)^N in operator norm. We achieve this by applying modern matrix martingale tools, which may of independent interest. In short (see  Sec.<ref> for a more thorough introduction), a martingale is a random process {B_k:  k=0,…,N} such that the distribution of B_k only depends on past elements B_k-1,…,B_1 (`causality') and also 𝔼[B_k|B_k-1⋯ B_1]=B_k-1 (`on average, tomorrow is the same as today').
  To control fluctuations, we introduce a martingale
  
  
  
  
  
  that interpolates between the extreme cases we need to compare:
  
      B_k =    (𝔼V)^N-kV_k ⋯ V_1   such that
      
      B_0 =    (𝔼V)^N   and  
       B_N= V_N ⋯ V_1.
  
  Note that adjacent elements of this process only differ in a single term: B_k arises from B_k-1 by replacing 𝔼V at position k (counted from the right) by a random realization V_k of V, and thus _k-1 B_k=B_k-1.  Moreover, 
  the current iterate B_k only depends on realizations in the past and the interpolation process {B_k: k=0,…,N} forms a martingale comprised of d × d matrices, where d=2^n is the Hilbert space dimension. The discrepancies are small: V_k - (𝔼 V)≤ 2t λ/N, because each realization of V is also very close to the identity matrix.
  
  
  
  
  
  
  
  
  Powerful tail bounds for matrix-valued martingales (which are relatively modern compared to their scalar counterparts) are available in the literature <cit.>.
  Adapting these results to the task at hand yields the bound
  
      Pr[  V_N ⋯ V_1 - (𝔼V)^N ≥ϵ/2 ] 
         ≤ 2d exp( - Nϵ^2/44t^2 λ^2),
  
  see Proposition <ref> below.
  In words, the product V_N ⋯ V_1 will concentrate around its expectation once N is sufficiently large. 
  Similar to more conventional random walks on integer lattices, the error is subgaussian with variance proportional to N · (λ^2 t^2 / N^2).  There is an extra dimensional factor d=2^n that arises because the martingale is matrix-valued. It converts to the number of qubits n=log(d) in the gate count N. For error parameters ϵ,δ∈ (0,1), 
  a gate count obeying
  
      N ≥ 44  t^2 λ^2/ϵ^2log (2d/δ)   implies
          V_N ⋯ V_1 - (𝔼 V)^N ≤ϵ2 with probability > 1-δ.
  
  
  
  Theorem <ref> can be derived by combining the bound (<ref>) for the deterministic bias with the bound (<ref>)  for random fluctuations. This provides a high probability error bound in the diamond distance when N ≥Ω(n t^2 λ^2 / ϵ^2).
  
  Corollary <ref> is derived by integrating the tail bounds of Theorem <ref>.
  This produces
  
  
  
      𝔼V_N ⋯ V_1 - (𝔼 V)^N ≲√(t^2 λ^2/Nlog_2 (d)),
  
  where the symbol ≲ suppresses a modest multiplicative constant. We refer to Section <ref> for details.
  
  With fixed input(Theorem <ref>), we
  use convexity to restrict our attention to a fixed, pure input state |ψ⟩.
  We then need to compare U|ψ⟩ with V_N ⋯ V_1 |ψ⟩ which can be achieved by constructing an interpolating random process that describes a vector-valued martingale.
  We then call another concentration inequality (Lemma <ref>)
  
      Pr[ (V_N ⋯ V_1 - 𝔼[V]^N)| ψ⟩_ℓ_2
          > ϵ]
          ≤exp( -ϵ^2 N/8e t^2 λ^2),
  
  which, in contrast, does not contain a dimensional factor. This is the reason why Theorem <ref> and Corollary <ref> do not depend on system size at all.
  
  
  
  
  § DISCUSSION AND OUTLOOK
  
  
  
  This work shows interesting characteristics of randomization that might help to further improve quantum simulation. (a)
  By studying typical unitary instances of qDRFIT, we have shown that L-independence of qDRIFT attributes to randomly sampling terms; mixing different realizations is not essential.
  (b)
  Gate complexities can be reduced substantially by restricting attention to a particular input state or/and target observable.
  Simple randomized compilation procedures, like qDRIFT, do not utilize extra information. We believe that more specialized algorithms might be able to exploit additional structure.
  (c) We have shown that channel averages can be much closer to the ideal evolution than any individual product formula, however, in the case of qDRIFT it is only a quadratic improvement 1/ϵ^2→ 1/ϵ. This can be seen as a manifestation of Jensen's inequality for convex functions (distance to ideal evolution) and we may also see such behavior in other quantum simulation algorithms.
  Up to our knowledge, we also provide the first matrix concentration analysis for a randomly sampled product formulas and – more generally – Hamiltonian simulation. Similar proof techniques readily apply to any random product formula, e.g. randomly permuted Suzuki formulas <cit.>, and in fact any (causal) random unitary product, e.g. randomized compiling <cit.>.
  
  Beyond random products, we expect the developed matrix concentration tools to be useful for controlling stochastic errors in other quantum computing applications, as well as analyzing properties of random Hamiltonians <cit.>.
  
  
   §.§ Acknowledgments:
   The authors want to thank John Preskill and Yuan Su for valuable inputs and inspiring discussions. Earl Campbell and Nathan Wiebe provided insightful comments, as well as encouraging feedback. Finally, we would like to thank the anonymous reviewers for in-depth
  comments and suggestions.
  CC is thankful for Physics TA Relief Fellowship at Caltech.
  HH is supported by the Kortschak Scholars Program.
  RK acknowledges funding from ONR Award N00014-17-1-2146 and ARO Award W911NF121054).
  JAT gratefully acknowledges funding from the ONR Awards N00014-17-1-2146 and N00014-18-1-2363
  and from NSF Award 1952777.
  
  
  
  
  
  
  
  
  
  
  § MATRIX AND VECTOR VALUED MARTINGALES.
  
  
  Let X_1,…,X_N be independent and identically distributed (i.i.d.) random variables.
  Then, the strong law of large numbers (LLN) implies that the sample mean μ̂ =(1/N)∑_k=1^N X_k converges to the actual mean μ = 𝔼X_k (μ̂≈μ as N →∞).  In fact, even for finite N, the sample average is likely to be close to the expectation.  This behavior is called concentration. It turns out that concentration is a surprisingly generic phenomenon, and it kicks in earlier than one might expect.
  Asymptotically large sample sizes (N →∞), which are essential for the law of large numbers and the central limit theorem, are not required establish concentration. An example is Bernstein's inequality for sums of bounded random numbers.
  [Bernstein's inequality]
  Let X_1,…,X_N be independent random variables that obey X_i=0 and |X_i| ≤ R almost surely. Then, for ϵ >0,
  
      Pr[ | ∑^N_k=1 X_k | ≥ϵ] ≤ 2exp(ϵ^2/2/v^2+Rϵ/3)   where   v=∑^N_k=1 X_k^2.
  
  
  
  Bernstein's inequality equips random sums with a LLN-type concentration behavior already for non-asymptotic sample sizes N ≳max{v^2/ϵ^2,R/ϵ}. However, it requires strong assumptions on the random variables involved. They need to be independent, bounded scalars.
  In fact, random processes may concentrate under much weaker conditions.
  
  Martingales form a richer family of processes that capture more realistic random processes and that still enjoy powerful concentration inequalities. Let us offer a minimal technical introduction following Tropp <cit.>.  Consider a filtration of the master sigma algebra _0⊂_1 ⊂_2 ⋯⊂_k ⊂⋯, where we denote the conditional expectation with respect to _k by the symbol _k. A martingale is a sequence { B_0, B_1, B_2, …} of random variables that satisfies
  
      σ(B_k)    ⊂_k    (causality), 
      _k-1 B_k    := [ B_k | B_k-1,…,B_0 ]= B_k-1   (status quo).
  
  Heuristically, we can think of k as a `time' index, and _k contains all possible events that are determined by the history up to time k. The present, aka B_k
  may depend on the past (i.e., B_0,…,B_k-1), but it cannot depend on the future (`causality').
  The second condition formalizes the requirement that, on average, tomorrow (B_k) is the same as today (B_k-1).
  
  A martingale sequence may be composed of scalars, e.g., a one-dimensional random walk, but we can also consider vector- or matrix-valued martingales. Analyzing concentration for vector- and matrix-valued martingales is an ongoing field of research; for example, see <cit.>,
  but many powerful concentration inequalities already exist.  In this work, we will use the matrix generalization of Freedman's inequality (due to one of the authors). Let 𝕄_d × d denote the space of complex-valued d × d matrices.
  
  [Matrix Freedman <cit.>] 
  Let {B_k:k=0,…,ℓ, ...}⊂𝕄_d × d be a  matrix martingale.
  Assume that the associated difference sequence C_k = B_k - B_k-1
  obeys C_k≤ R almost surely.  Define the random variable
  
      w_ℓ :=max( ∑_k=1^ℓ_k-1 C_k^† C_k
      , ∑_k=1^ℓ_k-1 C_k C_k^†).
  
  Then, for any τ >0
  
      Pr[ ∃ℓ : B_ℓ - B_0 ≥τ, w_ℓ≤ v ]
      ≤ 2 d exp( -τ^2/2/v +Rτ/3).
  
  
  
  
  
  
  
  
  
  
  
  This bound exponentially suppresses the probability of the following undesirable event:  the conditional variance w is small while the actual deviation is large. The intricacy of this bound is that the conditional variance w_ℓ is itself a random variable. However, for this work we will use a weaker but more transparent version.
  
   
  Let {B_k:k=0,…,N}⊂𝕄_d × d be a matrix martingale.
  Assume that the associated difference sequence C_k = B_k - B_k-1
  obeys C_k≤ R and its conditional variance obeys ∑_k=1^N _k-1 C_k C_k^†≤ v almost surely.
  Then
  
      Pr[  B_N - B_0 ≥τ]
      ≤ 2 d exp( -τ^2/2/v +Rτ/3)
        for any τ >0.
  
  
  
  
  
  
  
  
  
  
  To arrive at this, ignore the events for ℓ<N and use that w_ℓ≤ v holds almost surely.
  This simplified Matrix Freedman inequality closely resembles Bernstein's inequality. Actually, Fact <ref> can be viewed as a special case of Corollary <ref> where d=1 and B_k = ∑_k'=0^k X_k.
  But for matrix-valued martingales, the tail bound must depend on the dimension d.
  
  Recapitulating the proof of the general Matrix Freedman inequality would go beyond the scope of this work. It makes full use of Lieb's concavity theorem, stopping times, and Burkholder functions.  There is a slightly weaker result that follows from the Golden–Thompson inequality <cit.>; see also <cit.>.  
  
  
  
  
  
  Similar concentration inequalities are valid for martingales on the complex vector space ℂ^d.   Remarkably, these results are independent of the ambient dimension d.
   
  Let {B_k:k=0,…,N}⊂ℂ^d be a vector martingale taking valules in the inner product space ·_ℓ_2.  Assume that the associated difference sequence C_k = B_k - B_k-1
  obeys ∑_k-1C_k _ℓ_2^2 ≤ v and C_k_ℓ_2≤ R almost surely. Then, for any τ >0
  
      Pr[ B_N - B_0 _ℓ_2≥τ]
      ≤ 2 exp( -τ^2/2/v+Rτ/3).
  
  
  
  
  
  
  
  
  
  
  
  This is simplified version of <cit.>. To prove the above version, start with <cit.> for general martingales on Banach spaces. In our case, vectors are equipped with the standard ℓ_2 norm, set the smoothness constant D=1 (in the notation of <cit.>), and optimize for a λ like in Bernstein's inequality.
  
  Our concentration analysis of random product formulas with fixed inputs relies on another tool: uniform smoothness of the underlying space <cit.>.
  
  Given a vector martingale, uniform smoothness supplies a recursive/local bound on moment growth. Decomposing B_k=C_k+B_k-1, all we need is [C_k|B_k-1]=0, to apply the following result.
  
   
  Let x,y ∈ℂ^d be two random vectors that obey 𝔼[y | x ] =0. When q ≥ 2
  
      ( 𝔼x+y _ℓ_2^q)^2 / q≤( 𝔼x _ℓ_2^q )^2 / q + (q-1)
      ( 𝔼y _ℓ_2^q )^2 / q
  
  
  
  
  A slightly weaker version of this classic fact follows from a short argument.
  
  
  Start with Bonami's inequality <cit.> for normed vector spaces (denoting ·_ℓ_2 as ·_2):
  
      ( x+y _2^q +  x-y _2^q/2)^2/q≤1/2( x+√(q-1)y_2^2+ x-√(q-1)y_2^2) =  x_2^2+(q-1) y_2^2.
  
  
  This formula can be converted to the desired statement (Proposition <ref>) via elementary manipulations. We follow  <cit.>: take expectation and use triangle inequality for L_q/2 norm
  
      x+y _2^q +  x-y _2^q/2≤( x_2^2+(q-1) y_2^2)^q/2≤(( x_2^q)^2/q+(q-1)( y_2^q)^2/q)^q/2.
  
  Next, following <cit.>, by Jensen's in the first and Lyapunov's in the second  inequality,
  
      ( x+y _2^q)^2/q + ( x_2^q)^2/q/2   ≤( x+y _2^q)^2/q + ( x-y_2^q)^2/q/2
         ≤( x+y _2^q + x-y_2^q/2)^2/q≤
           ( x_2^q)^2/q+(q-1)( y_2^q)^2/q.
  
  By rearranging terms, we obtain the result with the constant 2(q-1), which is off by a factor of 2. The advertised constant (Proposition <ref>) can be obtained via a more involved trick <cit.>.
  Geometrically, this result expresses the uniform smoothness of the space (·_ℓ_2^q)^1/q.
  
  
  We refer to <cit.> for further exposition on uniform smoothness for general Schatten p-norm. Recently, this method has been applied to dynamic properties of random Hamiltonians <cit.>.
  For the task at hand, we can use Markov's inequality to convert such bounds on moment growth into a strong tail bound, similar to Fact <ref>. This is the content of Lemma <ref> below.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § TECHNICAL DETAILS AND PROOFS
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ 
  Proof of Theorem <ref> and Corollary <ref>: Approximation error under the worst-possible input
  
  
  
  
  The proofs of Theorem <ref> and Corollary <ref> were sketched in Section <ref>.
  This section contains the details. In Section <ref>, we first relate the diamond distance to the operator norm. This allows us to work with the operator norm, which is mathematically simpler. Then we bound the two error contributions arising from the deterministic bias (in Section <ref>), as well as random fluctuations (in Section <ref>). Finally, we combine the two
  bounds to obtain a convergence guarantee for randomly sampled product formulas. This is the content of Section <ref>.
  
  
  
    §.§.§ Conversion from diamond distance into operator norm
  
  
  
  The diamond distance is a rather intricate object. Although it can be phrased implicitly as a semidefinite program, analytical formulas are rare and far between. It is, however, sometimes possible to relate diamond distances to other figures of merit that are easier to access,
  The mixing lemma by Campbell <cit.> and Hastings <cit.> is one very useful example.
  
  
  
  Let U be a fixed unitary and V be a random approximation thereof that obeys
  
  V-U ≤ a almost surely. Then, the associated channels 𝒰(X)=UXU^† and 𝒱(X) = V X V^† obey
  
      12𝒱 -𝒰_♢≤ V - U +12a^2.
  
  
  
  More recent insights, in particular <cit.>, allow for sharpening this bound. In particular, we can completely remove the a^2/2-term on the r.h.s.
  
   
  Let 𝒰(ρ) = U ρ U^† and 𝒱(ρ) =V ρ V^† be unitary channels. Then,
  12𝒰(|ψ⟩⟨ψ|) - 𝒱 (|ψ⟩⟨ψ|) _1 ≤ (U-V) | ψ⟩_ℓ_2 for any pure state |ψ⟩.
  In turn, 12𝒰 - 𝒱_♢≤U - V.
  The latter relation generalizes to
  averages of random unitary channels:
  
  
      12𝒰 - 𝔼[𝒱] _♢≤U - 𝔼[V] .
  
  
  
  The first claim follows from the fact that stabilization is not required to compute the diamond distance between two unitary channels <cit.>. The second claim improves the mixing lemma.
  
  
  Fix an input |ψ⟩ and
  denote the output state vectors by |u ⟩=U| ψ⟩ and |v ⟩ = V | ψ⟩, respectively.
  Normalization ensures that
  these state vectors obey |⟨ u,v ⟩ | ≤ 1.
  
  Apply the Fuchs–van de Graaf relations <cit.> to convert the output trace distance into a (pure) output fidelity and conclude
  
      12|u ⟩⟨ u| - |v ⟩⟨ v| _1 = √(1-| ⟨ u,v ⟩|^2)
      = √((1+| ⟨ u,v ⟩|) (1-| ⟨ u,v⟩|))≤√(2(1-Re(⟨ u,v ⟩))
      = |u ⟩ - |v ⟩_ℓ_2.
  
  The first diamond distance bound then is a direct consequence of this relation. Use the fact that stabilization is not necessary for computing the diamond distance of two unitary channels <cit.> to conclude
  
      12𝒰 - 𝒱_♢
      = max_|ψ⟩⟨ψ|12𝒰(|ψ⟩⟨ψ|) - 𝒱(|ψ⟩⟨ψ| ) _1
      ≤max_| ψ⟩(U-V)| ψ⟩_ℓ_2=U-V  .
  
  Here, we have also used the definition of the operator norm.
  In order to handle expectation values, we need an additional argument. Let (p_k,V_k) be an ensemble of unitaries with weights p_k≥ 0 that obey ∑_k p_k =1. Then, Cauchy–Schwarz asserts
  
      | ⟨ψ |U^†𝔼[V] | ψ⟩|^2 = | ∑_k √(p)_k √(p)_k ⟨ψ |U^† V_k | ψ⟩ |^2
      ≤(∑_k p_k) ∑_k p_k | ⟨ψ | U^† V_k | ψ⟩ |^2
      =∑_k p_k | ⟨ψ | U^† V_k | ψ⟩ |^2
  
  for any unitary U and state |ψ⟩. Combined with Fuchs–van de Graaf, this observation delivers
  
      12𝒰(|ψ⟩⟨ψ|) - 𝔼[𝒱(|ψ⟩⟨ψ|)]_1
      ≤( 1 - ∑_k p_k | ⟨ψ |U^† V_k | ψ⟩|^2 )^1/2≤(1- | ⟨ψ|U^†𝔼[V]|ψ⟩|^2)^1/2
  
  for any pure input state |ψ⟩.
  This is enough to conclude 12 U|ψ⟩⟨ψ| U^†- 𝔼[V|ψ⟩⟨ψ|V^†]_1 ≤(U-𝔼[V])|ψ⟩_ℓ_2, much as before.
  We emphasize that this relation is true for any fixed unitary U and any unitary ensemble (p_k,V_k). This flexibility is essential to deduce the diamond distance bound, because 𝔼[𝒱] is not unitary and stabilization must be taken into account:
  
      12𝒰 - 𝔼 [𝒱]_♢   = max_|ψ⟩⟨ψ|12𝒰⊗ℐ (|ψ⟩⟨ψ|) - 𝔼[ 𝒱⊗ℐ (|ψ⟩⟨ψ|)] _1 
         = max_|ψ⟩⟨ψ|12 (U ⊗𝕀) |ψ⟩⟨ψ| (U ⊗𝕀)^† - 𝔼[ (V ⊗𝕀)|ψ⟩⟨ψ| (V ⊗𝕀)^†] _1 
         ≤max_|ψ⟩(U ⊗𝕀 - 𝔼[V ⊗𝕀])|ψ⟩_ℓ_2=  (U - 𝔼[V]) ⊗𝕀  = U - 𝔼[V] .
  
  This is what we needed to show.
  
  
  
  
    §.§.§ Controlling the deterministic bias
   
  
  Next, we establish a bound on the deterministic bias between the averaged channel
  and the ideal unitary evolution.
  
   
  Consider the i.i.d. unitary product constructed by the qDRIFT protocol (<ref>) for
  simulating U=exp (-i t H) for evolution time t.  Define the total strength λ= ∑_jh_j. Then
  
      U - 𝔼 [V_N ⋯ V_1] ≤t^2 λ^2/N.
  
  
  
  Note that the improved mixing lemma, Lemma <ref> above, allows for converting this statement into a diamond distance bound for the associated channels:
  
      12𝒰-𝔼 [𝒱_N ∘⋯∘𝒱_1]_♢≤U - 𝔼[V_N ⋯ V_1 ] ≤t^2 λ^2/N.
  
  
  This is a slight improvement over the main existing technical result regarding qDRIFT <cit.>.
  Indeed, Campbell labels the total average qDRIFT channel ℰ = 𝔼[𝒱_N ∘⋯∘𝒱_1],
  and he establishes that 12𝒰-ℰ_♢≤ (t^2λ^2/N)e^2 tλ/N in <cit.>.
  Both assertions become very similar in the large N limit, but (<ref>) is
  always tighter and the discrepancy can be quite pronounced for small and intermediate values of N.
  
  The proof of Proposition <ref> is based on an extension of the numerical bounds |e^ix-1| ≤ |x| and
  |e^ix-ix - 1 | ≤ x^2/2 for all x ∈ℝ to Hermitian matrices.
  
   
  Let X be Hermitian.  Then we have the zero-order bound
  exp (iX)-𝕀≤X and
  the first-order bound
  exp (iX)-iX-𝕀≤12X^2.
  
  
  These observations can be converted into accurate operator-norm bounds for the expected error of individual qDRIFT steps.
  
   
  Fix a Hamiltonian H=∑_j=1^L h_j and parameters N, t.
  Set U^1/N=exp( - i (t/N) H ) and λ = ∑_j=1^L h_j.
  Then, the random matrix V defined in (<ref>) obeys
  
      V - (𝔼 V)≤2 tλ/N (almost surely)  and   (𝔼 V) - U^1/N≤t^2 λ^2/N^2
  
  
  
  Streamline the notation from Figure <ref> by absorbing the scaling factor (t/N) into the random Hermitian matrix X.  In particular,
  V=exp(-iX), 𝔼V=𝔼[exp(-iX)], U^1/N=exp(-i𝔼[X]) and X=(tλ)/N almost surely.
  Observe that
  
      V-(𝔼V)≤exp(-iX)-𝕀+ 𝕀-𝔼[exp (-iX)]≤exp(-iX) - 𝕀 + 𝔼𝕀-exp(-iX) ,
  
  where the last inequality is Jensen's.
  Fact <ref> and uniform normalization 
   then imply exp(-iX)-𝕀≤X= (tλ)/N for any instance of the random matrix X.
  This uniform bound also covers the expected norm difference and we conclude V-(𝔼V) ≤ 2 t λ/N.
  The (tighter) second claim can be derived in a similar fashion.
  A combination of Jensen's inequality, Fact <ref>, and uniform normalization delivers
  
      (𝔼V)-U^1/N   = 𝔼[ exp (-iX)]-𝕀 +iX ] + ( 𝕀-i𝔼[X]-exp(-i𝔼[X]) ) 
         ≤𝔼exp (-iX) -𝕀 +iX 
      + exp (-i𝔼[X])-𝕀 + i𝔼[X]
         ≤12𝔼X ^2 + 12𝔼[X]^2
      ≤𝔼X ^2
      = ( t λ/ N )^2.
  
  This is the advertised result.
  
  
  We also need a statement regarding error accumulation over several applications of similar, but not identical, linear operators.
  It is a rather intuitive consequence of operator norm sub-multiplicativity and the triangle inequality.
  See <cit.> for related results.
  
   
  Let 𝔼V and U^1/N be matrices with bounded operator norm, i.e. 𝔼V≤ 1 and U^1/N≤ 1. Then
  
      (𝔼V)^N-(U^1/N)^N ≤ N (𝔼V)-U^1/N.
  
  
  
  
  The triangle inequality and sub-multiplicativity imply
  
      A_1 A_2 - B_1 B_2    =  (A_1 - B_1) A_2 + B_1 (A_2-B_2) ≤A_2 A_1 - B_1  + B_1 A_2 - B_2
  
  for any matrix quadruple A_1,A_2,B_1,B_2 with compatible dimensions.
  Use the assumed operator norm bounds to iteratively apply this relation and deduce the statement:
  
      (𝔼V)^N-(U^1/N)^N    =  (𝔼 V) (𝔼V)^N-1- U^1/N(U^1/N)^N-1
         ≤ (𝔼V)-U^1/N +  (𝔼V)^N-1- (U^1/N)^N-1≤⋯≤  N  (𝔼V)-U^1/N.
  
  This is the stated result.
  
  
  
  The main result of this section immediately follows from combining Lemma <ref> and Lemma <ref>.
  Decompose U=exp (- tH) into N steps U^1/N=exp (-i (t / N) H) and conclude
  
      U - (𝔼 V)^N 
      = (U^1/N)^N - (𝔼V)^N ≤ N  U^1/N-(𝔼 V)≤t^2 λ^2/N.
  
  This is what we had to show.
  
  
  
  
    §.§.§ Controlling random fluctuations
   
  
  In the previous subsection we have essentially recapitulated the state of the art regarding qDRIFT: the algorithm provides an accurate approximation in expectation over all possible random choices (deterministic bias).
  In this section, things start to get interesting.
  We want to show that a single realization of qDRIFT is likely to provide a good approximation, provided that the number of steps N is sufficiently large. In order to achieve this goal, we need to show that concrete fluctuations around the (accurate) expected behavior remain small:
  
      V_N ⋯ V_1 ≈𝔼[V_N ⋯ V_1] =(𝔼V)^N  with high probability.
  
  In words, we need to show that a product of i.i.d. random matrices concentrates sharply around its expectation value.
  This is an interesting and nontrivial problem, even in the (asymptotic) large N limit. While sharp concentration bounds for sums of i.i.d. random matrices have been available for more than a decade now <cit.>, our understanding of concentration for random matrix products is more limited; see <cit.> and references therein. There is a lot of math literature on random walks on Lie groups, but the focus is usually on asymptotic convergence and the machinery is different; see <cit.> and references therein.  The small-step regime has seen less development, although there are some asymptotic bounds <cit.>.
  
  Fortunately, the qDRIFT construction has several appealing features: the random unitaries V_N,…,V_1 are i.i.d. unit-norm matrices that are close to the identity matrix (V - 𝕀≤ tλ/N almost surely) and close to their expectation (V-(𝔼V)≤ 2 t λ/N almost surely).
  These properties allow us to use the matrix martingale formalism to derive a strong, nonasymptotic result on the quality of the approximation.
  
   
  Consider a Hamiltonian H=∑_j=1^L h_j with interaction strength λ = ∑_j=1^L h_j, and fix parameters N, t. Suppose that V_N,…,V_1 are i.i.d. instances of the random unitary d × d matrix V constructed by the qDRIFT protocol (<ref>). Then
  
  
      Pr[  V_N ⋯ V_1 - 𝔼[V_N ⋯ V_1]≥ϵ/2]
      ≤ 2 d exp( - N ϵ^2/44 t^2 λ^2)   for any ϵ∈ [0,4t λ].
  
  
  
  
  
  
  
  
  
  
  
  
  
  In particular, N ≥ (44t^2 λ^2/ϵ^2) log (2d / δ) implies that V_N ⋯ V_1 - 𝔼[V_N ⋯ V_1] ≤ϵ/2 with probability at least 1-δ.
  
  
  This statement provides a strong tail bound for random fluctuations in the small-error regime
  ϵ≤ 4t λ.
  As N increases, the probability of incurring (at least) error ϵ/2 diminishes exponentially.
  For ϵ > 4t λ, we have instead a subexponential tail bound:
  Pr[  V_n ⋯ V_1 - 𝔼[V_N ⋯ V_1] ≥τ]≤ 2d exp (- N ϵ/ 6t λ ).
  We refer to (<ref>) for a unified statement that covers both regimes.
  
  The proof technique deserves some exposition, as it is rather general and may be of independent interest.
  It heavily utilizes the martingale concentration tools introduced in Appendix <ref>.
  For fixed N,
  we interpolate between both sides of Rel. (<ref>) by means of a random process {B_k: k=0,…,N}:
  
      B_k = (𝔼 V )^N-kV_k ⋯ V_1   such that   B_0 = (𝔼V)^N  and  B_N = V_N ⋯ V_1.
  
  The increments of this random process are certainly not independent. For instance, B_k depends on the (random) choice of V_k and all previous choices V_k-1, …, V_1.
  Fortunately, we can recognize it as a (matrix-valued) martingale 
  satisfying the defining properties:
  
  
  
    * Causality: Each B_k is completely determined by the information we have collected up to step k. That is, the (random) choices of V_k,…,V_1.
  
    * Status quo: Conditioned on previous choices, the expectation of B_k+1 equals B_k: for 1 ≤ k ≤ N
  
      𝔼[B_k+1|V_k ⋯ V_1 ]
      = (𝔼 V)^N-(k+1)𝔼_V_k+1[V_k+1] V_k ⋯ V_1 = (𝔼 V)^N-kV_k ⋯ V_1 = B_k.
  
  This feature underscores similarities to an unbiased random walk.  On average, “tomorrow” (B_k+1) is the same as “today” (B_k).
  
  
  With this matrix martingale reformulation at hand, we can prove Proposition <ref> using a concentration inequality for matrix martingales, namely Corollary <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  We have already established that the random process B_k = (𝔼  V )^N-k V_k ⋯ V_1 forms a matrix martingale that interpolates between B_0 = 𝔼[V_N ⋯ V_1] = (𝔼 V)^N and V_N = V_N ⋯ V_1.
  The elements of the associated difference sequence are
  
      C_k    = B_k - B_k-1
      = (𝔼V)^N-k( V_k - (𝔼 V_k)) V_k-1⋯ V_1
        with   k=1,…,N.
  
  
  Recall that V_k = exp(-i X_j) for some Hermitian matrices X_j = tNλh_jh_j with index 1 ≤ j ≤ L. Boundedness (𝔼 V  ,  V_k ≤ 1) and Fact <ref> (exp(-i X)-𝕀≤X for X Hermitian) ensure
  
      C_k    ≤𝔼 V^N-k V_k - (𝔼V_k) V_k-1⋯ V_1 ≤ V_k - (𝔼V_k)
      =  V_k - 𝕀 - 𝔼[V_k - 𝕀]
         ≤ V_k - 𝕀 + 𝔼[V_k - 𝕀] ≤ 2 max_j exp( - i X_j ) - 𝕀≤ 2 max_j  X_j 
      = 2t λ/N
  
  almost surely.
  Set R=2 t λ/N, and
  invoke Corollary <ref>
  
  to conclude that
  
      Pr[ B_N - B_0 ≥τ]
      
      
      ≤ 2d exp( -
      Nτ^2/8 (t λ)^2+4(t λ) τ/3).
  
  The statement follows from bounding the somewhat complicated exponential by either exp(-3 τ^2/(8NR^2)) for τ≤ 2 λ t or by exp( - 3 τ^2/(8R) ) for τ≥ 2 λ t.  Last, we substitute τ = ϵ/2.
  
  
  
  
  In fact, the same proof works for any adapted small-step random walks on the unitary group. Such a generalization results in Theorem <ref> and we refer to Appendix <ref> for details.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    §.§.§ A bound for expected errors
   
  
  In the previous subsection, we established that a sufficiently long qDRIFT random product formula
  concentrates sharply around its expectation.  
  We can translate this statement into a bound on the expected fluctuation
  around the true evolution.
  
  
  
   
  Consider an n-qubit Hamiltonian H=∑_j=1^L h_j with total strength λ = ∑_j=1^L h_j.
  Fix parameters N, t, and assume that N ≥ n.
  Set 𝒰=U X U^† with U = exp( -i t H ), and suppose that 𝒱_N,…,𝒱_1 ∼𝒱
  are i.i.d. realizations of the qDRIFT protocol.  That is, 𝒱(X) = V X V^†,
  where V is defined by (<ref>). Then
  
      𝔼[ 12𝒰-𝒱_N ∘⋯∘𝒱_1 _♢] ≤t^2 λ^2/N+ C nt λ/N+C √(n t^2 λ^2/N)≈ C √(n t^2 λ^2/N),
  
  where C >0 is a (modest) numerical constant.
  The symbol ≈ denotes an accurate approximation in the large-N regime.
  
  
  It is instructive to compare this assertion to the original qDRIFT result <cit.>
  and the improvement in (<ref>):
  
      12𝒰 - 𝔼[ 𝒱_N ∘⋯∘𝒱_1 ]_♢≤t^2 λ^2/N.
  
  Note that the expectation over all possible realizations of all N unitary channels appears inside the diamond distance. This implies that qDRIFT performs well on average over many random realizations, provided that the number N of steps exceeds t^2 λ^2/ϵ. In contrast, (<ref>) has the expectation outside the diamond distance.
  
  Our result gives a much stronger conclusion: An individual realization of the randomized qDRIFT protocol does not deviate much from the target evolution, for any input states and observables, with very high probability.  The price for such an improvement is a larger number of steps that depends on the system size.  For n qubits, the gate complexity N ≥ C n t^2λ^2/ϵ^2 is sufficient to ensure ϵ-closeness on average. The quadratic scaling in the accuracy parameter ϵ is necessary (for large N) because of the central limit theorem for martingales.
  
  The appearance of the number n of qubits is a consequence of measuring closeness in diamond distance.
  To obtain
  
      ϵ≥𝔼[ 12𝒰-𝒱_N ∘⋯∘𝒱_1 _♢]
          = 𝔼[ 12max_ρ state𝒰_ρ - 𝒱_N ∘⋯∘𝒱_1 (ρ) _1 ],
  
  we need the random product formula to behave for all possible n-qubit input states ρ simultaneously.
  If we restrict our attention to any fixed input state, we can obtain a gate complexity that does not
  depend on n.  This is the topic of the next section.
  
  
  First, we relate the expected diamond distance to an expected operator norm distance and split it up into deterministic bias and (expected) fluctuations:
  
      𝔼[ 12𝒰 - 𝒱_N ∘⋯∘𝒱_1 _♢]
      
      
      ≤U - (𝔼V)^N  + 𝔼V_N ⋯ V_1 - (𝔼V)^N
  
  The first term is deterministic and controlled by Proposition <ref>: U -(𝔼 V)^N ≤ t^2 λ^2/N.
  The second term can be bounded by integrating the tail bound in Proposition <ref>,  or rather the tighter bound presented (<ref>); see <cit.>.
  For n qubits, we have d=2^n and conclude
  
      𝔼 V_N ⋯ V_1 - (𝔼 V)^N    = ∫_0^∞Pr[  V_N ⋯ V_1 - (𝔼 V)^N ≥τ] dτ
         ≤∫_0^∞min(1,2 · 2^n exp( - Nτ^2/2/4 λ^2 t^2 +2 λ t τ/3)) dτ
         ≤ Cmax{√(nt^2 λ^2/N), nt λ/N}≤
      2C( n t λ/N+ √(nt^2 λ^2/N)).
  
  Here, we have evaluated the integral by two segments: the integrand is at most 1 until τ∼max( √(nt^2λ^2/N), n t λ/N ); for larger τ, the expression decays exponentially and
  integrating it only produces a contribution of order 𝒪(max( √(t^2λ^2/N), t λ/N )). Also, 2C absorbs all constants.
  
  
  
  
  
  
  
   §.§ 
  Proof of Theorem <ref>: Approximation error under a single  input
  
  
  
  
  Proposition <ref> asserts that a single, random realization of the qDRIFT protocol (<ref>) accurately approximates a unitary target evolution with respect to the diamond norm: 
  
      𝔼[ 12𝒰 - 𝒱_N ∘⋯∘𝒱_1 _♢]
          = 𝔼[ 12max_ρ state𝒰(ρ) - 𝒱_N ∘⋯∘𝒱_N (ρ) _1 ] ≲ C √(n t^2 λ^2/N).
  
  Here, ≲ denotes an accurate approximation of the true bound in the large N regime. This bound scales linearly in the (qubit) system size n.
  The dependence on n should not come as a surprise, since the diamond norm produces a very stringent worst-case distance measure. As emphasized by the above reformulation, the approximation must be accurate even when we optimize to find the worst possible input state ρ.
  
  In Hamiltonian simulation, demanding such a stringent worst-case promise may be excessive. In most practical applications, the input state ρ is fixed and simple, e.g., a product state.
  In this more practical setting, we can obtain a gate complexity N that does not depend on the system size n.
  The main result of this section asserts
  
      max_ρ state𝔼[ 12𝒰(ρ) - 𝒱_N ∘⋯∘𝒱_1 (ρ) _1 ] ≤ C √(t^2 λ^2/N).
  
  In other words, fixing an arbitrary input state ρ  helps a lot. A total number of N = 4 ( t λ/ϵ)^2 steps ensures that qDRIFT produces an ϵ-accurate output state, with respect to trace distance.
  
  
  
  The proof is similar in spirit to the argument behind Proposition <ref>.   We construct a vector martingale that describes the evolution of the state.  
  
  We control the behavior of this martingale using the uniform smoothness of the L_q(ℓ_2) norm.
  This argument is inspired by the work <cit.> on concentration of random matrix products.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    §.§.§ Approximation error for a fixed state
  
  
  In this section, we state and prove our main technical result
  on the action of the qDRIFT protocal on a fixed input state.
  
   
  Consider a Hamiltonian H=∑_j=1^L h_j with total strength λ = ∑_j=1^L h_j.
  Fix evolution time t and a gate count N that obeys N ≥ (t λ)^2. Let 𝒱_1, …, 𝒱_N be the i.i.d. random unitary evolution
  operators constructed by the qDRIFT protocol (<ref>).
  
  Then,
  
      max_ρ state𝔼[ 12𝒰(ρ) - 𝒱_N ∘⋯𝒱_1(ρ) _1 ]
          ≤ 4 √(t^2λ^2/N).
  
  Moreover, for ϵ > 0,
  
      max_ρ statePr[ 12𝒰(ρ) - 𝒱_N ∘⋯𝒱_1(ρ) _1
          > ϵ]
          ≤exp( -ϵ^2 N/32e t^2 λ^2).
  
  
  
  
  First, we reduce the problem to a question about pure states.
  For any q ≥ 2, Markov's inequality implies that
  
      Pr[ 12𝒰(ρ) - 𝒱_N ∘⋯𝒱_1(ρ) _1
          > ϵ]
             ≤ϵ^-q𝔼[ 2^-q𝒰(ρ) - 𝒱_N ∘⋯𝒱_1(ρ) _1^q ].
  
  The right-hand side of this equation is a convex function of the state ρ.
  Thus, the maximum 
  over all states is attained
  at a pure state.
  As a consequence, we can establish both claims in the proposition by limiting
  our attention to an (unknown) pure state  ρ = |ψ⟩⟨ψ|
  that does not depend on the random unitary channels 𝒱_k (X) = V X V^†.
  
  Next, we convert the trace distance of the output states into a Euclidean
  distance on the state vectors themselves.  The power q ≥ 2 will remain
  fixed until the last step of the argument.  Lemma <ref> implies
  
      ( 𝔼[ 2^-q𝒰(|ψ⟩⟨ψ|)-𝒱_N ∘⋯∘𝒱_1 (|ψ⟩⟨ψ|) _1^q  ] )^1/q≤( 𝔼 (  V_N ⋯ V_1 -U )| ψ⟩_ℓ_2^q )^1/q
                 ≤ 2 max{deterministic bias |ψ_bias⟩( 𝔼[ V_N ⋯ V_1 ] - U )|ψ⟩_ℓ_2,  ( 𝔼random fluctuation |ψ_rand⟩( V_N ⋯ V_1 - 𝔼[ V_N ⋯ V_1 ] )| ψ⟩_ℓ_2^q )^1/q}.
  
  The last bound follows from the triangle inequality
  and (a+b)^q ≤ 2^q max{ a^q, b^q } for a, b ≥ 0.
  
  We have split up the difference into two components, a deterministic bias
  and a random fluctuation.
  To control the deterministic bias, we simply apply Proposition <ref>:
  
      (𝔼[ V_N ⋯ V_1 ] - U) |ψ⟩_ℓ_2
      =  ((𝔼V)^N-U)| ψ⟩_ℓ_2≤ (𝔼V)^N-U ≤(t λ)^2/N.
  
  We will see that the bias is always negligible in comparison with the fluctuation.
  To control the second term, we need the following lemma.
  
   
  Let V_N,…,V_1 be i.i.d. unitaries that implement the qDRIFT protocol (<ref>) with parameters t and λ.
  Then, for any q ≥ 2,
  
      (𝔼 (V_N ⋯ V_1 - 𝔼[V_N ⋯ V_1])| ψ⟩_ℓ_2^q )^1/q≤ 2 √((q-1)(t λ)^2/N).
  
  
  
  We will establish this lemma below.
  The basic idea behind the proof is to express the random vector
  using a martingale sequence, similar to the matrix case. We could have already call vector martingale tail bounds (Fact <ref>) to arrive at the desired results. However, we will demonstrate the same results via markov's inequality and moment bounds (uniform smoothness, Proposition <ref>) which we introduced in Appendix <ref>.
  
  Introduce the inequalities from (<ref>) and Lemma <ref>
  into the bound (<ref>).  We obtain
  
      ( 𝔼[ 2^-q𝒰(|ψ⟩⟨ψ|)-𝒱_N ∘⋯∘𝒱_1 (|ψ⟩⟨ψ|) _1^q ] )^1/q
          ≤ 4 √((q-1) (tλ)^2/N).
  
  We have used the assumption that N ≥ (tλ)^2 to see that the second branch
  of the maximum always dominates the first.
  
  We may now complete the proof.
  To obtain the expectation bound, we set q = 2 in (<ref>)
  and apply Lyapunov's inequality.
  To obtain the probability bound, we combine (<ref>) and (<ref>)
  to arrive at
  
      Pr[ 12𝒰(ρ) - 𝒱_N ∘⋯𝒱_1(ρ) _1
          > ϵ]
          ≤( 16q (tλ)^2/ϵ^2 N)^q/2.
  
  Select q = (ϵ^2 N)/(16et^2 λ^2) to obtain the stated result.
  The resulting probability bound is vacuous unless q ≥ 2.
  
  
  
  
  
    §.§.§ Proof of Lemma <ref>
  
  
  In this section, we establish the bound on the size of the fluctuations.
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Fix a vector | ψ⟩, and introduce a sequence of random vectors: |ψ_k ⟩ = ∏_i=1^k V_i | ψ⟩ for 1≤ k ≤ N.  As a consequence, (V_N ⋯ V_1 - 𝔼[V_N ⋯ V_1])|ψ⟩ = |ψ_N ⟩ - 𝔼[| ψ_N⟩].
  We can recast this difference as a sum of two random vectors that are conditionally orthogonal in expectation:
  
  
      𝔼 | ψ_N⟩ - 𝔼[ | ψ_N ⟩] _ℓ_2^q
         = 𝔼 (V_N - 𝔼[V_N])|ψ_N-1⟩ + 𝔼[V_N]( |ψ_N-1⟩ - 𝔼[ | ψ_N-1⟩] ) _ℓ_2^q =: 𝔼y + x _ℓ_2^q.
  
  Indeed, 𝔼[y | x]=𝔼[V_N- (𝔼 V_N)] |ψ_N-1⟩ = 0.
  We can apply uniform smoothness (Proposition <ref>) to split up the contributions:
  
      (𝔼 | ψ_N⟩ - 𝔼[ | ψ_N ⟩] _ℓ_2^q)^2/q
         ≤ (q-1) (𝔼 (V_N - (𝔼V_N))|ψ_N-1⟩_ℓ_2^q)^2/q
              + (𝔼 (𝔼 V) (|ψ_N-1⟩ - 𝔼[ | ψ_N-1⟩]) _ℓ_2^q)^2/q
         ≤ (q-1) (𝔼 V_N - 𝔼[V_N]^q)^2/q + (𝔼 | ψ_N-1⟩ - 𝔼[ | ψ_N-1⟩]_ℓ_2^q)^2/q.
  
  We can now iterate this argument to conclude that 
  
      (𝔼 | ψ_N⟩ - 𝔼[ | ψ_N ⟩] _ℓ_2^q)^2/q
      
      ≤ (q-1) ∑_k=1^N (𝔼V_k - 𝔼[V_k]^q)^2/q
      = (q-1) N  (𝔼V - 𝔼[V] ^q)^2/q
  
  Invoke Lemma <ref>, using the properties of the random unitaries constructed by qDRIFT:
  
      (𝔼V - 𝔼[V] ^q)^2/q≤(2 t λ/N)^2.
  
  Combine the last two displays to reach the stated result.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Proof of Theorem <ref> and Corollary <ref>
  
  
  The proof techniques for establishing Theorem <ref> and Theorem <ref> are remarkably general and we condense them into Theorem <ref>. Let us reiterate the premise. Consider approximating a target unitary product U=U_N⋯ U_1 by a random unitary V=V_N⋯ V_1 such that { V_k} satisfy:(i) Causality: for 1 ≤ k ≤ N the random selection of V_k can only depend on previous choices for V_1,…,V_k-1:
  
      Pr[ V_k |V_N,…,V_k+1 ,V_k-1,…, V_1 ]
      =   Pr[ V_k| V_k-1,…,V_1 ]
  
  
  
  
  
  
  (ii) Accurate approximation: each realization of V_k (and their conditional expectation) must be close to the ideal unitary U_k. More precisely, let R,b_k >0 be constants such that 
      V_k-_k-1 V_k≤ R   and  𝔼_k-1 V_k - U_k ≤ b_k, where    𝔼_k-1V_k:=𝔼[ V_k| V_k-1,…,V_1 ] .
  
  almost surely for all 1 ≤ k≤ N.
  Note this is more general then we need to prove Theorem <ref>; this is to take into account the cases when the conditional variance may be much smaller than NR^2.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Recall the decomposition of approximation error into a deterministic bias and a random fluctuation:
  
      V_N ⋯ V_1 - U_N⋯ U_1    ≤deterministic bias (𝔼 V )^N - U_N⋯ U_1 
      +
      random fluctuation V_N ⋯ V_1 - 𝔼[ V_N ⋯ V_1 ] .
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  The deterministic bias can be once again be controlled by a telescoping sum,
  
      (𝔼 V )^N - U 
        
           ≤∑_k=1^N b_k.
  
  Note that this also controls the performance of channel - _♢ by Lemma <ref>.
   For the random fluctuations, tweaking the proofs of Theorem <ref> and Theorem <ref> implies the following general result.
  
  
  
  
  
  
  Let {V_1, V_2, ⋯, V_N }⊂ U(d) an adapted(causal) random unitary matrices that obey
  
      ∑^N_k=1_k-1‖(V_k -  V_k)(V^†_k -  V^†_k)‖≤σ^2
          and  ‖ V_k -  V_k ‖≤ R. (almost surely.)
  
  for some constants v, R ≥ 0.
  Then,  the product of N random unitaries satisfies a concentration inequality:
  
      Pr(‖ V_N⋯ V_1 - [ V_N ⋯ V_1] ‖≥ϵ  ) ≤ 2d exp( -ϵ^2/2/ v+ Rϵ/3)   for ϵ >0.
  
  For a fixed, but arbitrary, input state ρ, concentration is independent of the ambient dimension d:
  
      max_ρ statePr( 12[𝒱_N ∘⋯𝒱_1(ρ)] - 𝒱_N ∘⋯𝒱_1(ρ) _1
          > ϵ)
          ≤ 2exp( -ϵ^2/2 /v+Rϵ/3)   for ϵ >0.
  
  
  
  
  
  For a fixed input state, we would call the vector Freedman inequality (Fact <ref>) instead of uniform smoothness (Proposition <ref>) in Theorem <ref>.
  There are several recent independent papers that also use matrix martingale tools to study products of random matrices that are close to the identity.  The work <cit.> addresses the problem using uniform smoothness tools.  The paper <cit.> uses the matrix Freedman inequality; their proof is quite similar to ours.  In contrast, we are interested in unitary products, which allows for additional simplifications.  For more background on matrix martingales, see <cit.>.
  
  It is instructive to illustrate these improvements by example. In qDRIFT, all steps have a uniform bound R, but in the fully general statement the variance v can differ from the crude uniform bound NR^2. In such a regime, the sub-exponential tail of size ^-3ϵ/2R can start playing a role.
  
  
  Lastly, for illustration in Theorem <ref> we give loose estimates on the variance to avoid complication with the heavy tail effects. Plugging in the parameters v = ra^2, R = 2a as V- V ≤ V -V'≤ 2  U-V = 2a translates to a typical fluctuation ϵ^2 ∼ n ra^2(and ϵ^2 ∼ ra^2 for fixed input). We conclude the proof by combining the bound for the deterministic bias and random fluctuation.
  
  
  
  Consider randomly permuting the 2k-th order Suzuki formulas as
  V_k = S^σ_2k(t/r) with uniform probability 1/L!. Then by direct calculation <cit.>:
  
      V-U   ≤ O((tΛ L)^2k+1/r^2k+1)=a ,
       V -U   ≤  O(tΛ/r)^2k+1L^2k)=b
  
  by Theorem <ref>,
  
      ϵ_det   ≤ ra =O( (tΛ L)^2k+1/r^2k )
      ϵ_typ   = O(√(nr)a)+2rb = O(√(n)(tΛ L)^2k+1/r^2k+1/2 +(tΛ)^2k+1(L/r)^2k)
      ϵ_fix   = O(√(r)a)+2rb = O( (tΛ L)^2k+1/r^2k+1/2+(tΛ)^2k+1(L/r)^2k) 
      ϵ_avg   ≤ 2rb  =O(  (tΛ)^2k+1(L/r)^2k)
  
  Which translates to the sufficient gate counts N =rL
      
      N_det   ≳ tΛ L^2 (tΛ L/ϵ)^1/2k
      
          N_typ   ≳ tΛ L^2(nΛ Lt/ϵ^2)^1/(4k+1)+t Λ L^2 (tΛ/ϵ)^1/2k
      
          N_fix   ≳ tΛ L^2(Λ Lt/ϵ^2)^1/(4k+1)+t Λ L^2 (tΛ/ϵ)^1/2k
      
          N_avg   ≳  t Λ L^2 (tΛ/ϵ)^1/2k
  
  
  
  
  
  § ASYMPTOTIC TIGHTNESS
   
  
  
  
  
  
  It is natural to wonder whether the bound (<ref>) is tight for some Hamiltonian H=∑_j h_j, i.e., whether N = Ω( n λ^2 t^2/ϵ^2) is also necessary to achieve concentration.
  More precisely, we want to understand whether the dependences on system size n = log_2(d), evolution time t and interaction strength λ=∑_j h_j are also necessary to control the typical deviation of the unitary random walk we considered.
  
  In the context of matrix concentration inequalities, this question has been thoroughly addressed <cit.>. The answer is affirmative for sums of bounded matrices: concentration inequalities are tight and saturated for collections of commuting matrices. Although in this work we consider products of random matrices,  we are still using a telescoping sum in the small step regime and expect an analogy.
  
  This observation motivates us to look at artificial Hamiltonians whose associated unitary evolution saturates the upper bounds put forth in this work.
  The cases we can handle lie at the two extremes: either the sum of single-site Pauli Zs or the sum of all 2^n many-body Pauli Zs.  We will see the presence of
  the system size factor n = log_2(d) at both extremes, so one may believe the same to hold for the intermediate q-local cases.
  However, this factor arises for very different reasons.
  It arises in the single-site case, because the operator norm completely factorizes into n constituents (one for each term).
  For Hamiltonians that encompass all 2^n many-body Zs, it comes from the fact that diagonal entries are nearly independent, so the union bound we used in Section <ref> is tight.  Independence of entries requires the presence of all many-body terms, and does not extend to the few-body case.
  
  The multivariate central limit Theorem will be crucial for analyzing both cases, as it greatly simplifies the analysis in large N limit.
  
  [CLT for the multinomial distribution]
  
  The multinomial distribution m=(m_1,…,m_K)∼Mult(N,(1/K,…,1/K)) (roll a fair K-sided dice N times) obeys a 
  central limit theorem (CLT):
  
      1√(N)(m-𝔼m) ∼𝒩(0,Σ)   in distribution as N →∞.
  
  The covariance matrix is
  Σ = 1K(𝕀-1KJ), where J denotes the K × K matrix of ones.
  
  
  
  
   §.§ Sum of single site Pauli-Z operators
  
  
  
  
  This example demonstrates the saturation of our martingale bounds for single site Hamiltonians that factorize completely.
  To this end, we revisit a variant of the n-qubit example Hamiltonian discussed in Section <ref>:
  
      H= ∑_k=1^n Z_k   where   Z_k = (k-1)-times𝕀⊗⋯⊗𝕀⊗ Z ⊗(n-k)-times𝕀⊗⋯⊗𝕀  for 1 ≤ k ≤ n.
  
  
   
  Suppose that we wish to obtain an N-term approximation of the time evolution U = exp(-itH) associated with the n-qubit Hamiltonian (<ref>) for evolution time t. In the large N limit (CLT), the qDRIFT approximation (<ref>) incurs an operator norm error that matches the (upper) bound from Corollary <ref> (and indirectly Theorem <ref>) up to a constant factor:
  
      𝔼U-V_N ⋯ V_1≥√(2/π)√((n-1) (tλ)^2/N)- 1/2(n-1) (t λ)^2/N.
  
  
  
  We have chosen to state this result directly in terms of operator norm deviation. A conversion into diamond distance is also possible: 12𝒰-𝒱_♢≥12U-V for any pair of unitary channels.
  This conversion rule readily follows from the geometric characterization of 12𝒰-𝒱_♢ provided in <cit.>.
  
  
  Each of the n terms in the Hamiltonian (<ref>) has unit operator norm (Z_k=1) and the total strength is λ=∑_k=1^n Z_k=n. For fixed N and t, each short-time approximation (<ref>) has the form V_k= exp( -itnNZ_k(i)), where each k(i) is an index chosen uniformly from the set {1,…,n} (multinomial distribution). Since all Z_ks commute, we can rewrite the entire product formula as
  
      V_N ⋯ V_1 = exp( - itλ/N∑_i=1^N Z_k(i))
      = exp( - itλ/N∑_k=1^n m_k Z_k ).
  
  Here, we have introduced the count statistics m_k for each site label k, that is the number of times location k has been selected throughout N independent selection rounds, to rearrange the sum.
  This count statistics obeys m̅_k = 𝔼 m_k = N/n=N/λ for each 1 ≤ k ≤ n. We can use this observation to re-express the target unitary U in a compatible fashion:
  
      U = exp( - it ∑_k=1^n Z_k ) = exp( - itλ/N∑_k=1^n m̅_k Z_k ).
  
  Unitary invariance then implies that the operator norm difference between both unitaries becomes
  
      V_N ⋯ V_1 - U = exp( - itλ√(N)∑_k=1^N m_k - m̅_k/√(N) Z_k) - 𝕀.
  
  This is a promising starting point. The multinomial CLT (Fact <ref>) ensures that the n centered and normalized random variables s_k = (m_k-m̅_k)/√(N) approach the coefficients of a Gaussian vector s ∈ℝ^n with covariance matrix Σ = 1n( 𝕀 - 1nJ ).
  This, in particular implies 𝔼 s_k = 0 and 𝔼s_k^2 = 1n(1-1n)=σ^2 for all 1 ≤ k ≤ n. We can capitalize on this observation by simplifying (<ref>) via a second-order Taylor expansion. Set X = -tλ√(N)∑_k s_k Z_k for brevity and apply Fact <ref> to obtain
  
      V_N ⋯ V_1 - U = exp(iX)-𝕀≥iX - exp(iX)-iX - 𝕀≥X - 12X^2.
  
  This relation is preserved under expectations and we obtain
  
      𝔼V_N ⋯ V_1 - U≥tλ/√(N)𝔼∑_k s_k Z_k - 1/2( tλ/√(N))^2 𝔼∑_k s_k Z_k^2.
  
  Let us focus on the leading order term first. The particular structure of the Hamiltonian (<ref>) – each Z_k is the tensor product of a single Pauli-Z matrix at location k with (n-1) identity matrices – ensures that the operator norm factorizes nicely. Use X ⊗𝕀 + 𝕀⊗ Y = X+Y iteratively to conclude
  
      𝔼∑_k s_k Z_k = 𝔼∑_k=1^n s_k Z_k = ∑_k=1^n |s_k| N →∞= n√(2/π1/n(1-1/n))=√(2/π(n-1)),
  
  because the CLT asserts that each |s_k| approaches a half-normal random variable with σ^2 = 1n(1-1n).
  
  To bound the quadratic term, we combine the factorization trick from above with a well-known relation among ℓ_p-norms in ℝ^n:
  
      𝔼∑_k=1^n s_k Z_k^2 = 𝔼( ∑_k=1^n |s_k| )^2
      = 𝔼s_ℓ_1^2 ≤ n 𝔼s_ℓ_2^2 = n ∑_k=1^n 𝔼 s_k^2 = n^2 σ^2 = (n-1).
  
  No CLT is required for this argument.
  Inserting both bounds into Eq. (<ref>) completes the argument.
  
  
  
  
  
  
  
  
  
   §.§ Sum of many-body Pauli-Z operators
  
  
  
  Let us revisit the example Hamiltonian from Sec. <ref>, albeit without additional sign factors. Recall the multi-indices p =(p_1,…,p_n) ∈{0,1}^n  and set
  
      H = ∑_𝐩∈{0,1}^n Z_𝐩
      = ∑_𝐩∈{0,1}^n Z^p_1⊗⋯⊗ Z^p_n,
  
  where we use the conventions Z^1=Z and Z^0 = 𝕀.
  This Hamiltonian is not local.  All constituents commute and have the same operator norm: Z_𝐩=1 for all 𝐩∈{0,1}^n.
  This in turn implies that the total strength λ=∑_𝐩 Z_𝐩=2^n equals the Hilbert space dimension. It is also worthwhile to point out that each term is diagonal in the computational basis |𝐛⟩ = |b_1,…,b_n ⟩ with 𝐛=(b_1,…,b_n)∈{0,1}^n. Overlaps of the Hamiltonian terms with computational basis states are given by
  
      ⟨𝐛 | Z_𝐩 | 𝐛⟩ = (-1)^⟨𝐛,𝐩⟩ = (-1)^∑_i b_i p_i∈{± 1 }.
  
  The following claim highlights that our findings are tight for asymptotically large step sizes N. This complements the example upper bound derived in Sec. <ref>, as well as Theorem <ref>.
  
   
  Suppose that we wish to obtain an N-term approximation of the time evolution U = exp( -it H) associated with the n-qubit Hamiltonian (<ref>) for evolution time t.
  In the large N limit (CLT) the qDRIFT approximation
  (<ref>) incurs an operator norm error that matches the (upper) bound from Corollary <ref> (and indirectly Theorem <ref>) up to a constant factor:
  
      𝔼U - V_N ⋯ V_1 ≥1/2√(n (tλ)^2/N) -  2 (n+12) (t λ)^2/N
  
  
  
  The conversion rule 𝒰-𝒱_♢≥U-V (for unitary channels) <cit.> once more allows for addressing the expected diamond distance as well.
  
  
  
  Each of the 2^n terms in the Hamiltonian (<ref>)  has unit operator norm (Z_p=1 for all p∈{0,1}^n) and the strength is λ = ∑_pZ_p=2^n.
  For fixed N and t,
   each short-time approximation (<ref>) has the form V_k = exp (-it λN Z_p(i)), where p(i) is a string chosen uniformly at random from all 2^n possibilities (multinomial distribution). Since all Z_ps commute,
   we can rephrase and simplify the expected operator norm difference in a fashion analogous to the proof of Proposition <ref>:
  
      𝔼V_N ⋯ V_1 - U≥tλ/√(N)𝔼∑_ps_pZ_p-1/2( t λ/√(N))^2 𝔼∑_ps_pZ_p^2.
  
  Here, s_p=(m_p-m̅_p)/√(N) is the centered and normalized variant of the count statistics m_p associated with bit string p∈{0,1}^n, that is the number of times the Hamiltonian term Z_p has been selected throughout N independent selection rounds.
  The multinomial CLT (Fact <ref>) asserts that the 2^n centered and normalized random variables 
  s_p approach distinct coefficients of a 2^n-dimensional Gaussian vector with covariance matrix
   Σ = 12^n(𝕀-12^n J )=12^n(𝕀-|1⟩⟨1|), where |1⟩ = 12^n∑_b∈{0,1}^n|b⟩ (the normalized all-ones vector in ℝ^2^n).
   In contrast to before, the individual contributions to this operator norm don't factor nicely anymore. Establishing tight bounds requires additional analysis.
  
  Let us focus on the (leading) first-order term for now. All matrix summands in the expression commute and are diagonal in the computational basis |b⟩ with b=(b_1,…,b_n) ∈{0,1}^n. This ensures that the operator norm is attained at a computational basis state:
  
      ∑_p Z_ps_p = max_b∈{0,1}^n| ∑_p s_p⟨b|Z_p| b⟩| = max_b∈{0,1}^n| ∑_p (-1)^⟨b,p⟩ s_p|,
  
  where the last equation is due to Rel. (<ref>).
  This expression is proportional to the largest entry (in modulus) of the Walsh-Hadamard transform of the 2^n-dimensional vector s with entries s_p for p∈{0,1}^n. More precisely,
  
      max_b∈{0,1}^n| ∑_p (-1)^⟨b,p⟩ s_p| =  2^n/2Had^⊗ n s _ℓ_∞=: 2^n/2ŝ_ℓ_∞  where  Had = 1√(2)(
      [  1  1;  1 -1 ]).
  
  We emphasize that the Walsh-Hadamard transform is an orthogonal transformation, which also applies to the limiting covariance matrix of ŝ=Had^⊗ ns (CLT):
  
      Σ̂ = 12^nHad^⊗ n(𝕀-|1⟩⟨1| ) Had^⊗ n
      =12^n( 𝕀 - |0,…,0 ⟩⟨ 0,…,0| ).
  
  Hence, the CLT asserts that the transformed vector ŝ approaches a standard Gaussian vector with 2^n-1 degrees of freedom: ŝ=(0,g_2,…,g_2^n)^T with g_i i.i.d.∼𝒩(0,2^-n) (one degree of freedom is erased by the normalization constraint ∑_pm_p=N of the count statistics).
  The bound on the expected leading order contribution now follows from invoking the well-known fact that the expected maximum of K standard Gaussian random variables with equal variance σ^2 is lower-bounded by 0.265√(log(K)σ^2), see e.g. <cit.>:
  
      tλ√(N)𝔼∑_ps_pZ_p
      =   t λ/√(N)2^n/2𝔼ŝ_ℓ_∞N →∞=t λ/√(N) 2^n/2𝔼max_2 ≤ i ≤ 2^n|g_i| 
      ≥    0.625 tλ/√(N)2^n/2√(log(2^n-1)2^-n)
      ≥1/2√(n (tλ)^2/N).
  
  Here, we have used the numerical bound 0.625√(log(2^n-1)/n)≥ 0.5 which is valid for any n ≥ 3 (for n=2 the ratio is slightly smaller).
  This completes the argument for the leading term in Eq. (<ref>).
  
  Moving on to the quadratic term in Eq. (<ref>), we employ a similar strategy. Observe
  
      𝔼∑_p Z_p^2
       =   𝔼∑_p,p' s_p s_p' Z_p Z_p'
      = 𝔼max_b∈{0,1}^n| ∑_p,p' s_p s_p'⟨b|Z_p Z_p' |b⟩| 
      
      =   𝔼max_b∈{0,1}^n| ∑_p (-1)^⟨b,p⟩ s_p∑_p' (-1)^⟨b,p'⟩ s_p'|,
  
  where the last equation follows from combining Eq. (<ref>) with the appealing group structure of the Z_p's: Z_pZ_p'=Z_p⊕p', where ⊕ denotes entry-wise addition modulo 2 (the set of all Z_p's form a maximal stabilizer group). We can now recognize two independent Walsh-Hadamard transforms of the 2^n-dimensional vector s in this expression:
  
      𝔼max_b∈{0,1}^n| ∑_p (-1)^⟨b,p⟩ s_p∑_p' (-1)^⟨b,p'⟩ s_p'|
      = 2^n 𝔼max_b∈{0,1}^n| ŝ_b|^2
  
  We already know from the CLT that the 2^n-dimensional Walsh-Hadamard transform of s approaches a standard Gaussian vector: ŝ=(0,g_2,…,g_2^n)^T with g_i i.i.d.∼𝒩(0,2^-n).
  In the large N limit (CLT), the r.h.s. of the above display becomes an expected maximum of K=2^n-1 squares of i.i.d. Gaussian variables with mean zero and variance σ^2=2^-n.
  Such expected maxima can be bounded using standard arguments, see e.g. <cit.>: 𝔼max_1 ≤ i ≤ K |g_i|^2 ≤ 4 σ^2 log (√(2)K) (the constants are chosen based on simplicity, not tightness).
  This allows us to conclude
  
      𝔼∑_p Z_p^2
      =2^n 𝔼max_b∈{0,1}^n| ŝ_b|^2 N →∞= 2^n 𝔼max_2 ≤ i ≤ N|g_i |^2
      ≤ 2^n 4 σ^2 log(√(2)(2^n-1))
      
      ≤ 4(n+1/2).
  
  Inserting linear and quadratic bound into Eq. (<ref>) completes the argument.
  
  
  
  
  
  10
  
  AKN98:Quantum-Circuits
  D. Aharonov, A. Kitaev, and N. Nisan.
  Quantum circuits with mixed states.
  In Proceedings of the Thirtieth Annual ACM Symposium on Theory
    of Computing, STOC ’98, page 20–30, New York, NY, USA, 1998. Association
    for Computing Machinery.
  
  AW02:Strong-Converse
  R. Ahlswede and A. Winter.
  Strong converse for identification via quantum channels.
  IEEE Trans. Inform. Theory, 48(3):569–579, 2002.
  
  an2019quantum
  D. An and L. Lin.
  Quantum linear system solver based on time-optimal adiabatic quantum
    computing and quantum approximate optimization algorithm.
  preprint arXiv:1909.05500, 2019.
  
  Babbush_2019
  R. Babbush, D. W. Berry, and H. Neven.
  Quantum simulation of the sachdev-ye-kitaev model by asymmetric
    qubitization.
  Phys. Rev. A, 99:040301, 2019.
  
  babbush2018low
  R. Babbush, N. Wiebe, J. McClean, J. McClain, H. Neven, and G. K.-L. Chan.
  Low-depth quantum simulation of materials.
  Phys. Rev. X, 8:011044, 2018.
  
  berry2015hamiltonian
  D. W. Berry, A. M. Childs, and R. Kothari.
  Hamiltonian simulation with nearly optimal dependence on all
    parameters.
  In 2015 IEEE 56th Annual Symposium on Foundations of
    Computer Science—FOCS 2015, pages 792–809. IEEE Computer Soc., Los
    Alamitos, CA, 2015.
  
  berry2020time
  D. W. Berry, A. M. Childs, Y. Su, X. Wang, and N. Wiebe.
  Time-dependent Hamiltonian simulation with L^1-norm scaling.
  Quantum, 4:254, 2020.
  
  campbell_mixing16
  E. Campbell.
  Shorter gate sequences for quantum computing by mixing unitaries.
  Phys. Rev. A, 95:042306, 2017.
  
  campbell2019random
  E. Campbell.
  Random compiler for fast hamiltonian simulation.
  Phys. Rev. Lett., 123:070503, 2019.
  
  chen2021concentration
  C.-F. Chen.
  Concentration of otoc and lieb-robinson velocity in random
    hamiltonians, 2021.
  
  chen2019operator
  C.-F. Chen and A. Lucas.
  Operator growth bounds from graph theory, 2019.
  
  chen2021optimal
  C.-F. Chen and A. Lucas.
  Optimal frobenius light cone in spin chains with power-law
    interactions, 2021.
  
  childs2019faster
  A. M. Childs, A. Ostrander, and Y. Su.
  Faster quantum simulation by randomization.
  Quantum, 3:182, 2019.
  
  childs2019theory
  A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu.
  A theory of trotter error.
  preprint arXiv:1912.08854, 2019.
  
  LCU
  A. M. Childs and N. Wiebe.
  Hamiltonian simulation using linear combinations of unitary
    operations.
  Quantum Inf. Comput., 12(11-12):901–924, 2012.
  
  Christofides2008martingales
  D. Christofides and K. Markström.
  Expansion properties of random Cayley graphs and vertex transitive
    graphs via matrix martingales.
  Random Structures Algorithms, 32(1):88–100, 2008.
  
  dawson2005kitaev
  C. M. Dawson and M. A. Nielsen.
  The Solovay-Kitaev algorithm.
  Quantum Info. Comput., 6:81, 2006.
  
  faehrmann2021randomized
  P. K. Faehrmann, M. Steudtner, R. Kueng, M. Kieferova, and J. Eisert.
  Randomizing multi-product formulas for improved Hamiltonian
    simulation.
  preprint arXiv:2101.07808, 2021.
  
  FR13:Compressed-Sensing
  S. Foucart and H. Rauhut.
  A Mathematical Introduction to Compressive Sensing.
  Applied and Numerical Harmonic Analysis. Birkhäuser, 2013.
  
  garling_2007
  D. J. H. Garling.
  Inequalities: a journey into linear analysis.
  Cambridge University Press, Cambridge, 2007.
  
  georgescu2014quantum
  I. M. Georgescu, S. Ashhab, and F. Nori.
  Quantum simulation.
  Rev. Mod. Phys., 86:153–185, 2014.
  
  Gross2011:Recovering
  D. Gross.
  Recovering low-rank matrices from few coefficients in any basis.
  IEEE Trans. Inform. Theory, 57(3):1548–1566, 2011.
  
  Haah_2021
  J. Haah, M. B. Hastings, R. Kothari, and G. H. Low.
  Quantum algorithm for simulating real time evolution of lattice
    hamiltonians.
  SIAM J. on Comput., page FOCS18–250–FOCS18–284, 2021.
  
  hastings2016turning
  M. B. Hastings.
  Turning gate synthesis errors into incoherent errors, 2016.
  
  HNTR20:Matrix-Product
  D. Huang, J. Niles-Weed, J. A. Tropp, and R. Ward.
  Matrix concentration for products.
  preprint arXiv:2003.05437, 2020.
  
  huang2019near
  H.-Y. Huang, K. Bharti, and P. Rebentrost.
  Near-term quantum algorithms for linear systems of equations.
  preprint arXiv:1909.07344, 2019.
  
  kathuria2020concentration
  T. Kathuria, S. Mukherjee, and N. Srivastava.
  On concentration inequalities for random matrix products.
  preprint arXiv:2003.06319, 2020.
  
  Lashkari_2013
  N. Lashkari, D. Stanford, M. Hastings, T. Osborne, and P. Hayden.
  Towards the fast scrambling conjecture.
  Journal of High Energy Physics, 2013(4), Apr 2013.
  
  google2020THC
  J. Lee, D. Berry, C. Gidney, W. J. Huggins, J. R. McClean, N. Wiebe, and
    R. Babbush.
  Even more efficient quantum computations of chemistry through tensor
    hypercontraction.
  preprint arXiv:2011.03494, 2020.
  
  lloyd1996universal
  S. Lloyd.
  Universal quantum simulators.
  Science, 273(5278):1073–1078, 1996.
  
  low2017optimal
  G. H. Low and I. L. Chuang.
  Optimal hamiltonian simulation by quantum signal processing.
  Phys. Rev. Lett., 118:010501, 2017.
  
  Low_2019_qubitize
  G. H. Low and I. L. Chuang.
  Hamiltonian simulation by qubitization.
  Quantum, 3:163, 2019.
  
  maldacena2016remarks
  J. Maldacena and D. Stanford.
  Remarks on the Sachdev-Ye-Kitaev model.
  Phys. Rev. D, 94:106002, 2016.
  
  martyn2021grand
  J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang.
  A grand unification of quantum algorithms.
  preprint arXiv:2105.02859, 2021.
  
  mcardle2020quantum
  S. McArdle, S. Endo, A. Aspuru-Guzik, S. C. Benjamin, and X. Yuan.
  Quantum computational chemistry.
  Rev. Mod. Phys., 92:015003, 2020.
  
  oliveira2010concentration
  R. I. Oliveira.
  Concentration of the adjacency matrix and of the laplacian in random
    graphs with independent edges, 2010.
  
  oliveira2009spectrum
  R. I. Oliveira.
  The spectrum of random k-lifts of large graphs (with possibly
    large k).
  J. Comb., 1(3-4):285–306, 2010.
  
  ouyang2020compilation
  Y. Ouyang, D. R. White, and E. T. Campbell.
  Compilation by stochastic Hamiltonian sparsification.
  Quantum, 4:235, 2020.
  
  pinelis2012optimum
  I. Pinelis.
  Correction: “Optimum bounds for the distributions of martingales
    in Banach spaces” [Ann. Probab. 22 (1994), no. 4, 1679–1706;
    MR1331198 (96b:60010)].
  Ann. Probab., 27(4):2119, 1999.
  
  polchinski2016spectrum
  J. Polchinski and V. Rosenhaus.
  The spectrum in the Sachdev-Ye-Kitaev model.
  J. High Energy Phys., 2016(4):001, front matter+24, 2016.
  
  sachdev1993gapless
  S. Sachdev and J. Ye.
  Gapless spin-fluid ground state in a random quantum heisenberg
    magnet.
  Phys. Rev. Lett., 70:3339–3342, 1993.
  
  low_energy2020
  B. Sahinoglu and R. D. Somma.
  Hamiltonian simulation in the low energy subspace.
  preprint arXiv:2006.02660, 2020.
  
  subacsi2019quantum
  Y. Subaşi, R. D. Somma, and D. Orsucci.
  Quantum algorithms for systems of linear equations inspired by
    adiabatic quantum computing.
  Phys. Rev. Lett., 122:060504, 2019.
  
  suzuki1985decomposition
  M. Suzuki.
  Decomposition formulas of exponential operators and Lie
    exponentials with some applications to quantum mechanics and statistical
    physics.
  J. Math. Phys., 26(4):601–612, 1985.
  
  suzuki1991general
  M. Suzuki.
  General theory of fractal path integrals with applications to
    many-body theories and statistical physics.
  J. Math. Phys., 32(2):400–407, 1991.
  
  tropp2011freedman
  J. A. Tropp.
  Freedman's inequality for matrix martingales.
  Electron. Commun. Probab., 16:262–270, 2011.
  
  Tro12:User-Friendly
  J. A. Tropp.
  User-friendly tail bounds for sums of random matrices.
  Found. Comput. Math., 12(4):389–434, 2012.
  
  tropp2015introduction
  J. A. Tropp.
  An introduction to matrix concentration inequalities, 2015.
  
  VH14:Probability-Lecture
  R. van Handel.
  Probability in high dimension.
  Technical report, PRINCETON UNIV NJ, 2014.
  
  varj2012random
  P. P. Varjú.
  Random walks in compact groups.
  Doc. Math., 18:1137–1175, 2013.
  
  versendaal2019large
  R. Versendaal.
  Large deviations for random walks on lie groups, 2019.
  
  watrous2018book
  J. Watrous.
  The Theory of Quantum Information.
  Cambridge University Press, 2018.
  
  
  
  