Assessing the quantumness of the annealing dynamics via Leggett Garg’s inequalities: a weak measurement approach

Adiabatic quantum computation (AQC) is a promising counterpart of universal quantum computation, based on the key concept of quantum annealing (QA). QA is claimed to be at the basis of commercial quantum computers and benefits from the fact that the detrimental role of decoherence and dephasing seems to have poor impact on the annealing towards the ground state. While many papers show interesting optimization results with a sizable number of qubits, a clear evidence of a full quantum coherent behavior during the whole annealing procedure is still lacking. In this paper we show that quantum non-demolition (weak) measurements of Leggett Garg inequalities can be used to efficiently assess the quantumness of the QA procedure. Numerical simulations based on a weak coupling Lindblad approach are compared with classical Langevin simulations to support our statements.

These properties define what has been called "classicity" or "macrorealism". Based on the assumptions above, Leggett and Garg derived Bell's-like inequalities that any system behaving classically should obey 2 . Violations of these inequalities provide evidence of quantum behavior of a system if one accepts that the alternative to classical probabilities is quantum mechanics. Therefore these violations can be interpreted as an indicator of the "quantumness" of a system. Following Ref. 3 , in this section we briefly introduce the Leggett-Garg's inequalities and discuss their properties as witness of "quantumness". Let us begin with the definition of a classical dichotomic variable Q which can assume value +1 or -1: Q(t i ) = Q i stands for the measurement value of the observable at time t i . We denote with P i (Q i ) the probability of obtaining the result Q i at time t i . The correlation function C i, j can be defined as follows: where the subscripts of P remind us of the times at which the measurements were performed. Assumption A, that is "Macrorealism per se", guarantees that P i j can be obtained as the marginal probability of P i jk (Q i , Q j , Q k ).
The assumption of "Non-invasive measurability" allows to drop the subscripts of P i jk and use the P(Q 3 , Q 2 , Q 1 ) alone (with ∑ Q 3 ,Q 2 ,Q 1 P(Q 3 , Q 2 , Q 1 ) = 1) to calculate the three correlation functions: C 1,2 ,C 2,3 ,C 1, 3 . We obtain  √ 2). a) Residual Energy in the absence of measurements; b) Residual Energy in the case of projective measurements at times t 2 = 0.3 t/t f and t 3 = 0.6 t/t f ; c) Fidelity in absence of measurements; d) Fidelity in the case of measurements at times t 2 = 0.3 t/t f and t 3 = 0.6 t/t f where ± stands for Q = ±. Next, we define The upper bound of K a 3 corresponds to P(+, −, +) = P(−, +, −) = 0, giving K a 3 = 1; the lower bound, instead, K a 3 ≥ −3 corresponds to P(+−, +) + P(−, +, −) = 1. Besides the inequality other inequalities exist, that can be found in the literature. Various symmetry properties can be used to derive other constrains on the correlations. The change Q → −Q in K a 3 gives rise to the following inequality: Finally, the last, different, third order inequality can be obtained from K a 3 , just changing a sign: These are the only three different inequalities that can be derived from correlations to third order. Higher order inequalities can also be constructed.

Measurement scheme
In this Section we sketch an idealized measurement approach which can be extended from projective to weak measurement, to reduce the invasiveness of the classical measurement process. Resorting to a weak measurement scheme is unavoidable to allow for a successful annealing. Indeed in Fig. 1S we show the residual energy and the ground state population during the annealing dynamics in the absence and in the presence of two projective measurements in order to demonstrate the necessity of weakening the measurement approach. In the panels a) and c) we show the residual energy and the ground state population in the absence of measurements. At the annealing time the latter is approximately 1 and the former is nearly 0 which reveals that the quantum annealing has been successful. On the other hand, in the panels b) and d) we calculate the residual energy and the ground state population, while measuring one of the possible C 2,3 necessary to build the K 3 s. Clearly, at the measurement times t 2 and t 3 , the (projective) measurement procedure suddenly alters the population of the ground state. The ground state population is very poor at the annealing time, and the residual energy considerably larger that zero, signalling that the annealing procedure has failed. Thus approaching to the calculation of the Leggett-Garg's functions with weak measurements is necessary to evaluate the coherence of the system during the quantum annealing with negligible effects on the dynamics. For the sake of concreteness, we consider a superconducting flux qubit, as the system S, and an hysteretic DC SQUID, as the ancilla/detector (Fig.2S ). In a superconductiong flux qubit, when polarized by an external flux φ close to φ 0 /2, where φ 0 is the flux quantum, the current can flow clockwise or anti-clockwise. A spin degree of freedom can be associated to the circulating current, e.g. the state of the system can be denoted by |↑ , if the current is clockwise, while it is |↓ , if the current is anti-clockwise. On this basis of eigenstates of σ z , the Hamiltonian of the flux qubit can be written as in the main text: where Γ x is the tunnel coupling between the current states and Γ z = 2I p ( φ 0 2 − φ ), where I p is the magnitude of the current flowing in the flux qubit. By introducing a linear time dependence in the same Hamiltonian one obtains: which is the Hamiltonian that describes a linear annealing protocol. Here s = t/t f goes from 0 to 1 and t f is the total annealing time. The system evolves along its dynamics and, to evaluate the possible violation of the Leggett-Garg's inequalities, it is necessary to measure correlations of system variables at different times. The quantum state of the system can be read out by exploiting the mutual inductance M of the qubit with the ancilla DC SQUID, represented by the interaction Hamiltonian Here J is the current circulating in the DC SQUID. As described in ref. 4 , the measurement can be performed considering that the ancilla, appropriately polarized with J close to the critical current I c , by means of a flux Φ ∼ nφ 0 (n integer), can be either in a superconducting state with zero voltage V = 0 or in a dissipative state with a finite voltageV . Let the circulating current in the ancilla be J < I c , prior to measurement. With a short current pulse, the DC SQUID can be biased very close to its critical current I c . During the pulse the ancilla has a certain probability of staying in the V = 0 state, or switching to the dissipative state depending on the state of the qubit |↑ , |↓ . Indeed, the circulating current I p induces an electromotive force in the ancilla loop, which increases or decreases I b , its bias current, respectively. After the pulse, I b is maintained stable at a value I b = I c /2 while the ancilla relaxes in one of its two possible states in order to measure the voltage output. If the relaxation time of the ancilla T r is much larger than the so-called discrimination time T V one can obtain meaningful information from a measurement and infer the qubit state 5 . If T r is not long enough compared to T V , such that a single measurement cannot provide the full information to evaluate the ancilla voltage and, consequently, the current state of flux qubit, the measurement becomes minimally invasive and weakly perturbs the quantum coherence of the evolution. An estimate of T V , can be given by requiring that the signal-to-noise ratio approaches unity . This occurs when the spectral density S V ( f ) of the output noise of the detector at frequency f can be approximated as To be more specific, let us map the values V = 0 and V =V onto a dimensionless variable x which assumes values ±1: V =V (1 + x)/2. The probability P(x), of reading a value x after the measurement is

3/9
where ρ Q is the reduced density matrix of the flux qubit given in Eq.(30) and P ± is the probability of having a value of x = ±, as the result of the measurement. P ± (x) can be viewed as two gaussian distributions centered around x = ±1 with the variance D = T r T V in analogy with Ref. 6 where a conceptually similar approach is investigated. The change of the time of measurement in the experiment amounts to tuning the width of the P ± (x) peaks. The Ansatz of a Gaussian distribution is due to the fact that a long measurement process gives V = 0 or V =V with probability ρ Q↓↓ or ρ Q↑↑ , while, by taking a short interaction time between qubit and ancilla, strange voltage values are not excluded 7 . Repeating the experiment many times, a bimodal distribution is expected with two unequal peaks centered at V = 0 or V =V , respectively. Of course, no matter how weak the measurement is, the density matrix ρ Q of the system turns out to be slightly modified, depending on the outcome of the measurement of the variable x. Following Ref. 8 , the transformation from ρ Q (before the measurement) to ρ Q (after the measurement) is defined as (time label omitted): This expression can be written in a more convenient form. From Eq.(12) we get Let us denote with γ the x/D ratio and get then ρ Q↓↓ = ρ Q↓↓ e γ ρ Q↓↓ e γ + ρ Q↑↑ e −γ .
(15) Therefore, one obtains the following quantum-map from ρ to ρ : The value of x is stored for evaluating the correlation functions C i j . By tuning D, we are able to weaken the measurement till the post-measurement update in the density matrix is negligible. This is crucial if the goal is of looking at the Leggett-Garg's correlations during an annealing dynamics, without spoiling the quantum coherence of the qubit too much. To sum up, the annealing protocol that we have realized in the simulation is the following. Firstly, one prepares the system in the ground state of the Hamiltonian H(0). The system evolves under U = e −iH(t/t f )t . Computationally this means solving the differential equation for the density matrix Eq.(30) with a fourth-order Runge-Kutta algorithm. At fixed times t 1 ,t 2 (or t 2 ,t 3 or t 1 ,t 3 ) one performs two weak measurements, by extracting values of x, which are used to evaluate C i j and to update the density matrix at the corresponding times by means of the probability distribution of Eq. (11). To gain sufficient statistics, the same evolution is repeated up to 10 6 times and the Leggett-Garg's correlation functions are evaluated as an average on the different runs. In this way, the Leggett-Garg's inequalities can be tested with weak measurements with minimal perturbation of the system during its dynamics. The idealized measurement approach described here hides a number of experimental challenges. For a study on the back action of the detector on the flux qubit, on the problems related to the Joule heating in the dissipative state and on the fidelity of the weakness of the measurement we refer to Ref. 5 and references therein.

Lindblad approach to the quantum dissipative environment
Here T is the time ordering operator in real time. The evolution of the system and the bath in the absence of interaction is governed by The interaction between system and bath is described, in full generality, by where g α are coupling constants. Adiabatic switching of the interaction at negative times is assumed. We define the total Hamiltonian for the joint system-bath universe, and the time-dependent density operator ρ(t), whose dynamics is expressed by the Von Neumann equation and the full system-bath evolution operator, Moving to the interaction picture, we definẽ where and A α (t), B α (t) are the time-evolved operators, U(t, 0) andρ(t) satisfy the following differential equations: Given the two point correlation function B αβ (τ) ≡ B α (τ)B β (0) , the spectral-density matrix of the bath is: where itsreal and imaginary part are

5/9
Lindblad theory eventually leads to master equation for the reduced density matrix (representing only the system variables), where the adiabatic dissipator D is and the Lamb shift Hamiltonian takes the form They are written in terms of the Lindblad operators L αω (t), which are defined as where { ε a (t) } are the instantaneous eigenvectors of the system Hamiltonian.
The universe we study here is a single spin coupled to a bath of bosonic harmonic oscillators described by the Hamiltonian where b † k and b k are, respectively raising and lowering operators for the k-th oscillator with frequency ω k and the frequency spectrum is assumed to be continuous as usually in the spin-boson model [9][10][11] . The bath is assumed to be in thermal equilibrium at inverse temperature β = 1/k B T , so that Its density operator is just ρ B = e −β H B /Z .

The interaction between the system and the bath is
. The Fourier transform of the bath correlation function is: where Θ(±ω) are Heaviside functions 9 . The model is fully determined once the explicit form of the function J(ω) is given. In this paper, we consider an Ohmic bath 10 , characterized by where ω c is a high-frequency cut-off that is the maximum phonon energy and η is a dimensional parameter with dimensions of time squared. ω c has been chosen 25 Γ x in the simulations. In conclusion, we define α = ηg 2 and explicit γ as follows:

Simulation of a two-level-system dissipative dynamics with classical Langevin dynamics
In this Section we consider the flux x of a flux qubit as a classical variable and we study its classical thermal evolution in a time dependent double well potential U(x,t) that mimicks the superconducting flux qubit evolution. Such double well potentials have been used in the analysis of experimental realization of flux qubits 12 . By setting up a classical model evolution which corresponds to the single qubit evolution studied in the main text, we are able to obtain the LG functions for the classical system and highlight the close correspondence between the classical and the quantum outcome, when the quantum evolution provides LG functions very close to the classical limit. Let us start defining the potential U(x,t) at t = 0. The time independent Hamiltonian of the flux variable is given by: where the mass is typically related to the device capacitance and the potential has two minima. The lowest energy wavefunctions of the potential, have to correspond to the two states |↑ and |↓ of the qubit of Eq. (8). Parameters Γ x and Γ z are given respectively by The case Γ z = 0 corresponds to a symmetrical double well, while Γ x = 0 corresponds to a very high potential barrier, so that the states do not overlap and the matrix element between them vanishes.
To construct the potential at t = 0, we first approximate the potential U(x) near the two minima with two harmonic wells. Next, we consider the lowest energy wavefunctions centered in each of the two harmonic wells and the corresponding states |L and |R . The two wave functions always overlap to some extent, as long as the potential barrier between the wells in U(x) is finite and their wavefunctions have to be orthogonalized.
The states |↑ and |↓ are superpositions of |L and |R with wavefunctions in which the variable x is mostly confined in each of the two wells of the potential.