A double-slit proposal for quantum annealing

TL;DR

This paper proposes a double-slit experiment analogy for quantum annealing, demonstrating interference effects in the ground state probability that reveal coherence and are affected by thermal environments, linking theory to practical hardware tests.

Contribution

It introduces a novel interference-based test for coherence in quantum annealing and simplifies the analysis of the generalized Landau-Zener problem.

Findings

Oscillations in ground state probability indicate quantum coherence.

Thermal bath coupling dampens interference oscillations.

Theoretical tools connect quantum speedup mechanisms to experimental hardware.

Abstract

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.

Equations

$$H_S(t)=-a(t)X-b(t)Z,$$

$$H_S(s)=-\frac{1}{2}{E}_0[A(s)Z+B(s)Y],$$

$$A(s)=\Omega (s)\cos \theta (s),B(s)=\Omega (s)\sin \theta (s),$$

$$\tau (s)\equiv {\int }_0^s\text{differential}{s}^{\prime }\Omega (s^{\prime })$$

$$i\text{derivative}\tau |\psi ⟩=-\frac{1}{2}{E}_0{t}_f[\cos \theta (\tau )Z+\sin \theta (\tau )Y]|\psi ⟩$$

$$i\text{derivative}\tau |{\psi }_{\mathrm{ad}}⟩=H_{\mathrm{ad}}|{\psi }_{\mathrm{ad}}⟩,H_{\mathrm{ad}}(\tau )\equiv \frac{1}{2}\left(\text{derivative}\theta \tau X-E_0{t}_{f}Z\right).$$

$$i\text{derivative}\tau |{\psi }_{\mathrm{I}}⟩=H_{\mathrm{I}}(\tau )|{\psi }_{\mathrm{I}}⟩,H_{\mathrm{I}}(\tau )\equiv \lambda (\tau )X_{\mathrm{I}}(\tau ),$$

$$K_1(\tau )={\int }_0^{\tau }\text{differential}{\tau }_1{H}_{\mathrm{I}}({\tau }_1)={\varphi }_{\tau }(E_0{t}_f)S_{+}+\text{h.c.},$$

$${\varphi }_{\tau }(\omega )\equiv \frac{1}{2}{\int }_0^{\tau }\text{differential}{\tau }_1\text{derivative}\theta {\tau }_1{e}^{-i\omega {\tau }_1}.$$

$$U_{\mathrm{I}}^{(N)}(\tau )=I\cos {\eta }^{(N)}(\tau )-i{\hat{n}}^{(N)}(\tau )\cdot \vec{\sigma }\sin {\eta }^{(N)}(\tau ).$$

$$p_{0\leftarrow 0}^{(N)}$$

$$=|\cos {\eta }^{(N)}({\tau }_f)-i{n}_Z^{(n)}({\tau }_f)\sin {\eta }^{(N)}({\tau }_f){|}^2,$$

$$p_{0\leftarrow 0}^{(1)}=|⟨0|e^{-i|\varphi |X}|0⟩{|}^2={\cos }^2(|\varphi |),\varphi \equiv {\varphi }_{\tau }(E_0{t}_f).$$

$$\text{derivative}\theta \tau =c\left(e^{-[\alpha (\tau -({\tau }_f/2+\mu )){]}^2}+e^{-[\alpha (\tau -({\tau }_f/2-\mu )){]}^2}\right).$$

$$p_{0\leftarrow 0}^{(1)}$$

$$t_{\mathrm{ad}}$$

$$\text{derivative}\theta \tau (s)=\frac{B^{\prime }(s)A(s)-A^{\prime }(s)B(s)}{\Omega (s{)}^3},$$

$$\text{derivative}\theta s=\Omega (s)\text{evaluated}{\text{derivative}\theta \tau }_{\tau =\tau (s)},$$

$$H_{\mathrm{tot}}(s)=\frac{1}{2}\dot{\theta }(s)X_{\mathrm{I}}(s)+g{t}_f\vec{\mu }(s)\cdot \vec{\sigma }\otimes \tilde{B}(s),$$

$${\dot{\rho }}_S(s)$$

$$-(g{t}_f{)}^2\text{commutator}\vec{\mu }(s)\cdot \vec{\sigma }\Lambda (s){\rho }_S(s)+\text{h.c.},$$

$$\Lambda (s)={\int }_0^s\text{differential}{s}^{\prime }C(s,s^{\prime })U_{\mathrm{I}}(s,s^{\prime })\vec{\mu }(s^{\prime })U_{\mathrm{I}}^{†}(s,s^{\prime })\cdot \vec{\sigma },$$

$${\dot{\rho }}_S$$

$$\displaystyle\qquad-g^{2}t_{f}\gamma_{d}\big{(}\rho_{ba}\outerproduct{b}{a}+\rho_{ab}\outerproduct{a}{b}\big{)}$$

$$\displaystyle\qquad+g^{2}t_{f}\gamma_{t}\pqty{\rho_{aa}-e^{-\beta\Delta}\rho_{bb}}\big{(}\outerproduct{b}{b}-\outerproduct{a}{a}\big{)}\ ,$$

$$|a(s)⟩$$

$${\gamma }_d(s)$$

$$H_{\mathrm{LS}}(s)=g^2{t}_f(S\text{pqty}\Delta (s)\text{outerproduct}bb+S\text{pqty}-\Delta (s)\text{outerproduct}aa).$$

$$P^{\prime}_{G}(t_{f})=\bigg{(}P_{G}(t_{f})-\frac{1}{2}\bigg{)}e^{-\bar{\gamma}_{d}t_{f}}+P_{E}(\beta)$$

$${\bar{\gamma }}_d=g^2{\int }_0^1\text{differential}{s}^{\prime }{\gamma }_d(s^{\prime })$$

$$P_E(\beta )\equiv \frac{e^{\beta E_0/2}}{Z},Z=2\cosh (\beta E_0/2)$$

$$\exp (-i\phi J_x)J_z\exp (i\phi J_x)$$

$$U_{\mathrm{I}}(\tau )$$

Full text

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” Feynman et al. (1963). Here we propose a double-slit experiment for quantum annealing (QA). In analogy to Feynman’s particle-wave double-slit, 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 Kadowaki and Nishimori (1998); Das and Chakrabarti (2008), which has led to major efforts to build physical quantum annealers for the purpose of solving optimization and sampling problems Johnson et al. (2011); Weber et al. (2017); Quintana et al. (2017); Novikov et al. (2018), and significant debate as to whether quantum effects are at play in the performance of such devices Boixo et al. (2014); Shin et al. (2014). 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 Denchev et al. (2016) and entanglement is often viewed as a necessary condition which must be demonstrated Lanting et al. (2014); Albash et al. (2015), 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 Somma et al. (2012), and other examples are known where QA provides a speedup over classical simulated annealing Farhi et al. (2002); Crosson and Deng (2014); Muthukrishnan et al. (2016); Kong and Crosson (2015); Brady and van Dam (2016); Jiang et al. (2017). 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 del Campo (2013). When these transitions result in a coherent recombination of the ground state amplitude (a phenomenon known as a diabatic cascade Muthukrishnan et al. (2016); Brady and van Dam (2017)), the result is a wave-like interference pattern in the ground state probability as the anneal time is varied Wiebe and Babcock (2012); Wecker et al. (2016); Brady and van Dam (2018). 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 non-simple quantum problem” M.V. Berry (1995), a driven TLS near level crossings Grifoni and Hänggi (1998). The two-level paradigm was introduced long ago by Landau and Zener (LZ) L. D. Landau (1932); C. Zener (1932). The corresponding Hamiltonian for the generalized LZ problem is

[TABLE]

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 $-\infty$ to $\infty$. The problem has since been studied under numerous variations, including Landau-Zener-Stueckelberg interferometry where $b(t)$ is periodic Stueckelberg (1932); Ashhab (2017), the subject of various experiments Oliver et al. (2005); Sillanpää et al. (2006); Petta et al. (2010). Complete analytical solutions were limited until recently to certain particular functional forms of $b(t)$ with constant $a(t)$ Bambini and Berman (1981), a finite-range linear schedule for both $a(t)$ and $b(t)$ Vitanov and Garraway (1996), and periodic $a(t)$ and $b(t)$ Bezvershenko and Holod (2011). An analytical solution for general $b(t)$ but constant $a(t)$ was found in Ref. Barnes and Das Sarma (2012), which was then extended to general (but implicitly specified) $a(t)$ as well Barnes (2013); Messina and Nakazato (2014). 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 double-slit 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

II.1 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)\geq 0$ respectively decrease/increase to/from [math] with time $t\in[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” Perdomo-Ortiz et al. (2011); Chancellor (2017); King et al. (2018); Ohkuwa et al. (2018); Venturelli and Kondratyev (2018) 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

[TABLE]

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 $\ket{0}$ and $\ket{-i}$, respectively. Second, we parametrize the annealing schedules in the angular form

[TABLE]

where $\theta(0)=0$ and $\theta(1)=\pi/2$. Under this parametrization the eigenvalues of $H_{S}(s)$ are $\pm E_{0}\Omega(s)/2$, so the gap is $\Delta(s)=E_{0}\Omega(s)$. Thus, any non-trivial time-dependence of the gap is encoded in the time-dependence of $\Omega(s)$, which we refer to as the dimensionless gap. The quantity

[TABLE]

is the cumulative dimensionless gap. Third, changing variables from $s$ to $\tau$ to absorb $\Omega(s)$, the system satisfies the Schrödinger equation

[TABLE]

(we work in $\hbar=1$ units throughout). The Hamiltonian is diagonalized at each instant by the rotation $R_{X}(\theta)=e^{-i\theta X/2}$. Thus, fourth, we change into the adiabatic frame Klarsfeld and Oteo (1992); Nalbach (2014) with $\ket{\psi_{\mathrm{ad}}}=R_{X}(\theta)\ket{\psi}$, yielding:

[TABLE]

We call $\derivative{\theta}{\tau}$ 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}(\tau)=e^{-iH_{0}\tau}$. Letting $S_{\pm}=(X\pm iY)/2$ denote the spin raising and lowering operators we have $X_{\mathrm{I}}(\tau)=U_{0}^{\dagger}(\tau)XU_{0}(\tau)=e^{-iE_{0}t_{f}\tau}S_{+}+\textrm{h.c.}$, and obtain

[TABLE]

where $\ket{\psi_{\mathrm{I}}}=U_{0}^{{\dagger}}\ket{\psi_{\mathrm{ad}}}$ and $\lambda(\tau)=\frac{1}{2}\derivative{\theta}{\tau}$. Therefore, we see that in this adiabatic interaction picture the dynamics of the annealed TLS is a rotation about the time-dependent $X_{\mathrm{I}}$ axis with a rate equal to the angular progression.

II.2 Magnus expansion

The corresponding time-ordered propagator $U_{\mathrm{I}}(\tau)=T_{+}e^{-i\int_{0}^{\tau}\differential\tau^{\prime}H_{\mathrm{I}}(\tau^{\prime})}$ can be calculated in time-dependent perturbation theory using the Magnus expansion (reviewed in Appendix A) for the Hermitian operator $\mathcal{K}^{(N)}(\tau)=\sum_{n=1}^{N}K_{n}(\tau)$. The resulting $U_{\mathrm{I}}^{(N)}(\tau)=\exp[-i\mathcal{K}^{(N)}(\tau)]$ converges to $U_{\mathrm{I}}(\tau)$ uniformly with growing $N$, and is unitary at all orders Blanes et al. (2009). To first order:

[TABLE]

where

[TABLE]

To systematically go beyond first order we note that the $K_{n}(\tau)$ are $n$th order nested commutators, and hence closure of the $su(2)$ Lie algebra guarantees that at all orders $\mathcal{K}^{(N)}(\tau)=\eta^{(N)}(\tau)\hat{n}^{(N)}(\tau)\cdot\vec{\sigma}$, where $\eta^{(N)}(\tau)>0$, $\hat{n}^{(N)}(\tau)$ is a unit vector, and $\vec{\sigma}=(X,Y,Z)$. It thus follows that

[TABLE]

We will be concerned primarily with the probability of remaining in the ground state at the final time, denoted $p_{0\leftarrow 0}$. Since $\ket{\psi_{\mathrm{I}}(s)}=U_{0}^{{\dagger}}(\tau(s))R_{X}(\theta(s))\ket{\psi(s)}$, we have $\ket{\psi_{\mathrm{I}}(0)}=\ket{0}$ and $\ket{\psi_{\mathrm{I}}(1)}=-i\ket{0}$. Thus, to $N$th order:

[TABLE]

where the states $\ket{0}$ and $\ket{1}$ are the initial ground and excited states, and where $\tau_{f}\equiv\tau(1)$. To first order we find (see App. A for the explicit form of $U^{(1)}$):

[TABLE]

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. Brady and van Dam (2018) 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.

II.3 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 L. D. Landau (1932); C. Zener (1932): $A(s)=1-s$ and $B(s)=s$. To evaluate Eq. (9) we can change the integration variable to $s$ and approximate $\tau(s)\approx\tau_{f}s$ in the exponent, yielding $\phi_{\tau_{f}}(\omega)=\frac{1}{2}\int_{0}^{1}\differential s\frac{1}{s^{2}+(1-s)^{2}}e^{-i\omega\tau_{f}s}$ for the first-order 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 Vitanov and Garraway (1996). 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 motivates us to introduce new schedules with strong and controllable quantum interference.

II.4 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 $\tau_{f}/2\pm\mu$ (with $\mu<\tau_{f}/2$):

[TABLE]

Note that $\int_{0}^{\tau_{f}}d\tau\derivative{\theta}{\tau}=\theta(1)-\theta(0)=\frac{\pi}{2}$, which fixes $c$. If we assume that $\alpha\gg 1$ then we may approximate $\int_{0}^{\tau_{f}}$ by $\int_{-\infty}^{\infty}$ (we bound the approximation error in Appendix B). Thus $c=\alpha\sqrt{\pi}/4$ and Eq. (9) yields $\phi_{\tau_{f}}(\omega)=\frac{\pi}{4}e^{-i\omega\tau_{f}/2}e^{-[\omega/(2\alpha)]^{2}}\cos(\mu\omega)$. Using Eq. (12), to first order the ground state probability is then

[TABLE]

The ground state probability thus approaches its adiabatic limit of $1$ on a timescale of $t_{\mathrm{ad}}$ (set by the Gaussian width), while undergoing damped oscillations with a period of $t_{\mathrm{coh}}$. The oscillations are overdamped when $t_{\mathrm{ad}}<t_{\mathrm{coh}}$. In particular, a single Gaussian ($\mu=0$) can thus not give rise to oscillations.

We plot the ground state probability $p_{G}(t_{f})\equiv p_{0\leftarrow 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. 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 Barnes and Das Sarma (2012); Barnes (2013); Messina and Nakazato (2014). We give a bound on the first-order Magnus expansion approximation error in Appendix B.

II.5 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 $\pi/4$ beam-splitters: see Fig. 2. 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 $\xi=E_{0}t_{f}\int_{s_{-}}^{s_{+}}{\Omega(s^{\prime})}\differential{s^{\prime}}$, which results in oscillations. In Appendix C we derive this result via a simple interferometer-type model that predicts the curve marked DS $\Gamma=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.

II.6 Role of the angular progression

We emphasize that the angular progression

[TABLE]

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 $\Omega(s)$ that leaves $\derivative{\theta}{\tau}$ 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 purposes, 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 $\Omega(s)$ and the angular progression $\derivative{\theta}{\tau}$. We can determine $\tau(s)$ by solving the differential equation $\derivative{\tau}{s}=\Omega(s)$ subject to the boundary condition $\tau(0)=0$. Then $\theta(s)$ can be determined by solving the differential equation

[TABLE]

subject to appropriate boundary conditions. Together, $\Omega(s)$ and $\theta(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 $\theta(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 Somma et al. (2012). An example is shown in Fig. 2, 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 $\Omega(s)$ is constant, the main contribution to the angular progression is the suddenness of the schedule, i.e., a large $A^{\prime}(s)$ or $B^{\prime}(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 Hamiltonian is $H=H_{S}(t)+H_{B}+gY\otimes B$, where $B$ is the dimensionless bath operator in the system-bath 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 $\rho_{S}(0)\otimes\rho_{B}$, with $\rho_{B}=\exp(-\beta H_{B})/\mathrm{Z}$ the Gibbs state of the bath at inverse temperature $\beta$ and partition function $\mathrm{Z}=\mathrm{Tr}[\exp(-\beta H_{B})]$. We transform to the interaction picture with respect to $H_{B}$, so that $H\mapsto\tilde{H}(t)=H_{S}(t)+gY\otimes\tilde{B}(t)$, with $\tilde{B}(t)=U_{B}^{\dagger}(t)BU_{B}(t)$, and $U_{B}(t)=e^{-itH_{B}}$. The same series of transformations as those leading to Eq. (6) can be summarized as: $Y\otimes\tilde{B}(t)\mapsto t_{f}Y\otimes\tilde{B}(s)\mapsto t_{f}R_{X}(\theta)YR_{X}(-\theta)\otimes\tilde{B}(s)=t_{f}[\cos(\theta)Y+\sin(\theta)Z]\otimes\tilde{B}(s)$. After the final transformation to the $H_{0}$-interaction picture, the total Hamiltonian replacing $H_{\mathrm{I}}(\tau)$ in Eq. (7) becomes

[TABLE]

where $\vec{\mu}=(\sin\phi\cos\theta,\cos\phi\cos\theta,\sin\theta)$ is a unit vector in polar coordinates, with $\phi(s)\equiv-E_{0}t_{f}\tau(s)$, and henceforth the dot denotes $\derivative{s}$.

III.1 Redfield master equation in the adiabatic interaction picture

The time-convolutionless (TCL) expansion Breuer and Petruccione (2002) 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:

[TABLE]

where

[TABLE]

and $C(s,s^{\prime})=\mathrm{Tr}[\tilde{B}(s)\tilde{B}\pqty{s^{\prime}}\rho_{B}]=C^{*}(s^{\prime},s)$ is the bath correlation function. We assume that the bath is Ohmic with spectral density $J\pqty{\omega}=\eta\omega e^{-{\omega}/{\omega_{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}\ll\frac{\beta}{g^{2}\eta}$, which is always satisfied in our simulations.

III.2 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 Gaspard and Nagaoka (1999); Whitney (2008). Although this is not necessary for complete positivity Majenz et al. (2013), 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 to

[TABLE]

where $\rho_{ab}=\matrixelement{a}{\rho_{S}}{b}$, all quantities except $g$, $t_{f}$ and $\beta$ are $s$-dependent, and the effective dephasing and thermalization rates $\gamma_{d}$ and $\gamma_{t}$, respectively, and the basis $\{\ket{a},\ket{b}\}$, are given by

[TABLE]

Here $\ket{\epsilon_{\pm}(s)}=U_{0}^{\dagger}(s)\ket{\pm}$ are the instantaneous eigenvectors of $H_{\mathrm{I}}(s)$. The Lamb shift is:

[TABLE]

The functions $\gamma\pqty{\omega}/2$ and $S\pqty{\omega}$ are the real and imaginary parts of the one-sided Fourier transform of the bath correlation function, and are implicitly $\beta$-dependent (see Appendix G, where we also discuss the validity conditions for the RWA).

III.3 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. 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:

[TABLE]

where $P^{\prime}_{G}\pqty{t_{f}}$ and $P_{G}\pqty{t_{f}}$ denote the open and closed system success probabilities, respectively, where

[TABLE]

is the average thermalization rate, and where

[TABLE]

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}(\beta)$ and fit $\beta$; we find that the fitted $\beta$ is consistently slightly lower than the actual $\beta$ 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\leftarrow 0}(t_{f})$ as given in Eq. (14a), the three timescales $t_{\mathrm{coh}}$, $t_{\mathrm{ad}}$, and $1/\bar{\gamma}_{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 “double-slit” proposal will reveal the predicted interference oscillations for qubits that are sufficiently coherent, such as aluminum-based flux qubits Weber et al. (2017); Quintana et al. (2017); Novikov et al. (2018), Rydberg atoms Glaetzle et al. (2017); Pichler et al. (2018), or trapped ions Graß et al. (2016); Zhang et al. (2017). 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 Somma et al. (2012), 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 Hastings and Freedman (2013); Jarret et al. (2016)), due to the absence of a sign problem Albash and Lidar (2018); Crosson and Harrow (2018). 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.

Acknowledgements.

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 Bezanson et al. (2017) and the DifferentialEquations.jl package Rackauckas and Nie (2017) 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.

Appendix A Dyson and Magnus series

We repeatedly use the following elementary identity for $su(2)$ angular momentum operators:

[TABLE]

Note that the Pauli matrices are related via $J_{i}=\sigma_{i}/2$, $i\in\{x,y,z\}$.

Let us denote the solution of the adiabatic frame Hamiltonian given in Eq. (6) by $U_{\mathrm{ad}}(\tau)$. The adiabatic interaction picture propagator,

[TABLE]

the solution of Eq. (7), can be computed using the Dyson series expansion:

[TABLE]

Note that each term in the Dyson series contributes to the ground state amplitude if and only if it is an even power, and likewise to the excitation amplitude if and only if it is an odd power. Consequently, the amplitudes calculated from the Dyson series may not be unitary to a desired precision until the terms are calculated to a high enough order. For this reason we prefer the Magnus expansion Blanes et al. (2009), for which

[TABLE]

The first few terms are given by

[TABLE]

Using $U_{\mathrm{I}}^{(N)}(\tau)=\exp[-i\mathcal{K}^{(N)}(\tau)]$ and Eq. (8) we thus find

[TABLE]

where we wrote $\phi$ as a shorthand for $\phi_{\tau}(E_{0}t_{f})$, and where $\varphi=\arg(\phi)$. This directly results in Eq. (12).

To compute the second order Magnus term we use $X_{\mathrm{I}}(\tau)=e^{-iE_{0}t_{f}\tau}S_{+}+\textrm{h.c.}$ for the commutation relation

[TABLE]

so that

[TABLE]

Appendix B Error Analysis of Gaussian Angular Progression Schedules

B.1 Extension into full Fourier integrals

We discuss the general Gaussian angular progression

[TABLE]

where $\psi$ 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 $\pm\infty$) to find

[TABLE]

and thus arrive at the first order Magnus term

[TABLE]

We now show that the assumption of a full Fourier transform results in an exponentially small error in $\alpha\tau^{*}$ where

[TABLE]

Let $I$ be the finite time integral

[TABLE]

with $0<\mu<\tau_{f}$, and let $F$ be the full Fourier integral

[TABLE]

and define the error

[TABLE]

where

[TABLE]

Thus, in terms of the standard normal cumulative density function $\Phi_{G}(x)=\frac{1}{2}\quantity(1+\erf[x/\sqrt{2}])$,

[TABLE]

where $\operatorname{erfc}(x)=1-\erf(x)$ is the complementary error function, and we have set $u_{1}\equiv\alpha\mu$, and $u_{2}\equiv\alpha(\tau_{f}-\mu)$. The complementary error function is known to have exponentially small bounds Chiani et al. (2003). We can quickly derive an even tighter bound by writing the error in terms of the Faddeeva function $w(z)$

[TABLE]

where

[TABLE]

which is real and positive for imaginary $z$. If $\imaginary z>0$, the Faddeeva function has the integral representation (nis, , Eq. 7.7.2)

[TABLE]

from which we note

[TABLE]

and as $1/\absolutevalue{z-t}\leq 1/\imaginary z$,

[TABLE]

With this bound on the Faddeeva function, it is straightforward to obtain the error bound

[TABLE]

If $\alpha\tau^{*}>\sqrt{\pi}\approx 1.77$, then $\epsilon\leq 0.014$. If $\alpha\tau^{*}>2\sqrt{\pi}\approx 3.5$, then with $\epsilon\leq e^{-4\pi}/2\pi$, 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 $\alpha\tau^{*}\gtrsim 2$. In other words, if $\alpha^{2}=1/2\sigma^{2}$, where $\sigma^{2}$ is the variance of the normal distribution, the interval $[0,\tau_{f}]$ should contain the confidence interval of at least $2\sqrt{2}\sigma\approx 2.8\sigma$.

B.2 Second Order Term of the Magnus Expansion

Evaluating the second order of the Magnus expansion [Eq. (33)] yields:

[TABLE]

The integral is antisymmetric under the exchange of $\tau_{1}$ and $\tau_{2}$ due to the sine, so we can extend the time-ordered integrals into the whole square domain as

[TABLE]

For large $\alpha$ 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,\tau_{f}]\times[0,\tau_{f}]$ contains the circle $C$ centered at $(\mu,\mu)$ with radius $\tau^{*}$, so the region $\mathbb{R}^{2}-C$ contains a probability mass of $e^{-(\alpha\tau^{*})^{2}}$ (from a 2D Gaussian distribution). Since $\absolutevalue{\sin(\omega\absolutevalue{\tau_{2}-\tau_{1}})}<1$, the error in extending the region of integration is therefore bounded by

[TABLE]

Thus, up to an error of $\epsilon_{2}$, we may write

[TABLE]

where $T_{2}$ and $T_{1}$ are independent Gaussian random variables with variance $1/(2\alpha^{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\mu$ and [math] respectively, and both with variance $1/\alpha^{2}$. Finally, the random variable $\absolutevalue{T_{-}}$ is known to be distributed according to the folded-normal distribution centered at [math] (i.e. the half-normal distribution). Thus, the expectation value is precisely the imaginary part of the folded-normal characteristic function Tsagris et al. (2014):

[TABLE]

where the parent normal distribution has a mean $\mu_{-}=0$. In this case, the characteristic function simplifies to the Faddeeva function

[TABLE]

where again $r=\omega/2\alpha$. From Eq. (44), we see that

[TABLE]

where $D(z)$ is the Dawson function

[TABLE]

Thus, as $\expectationvalue{\sin(\omega\absolutevalue{T_{-}})}=\imaginary(g_{t_{-}}(\omega))$, the second order term in the Magnus propagator is

[TABLE]

B.3 Magnus Expansion Convergence and Error Bounds

Let $S\geq\norm{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 Blanes et al. (1998)

[TABLE]

With a Gaussian angular progression as given in Eq. (36), and noting that $\norm{e^{-i\omega\mu}S_{+}+e^{i\omega\mu}S_{-}}\leq 1$, it is then sufficient that

[TABLE]

for the Magnus expansion to be convergent. This means that $\psi\leq 2\xi$, and in particular the physically relevant range $\psi\in[0,\pi/2]$ ($\pi/2$ represents a balanced beam-splitter, and $\psi>\pi/2$ is equivalent to $\pi-\psi$) is within the convergence radius.

If $\mathcal{K}^{(n)}$ is the $n$th order truncation of the Magnus expansion, the error in the truncation is given by

[TABLE]

where $\{b_{m}\}$ is a sequence defined in Ref. Blanes et al. (1998) via various recurrence relations. For $\psi=\pi/2$ and for $\omega/(2\alpha)=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)

[TABLE]

we see that the angular progression $\derivative{\theta}{\tau}$ 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.

C.1 Single Gaussian step

Let us first consider a single Gaussian step, which Eq. (13) reduces to when $\mu=0$, $c=\alpha\sqrt{\pi}/2$. Under the same assumptions as those leading to Eq. (14), we then find $\phi_{\tau_{f}}(\omega)=\frac{\pi}{4}e^{-i\omega\tau_{f}/2}e^{-(t_{f}/t_{\textrm{ad}})^{2}}$, with $\omega=E_{0}t_{f}$. Thus, Eq. (31) gives us the first order Magnus expansion propagator in the interaction picture with

[TABLE]

and $\varphi=E_{0}t_{f}\tau_{f}/2$. The $X$-rotation matrix in Eq. (31c) thus becomes:

[TABLE]

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: $\alpha\gg E_{0}t_{f}$, or equivalently $t_{\mathrm{ad}}\gg t_{f}$. The perturbation is then sudden relative to the adiabatic timescale, and acts like a beamsplitter in a Mach-Zehnder (MZ) interferometer Oliver et al. (2005). In this limit $|\phi|\approx\pi/4$ and Eq. (31c) gives

[TABLE]

Recall that in the adiabatic interaction picture $\ket{\psi_{\mathrm{I}}(0)}=\ket{0}$. Thus, the first phase factor $e^{-i\varphi Z}$ has no effect, and we can picture a process by which the ground state $\ket{0}$ is instantly split into an equal superposition $\frac{1}{\sqrt{2}}(\ket{0}-i\ket{1})$ by the “Mach-Zender” matrix $M^{\mathrm{G}}_{\pi/2}$. These two states are then propagated freely by $U^{\dagger}_{0}(\tau_{f})=e^{i(E_{0}t_{f}\tau_{f}/2)Z}$, so they accumulate a relative phase of $ie^{iE_{0}t_{f}\tau_{f}}$. For a single Gaussian, interference due to this phase difference is clearly not picked up via a $Z$ basis measurement.

C.2 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), $\phi_{\tau_{f}}(\omega)=\frac{\pi}{4}e^{-i\omega\tau_{f}/2}e^{-(t_{f}/t_{\textrm{ad}})^{2}}\cos(\mu\omega)$, with $\omega=E_{0}t_{f}$. Eq. (31) now gives us the first order Magnus expansion propagator in the interaction picture with $|\phi|=\frac{\pi}{4}|\cos(\mu E_{0}t_{f})|e^{-(t_{f}/t_{\textrm{ad}})^{2}}$ and again $\varphi=E_{0}t_{f}\tau_{f}/2$.111Note that without the exponential decay factor $e^{-(t_{f}/t_{\textrm{ad}})^{2}}=e^{-(t_{f}/t_{\mathrm{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.

Let us now derive an equivalent MZ interferometer model. On the one hand, we already know from Eq. (12) that $p^{(1)}_{0\leftarrow 0}=\cos^{2}(|\phi|)$, i.e.

[TABLE]

This function has a quasiperiod (the distance between consecutive maxima) of $\pi/(\mu E_{0})$, a minimum of $\cos^{2}(\pi/4)=1/2$ at $t_{f}=0$, and a maximum of $1$. On the other hand, we may model the two-step narrow ($\alpha\gg E_{0}t_{f}$) Gaussian schedule as two consecutive, localized (at $\tau_{f}/2\pm\mu$) and non-overlapping ($\alpha\gg 1/\mu$) “beam-splitter” steps, separated by a dimensionless time interval of $2\mu$. Each beam-splitter is of the form given in Eq. (64), the only difference being that the first acts at $\tau_{f}/2-\mu$ (preceded by free evolution) and the second acts at $\tau_{f}/2+\mu$ (followed by free evolution). In between the beam-splitter action there is free evolution of duration $2\mu$. Ignoring the initial and final free evolutions (since the initial and final state we are interested are both $\ket{0}$, which is invariant under $U_{0}$) we expect to be able to write the propagator as the following ansatz:

[TABLE]

where we left the angle $\psi$ in the beam splitter matrix (63) unspecified in order to determine it by matching to the properties of $p^{(1)}_{0\leftarrow 0}=\cos^{2}(|\phi|)$. Carrying out the matrix multiplication and computing the expectation value, we find

[TABLE]

In order for this to match Eq. (65), we require a quasiperiod of $\pi/(\mu 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 $\psi=\pi/4$.

Therefore, considering Eq. (66), we have shown that the two-step Gaussian model is equivalent (in the large $\alpha$ limit) to a MZ interferometer with two unbalanced beamsplitters, separated by free propagation of duration $2\mu$ (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. (67). Namely, we use the interference model given in Eq. (67), with $\psi=\pi/4$. To calculate the interference fringe, the position of each of the two Gaussians is given by $s_{\pm}=(\tau_{f}/2\pm\mu)/\tau$. The phase factor $\mu E_{0}t_{f}$, which only holds in the large $\alpha$ limit, is replaced by $E_{0}t_{f}[\tau(s_{+})-\tau(s_{-})]=E_{0}t_{f}\int_{s_{-}}^{s_{+}}\differential{s^{\prime}}\Omega(s^{\prime})$, where $\tau(s)$ is the cumulative dimensionless gap [Eq. (4)]. The reason for this replacement is given in the following, alternative and more direct derivation of the interferometer model.

C.3 Two Gaussian steps: direct derivation of the interferometer model

Given the two-step Gaussian schedule, Eq. (13),

[TABLE]

where $\tau_{\pm}=\tau_{f}/2\pm\mu$, we can split the unitary generated by the adiabatic frame Hamiltonian, Eq. (61), into two parts:

[TABLE]

We now wish to apply the Magnus expansion separately to each of the unitaries $U_{\mathrm{ad}}(\frac{\tau_{f}}{2},0)$ and $U_{\mathrm{ad}}(\tau_{f},\frac{\tau_{f}}{2})$. Consider $U_{\mathrm{ad}}(\frac{\tau_{f}}{2},0)$. Inverting Eq. (27a), the first order Magnus expansion [Eq. (31)] gives

[TABLE]

where, using Eq. (9), now

[TABLE]

For $\alpha\gg 1$ we may extend the limits of integration over the interval $[0,\tau_{f}/2]$ to $\pm\infty$ without considering the second Gaussian step:

[TABLE]

where we used $c=\alpha\sqrt{\pi}/4$ as we found in the derivation of Eq. (14). We may thus write the explicit form of the interaction picture unitary as

[TABLE]

and the adiabatic frame unitary becomes:

[TABLE]

Repeating this calculation for the second adiabatic frame unitary $U_{\mathrm{ad}}(\tau_{f},\frac{\tau_{f}}{2})$, we obtain

[TABLE]

Thus, Eq. (69) becomes

[TABLE]

which describes an interferometer composed of two unbalanced ($\pi/4$) double beam-splitters, interrupted by free propagation of duration $\tau_{+}-\tau_{-}$ (ignoring the initial and final phases).

The phase accumulated between $\ket{0}$ and $\ket{1}$ is solely determined by the free evolution in Eq. (76),

[TABLE]

whose value is given by

[TABLE]

where in the second equality we used Eq. (4).

Appendix D Interference oscillations in the double-slit experiment imply quantum coherence in the computational basis

Here we prove that coherence in the energy eigenbasis implies, in general, coherence in the computational basis.

Let $H(t)$ denote an arbitrary, time-dependent TLS Hamiltonian, with instantaneous energy eigenbasis $\{\ket{\epsilon_{i}(t)}\}$. The TLS density matrix can be written in this basis as

[TABLE]

Let us define “coherence” with respect to a given basis as the off-diagonal elements of the density matrix in the same basis. We can compute the coherence in the computational basis $\{\ket{0},\ket{1}\}$ via

[TABLE]

where $\epsilon_{ij}(t)=\outerproduct{\epsilon_{i}(t)}{\epsilon_{j}(t)}$. The two bases are related via a unitary rotation:

[TABLE]

so that Eq, (80) reduces to:

[TABLE]

where we used $\tilde{\rho}_{00}+\tilde{\rho}_{11}=1$ and $\tilde{\rho}_{01}=\tilde{\rho}_{10}^{*}$. Equation (82) can be further simplified using $(\tilde{\rho}_{00}-\frac{1}{2})\sin(2\theta)-\real(\tilde{\rho}_{10})\cos(2\theta)=C(\cos\varphi\sin(2\theta)-\sin\varphi\cos(2\theta))$, where

[TABLE]

Additionally, by making use of the trigonometric identity $\sin(2\theta-\varphi)=\sin 2\theta\cos\varphi-\sin\varphi\cos 2\theta$, Eq. (82) can be written as

[TABLE]

Since $C\sin(2\theta-\varphi)\in\mathds{R}$, it follows that $\imaginary(\tilde{\rho}_{10}(t))\neq 0$ implies $\matrixelement{0}{\rho(t)}{1}\neq 0$. Therefore we next establish that indeed, $\imaginary(\tilde{\rho}_{10}(t))\neq 0$ in our double-slit proposal.

Consider the the ground state just before the first beam-splitter,

[TABLE]

with $\varepsilon/(\tau_{+}-\tau_{-})\ll 1$. This state evolves through the double-beam-splitter region [recall Eq. (76)]:

[TABLE]

where $U_{0}$ is given in Eq. (77) and $M_{|\phi|}$ is given in Eq. (31d).

After passing through the first beam-splitter, the system density matrix in the energy eigenbasis becomes

[TABLE]

It is useful to include a simple model of decoherence between energy eigenstates during the time interval $[\tau_{-},\tau_{+}]$, complementary to our master equation treatment. We can do so by introducing a continuous dephasing channel. This damps the phases by the factor $e^{-\Gamma\Delta\tau}$, where $\Delta\tau=\tau_{+}-\tau_{-}=2\mu$, and $\Gamma>0$ is the dephasing rate. Right before the second beam-splitter, the system density matrix is then:

[TABLE]

After passing through the second beam-splitter, the state becomes $\rho(\tau_{+}+\varepsilon)=M_{|\phi|}\rho(\tau_{+}-\varepsilon)M^{\dagger}_{|\phi|}$. We find, after some algebra:

[TABLE]

We now note from Eq. (72b) that $\absolutevalue{\phi}=\frac{\pi}{8}e^{-(t_{f}/t_{\mathrm{ad}})^{2}}$. Therefore we may conclude that $\imaginary(\tilde{\rho}_{10}(t_{f}))>0$, and $\imaginary(\tilde{\rho}_{10})\rightarrow 0$ only in the adiabatic limit ($t_{f}\gg t_{\mathrm{ad}}$, which implies $\absolutevalue{\phi}\rightarrow 0$). Note that Eq. (89a) generalizes Eq. (67) by including the effect of dephasing in the energy eigenbasis.

It is clear from Eq. (89) that oscillations in the ground state probability $P_{G}(t_{f})$, which are present for finite $\Gamma$, imply a non-vanishing $\imaginary(\tilde{\rho}_{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 $\Gamma$, such coherence vanishes only in the adiabatic limit.

Appendix E Derivation of the adiabatic-frame TCL2/Redfield master equation

We start from the Hamiltonian given in Eq. (17), which we write as

[TABLE]

where $\kappa\equiv gt_{f}$. Our goal is to derive a master equation for the system evolution. It is convenient to do so using the time-convolutionless (TCL) approach Breuer and Petruccione (2002). To do so we must first perform yet another interaction picture transformation, defined by $H_{\mathrm{I}}(s)$, with the associated unitary $U_{\mathrm{I}}(s,s^{\prime})=T_{+}\exp[-i\int^{s}_{s^{\prime}}H_{\mathrm{I}}(s^{\prime\prime})ds^{\prime\prime}]$, where $T_{+}$ denotes forward time-ordering. In this frame the total Hamiltonian $H_{\mathrm{tot}}(s)$ becomes

[TABLE]

We can now calculate the TCL expansion generated by the superoperator

[TABLE]

whereupon

[TABLE]

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:

[TABLE]

where $\rho_{B}$ is the initial state of the bath, and the joint initial state is assumed to be in the factorized form $\rho_{S}\otimes\rho_{B}$. Note that the TCL2 approximation coincides with the Redfield master equation Breuer and Petruccione (2002).

Let

[TABLE]

denote the bath correlation function. By explicitly tracing out the bath, $\mathcal{K}_{2}(s)$ can be written as

[TABLE]

where

[TABLE]

After transforming back to the Schrödinger frame with respect to $H_{\mathrm{I}}(s)$ we obtain:

[TABLE]

where

[TABLE]

Appendix F Necessary convergence criterion

F.1 General criterion

Assuming $[\rho_{B},H_{B}]=0$ the correlation function becomes homogeneous in time, so we use the shorthand notation $C\pqty{x-y}=C\pqty{x-y,0}=C\pqty{x,y}$. We define the following quantities to bound the error of the expansion:

[TABLE]

and denote $\tau_{B}^{(1)}\equiv\tau_{B}$, which has a natural interpretation as the bath correlation time Albash et al. (2012).

Note that

[TABLE]

where $c_{2}=O(1)$ is a constant arising from the number of terms in the TCL2 double commutator expression (94). We can similarly estimate the magnitude of the TCL4 terms:

[TABLE]

where $c_{4},c_{4}^{\prime}=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. (102) as follows. Considering the first expression, we first make a change of variables as

[TABLE]

Because $1\geq s\geq s_{1}\geq s_{2}\geq s_{3}$, the new integration limits can be obtained as

[TABLE]

which is $s\geq x_{1}\geq x_{2}\geq x_{3}\geq 0$. The Jacobian $\absolutevalue{\det J}=1$ in this case and the new integral becomes

[TABLE]

Now we make another change of variables, with

[TABLE]

The new integration limits can be obtained by

[TABLE]

The first line means that $2s-v\geq u\geq v$, while the second line gives $0\leq v\leq s$ since $0\leq x_{1}\leq s$. Thus:

[TABLE]

where in the last inequality we used $s\leq 1$.

The same can be done for the second integral in Eq. (102):

[TABLE]

Combining these two bounds thus finally yields:

[TABLE]

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.:

[TABLE]

F.2 Ohmic bath case

Let us assume a spin-boson noise model, for which

[TABLE]

where $b_{k}$ is a bosonic annihilation operator for mode $k$ with frequency $\omega_{k}$, and $g_{k}=g\xi_{k}$ is the associated system-bath coupling strength, where $\xi_{k}$ is dimensionless and $g$ has units of energy. A standard approach is to introduce a spectral density such that $\absolutevalue{g_{k}}^{2}\mapsto J\pqty{\omega}d\omega$. For an Ohmic bath we have

[TABLE]

where $\eta$ 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 at $\tilde{H}_{SB}(s)$, the bath correlation function for the Ohmic spectral density is

[TABLE]

an integral which may be evaluated explicitly in terms of the Polygamma function Albash et al. (2012). In particular, for large $\beta\omega_{c}$ and $t_{f}/\beta$, the correlation function can be expanded as

[TABLE]

This form indicates a transition from a Markovian regime of purely exponential decay with a timescale of $\tau_{B}\overset{\omega_{c}\rightarrow\infty}{\rightarrow}{\beta}/({2\pi})$, followed by a non-Markovian regime of power-law decay with a timescale of $\tau_{M}=\sqrt{2\beta/\omega_{c}}$. The transition occurs at a time $\tau_{\mathrm{tr}}\approx\beta\ln(\beta\omega_{c})$ Albash et al. (2012). In the Markovian limit $\omega_{c}\rightarrow\infty$ we may thus replace Eq. (115) by

[TABLE]

and hence the correlation function integral of Eq. (100) becomes

[TABLE]

which replaces every factor of $\tau_{B}/t_{f}$ arising from the same integral in the bounds in the previous subsection. In particular, we now have the necessary condition $\norm{\mathcal{K}_{4}(s)[\tilde{\rho}_{S}\otimes\rho_{B}]}\leq\pqty{c_{4}+c_{4}^{\prime}}\alpha^{4}\left(\eta/({2\pi t_{f}\beta})\right)^{2}<\norm{\mathcal{K}_{2}(s)[\tilde{\rho}_{S}\otimes\rho_{B}]}\leq c_{2}\alpha^{2}\eta/({2\pi t_{f}\beta})$. Eq. (111) can thus be rewritten in the Markovian Ohmic case as

[TABLE]

For finite $\omega_{c}$, one can refine this bound by replacing Eq. (100) with

[TABLE]

where $s_{\mathrm{tr}}=\tau_{\mathrm{tr}}/t_{f}$. For our purposes the bound (118) suffices and is satisfied in all the numerical results presented in the main text. Namely, we have $\frac{g^{2}\eta t_{f}}{\beta}\leq 0.16$.

Appendix G Rotating wave approximation

Let

[TABLE]

be the one-sided Fourier transform of the bath correlation function, where

[TABLE]

and where $\gamma_{s}\pqty{\omega}/2$ and $S_{s}\pqty{\omega}$ are the real and imaginary parts of $\Gamma_{s}\pqty{\omega}$. Explicitly Breuer and Petruccione (2002):

[TABLE]

Here $\mathcal{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_{\mathrm{I}}(s)$ as

[TABLE]

where $\ket{\epsilon(s)}$ is an eigenstate of $H_{\mathrm{I}}(s)$ with instantaneous energy $\epsilon(s)$. We can then define the operator

[TABLE]

where

[TABLE]

is the dimensionless Bohr frequency, and the sum is over all pairs $\epsilon(s),\epsilon^{\prime}(s)$ subject to the constraint $\epsilon^{\prime}(s)-\epsilon(s)=\omega(s)$. The interaction picture master equation (93) can then be written to second order, with the TCL2 generator (94) as

[TABLE]

To obtain this master equation, we apply the standard Markovian approximation: change the integration variable $s^{\prime}\mapsto s-s^{\prime}$ and replace the upper limit with $\infty$. The RWA consists of neglecting terms in Eq. (126) for which $\omega^{\prime}\neq\omega$. A necessary condition for the validity of the RWA is Lidar (2019):

[TABLE]

which, unfortunately, is not always satisfied for the two-step Gaussian schedule (13) because [recall Eq. (125)]

[TABLE]

for $s$ outside the Gaussian pulse region.

Nevertheless, the RWA results in the interaction picture adiabatic Markovian master equation in Lindblad form Albash et al. (2012):

[TABLE]

where

[TABLE]

is the Lamb shift, and

[TABLE]

is the dissipator.

We can explicitly calculate $A\pqty{\omega(s)}$. First, recalling that $H_{\mathrm{I}}(\tau)=\frac{1}{2}\derivative{\theta}{\tau}U_{0}^{\dagger}(\tau)XU_{0}(\tau)$ [Eq. (7)], we realize that the eigenvalues and eigenvectors of $H_{\mathrm{I}}(s)$ can be written as

[TABLE]

Also, from the sequence of transformations leading to Eq. (17), the interaction terms have the form

[TABLE]

Substituting these expressions back into Eq. (124), we obtain

[TABLE]

After undoing the interaction picture transformation with respect to $H_{\mathrm{I}}(s)$ and ignoring the phase factors in the $A\pqty{\omega}$ operators, we obtain the Schrödinger picture master equation, namely Eqs. (20)-(22) given in the main text. In deriving this result we made use of the Kubo-Martin-Schwinger (KMS) condition Breuer and Petruccione (2002)

[TABLE]

where $\Delta$ is the dimensionless Bohr frequency in units of $1/t_{f}$:

[TABLE]

Appendix H Derivation of the semi-empirical Eq. (23)

The semi-empirical formula (23) can be derived directly from Eq. (129). Let us first write Eq. (129) in terms of the quantities defined in Eq. (21b):

[TABLE]

We now follow the steps in Ref. Albash and Lidar (2015) to obtain the solution in this interaction picture. Eq. (H) can be split into two decoupled ordinary differential equations:

[TABLE]

where

[TABLE]

and

[TABLE]

Additionally, the KMS condition allows us to write $\gamma_{d}(s)$ in terms of $\mathcal{F}_{+}(s)$

[TABLE]

The solution of Eqs. (138) is given by:

[TABLE]

where the initial conditions are:

[TABLE]

The next step is to move back to Schrödinger picture

[TABLE]

and write the open system ground state probability in terms of $\tilde{\rho}_{S}$:

[TABLE]

where

[TABLE]

For simplicity, we further denote $U^{a}(t)=U_{\mathrm{I}}(t)U_{0}^{\dagger}(t)$, whose elements can be related to those of $U_{\mathrm{I}}(t)$ in the $\{\ket{0},\ket{1}\}$ basis:

[TABLE]

where $k,l\in\{0,1\}$ and $\phi(t)=-E_{0}t/2$. Then:

[TABLE]

Because $U_{\mathrm{I}}(t)$ is the closed system unitary, we have

[TABLE]

and

[TABLE]

Eq. (145) becomes:

[TABLE]

This result is exact and corresponds to the numerical solution in the TCL2 case shown in Fig. 3.

We now make two additional approximations in order to arrive at a simpler expression. First, we ignore the Lamb shift term $\Omega(s)$ in Eqs. (142), which leads to:

[TABLE]

Second, we substitute the solution given in Eqs. (142) into line (151b):

[TABLE]

where in the last line we used the weak coupling assumption, $g^{2}t_{f}\ll 1$.

With these two approximations, Eq. (151) becomes the semi-empirical formula (23) with $P_{E}(0)=1/2$. We note that it is well known that for time-independent Lindbladians the RWA master equation has the Gibbs state as its steady state Breuer and Petruccione (2002). We do not recover this result for the time-dependent case. Rather, we find that the time-dependent Redfield master equation (TCL2) converges to the Gibbs state $P_{E}(\beta)=\frac{e^{\beta E_{0}/2}}{Z}$, but with a temperature that differs from that of the bath state, as illustrated in Fig. 3.