# Complete physical characterization of QND measurements via tomography

L. Pereira<EMAIL_ADDRESS>Instituto de Física Fundamental IFF-CSIC,
Calle Serrano 113b, Madrid 28006, Spain  J.J. García-Ripoll  Instituto de
Física Fundamental IFF-CSIC, Calle Serrano 113b, Madrid 28006, Spain  T. Ramos
<EMAIL_ADDRESS>Instituto de Física Fundamental IFF-CSIC, Calle
Serrano 113b, Madrid 28006, Spain

###### Abstract

We introduce a self-consistent tomography for arbitrary quantum non-demolition
(QND) detectors. Based on this, we build a complete physical characterization
of the detector, including the measurement processes and a quantification of
the fidelity, ideality, and back-action of the measurement. This framework is
a diagnostic tool for the dynamics of QND detectors, allowing us to identify
errors, and to improve their calibration and design. We illustrate this on a
realistic Jaynes-Cummings simulation of superconducting qubit readout. We
characterize non-dispersive errors, quantify the back-action introduced by the
readout cavity, and calibrate the optimal measurement point.

#### Introduction.-

Quantum non-demolition (QND) detectors measure an observable preserving its
expectation value [1, 2]. This property is essential in back-action-free
quantum metrology [3, 4, 5, 6, 7], and in quantum protocols with feedback,
e.g. fault tolerant computation [8, 9, 10, 11, 12]. QND measurements are
typically implemented indirectly by monitoring a coupled subsystem, as
demonstrated in various AMO [13, 14, 15, 16, 17] and solid-state [18, 19, 20,
21, 22, 23] platforms. In superconducting circuits, the standard qubit
measurement is a dispersive readout [24] mediated by frequency shifts in an
off-resonant cavity. This is a near-QND process that approximately preserves
the qubit’s polarization [cf. Fig 1(b)], and enables rapid high-fidelity
single-shot measurements [25, 26], protection by Purcell filters [27, 28, 29,
30], fast reset [31], and simultaneous readout through frequency multiplexing
[32, 28, 30].

State-of-the-art QND detectors still face experimental challenges. A critical
problem is the exponential accumulation of back-action errors from repeated
applications of the detector, which limits the scaling of quantum
technologies. In superconducting qubit readout, such errors originate in
deviations from the dispersive limit in practical devices [33, 34, 35]. This
has motivated more complex QND measurement schemes [36, 37, 38, 39, 40, 41,
42, 43, 44], which also introduce their own sources of imperfection.

In order to make progress in the design and operation of QND measurements, we
need a complete and self-consistent diagnostic tool, which helps both with the
calibration of the detector and with describing its real dynamics. Many
experiments have focused on optimizing simple quantities such as the readout
fidelity and the QND-ness [43, 40]. However, these fidelities do not quantify
the QND nature of a measurement, but rather its projectivity and ideality as
shown below. Another standard approach is detector tomography [45, 46, 47, 48,
49]. This method characterizes destructive measurements via positive operator-
valued measurements (POVMs), but ignores the post-measurement state, and
therefore a description of the measurement back-action.

Figure 1: (a) QND detector tomography for generic measurements. A self-
consistent calibration requires sampling over input states $\rho_{k}$ and two
consecutive QND measurements, interleaved by a unitary operation $U_{j}$. This
allows us to generate the measurement processes $\mathcal{E}_{n}$, for each
possible outcome $n=1,\dots,N$, and to reconstruct them tomographically from
the conditional probabilities $p(n|k)$ and $p(mn|jk)$. (b) Setup for
tomographic characterization of QND qubit readout. It requires pulse control
on qubit $\Omega_{q}(t)$ and cavity $\Omega_{c}(t)$, as well as continuous
homodyne detection $\langle a+a^{\dagger}\rangle_{c}$. The detector can have
arbitrary qubit-cavity coupling $g$ and any imperfection such as qubit decay
$\gamma$ and dephasing $\gamma_{\phi}$.

In this work, we develop a complete physical characterization of QND
measurements and their back-action via quantum tomography. The protocol,
without pre-calibration of the QND detector, estimates both the POVM elements
and the quantum process operators associated to each measurement outcome. As
seen in Fig. 1(a), this requires two consecutive applications of the detector,
interleaved by unitary operations, and repeated over a set of input states.
The information contained in the process operators can be used to identify
errors, calibrate, and optimize the design of QND detectors. This can be done
either directly, or through the analysis of simple metrics such as readout
fidelity, QND-ness, and destructiveness, a precise bound of the measurement
back-action that we introduce below. The method can be applied to any
detector, but we illustrate its power simulating realistically the calibration
of superconducting qubit readout beyond the dispersive approximation. Our
study shows that state-of-the-art dispersive readout is a near-ideal
measurement at the optimum of QND-ness, but there are other regimes where it
becomes maximally QND with minimal back-action error, as revealed by the
destructiveness. Other tomographic approaches to non-destructive detectors
focus on near-ideal measurements only [50, 51], and thus do not provide a
characterization of their real QND nature.

#### General description of QND measurements.-

A non-destructive quantum measurement with $N$ outcomes is represented by $N$
completely positive maps $\mathcal{E}_{n}$, which add up to a trace-preserving
map $\mathcal{E}=\sum_{n}\mathcal{E}_{n}$. The $n$-th measurement outcome is
obtained with probability $p(n)=\operatorname{Tr}\\{\mathcal{E}_{n}(\rho)\\}$,
leaving the system in the post-measurement state
$\rho_{n}=\mathcal{E}_{n}(\rho)/p(n)$. Each $\mathcal{E}_{n}$ is unambiguously
represented by the Choi matrix $\Upsilon_{n}$, whose $d^{4}$ elements
$\Upsilon_{n}^{ijkl}$ read

$\displaystyle\Upsilon_{n}^{ijkl}=\bra{ij}\Upsilon_{n}\ket{kl}=\bra{i}\mathcal{E}_{n}(\ket{k}\\!\bra{l})\ket{j},$
(1)

with $\\{\ket{i}\\}$ a basis of the measured system with dimension $d$. The
Choi matrices give the post-measurement states
$\mathcal{E}_{n}(\rho)=\sum_{ijkl}\Upsilon_{n}^{ijkl}\rho^{kl}\ket{i}\\!\bra{j}$
111In this work, the Choi matrix $\Upsilon_{n}^{ijkl}$ encodes the
transformation of the density matrix [107], while
$\tilde{\Upsilon}_{n}^{ijkl}=\Upsilon_{n}^{ikjl}$ is a transposed version that
is Hermitian and non-negative. and also the measurement statistics
$p(n)=\operatorname{Tr}\\{\Pi_{n}\rho\\}$, with the POVM elements,
$\Pi_{n}=\sum_{ijk}\Upsilon_{n}^{kjki}|i\rangle\langle j|$. Conservation of
probability requires the completeness relation, $\sum_{n}\Pi_{n}=\openone$,
which imposes $\sum_{nk}\Upsilon_{n}^{kjki}=\delta_{ij}$ on the Choi
components.

A QND measurement of observable $O$ [2, 53, 54] is one where the unconditional
process $\mathcal{E}=\sum_{n}\mathcal{E}_{n}$ conserves the probability
distribution $p(n)$ over repeated measurements. Equivalently, $\mathcal{E}$
preserves the average of all compatible observables $O_{c}$

$\displaystyle\langle
O_{c}\rangle=\operatorname{Tr}\\{O_{c}\rho\\}=\operatorname{Tr}\\{O_{c}\mathcal{E}(\rho)\\},\qquad\forall[O_{c},O]=0.$
(2)

Ideal measurements are well known QND measurements where consecutive
detections project the system onto the same eigenstate of $O$. This requires
all Choi matrices to be projectors $\Upsilon_{n}=(\Upsilon_{n})^{2}$, a
sufficient condition to satisfy Eq. (2). However, as shown below, general QND
measurements are not ideal, and allow the state on consecutive detections to
change, provided the averages $\braket{O_{c}}$ remain constant.

#### Tomographic reconstruction of measurement processes.-

We have developed a self-consistent characterization of QND measurements,
based on two consecutive applications of the detector, interleaved with
unitary operations from a universal set of gates $U_{j}$ [cf. Fig.1(a)]. The
first measurement induces the processes $\mathcal{E}_{n}$, conditioned on the
detected outcome, while the unitary $U_{j}$ and the second measurement are
used for a process tomography of the detector itself. By repeating the
protocol over a set of input states $\rho_{k}$, we obtain the conditional
probabilities $p(n|k)=\operatorname{Tr}(\Pi_{n}\rho_{k})$ and
$p(mn|jk)=\operatorname{Tr}(\Pi_{m}U_{j}\mathcal{E}_{n}(\rho_{k})U_{j}^{\dagger})$
after the first and second measurements, respectively. From these
distributions we reconstruct the POVMs $\Pi_{n}$ and Choi matrices
$\Upsilon_{n}$ that best approximate the measurement, without a pre-
calibration of the detector.

To recover matrices $\Pi_{n}$ and $\Upsilon_{n}$ that are meaningful and
satisfy all the physical constraints of a measurement, we use maximum
likelihood estimation (MLE) [48, 55, 56] in a two-step strategy. First, we
reconstruct the POVMs by minimizing the log-likelihood function
$f(\\{\Pi_{j}\\})=\sum_{n,k}\hat{p}(n|k)\log[\operatorname{Tr}(\Pi_{n}\rho_{k})]$,
which compares the experimental probabilities $\hat{p}(n|k)$ to the set of
feasible matrices $\\{\Pi_{n}\\}$ satisfying $\Pi_{n}\geq 0$ and
$\sum_{n}\Pi_{n}=\openone$. Finally, we estimate the Choi matrices
$\Upsilon_{n}$, minimizing a function
$f_{n}(\Upsilon_{n})=\sum_{mjk}\hat{p}(mn|jk)\log\operatorname{Tr}[(U^{\dagger}_{j}\Pi_{m}U_{j}\otimes\rho_{k}^{T})\tilde{\Upsilon}_{n}]$
which compares the experimental probabilities $\hat{p}(mn|jk)$ to a
parametrization of the Choi matrix $\Upsilon_{n}$ satisfying
$\tilde{\Upsilon}_{n}\geq 0$
($\tilde{\Upsilon}_{n}^{ijkl}=\Upsilon_{n}^{ikjl}$) and the POVM constraint
$\Pi_{n}=\sum_{ijk}\Upsilon_{n}^{kjki}|i\rangle\langle j|$. In total, QND
detector tomography solves $N+1$ optimization problems: one for POVMs of size
$d^{2}$, and $N$ for Chois of size $d^{4}$ 222We solve each minimization
problem using sequential least squares programming, satisfying the positivity
of operators via the Cholesky decomposition, and the completeness constraints
via Lagrange multipliers..

#### QND measurement quantifiers via tomography.-

We use the reconstructed Choi matrices $\Upsilon_{n}$ to quantify the
performance and QND nature of a measurement. Standard benchmarks for QND
detectors are readout fidelity $F=\sum_{n}p(n|n)/N$ and QND-ness
$Q=\sum_{n}p(nn|n)/N$, defined as the average probability that an observable’s
eigenstate $\ket{n}$ remains unchanged after one or two measurements,
respectively. These quantities are related to tomography via

$\displaystyle F={}$
$\displaystyle\frac{1}{N}\sum_{n}\braket{n}{\Pi_{n}}{n}=\frac{1}{N}\sum_{nj}\Upsilon_{n}^{njnj},$
(3) $\displaystyle Q={}$
$\displaystyle\frac{1}{N}\sum_{n}\braket{nn}{\Upsilon_{n}}{nn}=\frac{1}{N}\sum_{n}\Upsilon_{n}^{nnnn}.$
(4)

The readout fidelity $F$ quantifies how close the measurement is to a
projective one $\Pi_{n}=(\Pi_{n})^{2}$. QND-ness $Q$ and other similar
fidelities [51] quantify the overlap with an ideal measurement, satisfying
$\Upsilon_{n}=(\Upsilon_{n})^{2}$. Both are important measurement properties,
but none of them assess the QND nature of the detector and its back-action on
the observables (2). For instance, a maximum value $Q=1$ characterizes an
ideal measurement, but there are non-ideal measurements $Q\neq 1$ which are
close to QND. Non-ideal QND measurements are useful in quantum tasks that only
require evaluating observable averages such as variational algorithms or
sensing [58, 59, 60, 61, 62, 63].

We introduce the destructiveness $D$ as a precise quantifier of the QND nature
of a detector, regardless of how ideal it is. This new quantifier bounds the
change or back-action suffered by any observable compatible with the
measurement of $O$, in accordance to property (2):

$\displaystyle D={}$
$\displaystyle\frac{1}{2}\max_{||O_{c}||=1}||O_{c}-\mathcal{E}^{\dagger}(O_{c})||,\qquad[O,O_{c}]=0.$
(5)

Evaluating $D$ requires the tomographic reconstruction of the complete
measurement process,
$\mathcal{E}^{\dagger}(O_{c})=\sum_{ijkln}\left[\Upsilon_{n}^{klij}\right]^{\ast}O_{c}^{kl}|i\rangle\langle
j|$, and a maximization over all compatible operators of unit norm
$||O_{c}||=1$, with $||O||=\sqrt{\operatorname{Tr}(O^{\dagger}O)}$. In
practice, this maximization is done by finding the maximum eigenvalue of a
positive matrix [58]. For instance, in the case of the qubit observable
$O=\sigma_{z}$, the destructiveness reduces to
$D=||\sigma_{z}-\mathcal{E}^{\dagger}(\sigma_{z})||/\sqrt{8}$. Since $D$
verifies the general QND condition (2), it also bounds the change of the
probability distribution $p(n)$ over repeated measurements.

As shown below, the three quantities $F$, $Q$, and $D$ describe the most
important aspects of QND detectors, but there is further information to
extract from $\Upsilon_{n}$, which are the most general objects.

#### Calibration of qubit readout beyond dispersive approximation.-

In standard superconducting qubit readout [24], the qubit
$\\{\ket{g},\ket{e}\\}$ couples to an off-resonant cavity mode $a$ with
detuning $\Delta$ and coupling $g$ [cf. Fig. 1(b)]. For a highly anharmonic
qubit, such as the flux qubit, this interaction is described by a Jaynes-
Cummings (JC) Hamiltonian [64],

$\displaystyle
H_{\text{JC}}=\frac{\Delta}{2}\sigma_{z}+g(\sigma_{+}a+a^{\dagger}\sigma_{-})+\Omega_{c}(t)(a+a^{\dagger}),$
(6)

with Pauli operators $\sigma_{z}=\ket{e}\\!\bra{e}-\ket{g}\\!\bra{g}$,
$\sigma_{-}=\sigma_{+}^{\dagger}=\ket{g}\\!\bra{e}$, and a resonant drive
$\Omega_{c}(t)$ on the cavity. In the dispersive limit $\Delta\gg g$,
$H_{\text{JC}}$ approximates a dispersive model
$H_{\text{d}}=\frac{1}{2}(\Delta+\chi)\sigma_{z}+\chi\sigma_{z}a^{\dagger}a+\Omega_{c}(t)(a+a^{\dagger})$
that predicts a qubit-dependent displacement $\chi=g^{2}/\Delta$ on the cavity
resonance. In theory, by continuous homodyne detection
$\braket{a+a^{\dagger}}_{c}$ on the cavity, we can discriminate the qubit
state without destroying it. In practice, the non-dispersive corrections
implicit in $H_{\text{JC}}$ can slightly degrade the QND nature of the
measurement [33, 65].

To realistically quantify the performance and measurement back-action of
dispersive readout, we describe the dynamics with the full $H_{\text{JC}}$ and
a stochastic master equation (SME) [66, 67, 68, 69]. This formalism accounts
for the back-action of the continuous homodyne detection onto the qubit state,
as well as cavity decay $\kappa$, qubit decay $\gamma$, and qubit dephasing
$\gamma_{\phi}$ [58]. We simulate numerically the tomographic procedure,
solving the SME over many realizations of the experiment [70]. On each
trajectory, the qubit is prepared in one of the six states
$\rho_{k}\in\left\\{\ket{g},\ket{e},(\ket{g}\pm\ket{e})/\sqrt{2},(\ket{g}\pm
i\ket{e})/\sqrt{2}\right\\}$. We perform two single-shot measurements,
interleaved by a cavity reset time, and one of the three qubit gates
$U_{j}\in\left\\{\openone,\exp(-i\pi\sigma_{y}/2),\exp(-i\pi\sigma_{x}/2)\right\\}$.
In Fig. 2 we show a representative trajectory of the protocol, where the
outcome of each single-shot measurement is discriminated as
$\langle\sigma_{z}\rangle_{c}=-{\rm sign}(J)$ with
$J=\sqrt{\kappa}\int_{0}^{T}dt\langle a+a^{\dagger}\rangle_{c}$ the homodyne
current integrated over the duration $T$ of the readout pulse $\Omega_{c}(t)$
[58]. For simplicity of the simulation, we neglect imperfections in the qubit
state preparation and gates, performed with a local control $\Omega_{q}(t)$.
In a real experiment, these imperfections can be self-consistently separated
from intrinsic measurement errors by using standard gate set tomography [71,
72, 73]. Simulations consider state-of-the-art parameters of superconducting
qubit readout [74, 75, 26]: $g/2\pi=200{\rm MHz}$, $\kappa=0.2g$,
$\gamma=\gamma_{\phi}=10^{-4}g$, $T=8/\kappa\approx 32{\rm ns}$, and
$|\Omega_{c}|=0.173g$, corresponding to $\langle a^{\dagger}a\rangle\sim 1.5$
photons on cavity for $2\chi/\kappa=1$.

We calibrate the measurement by tuning $\Delta/g$ and computing via tomography
the quantifiers $1-F$, $1-Q$, and $D$ [cf. Fig. 3(a)]. We show predictions
using the realistic $H_{\text{JC}}$ interaction (solid), as well as the
dispersive model $H_{\text{d}}$ (dashed) to benchmark the results. We identify
three qualitatively different points of operation (i)-(iii) as indicated by
vertical lines in Fig. 3(a). For each of them, we display in Figs. 3(b)-(d)
the Choi matrices $|\Upsilon_{n}|$ for both measurement outcomes $n=e,g$,
where blue (orange) columns correspond to the JC (dispersive) predictions and
the upper color corresponds to the higher values.

Figure 2: Simulation of QND detector tomography for dispersive qubit readout.
(a) Pulse scheme on qubit (blue) and cavity (red) to implement state
preparation, gates, and homodyne measurements. (b) Cavity quadrature $\langle
a+a^{\dagger}\rangle_{c}$ conditioned on a single trajectory. (c) Average of
Pauli operators $\langle\sigma_{z}\rangle_{c}$ (blue) and
$\langle\sigma_{x}\rangle_{c}$ (light blue) conditioned on the same
trajectory. This realization corresponds to an initial state
$|+\rangle=(|g\rangle+|e\rangle)/\sqrt{2}$ on the qubit, a first measurement
with outcome $\ket{e}$, a cavity reset time, the use of gate
$\exp(-i\pi\sigma_{y}/2)$, and a second measurement with outcome $\ket{g}$.
Repeating this procedure over many trajectories with different inputs and
gates allows us to reconstruct the Choi matrices $\Upsilon_{g}$ and
$\Upsilon_{e}$. Figure 3: Measurement quantifiers and reconstructed Choi
matrices for QND qubit readout. (a) Readout infidelity $1-F$ (green circles),
QND-ness infidelity $1-Q$ (orange crosses), and Destructiveness $D$ (blue
squares) as a function of $\Delta/g$. Solid (dashed) lines correspond to
predictions for JC (dispersive) models. (b)-(d) Choi matrices $|\Upsilon_{n}|$
for $n=e,g$ measurement outcomes at the three representative values
$\Delta/g=[7.7,19.2,40]$, indicated by vertical lines in (a). Blue (orange)
bars correspond to predictions from JC (dispersive) models. On each column,
the upper color corresponds to the part of the longest bar that does not
overlap with the shortest one (and thus a single color indicates full
overlap). The probabilities to reconstruct $\Upsilon_{n}$ are estimated from
$2\times 10^{4}$ trajectories for each initial state $\rho_{k}$ and gate
$U_{j}$. Error bars correspond to one standard deviation, obtained from
$10^{3}$ bootstrap simulations. Parameters are indicated in text.

(i) Around $\Delta/g=7.7$, the dispersive model reaches its optimum in
fidelity and QND-ness, and the measurement is nearly ideal with projective
Choi matrices $\Upsilon_{n}^{(d)}\approx\ket{nn}\\!\bra{nn}$ [cf. orange
columns in Fig. 3(d)]. This occurs near $2\chi/\kappa\sim 1$ as expected by
theory [24]. In contrast, the Choi matrices of the JC model show strong
deviations from an ideal measurement [cf. blue columns in Fig. 3(b)]. The
population transfer $\ket{e}\to\ket{g}$ during measurement—mainly due to
cavity-induced Purcell decay [33]—reduces $\Upsilon_{e}^{eeee}$ and increases
$\Upsilon_{n}^{ggee}$. In addition, the increase of $\Upsilon_{n}^{genn}$ and
$\Upsilon_{n}^{egnn}$ corresponds to the growth of qubit coherences during
measurement—e.g. due to cavity-mediated qubit driving [65]. These non-
dispersive effects increase the infidelity and destructiveness [cf. Fig. 3(a)]
and shift the optimal working point to larger $\Delta/g$.

(ii) Around $\Delta/g=19.2$, the non-dispersive corrections decrease and the
QND-ness reaches an optimum $Q\approx 0.97$ where the actual measurement is
most ideal. From $\Upsilon_{n}$ in Fig. 3(c), we observe that non-dispersive
effects are present but strongly suppressed. In contrast to the dispersive
model, this optimum of QND-ness does not coincide with the optimum of fidelity
[cf. Fig. 3(a)]. This is explained by noting that QND-ness is influenced by
back-action $D$, while fidelity has no information about post-measurement
states and cannot depend on it. In the JC model, $D$ decreases monotonically
with $\Delta/g$, and thus the condition of highest $Q$ shifts towards larger
$\Delta/g$ in order to minimize back-action 333Refining the discrimination
criteria between measurement outcomes does not change this behavior [58]. In
contrast, the back-action $D$ of the dispersive model is nearly constant with
$\Delta/g$, and the optima of $Q$ and $F$ coincide.

(iii) Around $\Delta/g=40$ the system is deep in the dispersive limit. The
predictions for both models almost coincide, up to residual cavity-mediated
qubit driving [65]. Here, $1-F$ and $1-Q$ get worse, but $D$ reaches its
minimum, meaning that the realistic measurement is maximally QND, but less
ideal than in (ii). The loss of ideality is clearly manifested by the large
diagonal terms $\Upsilon_{n}^{mmll}$ [cf. Fig. 3(d)]. This behavior is
explained theoretically [58] by the low distinguishability between outcomes
$n=g,e$ when the cavity displacement $\sim g^{2}/\Delta$ is too small compared
to the measurement uncertainty $\sim\kappa$. The minimum value of $D$ is
attributed to how decoherence $\gamma=\gamma_{\phi}=10^{-4}g$ breaks the QND
condition (2) even when non-dispersive effects are suppressed [58].

We see that a tomography-based calibration allows us to characterize different
regimes of the detector, and to identify error sources via the process
matrices $\Upsilon_{n}$. When simulating the measurement dynamics with
$H_{\text{JC}}$, we account for all non-dispersive effects and the back-action
appearing in realistic superconducting circuit experiments, in a unified way.
We also consider imperfections due to larger intrinsic qubit decay and
dephasing on the operation points (ii) and (iii). The effect of intrinsic
decay $\gamma$ is qualitatively similar to Purcell decay, whereas pure
dephasing $\gamma_{\phi}$ has a negligible effect on the Choi matrices for
long readout pulses $T\gg 1/\kappa$ [58]. Finally, our physical analysis can
be used to choose an optimal working point for the detector. If one is
interested in a near-ideal QND measurement, optimizing $Q$ gives a good
compromise between fidelity $F$ and back-action $D$, as illustrated in (ii).
This requires a larger detuning $\Delta/g$ than expected by the standard
dispersive prediction [24]. However, if one is interested in a maximally QND
measurement with minimal back-action, $D$ is a more suitable quantity to
optimize. This may require loosing ideality of the measurement as shown in
(iii).

#### Conclusions and Outlook.-

We developed a tomographic procedure to calibrate and characterize arbitrary
QND detectors via a reconstruction of the measurement processes
$\Upsilon_{n}$. We applied the method to a realistic simulation of
superconducting qubit readout and identified important discrepancies between
the JC and dispersive models. We expect even larger non-dispersive effects
when taking into account the multi-level character of low anharmonic qubits
such as transmons [35]. This is because the effective dispersive shift is
reduced [58], and one requires lower detunings and stronger cavity drives to
reach an optimal readout. While the tomography will not change, an optimized
simulation of QND readout for multi-level qubits is numerically more
challenging [58] and lies outside the scope of this work.

Experimentally, the presented tomography requires only standard control and
measurement tools, and therefore it can be immediately implemented to
systematically analyze relevant effects on qubit readout such as the strong
driving regime [77, 35], leakage to higher levels [78, 79], or cross-talk [80,
81]. This understanding may help improve the QND measurement performance and
guide the design of alternative schemes [36, 37, 38, 39, 40, 41, 42, 43, 44].
Moreover, the method can be directly applied to other platforms such as QND
detectors of microwave [82, 83, 84, 85, 86, 87] or optical photons [88, 89,
90, 91]. State preparation and gate errors occurring in experiments can be
self-consistently included in the protocol by incorporating standard gate set
tomography [71, 72, 73]. Furthermore, the reconstruction of high-dimensional
Choi matrices could be done more efficiently using compressed sensing [92, 93,
94, 95], matrix-product-states [96, 97], or other advanced techniques [98, 99,
100, 101, 102, 103, 104].

From a fundamental point of view, our technique introduces an accurate
procedure to quantify the back-action and real QND nature of a measurement via
the destructiveness $D$. This complements the standard analysis in terms of
readout fidelity and QND-ness, and allows us to identify regimes of minimal
back-action regardless the ideality of the measurement. Non-ideal but highly
QND measurements may be also exploited for quantum information tasks that
require precise evaluations of expectation values [81] since the measurement
outcomes can be corrected via error mitigation strategies [54, 105, 106].

###### Acknowledgements.

This work has been supported by funding from Spanish project
PGC2018-094792-B-I00 (MCIU/AEI/FEDER, UE) and CAM/FEDER Project No.
S2018/TCS-4342 (QUITEMAD-CM). L.P. was supported by ANID-PFCHA/DOCTORADO-
BECAS-CHILE/2019-772200275. T.R. further acknowledges support from the EU
Horizon 2020 program under the Marie Skłodowska-Curie grant agreement No.
798397, and from the Juan de la Cierva fellowship IJC2019-040260-I.

## References

  * Braginsky and Khalili [1992] V. B. Braginsky and F. Y. Khalili,  _Quantum measurement_ (Cambridge University Press, 1992).
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Colangelo  _et al._ [2017] G. Colangelo, F. M. Ciurana, L. C. Bianchet, R. J. Sewell, and M. W. Mitchell, Nature 543, 525 (2017).
  * Møller  _et al._ [2017] C. B. Møller, R. A. Thomas, G. Vasilakis, E. Zeuthen, Y. Tsaturyan, M. Balabas, K. Jensen, A. Schliesser, K. Hammerer, and E. S. Polzik, Nature 547, 191 (2017).
  * Danilishin and Khalili [2012] S. L. Danilishin and F. Y. Khalili, Living Rev. Relativ. 15, 5 (2012).
  * Kimble  _et al._ [2001] H. J. Kimble, Y. Levin, A. B. Matsko, K. S. Thorne, and S. P. Vyatchanin, Phys. Rev. D 65, 022002 (2001).
  * Rossi  _et al._ [2020] M. A. C. Rossi, F. Albarelli, D. Tamascelli, and M. G. Genoni, Phys. Rev. Lett. 125, 200505 (2020).
  * Campbell  _et al._ [2017] E. T. Campbell, B. M. Terhal, and C. Vuillot, Nature 549, 172 (2017).
  * Fowler  _et al._ [2012] A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Phys. Rev. A 86, 032324 (2012).
  * Bermudez  _et al._ [2019] A. Bermudez, X. Xu, M. Gutiérrez, S. C. Benjamin, and M. Müller, Phys. Rev. A 100, 062307 (2019).
  * Chamberland  _et al._ [2020] C. Chamberland, K. Noh, P. Arrangoiz-Arriola, E. T. Campbell, C. T. Hann, J. Iverson, H. Putterman, T. C. Bohdanowicz, S. T. Flammia, A. Keller, G. Refael, J. Preskill, L. Jiang, A. H. Safavi-Naeini, O. Painter, and F. G. S. L. Brandão, arXiv:2012.04108 (2020).
  * Chen and et al. [Google Quantum AI] Z. Chen and et al. (Google Quantum AI), arXiv:2102.06132 (2021).
  * Leibfried  _et al._ [2003] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Rev. Mod. Phys. 75, 281 (2003).
  * Raha  _et al._ [2020] M. Raha, S. Chen, C. M. Phenicie, S. Ourari, A. M. Dibos, and J. D. Thompson, Nature Comm. 11, 1605 (2020).
  * Volz  _et al._ [2011] J. Volz, R. Gehr, G. Dubois, J. Estève, and J. Reichel, Nature 475, 210 (2011).
  * Gleyzes  _et al._ [2007] S. Gleyzes, S. Kuhr, C. Guerlin, J. Bernu, S. Deléglise, U. Busk Hoff, M. Brune, J.-M. Raimond, and S. Haroche, Nature 446, 297 (2007).
  * Grangier  _et al._ [1998] P. Grangier, J. A. Levenson, and J.-P. Poizat, Nature 396, 537 (1998).
  * Neumann  _et al._ [2010] P. Neumann, J. Beck, M. Steiner, F. Rempp, H. Fedder, P. R. Hemmer, J. Wrachtrup, and F. Jelezko, Science 329, 542 (2010).
  * Robledo  _et al._ [2011] L. Robledo, L. Childress, H. Bernien, B. Hensen, P. F. A. Alkemade, and R. Hanson, Nature 477, 574 (2011).
  * Nakajima  _et al._ [2019] T. Nakajima, A. Noiri, J. Yoneda, M. R. Delbecq, P. Stano, T. Otsuka, K. Takeda, S. Amaha, G. Allison, K. Kawasaki, A. Ludwig, A. D. Wieck, D. Loss, and S. Tarucha, Nature Nanotechnology 14, 555 (2019).
  * Xue  _et al._ [2020] X. Xue, B. D’Anjou, T. F. Watson, D. R. Ward, D. E. Savage, M. G. Lagally, M. Friesen, S. N. Coppersmith, M. A. Eriksson, W. A. Coish, and L. M. K. Vandersypen, Phys. Rev. X 10, 021006 (2020).
  * Wallraff  _et al._ [2005] A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. Lett. 95, 060501 (2005).
  * Walter  _et al._ [2017a] T. Walter, P. Kurpiers, S. Gasparinetti, P. Magnard, A. Potočnik, Y. Salathé, M. Pechal, M. Mondal, M. Oppliger, C. Eichler, and A. Wallraff, Physical Review Applied 7, 054020 (2017a).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * Vijay  _et al._ [2011] R. Vijay, D. H. Slichter, and I. Siddiqi, Phys. Rev. Lett. 106, 110502 (2011).
  * Walter  _et al._ [2017b] T. Walter, P. Kurpiers, S. Gasparinetti, P. Magnard, A. Potočnik, Y. Salathé, M. Pechal, M. Mondal, M. Oppliger, C. Eichler, and A. Wallraff, Phys. Rev. Applied 7, 054020 (2017b).
  * Reed  _et al._ [2010] M. D. Reed, B. R. Johnson, A. A. Houck, L. DiCarlo, J. M. Chow, D. I. Schuster, L. Frunzio, and R. J. Schoelkopf, Appl. Phys. Lett. 96, 203110 (2010).
  * Jeffrey  _et al._ [2014] E. Jeffrey, D. Sank, J. Y. Mutus, T. C. White, J. Kelly, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Megrant, P. J. J. O’Malley, C. Neill, P. Roushan, A. Vainsencher, J. Wenner, A. N. Cleland, and J. M. Martinis, Phys. Rev. Lett. 112, 190504 (2014).
  * Sete  _et al._ [2015] E. A. Sete, J. M. Martinis, and A. N. Korotkov, Phys. Rev. A 92, 012325 (2015).
  * Heinsoo  _et al._ [2018] J. Heinsoo, C. K. Andersen, A. Remm, S. Krinner, T. Walter, Y. Salathé, S. Gasparinetti, J.-C. Besse, A. Potočnik, A. Wallraff, and C. Eichler, Phys. Rev. Applied 10, 034040 (2018).
  * McClure  _et al._ [2016] D. T. McClure, H. Paik, L. S. Bishop, M. Steffen, J. M. Chow, and J. M. Gambetta, Phys. Rev. Applied 5, 011001(R) (2016).
  * Jerger  _et al._ [2012] M. Jerger, S. Poletto, P. Macha, U. Hübner, E. Il’ichev, and A. V. Ustinov, Appl. Phys. Lett. 101, 042604 (2012).
  * Boissonneault  _et al._ [2009] M. Boissonneault, J. M. Gambetta, and A. Blais, Phys. Rev. A 79, 013819 (2009).
  * Slichter  _et al._ [2012] D. H. Slichter, R. Vijay, S. J. Weber, S. Boutin, M. Boissonneault, J. M. Gambetta, A. Blais, and I. Siddiqi, Phys. Rev. Lett. 109, 153601 (2012).
  * Sank  _et al._ [2016] D. Sank, Z. Chen, M. Khezri, J. Kelly, R. Barends, B. Campbell, Y. Chen, B. Chiaro, A. Dunsworth, A. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, T. White, J. Wenner, A. N. Korotkov, and J. M. Martinis, Phys. Rev. Lett. 117, 190503 (2016).
  * Opremcak  _et al._ [2018] A. Opremcak, I. V. Pechenezhskiy, C. Howington, B. G. Christensen, M. A. Beck, E. Leonard, J. Suttle, C. Wilen, K. N. Nesterov, G. J. Ribeill, T. Thorbeck, F. Schlenker, M. G. Vavilov, B. L. T. Plourde, and R. McDermott, Science 361, 1239 (2018).
  * Govia  _et al._ [2014] L. C. G. Govia, E. J. Pritchett, C. Xu, B. L. T. Plourde, M. G. Vavilov, F. K. Wilhelm, and R. McDermott, Phys. Rev. A 90, 062307 (2014).
  * Siddiqi  _et al._ [2006] I. Siddiqi, R. Vijay, M. Metcalfe, E. Boaknin, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, Phys. Rev. B 73, 054510 (2006).
  * Krantz  _et al._ [2016] P. Krantz, A. Bengtsson, M. Simoen, S. Gustavsson, V. Shumeiko, W. D. Oliver, C. M. Wilson, P. Delsing, and J. Bylander, Nature Communications 7, 11417 (2016).
  * Dassonneville  _et al._ [2020a] R. Dassonneville, T. Ramos, V. Milchakov, L. Planat, E. Dumur, F. Foroughi, J. Puertas, S. Leger, K. Bharadwaj, J. Delaforce, C. Naud, W. Hasch-Guichard, J. J. García-Ripoll, N. Roch, and O. Buisson, Phys. Rev. X 10, 011045 (2020a).
  * Wang  _et al._ [2019] X. Wang, A. Miranowicz, and F. Nori, Phys. Rev. Applied 12, 064037 (2019).
  * Didier  _et al._ [2015] N. Didier, J. Bourassa, and A. Blais, Phys. Rev. Lett. 115, 203601 (2015).
  * Touzard  _et al._ [2019] S. Touzard, A. Kou, N. E. Frattini, V. V. Sivak, S. Puri, A. Grimm, L. Frunzio, S. Shankar, and M. H. Devoret, Phys. Rev. Lett. 122, 080502 (2019).
  * Noh  _et al._ [2021] T. Noh, Z. Xiao, K. Cicak, X. Y. Jin, E. Doucet, J. Teufel, J. Aumentado, L. C. G. Govia, L. Ranzani, A. Kamal, and R. W. Simmonds, arXiv:2103.09277 (2021).
  * Lundeen  _et al._ [2009] J. S. Lundeen, A. Feito, H. Coldenstrodt-Ronge, K. L. Pregnell, C. Silberhorn, T. C. Ralph, J. Eisert, M. B. Plenio, and I. A. Walmsley, Nature Physics 5, 27 (2009).
  * Luis and Sánchez-Soto [1999] A. Luis and L. L. Sánchez-Soto, Phys. Rev. Lett. 83, 3573 (1999).
  * D’Ariano  _et al._ [2004] G. M. D’Ariano, L. Maccone, and P. Lo Presti, Phys. Rev. Lett. 93, 250407 (2004).
  * Fiurášek [2001] J. Fiurášek, Phys. Rev. A 64, 024102 (2001).
  * Chen  _et al._ [2019] Y. Chen, M. Farahzad, S. Yoo, and T.-C. Wei, Phys. Rev. A 100, 052315 (2019).
  * Rudinger  _et al._ [2021] K. Rudinger, G. J. Ribeill, L. C. G. Govia, M. Ware, E. Nielsen, K. Young, T. A. Ohki, R. Blume-Kohout, and T. Proctor, arXiv:2103.03008 (2021).
  * Blumoff  _et al._ [2016] J. Z. Blumoff, K. Chou, C. Shen, M. Reagor, C. Axline, R. T. Brierley, M. P. Silveri, C. Wang, B. Vlastakis, S. E. Nigg, L. Frunzio, M. H. Devoret, L. Jiang, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. X 6, 031041 (2016).
  * Note [1] In this work, the Choi matrix $\Upsilon_{n}^{ijkl}$ encodes the transformation of the density matrix [107], while $\tilde{\Upsilon}_{n}^{ijkl}=\Upsilon_{n}^{ikjl}$ is a transposed version that is Hermitian and non-negative.
  * Vourdas [1990] A. Vourdas, Optics Communications 76, 164 (1990).
  * Ban [1998] M. Ban, Phys. Lett. A 249, 167 (1998).
  * James  _et al._ [2001] D. F. V. James, P. G. Kwiat, W. J. Munro, and A. G. White, Phys. Rev. A 64, 052312 (2001).
  * Shang  _et al._ [2017] J. Shang, Z. Zhang, and H. K. Ng, Phys. Rev. A 95, 062336 (2017).
  * Note [2] We solve each minimization problem using sequential least squares programming, satisfying the positivity of operators via the Cholesky decomposition, and the completeness constraints via Lagrange multipliers.
  * [58] Supplementary material for details on (i) the calculation of destructiveness and applications of non-ideal qnd measurements, (ii) the simulation of qubit readout with the stochastic master equation, (iii) the explicit form of choi matrices in the dispersive model, (iv) the effect of qubit decoherence on choi matrices, (v) the effect of outcome discrimination criteria on choi matrices, and (vi) a discussion on the multi-level character of qubits in qnd readout.
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, Nature Communications 5, 4213 (2014).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Havlíček  _et al._ [2019] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Nature 567, 209 (2019).
  * Thew  _et al._ [2002] R. T. Thew, K. Nemoto, A. G. White, and W. J. Munro, Phys. Rev. A 66, 012303 (2002).
  * Degen  _et al._ [2017] C. Degen, F. Reinhard, and P. Cappellaro, Reviews of Modern Physics 89, 035002 (2017).
  * Blais  _et al._ [2004] A. Blais, R. S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Phys. Rev. A 69, 062320 (2004).
  * Govia and Wilhelm [2016] L. C. G. Govia and F. K. Wilhelm, Phys. Rev. A 93, 012316 (2016).
  * Wiseman and Milburn [2010] H. M. Wiseman and G. J. Milburn,  _Quantum Measurement and Control_ (Cambridge University Press, 2010).
  * Gambetta  _et al._ [2008] J. Gambetta, A. Blais, M. Boissonneault, A. A. Houck, D. I. Schuster, and S. M. Girvin, Phys. Rev. A 77, 012112 (2008).
  * Laflamme  _et al._ [2017] C. Laflamme, D. Yang, and P. Zoller, Phys. Rev. A 95, 043843 (2017).
  * Yang  _et al._ [2018] D. Yang, C. Laflamme, D. V. Vasilyev, M. A. Baranov, and P. Zoller, Phys. Rev. Lett. 120, 133601 (2018).
  * Jacobs [2010] K. Jacobs,  _Stochastic Processes for Physicists: Understanding Noisy Systems_ (Cambridge University Press, 2010).
  * Merkel  _et al._ [2013] S. T. Merkel, J. M. Gambetta, J. A. Smolin, S. Poletto, A. D. Córcoles, B. R. Johnson, C. A. Ryan, and M. Steffen, Phys. Rev. A 87, 062119 (2013).
  * Greenbaum [2015] D. Greenbaum, arXiv:1509.02921 (2015).
  * Dehollain  _et al._ [2016] J. P. Dehollain, J. T. Muhonen, R. Blume-Kohout, K. M. Rudinger, J. K. Gamble, E. Nielsen, A. Laucht, S. Simmons, R. Kalra, A. S. Dzurak, and A. Morello, New J. Phys. 18, 103018 (2016).
  * Yamamoto  _et al._ [2014] T. Yamamoto, K. Inomata, K. Koshino, P.-M. Billangeon, Y. Nakamura, and J. S. Tsai, New Journal of Physics 16, 015017 (2014).
  * Yan  _et al._ [2016] F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A. P. Sears, D. Hover, T. J. Gudmundsen, D. Rosenberg, G. Samach, S. Weber, J. L. Yoder, T. P. Orlando, J. Clarke, A. J. Kerman, and W. D. Oliver, Nature Communications 7, 12964 (2016).
  * Note [3] Refining the discrimination criteria between measurement outcomes does not change this behavior [58].
  * Boissonneault  _et al._ [2010] M. Boissonneault, J. M. Gambetta, and A. Blais, Phys. Rev. Lett. 105, 100504 (2010).
  * Motzoi  _et al._ [2009] F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm, Phys. Rev. Lett. 103, 110501 (2009).
  * Wang  _et al._ [2021] C. Wang, M.-C. Chen, C.-Y. Lu, and J.-W. Pan, Fundamental Research 1, 16 (2021).
  * Hamilton  _et al._ [2020] K. E. Hamilton, T. Kharazi, T. Morris, A. J. McCaskey, R. S. Bennink, and R. C. Pooser, arXiv:2006.01805 (2020).
  * Bravyi  _et al._ [2021] S. Bravyi, S. Sheldon, A. Kandala, D. C. Mckay, and J. M. Gambetta, Phys. Rev. A 103, 042605 (2021).
  * Kono  _et al._ [2018] S. Kono, K. Koshino, Y. Tabuchi, A. Noguchi, and Y. Nakamura, Nature Phys 14, 546 (2018).
  * Besse  _et al._ [2018] J.-C. Besse, S. Gasparinetti, M. C. Collodo, T. Walter, P. Kurpiers, M. Pechal, C. Eichler, and A. Wallraff, Phys. Rev. X 8, 021003 (2018).
  * Lescanne  _et al._ [2020] R. Lescanne, S. Deléglise, E. Albertinale, U. Réglade, T. Capelle, E. Ivanov, T. Jacqmin, Z. Leghtas, and E. Flurin, Phys. Rev. X 10, 021038 (2020).
  * Dassonneville  _et al._ [2020b] R. Dassonneville, R. Assouly, T. Peronnin, P. Rouchon, and B. Huard, Phys. Rev. Applied 14, 044022 (2020b).
  * Essig  _et al._ [2021] A. Essig, Q. Ficheux, T. Peronnin, N. Cottet, R. Lescanne, A. Sarlette, P. Rouchon, Z. Leghtas, and B. Huard, arXiv:2001.03217 (2021).
  * Curtis  _et al._ [2021] J. C. Curtis, C. T. Hann, S. S. Elder, C. S. Wang, L. Frunzio, L. Jiang, and R. J. Schoelkopf, Phys. Rev. A 103, 023705 (2021).
  * Lvovsky and Raymer [2009] A. I. Lvovsky and M. G. Raymer, Rev. Mod. Phys. 81, 299 (2009).
  * Lobino  _et al._ [2008] M. Lobino, D. Korystov, C. Kupchak, E. Figueroa, B. C. Sanders, and A. I. Lvovsky, Science 322, 563 (2008).
  * Zhang  _et al._ [2012] L. Zhang, H. B. Coldenstrodt-Ronge, A. Datta, G. Puentes, J. S. Lundeen, X.-M. Jin, B. J. Smith, M. B. Plenio, and I. A. Walmsley, Nat. Photonics 6, 364 (2012).
  * Ramos and García-Ripoll [2017] T. Ramos and J. J. García-Ripoll, Phys. Rev. Lett. 119, 153601 (2017).
  * Gross  _et al._ [2010] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert, Phys. Rev. Lett. 105, 150401 (2010).
  * Riofrío  _et al._ [2017] C. A. Riofrío, D. Gross, S. T. Flammia, T. Monz, D. Nigg, R. Blatt, and J. Eisert, Nature Communications 8, 15305 (2017).
  * Ahn  _et al._ [2019a] D. Ahn, Y. S. Teo, H. Jeong, F. Bouchard, F. Hufnagel, E. Karimi, D. Koutný, J. Řeháček, Z. Hradil, G. Leuchs, and L. L. Sánchez-Soto, Phys. Rev. Lett. 122, 100404 (2019a).
  * Ahn  _et al._ [2019b] D. Ahn, Y. S. Teo, H. Jeong, D. Koutný, J. Řeháček, Z. Hradil, G. Leuchs, and L. L. Sánchez-Soto, Phys. Rev. A 100, 012346 (2019b).
  * Cramer  _et al._ [2010] M. Cramer, M. B. Plenio, S. T. Flammia, R. Somma, D. Gross, S. D. Bartlett, O. Landon-Cardinal, D. Poulin, and Y.-K. Liu, Nature Communications 1, 149 (2010).
  * Lanyon  _et al._ [2017] B. P. Lanyon, C. Maier, M. Holzäpfel, T. Baumgratz, C. Hempel, P. Jurcevic, I. Dhand, A. S. Buyskikh, A. J. Daley, M. Cramer, M. B. Plenio, R. Blatt, and C. F. Roos, Nature Physics 13, 1158 (2017).
  * Tóth  _et al._ [2010] G. Tóth, W. Wieczorek, D. Gross, R. Krischek, C. Schwemmer, and H. Weinfurter, Phys. Rev. Lett. 105, 250403 (2010).
  * Goyeneche  _et al._ [2015] D. Goyeneche, G. Cañas, S. Etcheverry, E. S. Gómez, G. B. Xavier, G. Lima, and A. Delgado, Phys. Rev. Lett. 115, 090401 (2015).
  * Zambrano  _et al._ [2019] L. Zambrano, L. Pereira, and A. Delgado, Phys. Rev. A 100, 022340 (2019).
  * Carmeli  _et al._ [2016] C. Carmeli, T. Heinosaari, M. Kech, J. Schultz, and A. Toigo, EPL (Europhysics Letters) 115, 30001 (2016).
  * Zambrano  _et al._ [2020] L. Zambrano, L. Pereira, D. Martínez, G. Cañas, G. Lima, and A. Delgado, Phys. Rev. Applied 14, 064004 (2020).
  * Pereira  _et al._ [2021] L. Pereira, L. Zambrano, and A. Delgado, Scalable estimation of pure multi-qubit states (2021), arXiv:2107.05691 .
  * Farooq  _et al._ [2021] A. Farooq, J. ur Rehman, and H. Shin, Matrix-completion quantum state tomography (2021), arXiv:2111.11071 .
  * Botelho  _et al._ [2021] L. Botelho, A. Glos, A. Kundu, J. A. Miszczak, O. Salehi, and Z. Zimborás, arXiv:2108.10927 (2021).
  * Geller [2021] M. R. Geller, Phys. Rev. Lett. 127, 090502 (2021).
  * Milz  _et al._ [2017] S. Milz, F. A. Pollock, and K. Modi, Open Syst. Inf. Dyn. 24, 1740016 (2017).
  * Note [4] Notice that the calculations in the main text using the two-level model do not depend on the sign of $\Delta$, but the inclusion of the $|f\rangle$ state breaks this symmetry.

Supplemental Material for:
Complete physical characterization of QND measurements via tomography

Luciano Pereira1, Juan José García-Ripoll1, and Tomás Ramos1

1 Instituto de Física Fundamental IFF-CSIC, Calle Serrano 113b, Madrid 28006,
Spain

## Contents

  * •

I.- Additional properties of Destructiveness

    * –

I.A.- Practical calculation of Destructiveness

    * –

I.B.- Destructiveness in the low back-action limit of dispersive readout

    * –

I.C.- Applications of non-ideal QND measurements

  * •

II.- Stochastic master equation for qubit readout via homodyne detection on
cavity

  * •

III.- Choi matrices in the dispersive model

  * •

IV.- Effect of qubit decoherence on Choi matrices

  * •

V.- Effect of outcome discrimination criteria on Choi matrices

  * •

VI.- Multi-level character of qubits in QND readout

##  I I.- Additional properties of Destructiveness

In this section, we give a practical recipe to calculate the Destructiveness
$D$ via the maximum eigenvalue of a specific positive semi-definite matrix
[cf. Sec. I.A], a derivation of a close formula for $D$ in the case of a qubit
[cf. Sec. I.A], and a numerical check that $D$ vanishes for qubit readout deep
in the dispersive limit and for zero qubit decoherence [cf. Sec. I.B].

###  I.1 I.A.- Practical calculation of Destructiveness

We define the Destructiveness as a precise measure of the back-action and QND
nature of a detector, given by Eq. (5) of the main text as,

$\displaystyle
D=\frac{1}{2}\max_{||O_{c}||=1}||O_{c}-\mathcal{E}^{\dagger}(O_{c})||,\qquad[O,O_{c}]=0.$
(7)

This involves a maximization over the set of all normalized and compatible
operators with $O$, and here we provide a simple recipe to calculate this
quantity. Let us consider the eigenvector decomposition of
$O_{c}=\sum_{j}o_{j}P_{j}$, with $P_{j}^{2}=P_{j}$ and $o_{j}\in\mathbb{R}$.
Using this decomposition, the destructiveness can be rewritten as

$\displaystyle
D^{2}=\frac{1}{4}\max_{||\boldsymbol{o}||=1}\sum_{jk}o_{j}B_{jk}o_{k},$ (8)

with
$B_{jk}=\operatorname{Tr}([P_{j}-\mathcal{E}^{\dagger}(P_{j})][P_{k}-\mathcal{E}^{\dagger}(P_{k})])$.
From Eq. (8), we see that $D$ can also be calculated as the largest eigenvalue
of the positive semi-definite matrix $B_{jk}$.

In addition, it is possible to find a close expression of $D$ in the qubit
case ${\cal O}=\sigma_{z}$, which reads

$\displaystyle D=||\sigma_{z}-\mathcal{E}^{\dagger}(\sigma_{z})||/\sqrt{8}.$
(9)

This expression is obtaining by noticing that the matrix $B_{jk}$ has at least
one null eigenvalue with eigenvector $\boldsymbol{o}=(1/\sqrt{2},1/\sqrt{2})$
since the identity $O_{c}=I$ is compatible with any observable and this does
not change with any quantum process. Using this observation in the qubit case,
we find that the no-null eigenvalue of $B_{jk}$ is simply $Tr(B)$, and that
its eigenvector is $\boldsymbol{o}=(1/\sqrt{2},-1/\sqrt{2})$. This eigenvector
defines the compatible observable $O_{c}=\sigma_{z}/\sqrt{2}$, which solves
the optimization problem (8). Evaluating, we obtain Eq. (9).

Figure 4: (a) Destructiveness $D$ obtained from QND tomography of dispersive
readout, for different values of qubit decay $\gamma$. From top to bottom,
$\gamma=10^{-2}$ (blue), $5\times 10^{-3}$ (orange), $10^{-3}$ (green),
$5\times 10^{-4}$ (red), $10^{-4}$ (purple). (b) Average Destructiveness for
each value of decay. The solid orange line is a linear fit,
$D=\alpha\gamma+\beta$, with $\alpha=27.7/g$ and $\beta=8.5\times 10^{-3}$.

###  I.2 I.B.- Destructiveness in the low back-action limit of dispersive
readout

In Fig. 3(a) of the main text we showed that the destructiveness $D$ does not
vanish deep in the dispersive limit of qubit readout, $\Delta\gg g$, despite
non-dispersive effects are strongly suppressed. Here, we show that this
behavior comes from the presence of qubit decoherence in our realistic model,
which breaks the exact QND condition. This occurs even in the dispersive model
$H_{\text{d}}$, which is exempt of non-dispersive effects.

To do so, we numerically compute $D$ by simulating the tomography via the
stochastic master equation with dispersive Hamiltonian $H_{\text{d}}$, and for
five values of the qubit decay $\gamma/g=[10^{-4},5\times
10^{-4},10^{-3},5\times 10^{-3},10^{-2}]$. For simplicity, we consider no
qubit dephasing, $\gamma_{\phi}=0$, but the results are qualitatively similar
when varying dephasing. We also consider a weak driving case,
$\Omega_{c}=g/2$, since here we are mainly interested in the effect of qubit
decay on $D$. Other parameters of the simulations are $\kappa=g/5$, and
$T=10/\kappa$, and we vary the detuning $\Delta/g$. Fig. 4(a) shows $D$ as a
function of $\Delta/g$, for the five values of decay $\gamma$. These results
are averaged over $10^{3}$ trajectories for each initial state and
measurement, indicated in the main text, and the error bars correspond to $3$
standard deviations obtained from $10^{3}$ bootstrappings. We can see that $D$
is almost constant in $\Delta/g$ for all the decays studied. This is
consistent with the dispersive model as there are no terms depending on
$\Delta/g$ that may spoil the QND condition. More importantly, we see that $D$
decreases with $\gamma$, showing that $\gamma$ is responsible for at least
part of the finite back-action quantified by $D$. In the dispersive model,
decoherence $\gamma$ is the only quantity that breaks the QND condition, and
therefore, it should be the only quantity that contributes to $D$ as well. To
show this, we plot the average of $D$ as function of $\gamma$ in Fig. 4(b),
and we perform a linear fit $D(\gamma)=\alpha\gamma+\beta$, with
$\alpha=27.7/g$ and $\beta=8.5\times 10^{-3}$. We observe a good agremeent
between data and fit within the statistical error, which suggests that there
exists a linear relationship between qubit decay and Destructiveness.
Extrapolating $D$ to the decoherence-free limit, $\gamma\rightarrow 0$, we
obtain that the destructiveness is lower than the statistical noise $D\leq
O(1/\sqrt{N_{t}})$, where $N_{t}=10^{3}$ is the number of trajectories used
for each input state and gate in the QND detector tomography. Therefore,
within the statistical precision of the tomographic reconstruction, we
conclude the method predicts a vanishing destructiveness $D\rightarrow 0$,
deep in the dispersive limit, and for zero decoherence.

###  I.3 I.C.- Applications of non-ideal QND measurements

The most popular QND measurement is the ideal measurement of an observable
$O$, that is, projections onto the eigenspaces of $O$. They are fundamental to
carry out error correction codes. It is not so well-known that non-ideal QND
measurements are also useful for quantum tasks, and in particular, for
algorithms that require multiple and precise evaluations of expectation values
[81] such as variational algorithms [59, 60], quantum machine learning [61],
quantum tomography [55, 62], or quantum sensing [63]. Given that the QND
measurements do not change the expected value of the observable in consecutive
measurements, it can be estimated several times in order to increase the
sample of the observable, and then, reduce its uncertainty. However, to
retrieve the expectation value of the observable of a non-ideal QND
measurement it is not enough to employ the standard procedure since the POVM
do not agree with the basis of eigenvectors of the observable. Instead, we
have to employ a mitigation techniques [53, 54]. In the following section, we
explain this procedure and show a real example of using a non-ideal QND
measurement.

Let us consider a generalized measurement $\mathcal{E}_{n}$. Each of the
measurement processes can be represented by a set of Kraus operators $K_{nm}$,
such as $\mathcal{E}_{n}(\rho)=\sum_{m}K_{nm}\rho K_{nm}^{\dagger}$. In this
representation we have access to the POVM elements as
$\Pi_{n}=\sum_{m}K_{nm}^{\dagger}K_{nm}$ and the Choi matrices as
$\Upsilon_{n}=\sum_{m}K_{nm}\otimes K_{nm}^{*}$. For an ideal measurement both
POVM elements and Krauss operators are simply projectors
$K_{m}=\Pi_{m}=\ket{m}\\!\bra{m}$, with $\ket{m}$ the basis of eigenstates of
$O$. Non-ideal measurements are more general and if they are QND, they satisfy
the condition,

$\operatorname{Tr}\left(O_{c}\sum_{m}K_{nm}\rho
K_{nm}^{\dagger}\right)=\operatorname{Tr}(O_{c}\rho),$ (10)

for any observable $O_{c}$ compatible with $O$. Any measurement with
destructiveness $D=0$ fulfils Eq. (10). In terms of the Krauss operators, a
sufficient condition to satisfy the Eq. (10) is that $[K_{nm},O_{c}]=0$. This
kind of QND measurements also satisfy $[\Pi_{n},O_{c}]=0$, that is, each POVM
element is compatible with the observable $O_{c}$. This property allows us to
retrieve the expected value of $O$ from $\Pi_{n}$ thanks to mitigation
techniques even though it is not an ideal measurement of $O$ [81].

These mitigation techniques work as follows. First, if POVM elements $\Pi_{n}$
are compatible with $O$, we can write them as
$\Pi_{n}=\sum_{m}b_{nm}\ket{m}\\!\bra{m}$. Second, let us denote
$\boldsymbol{p}$ the probability distribution of an ideal measurement of $O$
and $\boldsymbol{\tilde{p}}$ of the non-ideal measurement $\Pi_{n}$, which are
related by the equation $\boldsymbol{\tilde{p}}=B\boldsymbol{p}$, with $B$ a
matrix with elements $b_{nm}$. The ideal measurement results can be obtained
from the non-ideal ones by inverting the matrix,

$\displaystyle\boldsymbol{p}=B^{-1}\boldsymbol{\tilde{p}},$ (11)

allowing us to obtain the expected value of $O$. Notice that this procedure
can also be applied to destructive measurements such as $[\Pi_{n},O]=0$,
allowing us to improve the measurement results in that case as well. However,
in the case of a QND measurements the mitigation procedure also guarantees
that the expectation value $\langle O\rangle$ is preserved in consecutive
measurements, allowing to obtain a better precision. In this sense, the
destructiveness $D$ —introduced in this work—gives an upper bound to the
change of the expected value of any observable $O_{c}$ compatible with $O$,

$\displaystyle|\operatorname{Tr}(O_{c}\mathcal{E}(\rho))-\operatorname{Tr}(O_{c}\rho)|\leq
D,$ (12)

and therefore it can be used to estimate the effectiveness of the mitigation
techniques on any specific non-ideal but highly QND measurement.

A simple example of a non-ideal QND measurement that satisfies the condition
$[K_{nm},O_{c}]=0$ is given by the following set of Krauss operators:
$K_{g0}=\sqrt{1-2\epsilon}\ket{g}\\!\bra{g}$,
$K_{g1}=\sqrt{\epsilon}\sigma_{z}$,
$K_{e0}=\sqrt{1-2\epsilon}\ket{e}\\!\bra{e}$, and
$K_{e1}=\sqrt{\epsilon}\sigma_{z}$, with $\epsilon$ a real parameter. This is
a non-ideal QND measurement for any value $\epsilon<1/2$, but in the special
case of $|\epsilon|\ll 1$, it becomes a near-ideal measurement
$\\{\ket{g}\\!\bra{g},\ket{e}\\!\bra{e}\\}$ affected by phase flip noise,
which is a standard noise model for a realistic detector. The POVM elements of
this general non-ideal measurement are
$\Pi_{g}=\ket{g}\\!\bra{g}+\epsilon\sigma_{z}$ and
$\Pi_{e}=\ket{e}\\!\bra{e}-\epsilon\sigma_{z}$, and its mitigation matrix
reads

$\displaystyle B=\begin{pmatrix}1-\epsilon&\epsilon\\\
\epsilon&1-\epsilon\end{pmatrix}.$ (13)

Therefore, performing the above non-ideal measurements and inverting the
linear system (11), we can determine $\braket{O}$ as if the measurement would
be ideal (for any $\epsilon<1/2$). The same mechanism applies to more complex
non-ideal QND measurements.

##  II II.- Stochastic master equation for qubit readout via homodyne
detection on cavity

In this section, we explain how to model the dynamics of a dispersive qubit
readout within the formalism of a stochastic Master equation (SME) [66, 67,
68].

A single-shot measurement is carried out by applying a resonant pulse
$\Omega_{c}(t)$ on the cavity along a time $T$. The coherent dynamics of the
measurement is given either by the Jaynes-Cummings $H=H_{\text{JC}}$ or the
dispersive $H=H_{\text{d}}$ Hamiltonians, given in the main text. The cavity
is continuously monitored via homodyne detection of transmitted photons, which
introduces a stochastic dynamics on the qubit-cavity system. This dynamics is
describe by the following SME [67, 68, 69]

$\displaystyle d\rho=$ $\displaystyle-i[H,\rho]dt+\kappa\mathcal{D}[a]\rho
dt+\sqrt{\kappa}\mathcal{M}[a]\rho dW$
$\displaystyle+\gamma\mathcal{D}[\sigma_{-}]\rho
dt+\frac{\gamma_{\phi}}{2}\mathcal{D}[\sigma_{z}]\rho dt,$ (14)

where the continuous homodyne measurement is described by the the Wiener
process $dW$ and the super-operator $\mathcal{M}[A]\rho=(A-\langle
A\rangle)\rho+\mathrm{H.c.}$. Markovian decay of the cavity $\kappa$, decay of
the qubit $\gamma$, and dephasing of the qubit $\gamma_{\phi}$ are described
by standard Lindblad super-operators, $\mathcal{D}[A]\rho=A\rho
A^{\dagger}-\\{A^{\dagger}A,\rho\\}/2$. We numerically integrate this equation
using an implicit Euler algorithm [70]. The continuous measurement gives us
access to the homodyne current over a single trajectory,

$\displaystyle J(t)=\sqrt{\kappa}\braket{a+a^{\dagger}}_{c}(t)+\xi(t),$ (15)

where $\braket{a+a^{\dagger}}_{c}$ is the intracavity quadrature conditioned
to the trajectory, and $\xi(t)=dW/dt$ is the vacuum shot-noise
$\braket{\braket{\xi(t)\xi(t^{\prime})}}=\delta(t-t^{\prime})$. The direction
of the qubit-dependent phase shift of the resonator can be determined from the
homodyne current, which allows us to discriminate the state of the qubit
without destroying it. To reduce the noise of the signal, we integrated $J(t)$
over the measurement time $T$, obtaining,
$J=\int_{0}^{T}dtJ(t)=\sqrt{\kappa}\int_{0}^{T}dt\braket{a+a^{\dagger}}_{c}$.
Since the qubit-dependent cavity displacements are in opposite directions, the
measurement outcome can be discriminated from the sign of the integrated
signal $J$ as $\langle\sigma_{z}\rangle_{c}=-{\rm sign}(J)$.

Figure 5: Comparison of three independent methods to estimate the readout
infidelity $1-F$ of dispersive qubit readout. Blue corresponds to the
analytical formula (II), orange to the direct definition via probablities
$F=[p(g|g)+p(e|e)]/2$, and green to the tomographic reconstruct and the use of
Eq. (3) of the main text. The error bars correspond to 5 standard deviations
obtained with $10^{3}$ bootstrap simulations.

In order to benchmark our code that solves Eq. (14), we numerically calculate
the readout fidelity $F$ via two methods. The first method consists in
computing the probabilities $p(g|g)$ and $p(e|e)$ by a direct simulation of
the experiment with a single qubit readout. The second method is applying our
tomographic protocol involving two measurements and then obtaining $F$ via Eq.
(3) of the main text. In addition, we compare these two numerical methods with
the analytical prediction in the case of the dispersive model, given in Ref.
[42] by

$\displaystyle F=1-\frac{1}{2}{\rm erfc}({\rm SNR}/2),$ (16)

where ${\rm erfc}$ is the complementary error function and ${\rm SNR}$ is the
signal-to-noise ratio,

$\displaystyle{\rm SNR}=\frac{4\Omega_{c}\sin(\varphi)T}{\sqrt{2\kappa
T}}\bigg{[}$ $\displaystyle 1-\frac{4}{\kappa T}\Big{(}\cos(\varphi/2)^{2}$
$\displaystyle-\frac{\sin(\chi T+\varphi)}{\sin(\varphi)}e^{-\kappa
T/2}\Big{)}\bigg{]}.$ (17)

Here, $\varphi=2\arctan(2\chi/\kappa)$. For the simulations, we set the
parameters on $g$ units as $\kappa/g=1/5$, $\gamma=\gamma_{\varphi}=0$,
$\Omega_{c}/g=1/10$, and $T=10/\kappa$. The average number of photons in the
resonator is $\braket{a^{\dagger}a}=0.5$. We consider here a weak driving case
since it requires a low number of trajectories to achieve a good precision in
the numerical calculations. For each initial state and gate needed for the
tomography (indicated in the main text), we perform $10^{3}$ trajectories. The
error bars are 5 times the standard deviation obtained with $10^{3}$ bootstrap
simulations. The results are shown in Fig. 5. We see that the three
independent estimations of $F$ —direct, tomography, and analytical— agree well
within an error of 5 standard deviations.

##  III III.- Choi matrices in the dispersive model

In this section, we study the general structure of Choi matrices
$\Upsilon_{n}$ predicted by the dispersive model. With this we show that the
diagonal form of $\Upsilon_{n}$ in Fig. 3(d) of the main text is
characteristic from a low measurement indistinguishably occurring deep in the
dispersive limit.

To do so, let us consider that the qubit is initially in the state
$\ket{\Psi}=\psi_{e}\ket{e}+\psi_{g}\ket{g}$ and the cavity is empty
$|0\rangle_{r}$. Modeling the system with the dispersive Hamiltonian
$H_{\text{d}}$, the coherent evolution of the joint system
$\ket{\Psi}=\ket{\Psi}_{q}|0\rangle_{r}$ after a time $t$ is given by [24]

$\displaystyle\ket{\Psi(t)}=\psi_{g}\ket{g}\ket{\alpha_{g}(t)}+e^{-i\Delta
t}\psi_{e}\ket{e}\ket{\alpha_{e}(t)}.$ (18)

Here, $|\alpha_{j}(t)\rangle$ is the resonator state whose phase is shifted
conditioned on the qubit state $j\in\\{g,e\\}$. Therefore, we can infer the
qubit state by measuring the cavity field. This evolution is described by the
unitary operation,

$\displaystyle U(t)=\ket{g}\\!\bra{g}\otimes D(\alpha_{g}(t))+e^{-i\Delta
t}\ket{e}\\!\bra{e}\otimes D(\alpha_{e}(t)),$ (19)

where $D(\alpha)$ is the cavity displacement operator. Reading the qubit state
requires continuous homodyne detection during a time $T$ of the cavity
quadrature operator $\tilde{Q}=\sqrt{\kappa}\left(a+a^{\dagger}\right)$. The
outcome of the qubit measurement is discriminated form the integrated homodyne
current, $J=\frac{1}{T}\int_{0}^{T}dt\langle\Psi(t)|\tilde{Q}\ket{\Psi(t)}$.
This can be obtained as the expected value of $\tilde{Q}$ in the time average
density matrix,
$\mathcal{E}_{qr}(\ket{\Psi}\\!\bra{\Psi})=\frac{1}{T}\int_{0}^{T}dt\ket{\Psi(t)}\\!\bra{\Psi(t)}$.
In addition, the reduced evolution for the qubit can be obtained by tracing
over the cavity
$\mathcal{E}_{q}(\ket{\Psi}\\!\bra{\Psi})=\operatorname{Tr}_{r}[\mathcal{E}_{qr}(\ket{\Psi}\\!\bra{\Psi})]$.
Considering the eigenvector basis of $\tilde{Q}$, $\\{|q\rangle\\}$, we have
that the Kraus representation of the qubit dynamics is

$\displaystyle\mathcal{E}_{q}(\ket{\Psi}\\!\bra{\Psi})=\int_{0}^{T}dt\int_{-\infty}^{\infty}dqK_{q}(t)\ket{\Psi}\\!\bra{\Psi}K_{q}(t)^{\dagger},$
(20)

with Kraus operators $K_{q}(t)={}_{r}\langle q|U(t)|0\rangle_{r}/\sqrt{T}$
[2]. Considering the wave function of the cavity states $\phi_{j}(q,t)=\langle
q|\alpha_{j}(t)\rangle$, the Kraus operators can be expressed as

$\displaystyle
K_{q}(t)=\frac{1}{\sqrt{T}}(\phi_{g}(q,t)\ket{g}\\!\bra{g}+\phi_{e}(q,t)e^{-i\Delta
t}\ket{e}\\!\bra{e}).$ (21)

We can divide the trajectories in order to determine the outcome of the
measurement. If the homodyne current is larger than a predefined value
$\delta$, we assign the result $e$, and otherwise we assign $g$. For instance,
in the special case of a dispersive model without decoherence both states
displace the cavity by the same amount, but in opposite directions, having
$\delta=0$. Using this separation in the general case, the components of the
Choi matrices are given by

$\displaystyle\Upsilon_{g}^{ijkl}=$
$\displaystyle\int_{0}^{T}dt\int_{-\infty}^{\delta}dq\bra{i}K_{q}(t)\ket{k}\\!\bra{l}K_{q}(t)^{\dagger}\ket{j},$
(22) $\displaystyle\Upsilon_{e}^{ijkl}=$
$\displaystyle\int_{0}^{T}dt\int_{\delta}^{\infty}dq\bra{i}K_{q}(t)\ket{k}\\!\bra{l}K_{q}(t)^{\dagger}\ket{j}.$
(23)

We can further simplify the Choi matrices as,

$\displaystyle\Upsilon_{g}=$
$\displaystyle\epsilon_{g}\ket{gg}\\!\bra{gg}+(1-\epsilon_{e})\ket{ee}\\!\bra{ee}$
$\displaystyle+\zeta_{g}\ket{ge}\\!\bra{ge}+\zeta_{g}^{*}\ket{eg}\\!\bra{eg},$
(24) $\displaystyle\Upsilon_{e}=$
$\displaystyle(1-\epsilon_{g})\ket{gg}\\!\bra{gg}+\epsilon_{e}\ket{ee}\\!\bra{ee}$
$\displaystyle+\zeta_{e}^{*}\ket{ge}\\!\bra{ge}+\zeta_{e}\ket{eg}\\!\bra{ge},$
(25)

where we have defined the four auxiliary quantities
$\epsilon_{g}=\frac{1}{T}\int_{0}^{T}dt\int_{-\infty}^{\delta}dq|\phi_{e}(q,t)|^{2}$,
$\epsilon_{e}=\frac{1}{T}\int_{0}^{T}dt\int_{\delta}^{\infty}dq|\phi_{g}(q,t)|^{2}$,
$\zeta_{g}=\frac{1}{T}\int_{0}^{T}dt\int_{-\infty}^{\delta}dqe^{i\Delta
t}\phi_{g}(q,t)\phi_{e}(q,t)^{*}$, and
$\zeta_{e}=\frac{1}{T}\int_{0}^{T}dt\int_{\delta}^{\infty}dqe^{-i\Delta
t}\phi_{e}(q,t)\phi_{g}(q,t)^{*}$.

The general form of $\Upsilon_{n}$ in Eqs. (III)-(III) deviate from projectors
and thus we conclude that the measurement is not ideal in general. In the case
without decoherence, the wave functions correspond to coherent states. In this
case, the outcomes of the measurement cannot be perfectly discriminated since
the coherent states are not orthogonal. Therefore, we have that $\epsilon<1$
and $|\zeta_{j}|>0$, and that the dispersive readout cannot be a perfectly
ideal measurement. When the cavity displacement $\sim g^{2}/\Delta$ is too
small compared to the measurement uncertainty $\sim\kappa$, the overlap
between the wave functions increase, implying that $\epsilon_{j}$ decrease and
$|\zeta_{j}|$ increase. We can see this effect in the Choi matrices of Fig.
3(d) of the main text, which exhibits large diagonal terms
$\Upsilon_{n}^{gege}$ and $\Upsilon_{n}^{egeg}$.

##  IV IV.- Effect of qubit decoherence on Choi matrices

In this section, we simulate the tomographic characterization of dispersive
readout as a function of qubit decay $\gamma$ and dephasing $\gamma_{\phi}$ to
study its effect on the Choi matrix components $\Upsilon_{g}$ and
$\Upsilon_{e}$.

Figure 6: Choi matrices of dispersive readout with JC model in presence of
decay $\gamma$ for (a) $\Delta/g=19.2$ and (b) $\Delta/g=40$. Columns of
different colors represent different values of decay, which increase in the
direction indicated by the arrow as $\gamma/g=10^{-4}$ (blue), $5\times
10^{-3}$ (green), $10^{-2}$(yellow). Figure 7: Choi matrices of dispersive
readout with JC model in presence of dephasing $\gamma_{\phi}$ for (a)
$\Delta/g=19.2$ and (b) $\Delta/g=40$. Columns of different colors represent
different values of dephasing, which increase in the direction indicated by
the arrow as $\gamma_{\phi}/g=10^{-4}$ (blue), $5\times 10^{-3}$ (green),
$10^{-2}$(yellow).

The simulations of QND detector tomography are done using the JC Hamiltonian,
and with $N_{t}=2\times 10^{3}$ trajectories for each initial state and gate
(indicated in the main text). Parameters are $\kappa/g=1/5$, $T=6/\kappa$, and
$\Omega_{c}/g=\sqrt{3}/10$. For $\Delta/g$, we consider the operating points
(ii) and (iii) of Fig. 3(a) of the main text, that is a near-ideal measurement
$\Delta/g=19.2$, and the deep dispersive limit $\Delta/g=40$, respectively.
For both cases, we compute the Choi matrices $|\Upsilon_{n}|$ for three values
of decay $\gamma/g$ as shown by columns of different colors in Fig. 6.
Similarly, we consider three values of dephasing $\gamma_{\phi}/g$, and the
corresponding $|\Upsilon_{n}|$ are shown in Fig. 7. Both operating points (ii)
and (iii) show similar behavior when increasing qubit and dephasing,
indicating that decoherence manifests in $\Upsilon_{n}$ independently of
$\Delta/g$. From Fig. 6 we see that the main effect of decay $\gamma$ is to
increase the Choi components $\Upsilon_{n}^{ggee}$ and decrease
$\Upsilon_{n}^{eeee}$. This behavior is consistent with the effect of cavity-
induced Purcell decay discussed in the main text. From Fig. 7 we see that the
dephasing has marginal effect on the Choi matrix because all the deviations in
its components are on the order of magnitude or lower than the statistical
error $O(1/\sqrt{N_{t}})$. This is reasonable since dephasing cannot
appreciably affect the dynamics of the qubit after it has been projected to
$\ket{g}$ or $\ket{e}$ by the measurement. Dephasing may be relevant for short
measurement times $T$, when the projection has not yet fully happened, but in
all our simulations we have considered $T\gtrsim 6/\kappa$.

##  V V.- Effect of outcome discrimination criteria on measurement quantifiers

The simulation of QND readout carried out in the main text uses a simple
discrimination strategy based on the sign of the integrated homodyne current
(cf. Sec. II of SM). In this section, we consider two ways of refining the
discrimination method between outcomes and we calculate its effect on the
quantifiers $F$, $Q$, and $D$ obtained from the QND detector tomography. A
first refining comes from noticing that the homodyne current at short times is
not enough to discriminate the outcome of the measurement, and therefore it is
better to employ an integrated homodyne current that ponderates more the
values at longer times instead of using a flat integral such as in Sec. II of
the SM. This can be done by introducing a weighting function $w(t)$ in the
determination of the integrated homodyne current [26] as

$\displaystyle J=\int_{0}^{T}dtw(t)J(t),$ (26)

with $J(t)=\langle a+a^{\dagger}\rangle_{c}$ the time-dependent homodyne
current of a single trajectory and $T$ the integration time. A second simple
refinement comes from noticing that when non-dispersive effects are included,
the displacement of the cavity by the ground and the excited states are not
symmetric. Specifically, the displacement of the excited state is lower than
the ground state. Because of this, one obtains a better discrimination of
outcomes when placing the discrimination threshold on the integrated homodyne
at a finite value $\delta$ instead of simply at zero, as done in Sec. II of
the SM and in the main text. Both the weighting function $w(t)$ and threshold
$\delta$ can be calibrated from homodyne currents $J_{g}(t)$ and $J_{e}(t)$
obtained from solving the exact master equation when the qubit is initially in
$\ket{g}$ or $\ket{e}$, respectively. In particular, they are given by
$w(t)=|J_{g}(t)-J_{e}(t)|$ and $\delta=(J_{g}(T)+J_{e}(T))/2$. In a real
experiment, it is not necessary to perform these simulations to calibrate
$w(t)$ and $\delta$ because they can also be determined from a set of
experimental trajectories [26].

In Fig. 8 we compare the three measurement quantifiers $F$, $Q$, and $D$
obtained from the tomography using either the simple discrimination of the
main text (solid) or the optimized discrimination explained above (dashed).
For the simulations we consider the full JC model with the same parameters of
the main text. Our main conclusion is that refining the discrimination
criteria improves the overall quality of the measurement since the readout and
QND-ness infidelities reduce, but the specific features of the measurement
remain unchanged. In particular, for both discrimination strategies, the
minimum of $1-Q$ occurs at a larger detuning $\Delta/g$ than the minimum of
$1-F$. The effect of improving the discrimination is that both minima move
towards larger detuning $\Delta/g$ but they keep their distance roughly the
same on order $\Delta/g\sim 2$ (the minimum of $1-F$ moves from $\Delta/g\sim
17$ to $\Delta/g\sim 19$, while the minimum of $1-Q$ moves from $\Delta/g\sim
19$ to $\Delta/g\sim 21$). As discussed in the main text, these difference of
the minima of $1-F$ and $1-Q$ has to do with the back-action introduced by the
cavity in the JC model, which reduces at larger $\Delta/g$ as quantified by
$D$. Finally, it is also interesting to notice that the destructiveness $D$ is
not appreciably changed when improving the discrimination criteria, which is
consistent for a quantifier of the physical back-action of a measurement since
this cannot depend on the post-processing strategy used to interpret the data.

Figure 8: Readout infidelity $1-F$ (green circles), QND-ness infidelity $1-Q$
(orange crosses), and Destructiveness $D$ (blue squares) as a function of
$\Delta/g$ for QND qubit readout in the case of two different discrimination
criteria: Solid lines correspond to predictions with the simple discrimination
method used in the main text, while dashed lines correspond to the optimized
discrimination method with nonzero threshold and weighting function. In both
cases, simulations are done with JC model and for the same parameters as in
the main text. Probabilities were estimated from $10^{4}$ trajectories for
each initial state and gate necessary for the tomography. Error bars
correspond to one standard deviation, obtained from $10^{3}$ bootstrap
simulations.

## VI.- Multi-level character of qubits in QND readout

Superconducting qubits are more precisely modeled as multi-level systems with
Hamiltonian $H_{q}=\sum_{k}\omega_{k}\ket{k}\bra{k}$, where $\omega_{k}$ is
the resonance frequency of each of their eigenstates $\ket{k}$
($k=0,1,2,\dots$) [24, 75, 35]. The relevance of the multi-level character of
qubits depends on the size of their anharmonicity
$\alpha_{q}=(\omega_{2}-\omega_{1})-(\omega_{1}-\omega_{0})$ compared to the
other external control fields applied on the qubit as we see below.

To include the multi-level character of qubits in the simulation of QND
readout, one must generalize the JC Hamiltonian of the main text to a multi-
level JC Hamiltonian [35],

$\displaystyle H_{\text{JC}}^{\rm{(ML)}}={}$
$\displaystyle\sum_{k}\omega_{k}\ket{k}\bra{k}+\omega_{c}a^{\dagger}a+\Omega_{c}(ae^{i\omega_{d}t}+a^{\dagger}e^{-i\omega_{d}t})$
$\displaystyle+\sum_{k\neq
k^{\prime}}g_{kk^{\prime}}\ket{k^{\prime}}\bra{k}(a+a^{\dagger}).$ (27)

Here, $\omega_{c}$ is the cavity frequency, $\omega_{d}$ the driving
frequency, $\Omega_{c}(t)$ the cavity driving strength, and $g_{kk^{\prime}}$
the qubit-cavity couplings for the qubit transition
$\ket{k}\rightarrow\ket{k}^{\prime}$ (satisfying
$g_{kk^{\prime}}=g_{k^{\prime}k}$). Notice that Eq. (27) describes qubit-
cavity interaction processes outside the rotating wave approximation (RWA)
which is reported to be relevant when the cavity is highly populated
$\braket{a^{\dagger}a}\gg n_{c}=(g/2\Delta)^{2}$ [35]. The dissipative part of
the dynamics is governed by the same stochastic master equation presented in
Sec. II of this SM, replacing the Lindblad terms of $\sigma_{-}$ and
$\sigma_{z}$ by the various transition operators $\ket{k^{\prime}}\bra{k}$
($k\neq k^{\prime}$) and projectors $\ket{k}\bra{k}$, respectively.

The Hamiltonian (27) describes the complete coherent dynamics of the standard
$\ket{g}=\ket{0}$ and $\ket{e}=\ket{1}$ states of the qubit, but also includes
all other higher excited states $\ket{k}$ with $k\geq 2$. However, if one is
only interested in a QND measurement of the two-level manifold
$\\{\ket{g},\ket{e}\\}$ and in the moderate driving regime
($\braket{a^{\dagger}a}\ll n_{c}$), the most relevant deviation from the naive
two-level approximation comes from the second excited level $\ket{f}=\ket{2}$.
In this case, we can neglect the non-RWA terms and the JC multi-level
Hamiltonian (27) can be approximated by a qutrit Hamiltonian, which in the
rotating frame with the driving frequency $\omega_{d}$ can be written as:

$\displaystyle H_{\text{JC}}^{\rm{(3)}}={}$
$\displaystyle\frac{\Delta_{q}}{2}\sigma_{z}+\frac{3}{2}\left(\Delta_{q}+\alpha_{q}\right)|f\rangle\langle
f|+\Delta_{c}a^{\dagger}a$
$\displaystyle+g(\sigma^{+}a+a^{\dagger}\sigma^{-})+\sqrt{2}g\left(\ket{f}\bra{e}a+a^{\dagger}\ket{e}\bra{f}\right)$
$\displaystyle+\Omega_{c}(t)(a+a^{\dagger}).$ (28)

Here, $\sigma_{-}=\ket{g}\bra{e}$ and
$\sigma_{z}=\ket{e}\bra{e}-\ket{g}\bra{g}$ are the Pauli operators,
$\Delta_{q}=(\omega_{1}-\omega_{0})-\omega_{d}$ is qubit-drive detuning,
$\Delta_{c}=\omega_{c}-\omega_{d}$ the cavity-drive detuning,
$\alpha_{q}=(\omega_{2}-\omega_{1})-(\omega_{1}-\omega_{0})$ the qubit
anharmonicity, and we simplified the notation for the couplings as
$g_{01}=g_{10}=g$ and $g_{12}=g_{21}=\sqrt{2}g$. To gain insight on the main
effects of the extra $\ket{f}$ level on the QND readout scheme, we apply
standard perturbation theory on Eq. (28) [24] and project onto the
$\\{\ket{g},\ket{e}\\}$ manifold, obtaining a modified dispersive Hamiltonian
$H_{\text{D}}^{\prime}$ for an approximated two-level system:

$\displaystyle H_{\text{D}}^{\prime}={}$
$\displaystyle\frac{1}{2}(\Delta+\Delta_{c}+\chi^{\prime})\sigma_{z}+\left(\Delta_{c}-\frac{g^{2}}{\Delta+\alpha_{q}}\right)a^{\dagger}a$
$\displaystyle+\chi^{\prime}\sigma_{z}a^{\dagger}a+\Omega_{c}(t)(a+a^{\dagger}).$
(29)

This two-level model is valid in the dispersive regime
$g\ll|\Delta|,|\Delta+\alpha_{q}|$, with
$\Delta=(\omega_{1}-\omega_{0})-\omega_{c}$ the qubit-cavity detuning. Setting
$\Delta_{c}=g^{2}/(\Delta+\alpha_{q})$ to compensate the cavity shift, Eq.
(29) has exactly the same form as the standard dispersive Hamiltonian used in
the main text, except that here the dispersive shift $\chi^{\prime}$ is
modified by the finite anharmonicity $\alpha_{q}$ of the qubit as [26]

$\displaystyle\chi^{\prime}=g^{2}\left(\frac{1}{\Delta}-\frac{1}{\Delta+\alpha_{q}}\right)=\frac{g^{2}}{\Delta}\frac{\alpha_{q}}{\Delta+\alpha_{q}}.$
(30)

Therefore, for a supeconducting qubit with anharmonicity much larger than the
detuning $|\alpha_{q}|\gg|\Delta|$, we recover the dispersive model of the
main text with $\chi=g^{2}/\Delta$ as the $\ket{f}\leftrightarrow\ket{e}$
transition is much more off-resonant than $\ket{e}\leftrightarrow\ket{g}$.
However, for low anharmonic qubits (such as transmons [26, 35]), one typically
has $|\alpha_{q}|\lesssim|\Delta|$ (with $\alpha_{q}<0$) and therefore the
main effect of the $\ket{f}$ state is to reduce the magnitude of the
dispersive shift by a factor $\alpha_{q}/(\alpha_{q}+\Delta)$ compared to the
naive two-level model. Importantly, one can reduce $\Delta$ to obtain the same
effective shifts as predicted by the naive two-level model. For instance, for
$\alpha_{q}<0$ and $\Delta<0$, if we choose
$\Delta^{\prime}=|\alpha_{q}|/2+\sqrt{|\alpha_{q}|^{2}/4+|\Delta||\alpha_{q}|}$
we get $|\chi^{\prime}(\Delta^{\prime})|=|\chi(\Delta)|$ (the sign of the
shift does not affect the measurement performance). This means that results
for the dispersive model of Fig. 3 of the main text (obtained in the range
$|\Delta|/g\in[3,40]$) can be re-scaled to account for the most important
multi-level effects of the qubit, but at detunings in the shorter range
$\Delta^{\prime}/g\in[3,9]$ 444Notice that the calculations in the main text
using the two-level model do not depend on the sign of $\Delta$, but the
inclusion of the $|f\rangle$ state breaks this symmetry. The above is valid
deep in the dispersive regime ($g\ll|\Delta|,|\Delta+\alpha_{q}|$) and at
moderate drive power ($\braket{a^{\dagger}a}\ll n_{c}$), but studying non-
dispersive effects in the multi-level model is much more complex than for a
two-level system [35].

In practice, when performing QND readout of low anharmonic qubits, we expect
that deviations from the dispersive model $H_{D}^{\prime}$ should become more
important than predicted by the two-level model in Fig. 3 of the main text.
This is because the multi-level character of the qubit reduces the effective
dispersive shift $\chi^{\prime}$ as indicated in Eq. (30), and therefore one
needs to perform the measurement at lower detunings $\Delta$ and larger cavity
drives $\Omega_{c}$ to obtain a similar resolution for discriminating the
outcomes $n=g,e$. For instance, for a transmon qubit with a typical
anharmonicity of $\alpha_{q}/2\pi\sim 300$ MHz, we estimate using Eq. (16)
that reaching an optimal measurement performance such as for point (ii) in
Fig. 3, we require to reduce the detuning $\Delta/g$ by a factor of $3$, as
well as to increase the product $\Omega_{c}T$ by a factor of $2.5$ (with $T$
the measurement time). These two changes in parameters make the dispersive
approximation ($g\ll|\Delta|,|\Delta+\alpha_{q}|,\braket{a^{\dagger}a}\ll
n_{c}$) less valid and therefore the non-dispersive effects should be more
pronounced.

Exactly quantifying these non-dispersive effects with a full tomography of
multi-qubit QND readout using Hamiltonians (27) or (28) is numerically very
challenging. This is not only because the qubit has a larger dimension, but
also because one needs to account for many more photons in the dynamics of the
cavity due to the increase of $\Omega_{c}T$ to compensate the reduction of
$\chi^{\prime}$. We estimate that increasing $\Omega_{c}T$ by a factor 2.5
with respect to what was used in Fig. 3, implies duplicating the integration
time from $T=8/\kappa$ to $T=16/\kappa$, as well as increasing the steady-
state cavity population $\braket{a^{\dagger}a}$ from 1.5 to 2.5. This change
of parameters seems small but has a great impact on the time to carry out our
simulations since solving the stochastic master equation with high precision
requires many trajectories. In particular, when using the qutrit model (28)
with the above new parameters, the dimension of the total Hilbert space
becomes $33$ ($3$ states in qutrit and $11$ photons of cavity) instead of the
dimension $16$ used in obtain Fig. 3 of the main text ($2$ states in qubit and
$8$ photons in cavity). This means that each trajectory simulated with the
stochastic master equation will take 120 seconds instead of 11 seconds, making
the whole simulation of the QND detector tomography to last 10 times longer.
Since doing our full QND detector tomography with the precision of Fig. 3 of
the main text requires 20.000 individual trajectories for each of the 18
combinations of initial states and gates, our simulations took around 10 days,
implying that the qutrit simulation will around 3 months. Optimizing such a
simulation is unfeasible with the standard numerical methods used in this
work, but it is possible to develop radically different and more efficient
simulation techniques to provide a systematic study of QND readout of multi-
level qubits in the near future. Independent of this, all these interesting
non-dispersive effects can be immediately experimentally explored by
implementing our QND detector tomography procedure to a low anharmonic qubit
such as a transmon.