Time-Dependent Hamiltonian Reconstruction using Continuous Weak Measurements

Reconstructing the Hamiltonian of a quantum system is an essential task for characterizing and certifying quantum processors and simulators. Existing techniques either rely on projective measurements of the system before and after coherent time evolution and do not explicitly reconstruct the full time-dependent Hamiltonian or interrupt evolution for tomography. Here, we experimentally demonstrate that an a priori unknown, time-dependent Hamiltonian can be reconstructed from continuous weak measurements concurrent with coherent time evolution in a system of two superconducting transmons coupled by a flux-tunable coupler. In contrast to previous work, our technique does not require interruptions, which would distort the recovered Hamiltonian. We introduce an algorithm which recovers the Hamiltonian and density matrix from an incomplete set of continuous measurements and demonstrate that it reliably extracts amplitudes of a variety of single qubit and entangling two qubit Hamiltonians. We further demonstrate how this technique reveals deviations from a theoretical control Hamiltonian which would otherwise be missed by conventional techniques. Our work opens up novel applications for continuous weak measurements, such as studying non-idealities in gates, certifying analog quantum simulators, and performing quantum metrology.


I. INTRODUCTION
A central challenge in building reliable quantum computers and simulators is the characterization and validation of the Hamiltonians which they generate.This problem is of interest in both the quantum gate characterization and quantum simulation communities.The former seeks to quantify the errors in the quantum gates used in various experimental platforms and understand what limits their performance.Techniques such as quantum process tomography [1][2][3], gate set tomography [4][5][6][7], and randomized benchmarking [8][9][10][11] have been developed and employed experimentally to characterize the type and level of errors present in sets of quantum gates in superconducting qubit, trapped ion, spin qubit, and quantum dot platforms [12][13][14].However, such techniques do not necessarily elucidate what causes a gate to deviate from its idealization at the level of the Hamiltonian.For example, in superconducting qubit platforms, gates are performed on qubits by applying RF pulses generated by room temperature electronics.However, distortions in the pulse shape [15][16][17][18][19] due to cryogenic microwave components, cross-talk between qubits, and deviations from approximations in the Hamiltonian cause the actual Hamiltonian as seen by the qubit to differ from the theoretical Hamiltonian generated by those pulses.We can then ask the following question: if we apply a pulse of duration t intended to generate a quantum gate U = T exp(−(i/ ) t 0 H(t )dt ) (where T is the time-ordering operator), how does the actual time-dependent Hamiltonian H(t) produced deviate from the ideal or theoretical one?Current gate characterization techniques do not answer this question, as they do not directly compute the Hamiltonian H(t).A complete description of H(t), though, would shed light on the sources of nonideality of a gate, such as driving of undesirable terms in the Hamiltonian, or could be used to improve optimal control schemes [20][21][22][23].
This question also arises naturally in the setting of analog quantum simulators, where the goal is to validate the Hamiltonians used in simulating the evolution of many-body quantum states [24][25][26][27][28][29][30].Here a system of qubits evolves under an interacting many-body Hamiltonian, such as in the Fermi-Hubbard model.For general many-body systems, the number of terms in the Hamiltonian can scale exponentially in the number of degrees of freedom.However, for systems with local interactions, this scaling is only polynomial.This observation has been used to propose several techniques for reconstructing Hamiltonians, for example by constructing the so-called "correlation matrix" [31][32][33].Still, such techniques rely on interrupting the evolution to perform projective measurements.As the system controls, qubit evolution, and measurement apparatus cannot change instantaneously, the reconstructed Hamiltonian differs from the uninterrupted Hamiltonian seen by the qubit.
In this work, we introduce an experimental technique and theoretical tool for reconstructing an a pri-ori unknown time-dependent Hamiltonian using continuous weak measurements [34][35][36][37][38][39][40][41].Instead of interrupting the Hamiltonian evolution for projective measurements, we concurrently measure the system and evolve it under the Hamiltonian, and our reconstruction technique accounts for the effect of the continuous measurement tone, thereby avoiding introducing spurious dynamics into the reconstructed Hamiltonian.We experimentally apply and validate this technique for a system of two superconducting qubits coupled by a parametric coupler, which can be modulated to produce a variety of two-qubit entangling Hamiltonians [42][43][44][45][46][47][48][49].By continuously measuring the field amplitude while applying pulses to the qubits and coupler, we can recover the qubits' z coordinates σ z (t) during the evolution.Next, to recover the Hamiltonian generated by those pulses, we introduce a reconstruction algorithm which takes as input the qubits' z coordinates as a function of time and outputs an estimate of the Hamiltonian expressed in the Pauli basis.
We validate our technique by applying it to various single-qubit and two-qubit Hamiltonians.We compare the reconstructed Hamiltonians to estimates made from independent calibrations and find good agreement.We also achieve an average of 97.5% (97.1%) fidelity between the final state as predicted by the reconstructed single (two) qubit Hamiltonians and tomography of the system at the end of the time evolution.
Compared to gate characterization techniques, which use only initial and final states as input, our method provides a complete time-dependent description of the dynamics produced by a pulse.Compared to previous timedependent Hamiltonian reconstruction techniques [31][32][33]50], our method does not rely on interrupting the Hamiltonian nor on having carefully calibrated Hamiltonian parameters [36,39,[51][52][53][54][55][56][57][58][59][60][61][62][63].The theoretical and experimental techniques introduced in this work therefore present an alternative approach for analyzing the dynamics of quantum systems, and our results suggest that this technique may prove useful in analyzing quantum gates, especially in emerging quantum platforms, certifying analog quantum simulators, and performing quantum metrology.
The paper is organized as follows.In Sec.II, we explain the experimental setup and how the continuous weak measurements are performed.We then formally state the problem, define the algorithm for reconstructing Hamiltonians for one and two qubits, and state the requirements.In Sec.III, we apply the Hamiltonian reconstruction algorithm to Hamiltonians generated by various single qubit pulses.In Sec.IV, we demonstrate the operation of the flux-tunable coupler to generate entangling two-qubit Hamiltonians.In Sec.V, we then experimentally validate the algorithm on two-qubit Hamiltonians and, as a novel application of this technique, estimate a dynamical coherent fidelity.Finally, in Sec.VI we comment on the limitations of the technique and how they could be overcome and propose how this technique could be applied in other settings.

II. EXPERIMENTAL PROCEDURE
In superconducting qubit experiments, Hamiltonians are applied to the device using RF pulses, which are generated by room temperature electronics and sent to the device via cryogenic wiring.However, the Hamiltonian as seen by the device may differ from an ideal or theoretical one if, for example, the pulse shape is distorted by the wiring.Cross-talk between qubits or drive lines and pulse miscalibration also drive undesirable terms in the Hamiltonian.Finally, the theoretical Hamiltonian describing the system generally only approximates the dynamics, and high power pulses cause violations of the assumptions of those approximations.Our goal, then, is to reconstruct the true, time-dependent system Hamiltonian, which we refer to as the "target" Hamiltonian H t (t).For example, a general single qubit Hamiltonian H t (t) can be expressed as where σ = (σ x , σ y , σ z ), a vector of the three Pauli operators, and we seek to reconstruct the drive amplitudes Ω(t).We now introduce our procedure for using continuous weak measurements to reconstruct target Hamiltonians for systems of one and two qubits.We first detail how measurements on our device are performed and then explain Hamiltonian reconstruction algorithm.

A. Experimental Setup
Our experiment consists of two superconducting transmon qubits [64], labeled Q1(Q2) in Fig. 1(a), each dispersively coupled to a coplanar waveguide resonator, labeled R1(R2).The two qubits are coupled through a high frequency flux-tunable coupler consisting of an asymmetric SQUID and a large U-shaped paddle to achieve high capacitive coupling to the transmons [42][43][44][45][46][47][48][49].The frequency of the coupler and qubits and the resulting static coupling between the two qubits can be tuned by applying a DC bias on the flux line, which changes the flux Φ ext through the SQUID.Throughout this work, we operate with Φ ext = 0.The frequencies of Q1(Q2) and R1(R2) at this point are ω q /2π = 5.319 GHz(5.271GHz) and ω res /2π = 6.631GHz(6.529GHz), where ω q is the |0 → |1 transition frequency of the transmon and ω res is the resonator frequency.In Sec.IV, we demonstrate parametric modulation of the coupler, which generates a two-qubit entangling Hamiltonian, and then in Sec.V, we apply our reconstruction technique to both qubits evolv- ing under such target Hamiltonians.
In the interaction picture, the static dispersive Hamiltonian H disp for the system is approximately where χ i is the full dispersive shift of the i-th resonator frequency conditioned on the state of the i-th qubit, a i (a † i ) is the annihilation (creation) operator for the i-th resonator, and σ is the Pauli z operator on the i-th qubit.In our device, χ 1 /2π(χ 2 /2π) is 0.64 MHz(0.75MHz).The resonators are each coupled through a Purcell filter [65][66][67] to a transmission line at a rate κ/2π = 11.78MHz (14.47 MHz).A complete diagram of the experimental setup is provided in App.A, and a full table of device parameters is provided in App.B. We continuously measure the qubits as they evolve under H t (t) by simultaneously applying a weak measurement (WM) tone on the transmission line at the midpoint frequency between the two qubit statedependent frequencies of the resonator (top of Fig. 1(a)).The tone populates the resonator with mean photon number n = a † a ≈ 0.94(0.82),and the reflected resonator field is then measured as follows [63].We amplify both quadratures of the reflected field using a traveling wave parametric amplifier (TWPA) [68], and after further amplification and subsequent demodulation, digitize both quadratures at 1 gigasample per second (GSa/s) to generate voltage records.
Each shot of the experiment produces a noisy voltage record (Fig. 1(b)) carrying a small amount of information about the qubit z coordinate.To estimate the ensemble average z(t) = Tr(ρ(t)σ z ) from the voltage records, we first average many records and filter and downsample them to 2 ns resolution (see App. E for more details).When χ/κ 1, the average resonator field is proportional to z(t), retarded by a delay τ ≈ 2/κ: where Ω wm is the amplitude of the measurement tone.
To extract z(t) for Q1(Q2), we shift the averaged, filtered trace by τ = 27 ns(22 ns).When the qubit state z(t) does not evolve too quickly compared to the resonator linewidth, the resonator field adiabatically follows the qubit [69].For example, if the qubit is driven around the x axis with amplitude Ω X , we require that 2Ω X κ (App.C).We therefore strongly couple the resonators to the readout port.
To reconstruct H t (t), we perform this procedure for various initial states of the qubits, resulting in a collection of ensemble averages shown in Fig. 1(c).To rescale the voltage traces, we additionally prepare the qubits in |0 and |1 and perform the continuous measurement described above but without any target Hamiltonian applied to the qubits.These correspond to the lower z = 1 and upper z = −1 limits for the resonator field, which we use to rescale the voltage records in Fig. 1(c).Finally, the extracted z(t) for each qubit is inputted to the Hamiltonian reconstruction algorithm.
Besides allowing us to measure z(t), the WM tone induces two additional effects on the qubits.First, it shifts the qubit frequency by χn ∝ Ω 2 wm according to Eq. ( 2).Resonant single qubit pulses must therefore be applied at this shifted frequency.Second, the qubit dephases more rapidly, with the measurement-induced dephasing rate Γ d ∝ n.Using Ramsey measurements with the simultaneous WM tone, we measure Γ d [37] and supply it as an input to the algorithm, which takes this rate into account so that the extra dephasing does not impact the estimate of H t (t).This is explained in more detail in Sec.II B.
Because the algorithm accounts for the additional dephasing, using higher measurement power appears optimal, as SNR ∝ n, and increasing SNR reduces the number of samples and acquisition time needed to estimate z(t) accurately.However, because Γ d ∝ n, distinct initial states become indistinguishable more rapidly at higher measurement power.We will see in Sec.II B how this poses a problem for the algorithm, but by considering the limit of projective measurement, this obstacle can be understood physically as a consequence of the quantum Zeno effect.In this limit, because the qubit state is pinned to the poles of the Bloch sphere, it no longer undergoes coherent dynamics.Therefore, we cannot recover those dynamics from the measurement record.Considering this tradeoff between acquisition time and qubit collapse, we calibrate our measurement strength to ensure 1/Γ d t p for all pulse durations t p that we study.

B. Reconstruction Algorithm
We now explain the algorithm (Fig. 1(c)), which takes as input the ensemble average z(t) for various initial states, along with Γ d and full tomography of the initial states {ρ(0)}, and estimates the time-dependent target Hamiltonian H t (t) (Fig. 1(d)).
We first concentrate on the setup for a single qubit system, state the problem, and set up the notation.Formally, we consider an a priori unknown time-dependent Hamiltonian H t (t), for a duration t p .
The qubit, described by the ensemble-averaged density matrix ρ, evolves according to where we have ignored Lindblad operators other than measurement-induced dephasing for now.We express H t (t) in a complete operator basis, such as the Pauli operator basis in Eq. ( 1), with time-dependent amplitudes {Ω i (t)} to be determined.We measure z(t) using continuous weak measurement simultaneous with H t (t) with time resolution ∆t and discretize time with index n so that t n+1 = t n + ∆t.
In the Heisenberg picture, the discretized equations of motion for a qubit at first order in ∆t are: where ijk is the Levi-Civita symbol, repeated indices are summed, and the indices take values in {1, 2, 3} which correspond to the three directions {x, y, z}.To recover the drive amplitudes Ω j (t n ), it is convenient to express the Eq. ( 6) more concisely in terms of matrices M , D, and Γ as We now show that if the entire vector of operators σ(t n ) is known for all t and at least S ≥ 2 linearly independent σ(t n ) are used, then there exists a unique solution for Ω(t n ) in Eq. (7).Let σ [s] (t n ) denote the s-th vector out of S at time t n , and let so that Σ(t n ) has 3S elements, D S is a block-diagonal matrix of dimensions 3S×3S, and M S (t n ) has dimensions 3S × 3. Then Eq. ( 7) becomes This can be inverted to solve for Ω(t n ) as long as M S (t n ) has a left inverse (M S (t n )) + , which is true if and only if it has rank 3 [70].This condition can be satisfied for S ≥ 2.
For example, M 2 (t n ) has rank 3 if it is constructed from the two initial states |0 and |+ = (|0 + |1 )/ √ 2, but this condition cannot be satisfied for S = 1.For S ≥ 2, then: Viewed in the Schrodinger picture, this yields a procedure by which the drive amplitudes Ω(t n ) are recovered from a time series of Bloch vectors σ(t n ) starting from S initial states.Moreover, the dephasing due to the weak measurement tone is taken into account, and the solution in Eq. ( 13) may be extended to account for other dissipative effects.If all three components of the Bloch vector are measured continuously, Eq. 13 can be used to reconstruct the Hamiltonian.However, in our experiment, we measure only σ z (t n ) from the resonator field as in Eq. 3, so we have only one component of the Bloch vector.In that case, the preceding analysis simplifies as follows.Let (14) be reduced versions of the Pauli vector, drive amplitude vector, and matrix M .Similarly to (11), let The matrix D S is not necessary in this case because the measurement-induced dephasing does not appear in the evolution of the measured coordinate σ z .Then the firstorder evolution of all S states can be summarized as Again, if S ≥ 2, then Ω (t n ) can be solved: Note that M ( σ (t n )) depends on σ x (t n ) and σ y (t n ), which are not measured.However, full tomography of ρ(0), the initial state, gives the entire Bloch vector σ(0), and σ(t n ) for t > 0 is computed by first solving for the amplitudes Ω (t n−1 ) using Eq. ( 17) and then integrating Eq. ( 4) for one timestep.Alternating between solving Eq. ( 17) to obtain Ω (t n ) and integrating Eq. (4) to obtain Σ S (t n ), we recover Ω X (t) and Ω Y (t).Although D S does not appear in Eq. ( 17), the measurement-induced dephasing still affects the reconstruction during the integration step.Dephasing also shrinks the Bloch vector and will tend to make M S become singular, so it is essential to keep 1/Γ d t p .The inability to recover Ω Z (t) is a limitation of applying a first-order update: the qubit z coordinate is insensitive at first order in time to drives around the z axis.This limitation can be overcome by going to a second order update scheme, which we detail further in App.D. In this work, however, we restrict our attention to the first-order update methods due to experimental limitations that we discuss later, taking care to ensure that Ω Z (t) ≈ 0. Alternatively, if Ω Z (t) is nonzero but known a priori, it is possible to account for it by including it in H t (t) when integrating Eq. ( 4).We refer to this process as "preconditioning" Ω Z (t).
The discussion above can be generalized to Hamiltonians acting on Q qubits with total Hilbert space dimension d = 2 Q , and the generalization of Eq. ( 17) then requires S ≥ d(d − 1)/Q initial states.The amplitudes associated to Pauli terms which are tensor products of the identity and σ z cannot be recovered, of which there will be d − 1 out of d 2 − 1 total.Just as for Ω Z (t) in the single qubit case, if these amplitudes are known, they can still be preconditioned.
We summarize the requirements to perform Hamiltonian reconstruction and refer the reader to App.D for implementation details.The algorithm takes as input the time series z coordinate data and the initial Bloch vector for various initial states, calibrated values for dissipative rates, and any preconditioned drive amplitudes.It is not necessary that the initial states be prepared perfectly along any particular axis, but that M S is full rank.The algorithm returns an estimate of the drive amplitudes Ω(t), except those for which it cannot solve using first-order updates, and an estimate of Bloch vector at all times for each initial state.

III. SINGLE QUBIT HAMILTONIAN RECONSTRUCTION
We first apply and validate the algorithm on a single qubit (Q1) by applying three different target Hamiltonians.The pulse sequence for these experiments is shown schematically in Fig. 2(a).The initial preparation pulse prepares the qubit in one of six initial states {|±X , |±Y , |±Z }, where, e.g.|+X corresponds to the +1 eigenstate of the Pauli operator σ x .The same target Hamiltonian is then applied to each initial state for t p = 250 ns and takes the general form The sequence concludes with full state tomography, which we use not as an input for the algorithm, but as validation of the algorithm's output.Full tomography of each initial state is also performed, as the complete density matrix at t = 0 is an input to the algorithm.
For each target Hamiltonian, we vary the initial state and average, filter, and deconvolve ∼ 5.4M trajectories with 2 ns time resolution to obtain a smooth estimate of FIG. 2. Single qubit Hamiltonian reconstruction.(a) Pulse sequence.We weakly measure Q1 while applying a qubit control pulse (duration tp).After the pulse we use state tomography as independent validation of the reconstructed Hamiltonian Ht(t).A 40 ns buffer is added between consecutive pulses to ensure the resonator has reached a steady state.(b) Reconstructed Hamiltonian (∆t = 2 ns) for three different pulses: (i) Calibrated π-pulse with flat top and cosine ramp.(ii) Same as in (i) plus additional out of phase waveform (iii) Same as in (i) plus a sinusoidal amplitude "error" pulse.(iv) Extracted π-pulse error estimated from the reconstructed Hamiltonian subtracted from calibrated estimate.(c) Post pulse reconstruction infidelity quantifying overlap between predicted final state and state tomography for the three pulses in (b).(d) Decoupling of resonator from qubit.For fast pulses, the qubit coordinate z(t) obtained from the resonator evolution (solid lines, transparency indicates evolution after the end of the π pulse) deviates from the QuTiP simulation (dashed lines).(e) Reconstruction infidelity averaged over all initial states for different gate durations.Error bars represent spread of 1 − Fr over different initial states used.The qubit simulation (dashed line) which includes the readout resonator captures the measured trend and shows that experimental value saturates the minimum for all pulse lengths.Low reconstruction fidelity for fast pulses is explained by the deviation observed in (d).
z(t) (see App. E), which are then inputted to the Hamiltonian reconstruction algorithm.
We validate the output of the algorithm in three ways.First, we compare the reconstructed pulse shape to the applied waveform and confirm that the drive axis agrees with the phase(s) of the applied pulse(s).Next, we compare the reconstructed amplitude to an estimate obtained by an amplitude calibration, i.e. by fitting the rate of Rabi oscillations vs. drive amplitude.Lastly, the algorithm reconstructs the full density matrix at each timestep, so we compare the final state predicted by the algorithm to the experimental tomography of the final state.We call the agreement between these two states the "reconstruction fidelity" F r , defined as follows: where ρ tomo is the density matrix as measured by full tomography of the final state and ρ rec is the final density matrix as predicted by the algorithm, computed using the tomography of the initial states {ρ(0)}, the reconstructed Hamiltonian amplitudes Ω X (t) and Ω Y (t), and the measurement-induced dephasing Γ d .We present the results for three target Hamiltonians in Fig. 2(bc).The first target Hamiltonian (Fig. 2(b)(i)) is generated by a t p = 250 ns π pulse around the x axis with waveform consisting of a flat-top and cosine ramp edges.The reconstructed pulse shape Ω X (t) qualitatively matches with that of the applied waveform and quantitatively agrees with the amplitude of the estimate.The algorithm also concludes that there is zero drive along the y axis, as expected.Finally, tomographic validation of the reconstructed final state in Fig. 2(c) shows that the algorithm achieves an average reconstruction fidelity of 97.2%.
Fig. 2(c) also shows that the reconstructed Hamiltonian produces different fidelities for different initial states.For example, the |±X states are insensitive to the shape of Ω X (t) in pulse (i) and achieves the highest fidelity of the three preparation axes.To assess the performance of the algorithm when the target Hamiltonian acts non-trivially on all six initial states, we next apply it to a target Hamiltonian generated by two simultaneous outof-phase π pulses with the same waveform as in (i).The reconstructed Hamiltonian (shown in Fig. 2(b)(ii)) again shows good agreement with the estimates of Ω X (t) and Ω Y (t) and achieves 98.4% reconstruction fidelity.This indicates that the high fidelity achieved in (i) is not an artifact of the choice of initial states.
Finally, we demonstrate that the reconstructed Hamiltonian gives information beyond what is provided by characterization techniques which use as input only the initial and final states.We apply our technique to a target Hamiltonian generated by a sinusoidal "error" pulse with period t p = 250 ns superimposed on the π pulse from (i) (Fig. 2(b)(iii)).We motivate this experiment by viewing the superimposed error as a deviation from a theoretical control Hamiltonian H c (t) [20,21,23], which we take to be a flat-top pulse with cosine ramp edges.Importantly, this deviation is undetectable in the final state tomography, as its effect integrates to zero.For visual clarity of the error pulse we use the estimate for Ω X (t) from (i) as H c (t) and plot it as the dashed curves in Fig. 2(b)(iii).The difference between the reconstructed amplitude and this estimate, plotted in Fig. 2(b)(iv), shows that the algorithm correctly reveals the sinusoidal deviation from H c (t).Thus although H t (t) and H c (t) produce the same unitary at t = t p , our technique reveals the difference between them.
Our reconstruction technique can, in principle, be applied to target Hamiltonians generated by faster pulses.However, when the rate of qubit evolution exceeds the resonator linewidth (Ω κ/2), the system no longer remains in the adiabatic regime, and the qubit z coordinate can no longer be extracted easily from the reflected resonator field [69].In Fig. 2(d), we show that the reconstructed qubit evolution for π pulses of varying length deviates from z(t) obtained in QuTiP simulations of the qubit and resonator [71,72] when the adiabatic condition is violated.This results from the well-studied decoupling of the resonator from a rapidly driven qubit [69,73].Poor estimation of the qubit z coordinate, in turn, leads to an overall lower reconstruction fidelity, shown averaged over the initial states in Fig. 2(e).The dashed black line indicates the error in the reconstruction algorithm when the simulated resonator field is used to estimate z(t) and then inputted to the algorithm.It represents the minimum reconstruction error one can expect from our protocol for each pulse duration.Fig. 2(e) shows that for fast pulses, the primary source of error in reconstruction is violation of the adiabatic condition.The remaining difference between the experimental and simulated values for F r may be explained by additional transience and dephasing of the qubit state during the tomography and projective measurement pulses.Our technique therefore performs reliably when the amplitudes Ω κ/2 = 5.89 MHz.
Having validated the algorithm for a single qubit and established the regime in which our technique performs well, we next validate it on a system of two qubits.We apply the algorithm to two-qubit entangling target Hamiltonians generated by parametrically modulating the coupler, which we discuss next.

IV. XY INTERACTION
The first-order update rule described in Sec.II B does not allow for the reconstruction of Ω IZ (t), Ω ZI (t), and Ω ZZ (t).The presence of such terms in the system can limit the performance of the reconstruction algorithm, so it is important to minimize and to characterize them.The first two unwanted terms correspond to off-resonant driving and can be controlled by tuning the drive frequency.The third term can arise due to ZZ interactions between qubits.Our device uses a tunable coupler architecture (Fig. 3(a)) in which the static ZZ coupling between the qubits can be tuned by biasing the coupler SQUID flux Φ ext [42][43][44][45][46][47][48][49].Further, by modulating Φ ext using an additional RF pulse on the flux line FF, it is possible to activate a variety of entangling two-qubit gates such as iSWAP and CZ.Here, we review the physics of the coupler and demonstrate operation and measurements of this system.
We first examine the static properties.The two qubits Q1 and Q2 are coupled directly through capacitance between transmon pads as well as indirectly through the tunable coupler C12.The full two-level Hamiltonian  30)) is first-order insensitive to the flux (inset).The static ZZ rate ζ is measured using a JAZZ measurement [74], and the data are fit using models described in App.H full of this system is + σ (2) where σ ± are the raising/lowering operators for the qubits, σ (c) z is the Pauli z operator for the coupler mode, σ (c) ± are the raising/lowering operators for the coupler mode, ω c is the coupler frequency, g 12 is the direct coupling between the qubits, and g ic is the coupling between qubit i and the coupler.The coupler frequency ω c tunes with the flux as where Φext = Φ ext /Φ 0 , Φ 0 = h/(2e), d is the SQUID asymmetry, and ω c,0 is the coupler frequency at zero flux.Denote the qubit-coupler detuning ∆ i (Φ ext ) ≡ ω . By applying the Schrieffer-Wolff transformation, the effective, static two-qubit Hamiltonian H eff to second order in |g ic /∆ i (Φ ext )| is approximately [44]: Hqq / = (J(Φ ext ) + g 12 ) (σ + σ (2) where g 12 is the direct coupling between the two qubits, g ic is the bare coupling between qubit i and the coupler, J(Φ ext ) is the indirect coupling between the qubits, and ω(i) q (Φ ext ) is the Lamb-shifted qubit frequency which also tunes with Φ ext : Because the coupler frequency is much higher than the qubit frequencies in this device (∆(Φ ext ) < 0), the two sources of coupling can have opposite sign and the static ZZ interaction can therefore be eliminated by biasing Φ ext .We use Ramsey measurements and JAZZ measurement [74] to measure the qubit frequencies and static ZZ rate ζ along the tuning curve, shown in Fig. 3(b).
Biasing Φ ext to 0.23Φ 0 achieves ζ = 0, but exposes the qubits to greater flux noise than at Φ ext = 0, as |∂ω We therefore operate at Φ ext = 0 and find ζ/2π = 28.1 kHz.However, since ζ is much smaller than the typical drive amplitudes for our t p = 250 ns pulses (Fig. 2(b)), we expect that the corrections will be negligible.
We next turn to the dynamical properties of the coupler (see App. F for details).Modulating the coupler flux Φ ext with frequency ω Φ near (ω q )/2 and amplitude ε generates the following effective two-qubit Hamiltonian [42]: + σ (2) + ).( 31) To measure the exchange interaction generated by Eq. (31) we prepare the two qubits in the state |10 (|Q 1 Q 2 ) and measure the populations for varying pulse lengths and modulation frequencies around (ω q )/2.The resulting population oscillations between the two states are shown in Fig. 3(c), and a linecut shows the populations of the |01 and |10 states oscillating out of phase with each other, taking 82 ns to share the excitation equally between the two qubits, generating a maximally entangled state.
By calibrating the phase and frequency of the coupler drive, we use Eq.(31) to generate an approximate entangling Hamiltonian H in the Pauli basis as where . By further calibrating the amplitude and pulse duration, we then generate the so-called XY gates [45] parameterized by two angles β and θ: for β ∈ [0, π) and θ ∈ [0, 2π), which correspond to the phase and amplitude, respectively.
In the next section, we apply the Hamiltonian reconstruction algorithm to the two qubits evolving under the exchange interaction generated by modulation of the coupler.For several choices of (β, θ), we evaluate the performance of the algorithm, connect the deviations in the recovered dynamics to the calibration, and evaluate the fidelity of the corresponding two-qubit gate.

V. TWO QUBIT HAMILTONIAN RECONSTRUCTION
To validate this technique on a system of two qubits, we apply target Hamiltonians to both Q1 and Q2 and apply separate weak measurement tones on R1 and R2.We explore both non-interacting and interacting Hamiltonians but reserve the results of the former for App.G.The pulse sequence for reconstructing interacting Hamiltonians is shown in Fig. 4(a).For each Hamiltonian, we prepare 16 initial states {|+X , |+Y , |±Z } ⊗2 .Note that each initial state is a product state, which can be prepared quickly and with high fidelity using single qubit pulses.Target Hamiltonians H t (t) are generated by modulating the coupler for a duration of t p = 250 ns.Full tomography of the final state is performed for comparison with the output of the reconstruction algorithm.For each target Hamiltonian and initial state, we average ∼ 2.7M trajectories to generate the input for the algorithm.
Modulating the coupler as described in Sec.IV ideally produces nontrivial action only in the single-excitation subspace as in Eq. (32) and can be summarized by the Bloch sphere with poles corresponding to {|01 , |10 } shown in Fig. 4(b).The three ideal paths shown correspond to pulses as follows: (i) a π/2 pulse around the x axis with flat top and cosine ramp edges, generating XY (0, −π/2), (ii) a -π/2 pulse around the y axis, generating XY (π/2, π/2), (iii) a -π pulse around the x axis, generating XY (0, π).
A general two-qubit Hamiltonian H t (t) takes the form and contains 15 nontrivial amplitudes of which we reconstruct 12.The four most relevant reconstructed amplitudes for the three pulses are shown in Fig. 4(c), as we find negligible signal in the others.For the full set of 12 terms, see App.G.We now discuss the results for each Hamiltonian in turn.Hamiltonian (i) (shown in blue in Fig. 4(b)) is generated by the term σ x σ (2) y .We therefore expect that Ω XX (t) = Ω Y Y (t) and that Ω XY (t) = Ω Y X (t) = 0. We estimate an expected amplitude profile by measuring the rate of |01 ↔ |10 oscillations at the drive power used, similar to Fig. 3(d).We note that Eqs. ( 31) and (32) suggest that Ω XX (t) ∝ ε(t) 2 , so the amplitudes will be approximately the square of the waveform.Fig. 4(c) confirms that the reconstructed Hamiltonian agrees well with this expected amplitude profile.We suspect that the small reconstructed amplitudes Ω XY (t) and Ω Y X (t) are due to a miscalibration in the phase of the coupler pulse.The reconstruction infidelity in Fig. 4(d) also confirms that the final state as predicted by the algorithm agrees with the experimental tomography with a fidelity of 97.6%.
We next investigate how this technique performs on XY and Y X terms by driving the coupler to generate the Hamiltonian (ii), which is proportional to σ x (shown in orange in Fig. 4(b)).Again, Fig. 4(c) shows that the amplitudes of the XY and Y X terms are in good qualitative agreement with the expected amplitude shapes and that the XX and Y Y terms are not driven.This suggests that, compared to the pulse used for Hamiltonian (i), the phase of the coupler pulse was more accurately calibrated.We find a reconstruction fidelity of 97.2%, providing a second validation of the algorithm.
Finally, we investigate the performance of this technique when Ω IZ and Ω ZI are non-negligible.We drive the coupler to try to generate Hamiltonian (iii) (shown in purple in Fig. 4(b)) and find that although the qualitative agreement in Fig. 4(c) is good, the quantitative agreement in Fig. 4(d) is only 95%, noticeably worse than for the slower pulses.To determine the source of the error, we precondition the algorithm with nonzero values for Ω IZ and Ω ZI .By optimizing the reconstruction fidelity over these values, we find that choosing Ω IZ /2π = 763 kHz and Ω ZI /2π = 392 kHz improves the reconstruction fidelity by 1.3%.On the other hand, performing this optimization for Hamiltonians (i) and (ii) improves the reconstruction fidelity by only 0.1% and estimates Ω IZ /2π = 191 kHz Ω ZI /2π = 112 kHz for case (i) and Ω IZ /2π = 120 kHz Ω ZI /2π = 60 kHz for case (ii).We can understand the source of these single qubit Pauli terms as follows.Modulating the coupler with amplitude ε induces dynamical single qubit frequency shifts (σ Z ) terms proportional to ε 2 (∂ 2 ω(i) q /∂Φ 2 ext ) (evaluated at Φ ext = 0) [42].When the coupler is driven resonantly, these shifts vanish in an appropriate frame, yielding Eq. ( 31), but when it is driven off-resonantly, their difference cannot be eliminated completely.The improved reconstruction fidelity for pulse (iii) suggests that the modulation frequency of the pulse was off-resonant by ((Ω IZ − Ω ZI )/4)/2π = 93 kHz, resulting in an imperfect XY (0, π) pulse.While the relative values of Ω IZ and Ω ZI stem from imperfect calibration of the coupler modulation frequency, their absolute value is due to imperfect calibration of virtual Z phases [75] following the pulse and the optimization of these amplitudes with respect to post-pulse tomography.Using the Hamiltonian reconstruction, we can thus isolate sources of coherent error at the Hamiltonian level and connect them to errors in the controls, such as phase and frequency miscalibration.This experiment also demonstrates how significant single qubit frequency shifts in (iii) lead to model error due to using a first-order algorithm, which can be remedied by preconditioning.
The reconstructed time-dependent Hamiltonian H r (t) also allows us to compare the coherent dynamics of the qubits to those generated by a theoretical control Hamiltonian H c (t), such as one found by an optimal control algorithm [20,21,23].To quantify how different the dynamics produced by H r (t) and H c (t) are, we introduce the dynamical coherent fidelity F(H r (t), H c (t)) as follows: where U r (t) is the unitary generated by the reconstructed Hamiltonian, U c (t) is the unitary generated by the theoretical control Hamiltonian, and T is the time-ordering operator.Operationally, it measures on average how different the actions of the two unitaries are when acting on a pure state, during the course of evolution.Note that F(H r (t), H c (t)) need not decrease monotonically with t, as multiple Hamiltonians can ultimately generate the same final unitary U (t p ).We plot the dynamical coherent infidelity 1 − F(H r (t), H c (t)) in Fig. 4(e) for theoretical control Hamiltonians which give the three ideal trajectories on the single excitation Bloch sphere in Fig. 4(b).We observe that pulses (i) and (ii) steadily diverge from their ideal Hamiltonians due to imperfect phase and frequency calibrations and the dynamical single qubit frequency shifts to achieve 98.4% and 98.7% final coherent gate fidelity, respectively.However, in these cases, the final coherent gate fidelities do not differ appreciably from their lowest values during the pulses.Instead, we demonstrate in the inset of Fig. 4(d) how F(H r (t), H c (t)) reveals a significant divergence from H c (t) which the final gate fidelity cannot detect by evaluating it for the single qubit flat-top π pulse with a superimposed sinusoidal pulse (Fig. 2(b)(iii)).Although the final coherent gate fidelity is 99.4%, the dynamical coherent fidelity falls as far as 94.8% at t = 142 ns due to the superimposed sine pulse.This "worst case" value quantifies how far the true Hamiltonian strayed from the theoretical control Hamiltonian and indicates suboptimal control despite ultimately achieving high fidelity.Therefore, the full dynamical coherent fidelity calculated from the reconstructed Hamiltonian reveals more information about how closely the experimental Hamiltonian matches the theoretical control Hamiltonian than can be obtained from only initial and final measurements and enables us to quantify the discrepancy.

VI. CONCLUSIONS
We have introduced a technique for reconstructing the time-dependent Hamiltonian of an interacting system using continuous weak measurements and have experimentally applied and validated it for systems of one and two superconducting transmon qubits.We have shown that we are able to reconstruct nontrivial two-qubit Hamiltonians with a time resolution of 4 ns and achieve qualitative agreement with expected values and quantitative agreement between the predicted final state and experimental tomography of the final state.
We now comment on the limitations of this technique and how they may be overcome.First, performance of the algorithm suffers if terms in the Hamiltonian which are not calculated by the algorithm are significant.In our experiments, we are most likely to suffer from this when there are significant dynamics in the IZ, ZI, and ZZ terms.We demonstrated in Sec.V that when the coupler was driven with higher power to perform XY (0, π), imperfect calibration of the modulation frequency produced significant IZ and ZI dynamics, and by preconditioning Ω IZ (t) and Ω ZI (t) to be nonzero, we could improve the reconstruction fidelity.This limitation stems from the fact that calculating amplitudes for these three terms in the algorithm requires an update equation which is second order in time, which amplifies shot noise.Second-order update may therefore be achievable with more repetitions or with faster pulses and stronger measurement, provided that the adiabatic condition still holds.Naively increasing the number of repetitions requires the experiment to remain stable for longer periods of time.However, recent advances in active qubit reset [76] can be used to increase the repetition rate, thus reducing the total time requirement.Finally, although we reconstruct the time-dependent Hamiltonian from ensemble averages in this work, developing a reconstruction technique which uses a small set of noisy trajectories is an interesting future direction, as it may prove more efficient in terms of required repetitions.
We also showed that when the adiabatic condition 2Ω κ is violated, the resonator field does not follow the qubit as it evolves under the target Hamiltonian.This results in a poor estimate of the z coordinates, which leads to poor performance of the reconstruction algorithm.To reconstruct faster dynamics, it may be possible to increase the coupling κ, but the resonator and Purcell filter must be designed carefully to preserve qubit lifetimes.Alternatively, the technique developed in [69] could be used to reconstruct the qubit state even when the adiabatic condition is violated.The reconstructed states could then be used for Hamiltonian reconstruction.
The most immediate application of Hamiltonian reconstruction is to study novel quantum gates, particularly in nascent quantum platforms which do not yet have standard gate sets.For example, classical optimization can be used to estimate an optimal control Hamiltonian for a gate, and our technique can be used to validate it experimentally.By extending this technique to qutrits, one could further investigate leakage errors during gates.Second, our technique may be used to validate analog quantum simulators with locally tunable interactions by isolating subsystems and performing Hamiltonian reconstruction.Finally, our technique is well-suited for quantum metrology in which the qubit behaves as a sensor and drive amplitudes come from the environment rather than user control.

A. EXPERIMENTAL SETUP
The experiment presented in this work was performed with a chip of three superconducting transmon qubits, one of which was not used.The chip was cooled to 25 mK in an Oxford Instruments dilution refrigerator.The cryogenic wiring and room temperature electronics used in the setup are shown in Fig. 5. Single qubit pulses are generated by upconversion of IF pulses generated using a Tektronix 5014c Arbitrary Waveform Generator (AWG) at 1 GSa/s through IQ mixing with a 5.5 GHz local oscillator (LO) tone sourced by a Keysight MXG N5183B Signal Generator (SG).Pulses are then amplified and low-pass filtered at room temperature.Readout pulses are also generated by the same AWG and similarly upconverted by mixing the signal from the AWG with a 6.8 GHz LO tone sourced by an Agilent E8257D SG and then amplified at room temperature.Finally, low-frequency tunable coupler modulation pulses are generated directly by the AWG.All pulses pass through a DC block, are then attenuated at each temperature stage, and finally are low-pass filtered at the base stage before being passed to the chip.DC bias for the flux-tunable coupler is sourced by a Yokogawa GS200.The experiment is performed at a bias of -0.0752 mA at which point the qubit frequencies are first-order insensitive to the the coupler flux and we observe the longest coherence times.
The reflected measurement signals are redirected through cryogenic circulators, amplified by a traveling wave parametric amplifier (TWPA) at 25 mK, a high electron mobility transistor (HEMT) at 4 K, and a lownoise room-temperature amplifier.The room temperature signal is then downconverted to I and Q components at the IF by mixing with the same 6.8 GHz LO tone used to generate the readout pulses.The IF signals are then filtered and digitized at 1 GSa/s by an Alazar ATS 9373 Analog-to-Digital Converter (ADC).Finally, the digitized signals are demodulated in software.
The TWPA is pumped by a 7.98 GHz tone generated by a Hewlett-Packard 83732B SG.The pump tone is lowpass filtered at the base stage before entering the TWPA.The pump amplitude and frequency are chosen to maximize gain in the reflected readout signals at the two frequencies used to measure the qubits.Most of the power of the pump tone is used as a nulling tone to cancel the pump tone in the reflected readout signal at room temperature.This is achieved by phase-shifting and attenuating the nulling tone and then coupling to the amplified, reflected readout signal before downconversion to IF.

B. DEVICE PARAMETERS
The sample consists of three fixed-frequency transmons of which the two used in this work are coupled directly as well as through a flux-tunable coupler.Each trans-mon is coupled to a λ/4 coplanar waveguide resonator (CPW) used for dispersive readout.The resonators are further coupled to a Purcell filter which is then coupled to the feedline used to measure the resonator fields.We define and report the most relevant parameters for the two qubits used in this work in Table I.The resonator frequency ω res and linewidth κ are measured by spectroscopy.T 1 is measured by repeatedly preparing a qubit in the excited state and measuring the qubit population at different times.T * 2 and ω q are measured by repeated Ramsey experiments on the |0 → |1 transition, and the anharmonicity α q is measured by Ramsey measurements on the |1 → |2 transition.
Acquiring a dataset for a given pulse requires repeatedly preparing qubits in the chosen initial states and performing the pulse sequence presented in the main text.For our single (two) qubit datasets, we average over ∼ 5.4M (∼ 2.7M) repetitions, taking 8 hours.This method is limited by residual shot noise in the data, which in theory can be reduced by more measurements.Practically, however, stability may limit the number of repetitions that can be carried out.We find that for our system, T 1 and T * 2 fluctuate by 13% and 28%, respectively, over the acquisition time but still remain much longer than the duration of each repetition (∼ 500 ns).

C. RECOVERING QUBIT COORDINATES FROM CONTINUOUS WEAK MEASUREMENTS
In this section we derive analytically how the resonator response encodes information about the z coordinate of a qubit.The starting point for our derivation is the driven, dispersive qubit-resonator Hamiltonian in the rotating frame of the qubit and resonator drives.For simplicity, we assume that both drives are on resonance: where H q is the single-qubit Hamiltonian, Ω x (t) and Ω y (t) are the drive envelopes, H r is the resonator Hamiltonian, Ω wm is the measurement strength, H qr is the dispersive coupling between the qubit and resonator, and χ is the full dispersive shift.
The resulting Heisenberg equation of motion is It can be checked that the following implicit equation is a solution to Eq. (C5): This recurrence relation has solution We can understand this equation as follows: the resonator state is a series of filtered, time-ordered averages of the qubit z coordinate.At lowest non-trivial order, this gives Expanding the integrand, σ z (t−τ ) = σ z (t)− σz (t)τ +. . .: where τ ≡ 2/κ, Γ(x, y) is the lower incomplete gamma function, and boundary terms proportional to e −t/τ have been dropped.This result reveals the following intuition.Provided that the qubit state does not evolve too rapidly compared to the linewidth of the resonator, τ , the real part of the resonator state is proportional to the qubit state.Finally, the qubit state which appears in Eq. ( C10) is shifted by τ .In this experiment, the resonator linewidth (κ/2π) qubit 1(2) is 11.78 MHz(14.47MHz) which corresponds to a shift of 27 ns(22 ns).When we use the resonator field to estimate the ensemble average z component of the Bloch vector z(t) as the input for the Hamiltonian reconstruction algorithm, we must rescale the voltage records by identifying the values that correspond to z = ±1 as well as shift the records by τ .The approximation in Eq. (C10) applies when ż(t) < τ = 2/κ.The fastest source of qubit evolution arises from H q , so | ż(t)| ≤ Ω max = max 0≤t≤tp max(|Ω X (t)|, |Ω Y (t)).This leads us to the "adiabatic condition", 2Ω max < κ which ensures that the resonator field follows the qubit evolution.We show in Fig. 2(d-e) that when the qubit is driven too rapidly compared to the resonator linewidth, the resonator field "averages" the qubit evolution too much to be able to recover the qubit state from it, and the algorithm consequently produces a poor estimate of the Hamiltonian.

D. RECONSTRUCTION PROCEDURE
In this Appendix, we first give a more explicit description of the first-order update algorithm for a single qubit, which calculates the drive amplitudes Ω X (t) and Ω Y (t).
We then outline all of the steps of the algorithm.Finally, we introduce a second-order algorithm which solves for all three drive amplitudes.

First-Order Update Algorithm
Eq. ( 1) describes the coherent evolution of a qubit expressed in the Pauli operator basis.Eq. ( 13) can be used to reconstruct the drive amplitudes when all three components of the Bloch vector are measured at all times, but Eq. ( 17) is more useful in our context, as only the qubit z coordinate is measured simultaneously with the target Hamiltonian.We therefore concentrate on detailing the derivation and algorithm for this case.Below we denote the time at discrete time step n as t n .The Heisenberg equation of motion for σ z is where Γ 1 = γ ↓ + γ ↑ = 1/T 1 is the sum of the rates of spontaneous relaxation and absorption, and Γ ∆ = γ ↓ − γ ↑ .Discretizing Eq. (D1) at lowest order in ∆t, Eq. (D2) is linear in the drive amplitudes and can be solved simply if we always have the Bloch vector (x(t n ), y(t n ), z(t n )) and z(t n+1 ) for at least two states such that (y(t n ), −x(t n )) are linearly independent.We achieve this by performing the continuous weak mea-surement with the target Hamiltonian on S ≥ 2 initial states.Then, labeling z s (t n ) as the z coordinate of the qubit evolved to time t n from the s-th initial state, we can rewrite Eq. (D2) in terms of the components of the Bloch vector for different initial states as The S × 2 matrix in the first term of the right hand side (RHS) is the matrix M S (t n ) in Eq. ( 16), and has a left inverse when it has rank 2, so Eq.(D3) can be solved for the vector of amplitudes (Ω X (t n ), Ω Y (t n )) T by standard linear algebra and is equivalent to least-squares regression.We now outline the steps of the algorithm.We assume we are given {ρ 1 (0), ρ 2 (0), . . .ρ S (0)}, the full density matrix of each initial state, the decoherence rates Γ 1 , Γ d , and z(t n ) for 1 ≤ n ≤ N timesteps (t 0 = 0).The algorithm then proceeds as follows, starting at n = 0: 1. Solve for the drive amplitudes Ω X (t n ) and Ω Y (t n ) by solving Eq. (D3).
2. Integrate the Lindblad master equation for the density matrix in Eq. ( 4) to obtain the full density matrix at timestep t n .This provides access to all of the components of the Bloch vector at t n for computing the drive amplitudes at the next time step.The approach given thus far suffers from the absence of Ω Z (t) in Eq. (D3).By proceeding to an update equation which is second order in ∆t, we can solve for Ω Z (t).We use the Heisenberg equations for σ x (t) and σ y (t): where Γ 2 = Γ d +Γ φ +Γ 1 /2 and Γ φ is the intrinsic dephasing rate.Eq. (D2) can be promoted to a second-order equation as follows: We now treat Ω X (t n ), Ω Y (t n ), and ρ(t n ) as "known" and solve for the three unknown components Ω X (t n+1 ), Ω Y (t n+1 ), and Ω Z (t n ).
To satisfy the base case, we use Eq.(D3) to solve for Ω X (t 0 ), Ω Y (t 0 ) and as before assume that we have tomography of the initial state ρ(t 0 ).However, Eq. (D6) is now nonlinear in the unknowns, containing terms like Ω X (t n+1 )Ω Z (t n ) and Ω Y (t n+1 )Ω Z (t n ), and requires a nonlinear solver, such as Gauss-Newton method, gradient descent, or other methods of nonlinear regression.This method requires S ≥ 3 equations to solve for the amplitudes at each time step.
One limitation of this method is that when coupling between σ z (t) and the other qubit coordinates vanishes, i.e.Ω X (t) = Ω Y (t) = 0 for some t, z(t) is totally unaf-fected by Ω Z (t) and cannot be recovered from Eq. (D6).However, if they are known to vanish for all t, then we can apply a known spectroscopy drive, i.e. a drive around the x or y axes with known nonzero amplitude, to reintroduce coupling to Ω Z (t).A second drawback is that this method is second order in ∆t, as we are effectively computing the second derivative of z(t).However, each derivative we take amplifies high frequency noise.In our experiment, we find that applying the second order update rule results in a poor estimate of the target Hamiltonian.However, this is not a fundamental limitation and can in theory be overcome by more averaging and longer acquisition, provided that the system remains stable.
We now comment on how this technique extends to a system of two qubits.At first order, Eq. (D2) for each qubit now contains eight terms, with four in common (σ (1) x σ (2) x , σ x σ (2) y , σ y σ (2) x , and σ (1) y σ (2) y ).The first order update requires a minimum of six states, as each yields two linearly independent equations.Terms which are tensor products of only identity and Pauli z operators, i.e. σ do not appear at first order but do appear at second order, similar to Eq. (D6).More generally, for a system of Q qubits (total Hilbert space dimension d = 2 Q ) each continuously measured along the z axis, a minimum of d(d − 1)/Q initial states to apply the first order update to recover d 2 − d out of d 2 − 1 amplitudes.By going to second order we can recover all of the drives using d 2 −1 initial states.As before, we can use a spectroscopy drive to introduce deliberate coupling if necessary.Moreover, this technique requires only single qubit drives to ensure coupling to all drive terms.
Finally, we note that the formalism introduced can be extended to the setting where multiple qubit coordinates are continuously measured.For example, if two coordinates of a qubit, such as y(t) and z(t), are known for all t, then a first-order update method suffices and all three amplitudes can be recovered, provided that S ≥ 2 initial states are used.It can also be extended to higher dimensional systems by expanding the Hamiltonian in a complete basis, such as the generalized Gell-Mann matrices.

Efficient Reconstruction for Fast Evolution
Eqs. (D3) and (D6) provide an efficient means for estimating the amplitudes of the target Hamiltonian at each time step.However, their accuracy suffers at large drive amplitudes (relative to ∆t).On the other hand, minimizing an objective function which calls a matrix exponentiation function is computationally intensive.We combine the two approaches by splitting the Hamiltonian into a "fast" initial guess H g (t) and a "slow" perturbation δH(t), where the speeds refer to Ω/(1/∆t), and optimize the perturbation δH(t).We then view the RHS of the update equations as estimating ż(t) by żest (t; H g (t), δH(t)), a function of H g (t) and δH(t), and minimize the squared sums (over initial states) of żest (t; H g (t), δH(t)) − ż(t) over δH(t).The update equations in the Heisenberg picture are now written as and we optimize over δH(t) expressed in the Pauli operator basis.Since the matrix exponential does not depend on δH(t), it is not evaluated repeatedly by the optimizer.We find that this approach yields more accurate reconstruction, particularly when the dynamics are fast compared to the timestep ∆t.In our experiment, we are limited by the resonator linewidth κ due to the requirement that the resonator field adiabatically follow the qubit state, and it is not strictly necessary to use this approach.However, if the resonator linewidth far exceeds the time resolution of the data after filtering, then this technique is essential to model and reconstruct the dynamics accurately.

E. DETAILS OF DATA PROCESSING AND FILTERING
The Hamiltonian reconstruction algorithm requires as input the full density matrix of each initial state used, {ρ(0)}, the qubits' z(t) coordinates during H t (t), and finally an estimate of the rates of measurement-induced dephasing {Γ m }.In this section, we detail the steps taken to calculate these z(t) from the resonator field.
The output of the experiment is an ensemble of voltage traces measured with 1 ns resolution, which are averaged to estimate the ensemble averaged resonator field.We showed in App.C that this amplitude is proportional to the qubit z coordinate retarded by a delay τ = 2/κ.However, the data must first be filtered, deconvolved, and rescaled.Filtering is necessary as weak measurements leave residual high frequency noise due to finite averaging, and it is particularly necessary for weak measurement of both qubits, as we use a single feedline to si-multaneously probe both qubits' resonators.To achieve good readout isolation, we have designed the resonator frequencies in this device to be 102 MHz apart and lowpass filter each qubit's data to suppress dynamics of the other qubit as well as other high frequency shot noise.We downsample to 2 ns time resolution and use a thirdorder low-pass Butterworth filter with a 50 MHz critical frequency to smooth the data.Then, to determine the overall scaling to convert voltages to z(t), we also prepare the qubits along the z axis and measure the field without applying the target Hamiltonian.We use the difference of the mean values of these two measurements to rescale the voltages to lie in [−1, 1].Values that extend beyond these bounds are clipped.This produces the z(t) input for the Hamiltonian reconstruction algorithm.The reconstructed amplitudes {Ω IX (t), Ω IY (t) . . .Ω ZY (t)} are also filtered using a 5th order low-pass Butterworth filter with critical frequency 50 MHz, as the update methods described in Sec.D amplify high frequency noise.Naturally, filtering the data more aggressively smooths the inputs but can distort the Hamiltonian.We demonstrate this effect in Fig. 6 for single qubit Hamiltonian reconstruction applied to a π pulse and two qubit Hamiltonian reconstruction applied to simultaneous π pulses.When the critical frequency is too low, the voltage data and consequently the input to the algorithm z(t) are distorted.The algorithm then poorly reconstructs the Hamiltonian and fails to predict the final state of the system accurately.It is therefore necessary to keep the critical frequency high and instead decrease shot noise through more averaging, stronger measurement, or more efficient measurement (e.g. through quantum-limited amplifiers).

F. XY INTERACTIONS THROUGH PARAMETRIC MODULATION
In this Appendix, we review how parametric modulation of the flux-tunable coupler generates an entangling exchange interaction between the two qubits and derive its particular form.Our starting point for the derivation is the effective two-qubit Hamiltonian in the limit where the coupler frequency far exceeds the qubit frequencies, i.e. ∆ i (Φ ext ) = ω (i)  q − ω c (Φ ext ) (F1) Σ i (Φ ext ) = ω (i) q + ω c (Φ ext ) (F2) where ω c (Φ ext ) is the coupler frequency, which tunes with the flux Φ ext , and g ic is the direct capacitive couplings between the i-th transmon and the coupler.The coupler is also assumed to remain in its ground state.The effective Hamiltonian H eff is then Hqq / = (J(Φ ext ) + g 12 ) σ (1)  x σ (2)   x (F6) where ω(i) q (Φ ext ) is the qubit frequency, which now tunes with the coupler flux Φ ext , g 12 is the direct capacitive coupling between the transmons, J(Φ ext ) is the indirect coupling mediated by the coupler, and 1/∆(Φ ext ) = 1/∆ 1 (Φ ext ) + 1/∆ 2 (Φ ext ) (Σ(Φ ext ) is defined similarly).For a derivation of H eff from circuit quantization, see [44,49].

G. EXTENDED DATA FOR TWO QUBIT HAMILTONIAN RECONSTRUCTION
In this Appendix, we provide data for the extended Hamiltonian reconstruction of two qubits when the Hamiltonian is generated by single qubit pulses and when it is generated by a pulse on the coupler flux line.
In Fig. 7, we use the two-qubit Hamiltonian reconstruction to measure whether crosstalk occurs in this device when applying a π pulse on one of the two qubits at a time.We reconstruct 12 of the 15 possible terms in a twoqubit Hamiltonian (the remaining terms being Ω IZ (t), Ω ZI (t), and Ω ZZ (t)) and find no appreciable crosstalk.We also show the result for when no pulse is applied on either qubit.We find some residual fluctuations which may be an artifact of filtering, so we use the RMS value of these curves to estimate an uncertainty for reconstructions of nontrivial pulses, shown as the shaded band extending above and below each curve.
In Fig. 8, we show the complete reconstruction data for the entangling two-qubit Hamiltonians generated by modulating the coupler flux, after preconditioning Ω IZ (t) and Ω ZI (t) as described in the main text.We find no significant amplitudes outside the main expected four terms.

FIG. 1 .
FIG. 1. Hamiltonian reconstruction method.(a) To reconstruct an unknown target Hamiltonian Ht(t) of a system composed of two fixed frequency transmons (Q1, Q2) coupled via flux-tunable coupler (C12), we continuously measure the qubits while applying pulses to generate Ht(t).Single qubit target Hamiltonians are generated by pulses of duration of tp on the charge lines.Entangling two qubit target Hamiltonians are generated by applying RF pulses on the flux line (FF).In the adiabatic limit, the noisy resonator output field (schematically depicted in (b)) follows the evolution of the qubit z coordinate during the pulse.(c) We average, filter and deconvolve voltage records V1,2 and record the resulting qubits' z evolution for S ≥ 2 different initial states while keeping the qubit/coupler pulses fixed.The resulting ensemble averages z1,2(t) along with the calibrated measurement-induced dephasing rate Γ d are inputted to the Hamiltonian reconstruction algorithm.(d) The output of the algorithm is an estimate of the time-dependent target Hamiltonian expanded in the Pauli operator basis.

FIG. 3 .
FIG. 3. Two qubit gate operation.(a) Two fixed-frequency transmon qubits (Q1, Q2) are coupled to a tunable coupler (SQUID).The flux line FF is used to DC bias and modulate the coupler flux Φext.(b) By DC biasing the coupler, the qubit and bus frequencies and couplings can be tuned.The coupler is biased at Φext = 0 (marked by red star) where the qubit frequency (fit to Eq. (30)) is first-order insensitive to the flux (inset).The static ZZ rate ζ is measured using a JAZZ measurement[74], and the data are fit using models described in App.F. (c) Population oscillations of the |10 state due to driving the flux line (FF) with an AC flux drive at frequency ωΦ near 1 2 (ω FIG. 3. Two qubit gate operation.(a) Two fixed-frequency transmon qubits (Q1, Q2) are coupled to a tunable coupler (SQUID).The flux line FF is used to DC bias and modulate the coupler flux Φext.(b) By DC biasing the coupler, the qubit and bus frequencies and couplings can be tuned.The coupler is biased at Φext = 0 (marked by red star) where the qubit frequency (fit to Eq. (30)) is first-order insensitive to the flux (inset).The static ZZ rate ζ is measured using a JAZZ measurement[74], and the data are fit using models described in App.F. (c) Population oscillations of the |10 state due to driving the flux line (FF) with an AC flux drive at frequency ωΦ near 1 2 (ω (1) q − ω (2) q ).(d) Populations of the |10 and |01 states after preparing the qubits in |10 and driving the coupler at ωΦ/2π = 23.7 MHz for various pulse durations.The populations oscillate out of phase, indicating Q1 and Q2 are swapping an excitation.A full swap of an excitation occurs when θ = π (as in Eq. (33)), taking 163 ns.

FIG. 4 .
FIG. 4. Two qubit Hamiltonian reconstruction.(a) Pulse sequence.The coupler RF pulse generates an entangling two-qubit Hamiltonian.(b) We study three different pulses XY (β, θ) whose ideal action in the {|01 , |10 } subspace is schematically depicted on the Bloch sphere.(c) Reconstructed Hamiltonians for XY (0, − π2 ) (blue), XY ( π 2 , π 2 ) (orange), and XY (0, π) (pink) (∆t = 4 ns).The shaded band around each Pauli term represents the confidence interval, obtained from the average Pauli amplitude when no coupler pulse was applied.Single qubit and additional two qubit Pauli terms have no appreciable signal (see App. G) (d) The reconstruction infidelity (averaged over initial states) computed with respect to tomography after the pulse is small for all three pulses.Preconditioning the Hamiltonian reconstruction algorithm with an IZ and ZI term (magnitude derived from QPT) reduces this error further (right, shaded bars), particularly for XY (0, π).(e) Dynamical gate infidelity during the three pulses showing the accuracy of the reconstructed dynamics with respect to the assumed pulse.Inset shows dynamical gate fidelity for pulse (iii) in Fig.2(b) and reveals large deviation from an "ideal" control Hamiltonian despite high final gate fidelity.

3 FIG. 5 .
FIG. 5. Experimental setup.Details and specification of larger electronics are given in Appendix A.

TABLE I .
Calibrated parameters of the two transmon qubits and their respective resonators.All parameters are measured in the absence of any measurement tone.Uncertainties reflect standard deviation of parameters measured during 8 hours of data acquisition for each experiment.
3. Increment n and proceed to Step 1.The output of this algorithm is the drive amplitudes at N timesteps t 0 , t 1 , . . .t N −1 as well as the density matrix FIG.6.Average reconstruction fidelity vs. critical frequency of low-pass Butterworth filter applied to the data.Dashed line is provided as a guide for the eye.The reconstruction fidelity, averaged over initial states, Fr falls when the time series voltage data is filtered too strongly and distorts the dynamics.
AC is an AC shift of the qubit frequency and all derivatives are evaluated at Φ ext = Φ DC .By similarly expanding J(Φ ext ), we obtain a modulated coupling: Φ t + φ)) σ