A double-slit proposal for quantum annealing

We formulate and analyze a double-slit proposal for quantum annealing, which involves observing the probability of finding a two-level system (TLS) undergoing evolution from a transverse to a longitudinal field in the ground state at the final time $t_f$. We demonstrate that for annealing schedules involving two consecutive diabatic transitions, an interference effect is generated akin to a double-slit experiment. The observation of oscillations in the ground state probability as a function of $t_f$ (before the adiabatic limit sets in) then constitutes a sensitive test of coherence between energy eigenstates. This is further illustrated by analyzing the effect of coupling the TLS to a thermal bath: increasing either the bath temperature or the coupling strength results in a damping of these oscillations. The theoretical tools we introduce significantly simplify the analysis of the generalized Landau-Zener problem. Furthermore, our analysis connects quantum annealing algorithms exhibiting speedups via the mechanism of coherent diabatic transitions to near-term experiments with quantum annealing hardware.


I. INTRODUCTION
Feynman famously wrote that the double-slit interference experiment "...has in it the heart of quantum mechanics. In reality, it contains the only mystery" [1]. Here we propose a double-slit experiment for quantum annealing (QA). In analogy to Feynman's particle-wave doubleslit, the proposed experiment can only be explained by the presence of interference and would break down upon either an intermediate measurement or strong decoherence. We are motivated by the recent resurgence of interest in quantum annealing using the transverse field Ising model [2,3], which has led to major efforts to build physical quantum annealers for the purpose of solving optimization and sampling problems [4][5][6][7], and significant debate as to whether quantum effects are at play in the performance of such devices [8,9]. The mechanisms by which QA might achieve a speedup over classical computing remain hotly contested, and while tunneling is often promoted as a key ingredient [10] and entanglement is often viewed as a necessary condition which must be demonstrated [11,12], a consensus has yet to emerge. Yet, an explicit example is known where QA theoretically provides an oracle-based exponential quantum speedup over all classical algorithms [13], and other examples are known where QA provides a speedup over classical simulated annealing [14][15][16][17][18][19]. An essential feature in all these cases are diabatic transitions which circumvent adiabatic ground state evolution to enable the speedup, in the spirit of the idea of shortcuts to adiabaticity [20]. When these transitions result in a coherent recombination of the ground state amplitude (a phenomenon known as a diabatic cascade [16,21]), the result is a wave-like interference pattern in the ground state probability as the anneal time is varied [22][23][24]. We thus conjecture that coherent recombination of ground state amplitudes after coherent evolution between diabatic transitions can play a critical role in enabling quantum speedups in QA. The double-slit proposal we formulate and analyze here is designed to test for the presence of quantum interference due to such coherent evolution.
Viewed from a different perspective, our double-slit proposal joins a family of protocols designed to probe the dynamics of what Berry called the "simplest nonsimple quantum problem" [25], a driven TLS near level crossings [26]. The two-level paradigm was introduced long ago by Landau and Zener (LZ) [27,28]. The corresponding Hamiltonian for the generalized LZ problem is where X, Y and Z are the Pauli matrices. In the original protocol which LZ solved analytically, a(t) is constant, b(t) is linear in t, and t runs from −∞ to ∞. The problem has since been studied under numerous variations, including Landau-Zener-Stueckelberg interferometry where b(t) is periodic [29,30], the subject of various experiments [31][32][33]. Complete analytical solutions were limited until recently to certain particular functional forms of b(t) with constant a(t) [34], a finite-range linear schedule for both a(t) and b(t) [35], and periodic a(t) and b(t) [36]. An analytical solution for general b(t) but constant a(t) was found in Ref. [37], which was then extended to general (but implicitly specified) a(t) as well [38,39].
Here we consider the case of general schedules a(t) and b(t), and develop a simple to interpret, yet surprisingly accurate, low-order time-dependent perturbation theory approach, that allows us to identify a class of schedules exhibiting "giant" (relative to linear schedules) interference oscillations of the ground state population as a function of the total annealing time. Our proposal should in principle be straightforward to implement using, e.g., flux qubits, and toward this end we also study the effects of coupling to a thermal environment. The structure of this paper is as follows. In Sec. II we analyze the TLS quantum annealing problem in the closed system limit. We first transform to an adiabatic interaction picture and perform a Magnus expansion, which allows us to give a simple expression for the ground state probability in terms of the Fourier transform of a key quantity we call the angular progression. We then analyze both the LZ problem (with a linear schedule) and a new "Gaussian angular progression" schedule which gives rise to large interference oscillations. We explain how these oscillations can be interpreted in terms of a doubleslit experiment generating interference between ground state amplitudes. In Sec. III we analyze the problem in the presence of coupling to a thermal environment. We consider the weak-coupling limit both without and with the rotating wave approximation, and find the range of coupling strengths and temperatures over which the interference oscillations are visible, using parameters relevant for superconducting flux qubits. We find a simple semi-empirical formula that accurately captures all our open-system simulation results in terms of three physically intuitive quantities: the oscillation period, rate of convergence to the adiabatic limit, and damping due to coupling to the thermal environment. We express all three are in terms of the input parameters of the theory. Conclusions and the implications of our results are discussed in Sec. IV. A variety of supporting technical calculations and bounds are provided in the Appendix.

II. CLOSED SYSTEM ANALYSIS AND RESULTS
A. Adiabatic interaction picture for two-level system quantum annealing We first consider the closed system setting. Consider a two-level system (TLS) quantum annealing Hamiltonian in the standard form (1), where the annealing schedules a(t), b(t) ≥ 0 respectively decrease/increase to/from 0 with time t ∈ [0, t f ], where t f is the duration of the anneal. The schedules need not be monotonic, and our analysis thus includes "reverse annealing" [40][41][42][43][44] as a special case. The TLS can be a single qubit or the two lowest energy levels of a multi-qubit system separated by a large gap from the rest of the spectrum. Key to our analysis is a series of transformations designed to arrive at a conveniently reparametrized interaction picture. First, we rewrite Eq. (1) in the form where A(s) = 2a(t)/E 0 and B(s) = 2b(t)/E 0 are dimensionless schedules parametrized by the dimensionless time s = t/t f , and E 0 > 0 is the energy scale of the Hamiltonian. We have cyclically permuted the Pauli matrices for later convenience. The ground states of H S (0) and H S (1) are |0 and |−i , respectively. Second, we parametrize the annealing schedules in the angular form where θ(0) = 0 and θ(1) = π/2. Under this parametrization the eigenvalues of H S (s) are ±E 0 Ω(s)/2, so the gap is ∆(s) = E 0 Ω(s). Thus, any non-trivial timedependence of the gap is encoded in the time-dependence of Ω(s), which we refer to as the dimensionless gap. The quantity is the cumulative dimensionless gap. Third, changing variables from s to τ to absorb Ω(s), the system satisfies the Schrödinger equation (we work in = 1 units throughout). The Hamiltonian is diagonalized at each instant by the rotation R X (θ) = e −iθX/2 . Thus, fourth, we change into the adiabatic frame [45,46] with |ψ ad = R X (θ) |ψ , yielding: (6) We call dθ dτ the angular progression of the anneal. Finally, we transform into the interaction picture with respect to the free Hamiltonian H 0 = −E 0 t f Z/2 and its propagator U 0 (τ ) = e −iH0τ . Letting S ± = (X ± iY )/2 denote the spin raising and lowering operators we have X I (τ ) = U † 0 (τ )XU 0 (τ ) = e −iE0t f τ S + + h.c., and obtain where |ψ I = U † 0 |ψ ad and λ(τ ) = 1 2 dθ dτ . Therefore, we see that in this adiabatic interaction picture the dynamics of the annealed TLS is a rotation about the time-dependent X I axis with a rate equal to the angular progression.

B. Magnus expansion
The corresponding time-ordered propagator U I (τ ) = T + e −i τ 0 dτ HI(τ ) can be calculated in time-dependent perturbation theory using the Magnus expansion (reviewed in Appendix A) for the Hermitian operator K (N ) (τ ) = N n=1 K n (τ ). The resulting U (N ) I (τ ) = exp −iK (N ) (τ ) converges to U I (τ ) uniformly with grow-ing N , and is unitary at all orders [47]. To first order: where To systematically go beyond first order we note that the K n (τ ) are nth order nested commutators, and hence closure of the su(2) Lie algebra guarantees that at all orders is a unit vector, and σ = (X, Y, Z). It thus follows that We will be concerned primarily with the probability of remaining in the ground state at the final time, denoted p 0←0 . Since |ψ I (s) = U † 0 (τ (s))R X (θ(s)) |ψ(s) , we have |ψ I (0) = |0 and |ψ I (1) = −i |0 . Thus, to N th order: where the states |0 and |1 are the initial ground and excited states, and where τ f ≡ τ (1). To first order we find (see App. A for the explicit form of U (1) ): This conceptually elegant result already indicates that quite generally one may expect the ground state probability to oscillate as a function of the anneal time t f , before the adiabatic limit sets in, a conclusion also reached in Ref. [24] on the basis of either a large-gap (near-adiabatic limit) or very small gap (stationary phase approximation) assumption. Our analysis applies for arbitrary gaps.

C. LZ problem (linear schedule)
Let us first consider the simplest annealing schedule, namely a linear interpolation of the type considered in the original LZ problem [27,28]: A(s) = 1 − s and B(s) = s. To evaluate Eq. (9) we can change the integration variable to s and approximate τ (s) ≈ τ f s in the exponent, yielding φ τ f (ω) = 1 2 1 0 ds 1 s 2 +(1−s) 2 e −iωτ f s for the firstorder Magnus expansion. We compare this to the numerically exact solution in Fig. 1, which shows remarkably good agreement. The simplicity of our Magnus expansion approach should be contrasted with the analytical solution for linear schedules in terms of parabolic cylinder functions [35]. Also notable is that while a quantum interference pattern is visible, the oscillations are very weak and not controllable (see the insert of Fig. 1). This Here and in other plots we use parameters compatible with quantum annealing using flux qubits [4][5][6][7]. Also shown is the prediction of a simplified double-slit type analysis (dashed, red). Both the latter and the first order Magnus expansion result are in excellent agreement with the numerically exact solution. The effect of strong dephasing in the instantaneous energy eigenbasis is shown as well (dashed, black), obtained using a phenomenological noise model with dephasing parameter Γ described in Appendix D. In this case the interference oscillations are strongly damped. motivates us to introduce new schedules with strong and controllable quantum interference.

D. Strong quantum interference pattern via Gaussian angular progression
Our goal is to identify a family of annealing schedules that generate strong interference between the paths leading to the final ground state, such that "giant" oscillations of the ground state probability can be observed. Therefore we now introduce Gaussian angular progressions.
Suppose that the angular progression is two-step Gaussian, namely, a sum of two Gaussians centered at τ f /2±µ (with µ < τ f /2): Note that τ f 0 dτ dθ dτ = θ(1) − θ(0) = π 2 , which fixes c. If we assume that α 1 then we may approximate to first order the ground state probability is then The ground state probability thus approaches its adiabatic limit of 1 on a timescale of t ad (set by the Gaussian width), while undergoing damped oscillations with a period of t coh . The oscillations are overdamped when t ad < t coh . In particular, a single Gaussian (µ = 0) can thus not give rise to oscillations.
We plot the ground state probability p G (t f ) ≡ p 0←0 in Fig. 1, for a two-step Gaussian progression with parameters chosen to represent the underdamped case; the associated annealing schedules are shown in Fig. 2(a). The amplitude of the resulting pre-adiabatic oscillations seen in Fig. 1 is, as desired, much larger than that associated with the linear schedule. The accuracy of the first-order Magnus expansion is again striking, especially given its simplicity compared to the analytical solution approaches [37][38][39]. We give a bound on the first-order Magnus expansion approximation error in Appendix B.

E. Physical origin of the oscillations
What is the origin of the oscillations? The answer is an interference effect between the two paths created by the two-step schedule, which enforces a double-slit or an unbalanced Mach-Zender interferometer scenario, with π/4 beam-splitters: see Fig. 2(b). The first step is a perturbation that generates amplitude in the excited state, while the second step allows for some of this amplitude to recombine with the ground state. The relative phase between the two paths is ξ = E 0 t f s+ s− Ω(s ) ds , which results in oscillations. In Appendix C we derive this result via a simple interferometer-type model that predicts the curve marked DS Γ = 0 in Fig. 1, which is in excellent agreement with the numerically exact result.
A natural question is whether the observation of interference oscillations as a function of t f implies the existence of quantum coherence in the computational basis at t f . We give a formal proof that the answer is affirmative in Appendix D. An illustration is given in Fig. 1, for the case of dephasing in the instantaneous energy eigenbasis, which is equivalent to performing a measurement in this basis between the two Gaussian steps. The final ground state probability is then the sum of classical conditional probabilities through each beam-splitter, and as expected, the oscillations disappear.

F. Role of the angular progression
We emphasize that the angular progression is the sole quantity needed to determine the ground state probability, per Eqs. (9) and (12). In particular, per Eq. (15), any transformation of A(s), B(s) and Ω(s) that leaves dθ dτ invariant will not affect P G in the closed-system setting.
Note, furthermore, that specifying the angular progression does not uniquely determine the annealing schedules A(s) and B(s). This is advantageous for practical pur-poses, since such schedules are typically implemented via arbitrary waveform generators (AWGs) with bandwidth constraints that can be incorporated into the schedule design process. To determine these schedules we need to specify the dimensionless gap Ω(s) and the angular progression dθ dτ . We can determine τ (s) by solving the differential equation dτ ds = Ω(s) subject to the boundary condition τ (0) = 0. Then θ(s) can be determined by solving the differential equation subject to appropriate boundary conditions. Together, Ω(s) and θ(s) determine the annealing schedules A(s) and B(s) via Eq. (3). In the two-step Gaussian case this means integrating Eq. (13), which, for a constant gap, yields θ(s) as a sum of erf functions. A particularly interesting example of a dimensionless gap schedule is one that represents the presence of two avoided level crossings, a significant feature of the glued trees problem [13]. An example is shown in Fig. 2(a), representing an example of the procedure outlined above for numerical determination of the schedule. It is clear from Eq. (15) that the main contribution to the angular progression is the near-vanishing of the gap. In contrast, when Ω(s) is constant, the main contribution to the angular progression is the suddenness of the schedule, i.e., a large A (s) or B (s).

III. OPEN SYSTEM ANALYSIS AND RESULTS
While a phenomenological model of dephasing in the instantaneous energy eigenbasis already shows clearly how the interference pattern disappears under decoherence ( Fig. 1 and Appendix D), this is not a realistic model of decoherence. We thus examine the effect of coupling the TLS to a thermal environment that corresponds more closely to experiments, e.g., with superconducting flux qubits.
We consider a dephasing model wherein the total system-bath where B is the dimensionless bath operator in the systembath interaction, H S (t) is given in Eq. (2), H B is the bath Hamiltonian, and g is the coupling strength with units of energy. We assume a separable initial state ρ S (0) ⊗ ρ B , with ρ B = exp(−βH B )/Z the Gibbs state of the bath at inverse temperature β and partition function Z = Tr[exp(−βH B )]. We transform to the interaction picture with respect to H B , so that H → The same series of transformations as those leading to Eq. (6) can be summarized as: After the final transformation to the H 0 -interaction picture, the total Hamiltonian replacing H I (τ ) in Eq. (7) becomes where µ = (sin φ cos θ, cos φ cos θ, sin θ) is a unit vector in polar coordinates, with φ(s) ≡ −E 0 t f τ (s), and henceforth the dot denotes d ds .
A. Redfield master equation in the adiabatic interaction picture The time-convolutionless (TCL) expansion [48] provides a convenient and systematic way to derive master equations (MEs) without requiring an adiabatic or Markovian approximation. With the detailed derivation given in Appendix E, the 2nd order TCL (TCL2) ME in the adiabatic-frame can be written as: and is the bath correlation function. We assume that the bath is Ohmic with spectral density J(ω) = ηωe −ω/ωc . To ensure the validity of the TCL2 approximation-which is also known as the Redfield ME-we derive a general error bound in Appendix F, and apply this bound to the Ohmic case. We find the condition t f β g 2 η , which is always satisfied in our simulations.

B. Rotating Wave Approximation (RWA)
In general, the Redfield ME (18) does not generate a completely positive map, which can result in non-sensical results such as negative probabilities [49,50]. Although this is not necessary for complete positivity [51], a further rotating wave approximation (RWA) is usually performed. The resulting Lindblad-type ME also lends itself to a simpler physical interpretation. As detailed in Appendix G, this leads tȯ where ρ ab = a|ρ S |b , all quantities except g, t f and β are s-dependent, and the effective dephasing and thermalization rates γ d and γ t , respectively, and the basis {|a , |b }, are given by Here | ± (s) = U † 0 (s) |± are the instantaneous eigenvectors of H I (s). The Lamb shift is: The functions γ(ω)/2 and S(ω) are the real and imaginary parts of the one-sided Fourier transform of the bath correlation function, and are implicitly β-dependent (see Appendix G, where we also discuss the validity conditions for the RWA).

C. Numerical results
The numerical solutions of Eqs. (18) and (20) are shown in Fig. 3 for the two-step Gaussian schedule with parameters as in Fig. 1 and for the gap schedule plotted in Fig. 2(a). The main message conveyed by this figure is that oscillations are visible over a wide range (an order of magnitude) of temperatures and system-bath coupling strengths. We also note that for these parameter values the Redfield ME produces physically valid solutions, despite the concerns about complete positivity mentioned above. The Redfield ME results in consistently higher ground state probabilities than the RWA.
These numerical results are accurately reproduced in terms of a simple semi-empirical formula, also shown in Fig. 3, and derived in Appendix H: where P G (t f ) and P G (t f ) denote the open and closed system success probabilities, respectively, wherē is the average thermalization rate, and where is the ground state probability in the adiabatic limit, given by the thermal equilibrium value associated with H S (1) [Eq. (2)]. As seen in Fig. 3, the agreement is excellent with both the RWA result when we use P E (0) = 1/2 (the infinite temperature limit), and with the TCL2 results when we use P E (β) and fit β; we find that the fitted β is consistently slightly lower than the actual β values used in our simulations.

IV. DISCUSSION AND CONCLUSIONS
We have proposed a double-slit approach to quantum annealing experiments, exhibiting "giant" interference patterns, motivated by the role of coherent diabatic evolution in enabling quantum speedups. Our analytical approach based on a simple time-dependent expansion in the adiabatic interaction picture accurately describes the associated dynamics. The experimental observation of such interference oscillations then becomes a clear and easily testable signature of coherence in the instantaneous energy eigenbasis. The test is simple in principle: it involves a quantum annealing protocol that employs the proposed schedules, with a measurement of only the ground state population as a function of the anneal time t f . When the relative phase between the upper and lower paths to the ground state is randomized, the interference effect is weakened.
To explain these results we proposed an effective model that accurately explains the interference oscillations in terms of a few simple parameters. Namely, upon replacing P G (t f ) in Eq. (23) by p (1) 0←0 (t f ) as given in Eq. (14a), the three timescales t coh , t ad , and 1/γ d respectively characterize the oscillation period, Gaussian damping due to approach to the adiabatic limit, and exponential damping due to coupling to the thermal bath. We expressed all three timescales in terms of the input physical parameters of the problem [Eqs. (14b) and (24)], and together they completely characterize the oscillations and their damping.
We expect that an experimental test of our "doubleslit" proposal will reveal the predicted interference oscillations for qubits that are sufficiently coherent, such as aluminum-based flux qubits [5][6][7], Rydberg atoms [52,53], or trapped ions [54,55]. Such an experiment can be viewed as a necessary condition for quantum annealing implementations of algorithms exhibiting a quantum speedup, e.g., the glued trees problem [13], which rely on coherence between energy eigenstates. It appears relevant (if not essential) to use such coherence in order to bypass the common objection that stoquastic quantum annealing or adiabatic quantum computing are subject to, which is that they can be efficiently simulated using the quantum Monte Carlo algorithm when restricted to ground-state evolution (with some known exceptions [56,57]), due to the absence of a sign problem [58,59]. Therefore an experimental observation of the quantum interference pattern predicted here will bolster our confidence in the abilities of coherent quantum annealers to one day deliver a quantum speedup.

ACKNOWLEDGMENTS
We are grateful to L. Campos-Venuti, L. Fry-Bouriaux, M. Khezri, J. Mozgunov, and P. Warburton for insightful comments and discussions. We used the Julia programming language [60] and the DifferentialEquations.jl package [61] for some of the numerical calculations reported in this work. The research is based upon work (partially) supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. We repeatedly use the following elementary identity for su(2) angular momentum operators: Note that the Pauli matrices are related via J i = σ i /2, i ∈ {x, y, z}.
Let us denote the solution of the adiabatic frame Hamiltonian given in Eq. (6) by U ad (τ ). The adiabatic interaction picture propagator, the solution of Eq. (7), can be computed using the Dyson series expansion: (A4) The first few terms are given by Using U where we wrote φ as a shorthand for φ τ (E 0 t f ), and where ϕ = arg(φ). This directly results in Eq. (12).
To compute the second order Magnus term we use X I (τ ) = e −iE0t f τ S + + h.c. for the commutation relation so that Appendix B: Error Analysis of Gaussian Angular Progression Schedules

Extension into full Fourier integrals
We discuss the general Gaussian angular progression where ψ is the Bloch sphere rotation angle. In the main text we assumed that we can perform a full Fourier transform (i.e., integration limits extended to ±∞) to find and thus arrive at the first order Magnus term We now show that the assumption of a full Fourier transform results in an exponentially small error in ατ * where Let I be the finite time integral with 0 < µ < τ f , and let F be the full Fourier integral and define the error where Thus, in terms of the standard normal cumulative density function Φ G (x) = 1 where erfc(x) = 1 − erf(x) is the complementary error function, and we have set u 1 ≡ αµ, and u 2 ≡ α(τ f − µ). The complementary error function is known to have exponentially small bounds [62]. We can quickly derive an even tighter bound by writing the error in terms of the Faddeeva function w(z) which is real and positive for imaginary z. If Im z > 0, the Faddeeva function has the integral representation [63, Eq. 7.7.2] and as 1/|z − t| ≤ 1/ Im z, With this bound on the Faddeeva function, it is straightforward to obtain the error bound If ατ * > √ π ≈ 1.77, then ≤ 0.014. If ατ * > 2 √ π ≈ 3.5, then with ≤ e −4π /2π, the percent error in assuming a full Fourier transform is under 0.6 parts per million. Extending the limits of integration will not result in an appreciable error if ατ * 2. In other words, if α 2 = 1/2σ 2 , where σ 2 is the variance of the normal distribution, the interval [0, τ f ] should contain the confidence interval of at least 2 √ 2σ ≈ 2.8σ.

Second Order Term of the Magnus Expansion
Evaluating the second order of the Magnus expansion [Eq. (A8)] yields: The integral is antisymmetric under the exchange of τ 1 and τ 2 due to the sine, so we can extend the time-ordered integrals into the whole square domain as For large α we can let A extend over the entire plane. The error bound due to extending the integration limits can be found straightforwardly: the square [0, τ f ]×[0, τ f ] contains the circle C centered at (µ, µ) with radius τ * , so the region R 2 − C contains a probability mass of e −(ατ * ) 2 (from a 2D Gaussian distribution). Since |sin(ω|τ 2 − τ 1 |)| < 1, the error in extending the region of integration is therefore bounded by Thus, up to an error of 2 , we may write where T 2 and T 1 are independent Gaussian random variables with variance 1/(2α 2 ). We can perform a change of variables into T + = T 2 + T 1 and T − = T 2 − T 1 , which are independent random variables with sum and difference means 2µ and 0 respectively, and both with variance 1/α 2 . Finally, the random variable |T − | is known to be distributed according to the folded-normal distribution centered at 0 (i.e. the half-normal distribution). Thus, the expectation value is precisely the imaginary part of the folded-normal characteristic function [64]: where the parent normal distribution has a mean µ − = 0. In this case, the characteristic function simplifies to the Faddeeva function where again r = ω/2α. From Eq. (B11), we see that where D(z) is the Dawson function Thus, as sin(ω|T − |) = Im(g t− (ω)), the second order term in the Magnus propagator is

Magnus Expansion Convergence and Error Bounds
Let S ≥ K 1 be a bound on the operator norm of the first order term in the Magnus expansion. A sufficient condition for the convergence of the Magnus expansion is that [65] S ≤ ξ = 1.08686870 . (B25) With a Gaussian angular progression as given in Eq. (B3), and noting that e −iωµ S + + e iωµ S − ≤ 1, it is then sufficient that for the Magnus expansion to be convergent. This means that ψ ≤ 2ξ, and in particular the physically relevant range ψ ∈ [0, π/2] (π/2 represents a balanced beamsplitter, and ψ > π/2 is equivalent to π − ψ) is within the convergence radius.
If K (n) is the nth order truncation of the Magnus expansion, the error in the truncation is given by where {b m } is a sequence defined in Ref. [65] via various recurrence relations. For ψ = π/2 and for ω/(2α) = 0, 0.5 and 1.0, the corresponding second order truncation errors are 0.25, 0.1, and 0.008.

Appendix C: Double-slit interpretation
Having derived the adiabatic frame Hamiltonian given in Eq. (6) H ad (τ ) = 1 2 we see that the angular progression dθ dτ of an annealing schedule is the perturbation that causes transitions between the two levels of the system. While this perturbation is steady and small in the case of a linear schedule, Gaussian schedules in which the perturbation is localized suggest an appealing physical picture similar to a double-slit or interferometer model.

Single Gaussian step
Let us first consider a single Gaussian step, which Eq. (13) reduces to when µ = 0, c = α √ π/2. Under the same assumptions as those leading to Eq. (14), we then find φ τ f (ω) = π 4 e −iωτ f /2 e −(t f /t ad ) 2 , with ω = E 0 t f . Thus, Eq. (A6) gives us the first order Magnus expansion propagator in the interaction picture with and ϕ = E 0 t f τ f /2. The X-rotation matrix in Eq. (A6c) thus becomes: (C3) with the superscript G serving as a reminder that this is the Gaussian step case. Now let us suppose that the Gaussian profile is narrow: α E 0 t f , or equivalently t ad t f . The perturbation is then sudden relative to the adiabatic timescale, and acts like a beamsplitter in a Mach-Zehnder (MZ) interferometer [31]. In this limit |φ| ≈ π/4 and Eq. (A6c) gives Recall that in the adiabatic interaction picture |ψ I (0) = |0 . Thus, the first phase factor e −iϕZ has no effect, and we can picture a process by which the ground state |0 is instantly split into an equal superposition 1 √ 2 (|0 − i |1 ) by the "Mach-Zender" matrix M G π/2 . These two states are then propagated freely by U † 0 (τ f ) = e i(E0t f τ f /2)Z , so they accumulate a relative phase of ie iE0t f τ f . For a single Gaussian, interference due to this phase difference is clearly not picked up via a Z basis measurement.

Two Gaussian steps: indirect derivation of the interferometer model in the narrow Gaussian limit
If instead we consider a two-step Gaussian schedule [Eq. (13)], then as we already found before Eq. (14), φ τ f (ω) = π 4 e −iωτ f /2 e −(t f /t ad ) 2 cos(µω), with ω = E 0 t f . Eq. (A6) now gives us the first order Magnus expansion propagator in the interaction picture with |φ| = π 4 | cos(µE 0 t f )|e −(t f /t ad ) 2 and again ϕ = E 0 t f τ f /2. 1 Let us now derive an equivalent MZ interferometer model. On the one hand, we already know from Eq. (12) that p This function has a quasiperiod (the distance between consecutive maxima) of π/(µE 0 ), a minimum of cos 2 (π/4) = 1/2 at t f = 0, and a maximum of 1.
On the other hand, we may model the two-step narrow (α E 0 t f ) Gaussian schedule as two consecutive, localized (at τ f /2 ± µ) and non-overlapping (α 1/µ) "beam-splitter" steps, separated by a dimensionless time interval of 2µ. Each beam-splitter is of the form given in Eq. (C4), the only difference being that the first acts at τ f /2 − µ (preceded by free evolution) and the second acts at τ f /2 + µ (followed by free evolution). In between the beam-splitter action there is free evolution of duration 2µ. Ignoring the initial and final free evolutions (since the initial and final state we are interested are both |0 , which is invariant under U 0 ) we expect to be able to write the propagator as the following ansatz: where we left the angle ψ in the beam splitter matrix (C3) unspecified in order to determine it by matching to the properties of p (1) 0←0 = cos 2 (|φ|). Carrying out the matrix multiplication and computing the expectation value, we find (C7) In order for this to match Eq. (C5), we require a quasiperiod of π/(µE 0 ) (which is already the case), a minimum of 1/2 at t f = 0, and a maximum of 1. The latter two conditions force ψ = π/4. Therefore, considering Eq. (C6), we have shown that the two-step Gaussian model is equivalent (in the large α limit) to a MZ interferometer with two unbalanced 1 Note that without the exponential decay factor e −(t f /t ad ) 2 = e −(t f /t ad ) 2 the oscillations are completely undamped and the adiabatic limit is never reached. Thus it is clear that the finite width of the Gaussian steps is solely responsible for the onset of adiabaticity. beamsplitters, separated by free propagation of duration 2µ (the separation between the two Gaussians).
The double-slit (or MZ interferometer model) is remarkably accurate in terms of predicting the ground state probability. This is shown in Fig. 1, where we compare the numerically exact result and the solution of the simple interferometer model given by Eq. (C7). Namely, we use the interference model given in Eq. (C7), with ψ = π/4. To calculate the interference fringe, the position of each of the two Gaussians is given by s ± = (τ f /2 ± µ)/τ . The phase factor µE 0 t f , which only holds in the large α limit, is replaced by The reason for this replacement is given in the following, alternative and more direct derivation of the interferometer model.

Two Gaussian steps: direct derivation of the interferometer model
Given the two-step Gaussian schedule, Eq. (13), where τ ± = τ f /2 ± µ, we can split the unitary generated by the adiabatic frame Hamiltonian, Eq. (C1), into two parts: We now wish to apply the Magnus expansion separately to each of the unitaries U ad ( where, using Eq. (9), now For α 1 we may extend the limits of integration over the interval [0, τ f /2] to ±∞ without considering the second Gaussian step: where we used c = α √ π/4 as we found in the derivation of Eq. (14). We may thus write the explicit form of the interaction picture unitary as and the adiabatic frame unitary becomes: Repeating this calculation for the second adiabatic frame unitary U ad (τ f , Thus, Eq. (C9) becomes which describes an interferometer composed of two unbalanced (π/4) double beam-splitters, interrupted by free propagation of duration τ + − τ − (ignoring the initial and final phases).
The phase accumulated between |0 and |1 is solely determined by the free evolution in Eq. (C16), whose value is given by where in the second equality we used Eq. (4).
It is clear from Eq. (D11) that oscillations in the ground state probability P G (t f ), which are present for finite Γ, imply a non-vanishing Im(ρ 10 (t f )). Therefore we may conclude that the observation of interference oscillations in our proposed double-slit experiment are also evidence of coherence in the computational basis at t f . For finite Γ, such coherence vanishes only in the adiabatic limit. The different orders are called TCL2, TCL4, etc. We give details on the convergence criteria of this expansion in Appendix F.
To second order the TCL generator is: where ρ B is the initial state of the bath, and the joint initial state is assumed to be in the factorized form ρ S ⊗ ρ B . Note that the TCL2 approximation coincides with the Redfield master equation [48]. After transforming back to the Schrödinger frame with respect to H I (s) we obtain: where Λ(s) = and denote τ B ≡ τ B , which has a natural interpretation as the bath correlation time [66].
Note that where c 2 = O(1) is a constant arising from the number of terms in the TCL2 double commutator expression (E5). We can similarly estimate the magnitude of the TCL4 terms: where c 4 , c 4 = O(1) are constants arising from the number of terms in the TCL4 sum over multiple commutators and triple integral. We can bound the two integrals in Eq. (F3) as follows. Considering the first expression, we first make a change of variables as Because 1 ≥ s ≥ s 1 ≥ s 2 ≥ s 3 , the new integration limits can be obtained as The Jacobian |det J| = 1 in this case and the new integral becomes Now we make another change of variables, with The new integration limits can be obtained by The first line means that 2s−v ≥ u ≥ v, while the second line gives 0 ≤ v ≤ s since 0 ≤ x 1 ≤ s. Thus: where in the last inequality we used s ≤ 1.
The same can be done for the second integral in Eq. (F3): Combining these two bounds thus finally yields: In particular, to ensure the validity of the TCL2 approximation it should be the case that the TCL4 term is much smaller than TCL2, i.e.:

Ohmic bath case
Let us assume a spin-boson noise model, for which where b k is a bosonic annihilation operator for mode k with frequency ω k , and g k = gξ k is the associated systembath coupling strength, where ξ k is dimensionless and g has units of energy. A standard approach is to introduce a spectral density such that |g k | 2 → J(ω)dω. For an Ohmic bath we have where η is a parameter with dimensions of time squared. After transforming to the bath interaction picture and replacing t by s = t/t f to arrive atH SB (s), the bath correlation function for the Ohmic spectral density is an integral which may be evaluated explicitly in terms of the Polygamma function [66]. In particular, for large βω c and t f /β, the correlation function can be expanded as (F16) This form indicates a transition from a Markovian regime of purely exponential decay with a timescale of τ B ωc→∞ → β/(2π), followed by a non-Markovian regime of powerlaw decay with a timescale of τ M = 2β/ω c . The transition occurs at a time τ tr ≈ β ln(βω c ) [66]. In the Marko-vian limit ω c → ∞ we may thus replace Eq. (F16) by and hence the correlation function integral of Eq. (F1) becomes which replaces every factor of τ B /t f arising from the same integral in the bounds in the previous subsection. In particular, we now have the necessary con- Eq. (F12) can thus be rewritten in the Markovian Ohmic case as or For finite ω c , one can refine this bound by replacing Eq. (F1) with where s tr = τ tr /t f . For our purposes the bound (F19) suffices and is satisfied in all the numerical results presented in the main text. Namely, we have and where γ s (ω)/2 and S s (ω) are the real and imaginary parts of Γ s (ω). Explicitly [48]: Here P denotes the Cauchy principal value, and the s subscript is a reminder that t f has been factored out.
To perform the rotating wave approximation, let us first define the eigenspace projection operator of H I (s) as Π( (s)) = | (s) (s)| , where | (s) is an eigenstate of H I (s) with instantaneous energy (s). We can then define the operator is the dimensionless Bohr frequency, and the sum is over all pairs (s), (s) subject to the constraint (s) − (s) = ω(s). The interaction picture master equation (E4) can then be written to second order, with the TCL2 generator (E5) aṡ To obtain this master equation, we apply the standard Markovian approximation: change the integration variable s → s − s and replace the upper limit with ∞. The RWA consists of neglecting terms in Eq. (G7) for which ω = ω. A necessary condition for the validity of the RWA is [67]: which, unfortunately, is not always satisfied for the twostep Gaussian schedule (13) for s outside the Gaussian pulse region.