# Multiplexed Photon Number Measurement

Antoine Essig, Quentin Ficheux, Théau Peronnin, Nathanaël Cottet, Raphaël Lescanne, Alain Sarlette, Pierre Rouchon, Zaki Leghtas, Benjamin Huard

► To cite this version:

Antoine Essig, Quentin Ficheux, Théau Peronnin, Nathanaël Cottet, Raphaël Lescanne, et al. Multiplexed Photon Number Measurement. Physical Review X, 2021, 11 (3), 10.1103/Phys-RevX.11.031045. hal-03443020

HAL Id: hal-03443020

https://minesparis-psl.hal.science/hal-03443020v1

Submitted on 29 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Multiplexed Photon Number Measurement

Antoine Essig,¹,* Quentin Ficheux,¹,*, Théau Peronnin,¹ Nathanaël Cottet,¹ Raphaël Lescanne,²,³ Alain Sarlette,²,³,⁴ Pierre Rouchon,⁵,³ Zaki Leghtas,⁵,²,³ and Benjamin Huard¹

¹Univ Lyon, ENS de Lyon, CNRS, Laboratoire de Physique, F-69342 Lyon, France

²Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, Paris, France

³QUANTIC team, INRIA de Paris, 2 Rue Simone Iff, 75012 Paris, France

⁴Department of Electronics and Information Systems, Ghent University, Belgium

⁵Centre Automatique et Systèmes, Mines-ParisTech, PSL Research University, 60, bd Saint-Michel, 75006 Paris, France

<figure><img src="image_2.png" /></figure>

(Received 5 February 2021; revised 28 May 2021; accepted 29 June 2021; published 27 August 2021)

When a two-level system—a qubit—is used as a probe of a larger system, it naturally leads to answering a single yes-no question about the system state. Here we propose a method where a single qubit is able to extract, not a single, but many bits of information about the photon number of a microwave resonator using continuous measurement. We realize a proof-of-principle experiment by recording the fluorescence emitted by a superconducting qubit reflecting a frequency comb, thus implementing multiplexed photon counting where the information about each Fock state—from 0 to 8—is simultaneously encoded in independent measurement channels. Direct Wigner tomography of the quantum state of the resonator evidences the backaction of the measurement as well as the optimal information extraction parameters. Our experiment unleashes the full potential of quantum meters by replacing sequential quantum measurements with simultaneous and continuous measurements separated in the frequency domain.

DOI: 10.1103/PhysRevX.11.031045

Subject Areas: Quantum Physics, Quantum Information

# I. INTRODUCTION

The most general measurement of a quantum system consists of using a quantum apparatus as a probe. The system interacts with the probe before the probe gets measured projectively. In the simplest case, the probe is a qubit whose readout answers a yes-no question about the system state. Identifying what is the state of a system thus comes down to playing a game of “guess who?” A series of binary questions are asked iteratively to refine our knowledge about the state. Unlike the classical game, each answer disturbs the state of the system. To give a concrete example, in order to determine how many photons are stored in a cavity, one may ask “is there an even number of photons?” or a series of binary questions such as “are there *n* photons?” for each integer *n* [Fig. 1(a)].

Such experiments have been implemented with Rydberg atoms or superconducting circuits probing a microwave cavity [1,2] with the possible refinement of choosing what binary question should be optimally asked conditioned on

the previous answers [3]—using a feedback loop [4,5] or advanced pulse shaping [6,7]. Determining an arbitrary number of photons in the cavity between 0 and $2^m - 1$ thus takes at least *m* consecutive probe measurements since each answer provides at most 1 bit of information about the system state. This limitation originates from the encoding of the extracted information into the quantum state of the qubit. But is it the best use of a qubit to determine an observable with many possible outcomes such as a photon number?

We propose to use a qubit as an *encoder* of information about the cavity state in the many propagating modes of a transmission line. Assuming an ideal detector, we show that photon counting can then be implemented in a time independent of the number of photons. We demonstrate the practicality of this approach in an experiment where information about nine possible photon numbers (more than 3 bits) in a microwave resonator is simultaneously extracted by a single superconducting qubit into nine propagating modes of a transmission line. Owing to dispersive interaction [8,9], each photon number corresponds to a unique qubit frequency. When driving the qubit at nine test frequencies by multiplexing, the qubit simultaneously emits nine microwave signals that each reveals information about the photon number ranging from 0 to 8. Daring an analogy with communication protocols [10], previous measurement schemes with time series of binary

<footnote>*A. E. and Q. F. contributed equally to this work.</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Multiplexed photon counting setup and protocol. (a) Time-domain division multiplexing proceeds one binary question at a time. Frequency division multiplexing simultaneously retrieves multiple binary answers. (b) Scheme of the device in coplanar waveguide architecture. The storage mode (green) is coupled to a transmon multiplexing qubit (orange), which is directly coupled to a transmission line (rainbow). A directional coupler and broadband traveling-wave parametric amplifier (TWPA) allow us to probe the qubit in reflection. An additional transmon yes-no qubit (blue) and its readout resonator (purple) are used as a reference photon counter and for Wigner tomography (see Appendix B 2). (c) The frequency (color) of the multiplexing qubit encodes the storage photon number. The reduction in the reflection amplitude of the qubit at one of the frequencies reveals the number of photons in the storage mode, e.g., here two photons. (d) The qubit is probed by a frequency comb of amplitude Ω. The reflected pulse is amplified and digitized before numerical demodulation at every frequency $f_{MP} - k\chi_{s,MP}$. This multiplexing-demultiplexing process leads to reflection coefficients $r_k$ that each encodes the probability that $k$ photons are stored.</figcaption></figure>

questions used time division multiplexing, while our experiment demonstrates the analog of frequency division multiplexing, where the qubit alone acts as the frequency multiplexing transducer [Fig. 1(a)]. This experiment directly benefits from the recent bandwidth improvements of near-quantum-limited amplifiers [11], which enable us to bring the measurement process in the frequency domain. Current limitations in the cavity lifetime and detector efficiency prevent us from reaching single-shot readout of the photon number in this proof-of-principle experiment, and hence from decoding 3 bits of information per experiment. However, unlike in sequential measurement schemes, a single run of our experiment does provide, in parallel, partial information about each bit of the photon number. In addition, we manage to observe the multiplexed measurement backaction on the resonator using direct Wigner tomography, which allows us to measure the decoherence rate of the resonator induced by the measurement. We evidence an optimal qubit drive amplitude for information extraction, which matches the expected dynamics of a qubit under a multifrequency drive.

The article is organized as follows: Sec. II demonstrates standard photon counting using a qubit dispersively coupled to a microwave resonator and reviews the practical limitations of this technique. It provides a calibration of the photon number for the rest of the article. Section III describes photon counting using the fluorescence of a second qubit. Section IV presents a gedanken experiment that shows how multiplexed photon number measurement can reach an outcome in a time that does not depend on the number of photons. Section V presents an actual experiment that implements a version of the gedanken experiment where the photodetectors are replaced by heterodyne detectors in each frequency band. In practice, it comes down to experimentally probing the qubit in fluorescence with a frequency comb, thus revealing how information about cavity photon numbers is routed by the qubit onto the transmission line modes. Section VI evidences the measurement backaction of the multiplexed photon counting using direct Wigner tomography of the quantum state and shows how to maximize information extraction. A detailed comparison between the efficiency of our approach
compared to other measurement schemes is given in Appendix E.

## II. STANDARD PHOTOCOUNTING

The experiment implements both the standard approach and the multiplexed one in order to count the number of photons in a resonator dubbed the *storage mode*, which resonates at $f_s = 4.558$ GHz. Two off-resonant transmon qubits are coupled to the storage mode [see Fig. 1(b)]. The *yes-no qubit* with a frequency $f_{YN} = 3.848$ GHz is used to ask standard binary questions about the photon number or to perform storage-mode tomography, while the *multiplexing qubit* with a frequency $f_{MP} = 4.238$ GHz is used for fluorescence photon counting (Sec. III) and frequency multiplexed photon counting (Secs. V and VI). Both qubits are dispersively coupled to the resonator so that their frequency redshifts, respectively, by $\chi_{s,YN} = 1.4$ MHz and $\chi_{s,MP} = 4.9$ MHz per additional photon in the storage mode.

In the standard approach [2,8], which probes whether there are $k$ photons, the probability to have $k$ photons is encoded as the probability $P_e$ to excite the yes-no qubit by driving it with a $\pi$ pulse at $f_{drive} = f_{YN} - k\chi_{s,YN}$. The state of the yes-no qubit is read out using a dedicated resonator [Fig. 1(b)]. To demonstrate this photon counting ability, we use a microwave tone at $f_s$ to prepare the storage mode in a coherent state $|\beta\rangle = e^{-|\beta|^2/2} \sum_{n=0}^{+\infty} (\beta^n / \sqrt{n!}) |n\rangle$, which is a superposition of all Fock states with mean photon number $\bar{n} = |\beta|^2$. The probability $P_e$ is then measured and shows resolved peaks as a function of $f_{drive}$ for every photon number up to about 7 [Figs. 2(a) and 2(b)]. For the rest of the paper, we use this measurement as a calibration of the photon number in the storage mode. The linear relation between $\beta$ and the amplitude $V_s$ of the tone at $f_s$ is extracted using a master-equation-based model (see Appendix Sec. B 1) reproducing the measured $P_e$ [solid lines in Fig. 2(b)].

With this approach, the choice of binary question asked to the resonator is set by the frequency of the drive used to perform the conditional $\pi$ pulse on the qubit. Finding the photon number $N$ by asking successively whether the resonator is in state $|k\rangle$ for $k = 0, ..., N$ takes $N + 1$ consecutive measurements and is highly sensitive to measurement errors (see Appendix F 2). Each step reveals at most 1 bit of information about the photon number in the resonator since the qubit state encodes the information. The duration of this procedure can be reduced by adaptive measurement [3]. It reaches a minimal number of measurement steps using binary decimation at the expense of using a feedback loop to adjust the pulse sequence in real time [4,5] or numerical optimal control techniques such as gradient ascent pulse engineering [6,7]. Therefore, the best strategy using a qubit state as the probe of the photon number requires advanced control techniques and at least of

<figure><img src="image_4.png" /><figcaption>FIG. 2. Standard photon counting. The storage mode is prepared in a coherent state with an average photon number $\bar{n}$ using a microwave pulse at storage frequency and amplitude $V_s$. (a),(b) Measured probability $P_e$ that the yes-no qubit gets excited by a $\pi$ pulse at a frequency $f_{drive}$. Peaks appear at $f_{YN} - k\chi_{s,YN}$ and indicate the probability to store $k$ photons. The dots in (b) are cuts along the dashed lines in (a) and match the master equation model (solid lines), hence providing a calibration of $\bar{n}$ as a function of the drive amplitude $V_s$.</figcaption></figure>

the order of $\log_2(N)$ measurement steps that each requires a time of the order of $1/\chi_{s,YN}$ at best (see Appendix F 2).

## III. FLUORESCENCE PHOTON COUNTING

The intrinsic limitation of the standard approach is that measuring the qubit state can at most reveal 1 bit of information per step. It is possible to avoid this constraint by observing the qubit frequency directly instead of measuring its state. The multiplexing qubit is coupled to the transmission line so that when there are $k$ photons in the storage mode, the qubit emits a fluorescence signal into the mode of the transmission line that is centered at the qubit frequency $f_{MP} - k\chi_{s,MP}$. This encoding ability can be observed by driving the multiplexing qubit with a single microwave drive in reflection through the transmission line [Fig. 1(b)] [12–15]. The measured real part $\text{Re}(r)$ of the reflection coefficient of a microwave pulse at frequency $f_{probe}$ is reduced when the probe resonates with the qubit, hence revealing the photon number $k$ [Fig. 1(c)] [9]. This reduction arises from the coherent emission by the qubit in phase opposition with the reflected drive [16]. Therefore, on average, the distribution of photon numbers in the
<figure><img src="image_5.png" /><figcaption>FIG. 3. Fluorescence photon counting. (a),(b) Measured emission coefficient $1 - \text{Re}(r)$ as a function of the probe frequency $f_{\text{probe}}$, and the mean photon number $\bar{n}$ in the storage mode. The emission coefficient exhibits a resolved peak for each photon number. The dots in (b) are cuts along the dashed lines in (a) and are captured by a master equation model (solid lines).</figcaption></figure>

storage mode can be deduced from the relative amplitudes of the reduction of $\text{Re}(r)$ at each frequency $f_{\text{MP}} - l\chi_{s,\text{MP}}$.

In Figs. 3(a) and 3(b), we show the measured qubit emission coefficient $1 - \text{Re}(r)$ as a function of a single probe frequency $f_{\text{probe}}$ and of the initial amplitude of the storage-mode coherent state $\sqrt{\bar{n}}$. The measurement is performed using a drive strength $\Omega = \chi_{s,\text{MP}}/4$ (expressed as the corresponding Rabi frequency) and pulse duration of $2~\mu\text{s}$, which is smaller than the storage lifetime of $3.8~\mu\text{s}$. Resolved peaks develop for every photon number up to at least 9. Using the former calibration of $\bar{n}$, a master-equation-based model enables us to reproduce the measurement results (see Appendix G 1).

The observation of resolved peaks is due to our choice of parameters. We design the relaxation rate of the multiplexing qubit $\Gamma_{1,\text{MP}} = (42~\text{ns})^{-1}$ so that the decoherence rate $\Gamma_{2,\text{MP}} = \Gamma_{1,\text{MP}}/2$ is smaller than the dispersive shift $2\pi\chi_{s,\text{MP}}$. When peaks are separated, probing the qubit at one of its resonance frequencies $f_{\text{MP}} - k\chi_{s,\text{MP}}$ opens a communication channel with a maximal bandwidth $\Gamma_{2,\text{MP}}$ carrying information only about Fock state $|k\rangle$. We maximize the bandwidth of each channel by designing $\Gamma_{1,\text{MP}}$ as large as possible by adjusting the direct coupling to the transmission line, under the constraint of keeping the peaks resolved (see Appendix A 4).

Therefore, we show that both the fluorescence photon counting and the standard photon counting (Figs. 2 and 3) allow us to ask questions of the kind, “are there $k$ photons?” The important difference between both techniques is that only the fluorescence photon counting can be multiplexed. Indeed, for the standard technique, one needs to read out and reset the qubit at the end of each step. The readout step cannot be multiplexed as it always occurs at the readout mode frequency. In contrast, with the fluorescence readout, information about a given photon number $k$ is constantly extracted through the frequency mode $f_{\text{MP}} - k\chi_{s,\text{MP}}$ of the transmission line. It thus enables the key ingredient of our approach: the multiplexing measurement of reflection at every frequency $f_{\text{MP}} - k\chi_{s,\text{MP}}$. The qubit thus acts as an encoder of the state of the storage mode into the many modes of the transmission line at frequencies $\{f_{\text{MP}} - l\chi_{s,\text{MP}}\}_l$, which can collectively host much more than a single bit of information.

IV. GEDANKEN MULTIPLEXED EXPERIMENT

In this section, we show how, despite using a single qubit as well, multiplexed measurements are able to determine the photon number in a constant time in contrast with the standard approach. We consider an ideal detector for the propagating modes in order to better illustrate the power of multiplexing. The ideal detector is made of a frequency multiplexer followed by a perfect photodetector on each of its outputs (Fig. 4). The multiplexer is made of a parallel ensemble of bandpass filters that are each centered on the frequency $f_{\text{MP}} - k\chi_{s,\text{MP}}$ with a bandwidth $\chi_{s,\text{MP}}$. The protocol proceeds in three steps to count the number of photons in the storage mode starting in state $|\psi\rangle_s$, as detailed in Fig. 4. First, the multiplexing qubit is excited with a $\pi$ pulse that is short enough so that it prepares the qubit in the excited state irrespective on the number of photons. Second, the qubit decays in the transmission line converting its excitation into a single photon contained in a propagating wave packet whose envelope decays at a rate $\Gamma_{1,\text{MP}}$. In the limit where $\Gamma_{1,\text{MP}} \ll \chi_{s,\text{MP}}$, and without pure dephasing of the qubit, the photon emission produces an entangled state between the storage mode and the propagating modes of the line

$$ \sum_k \langle k | \psi \rangle_s \otimes_j | \delta_{k,j} \rangle_j \otimes | k \rangle_s, $$

where $|\cdot\rangle_j$ represents the quantum state of the propagating mode going through the multiplexer on branch $j$ corresponding to frequencies in the band $[f_{\text{MP}} - (j + 1/2)\chi_{s,\text{MP}}, f_{\text{MP}} - (j - 1/2)\chi_{s,\text{MP}}]$ [see Fig. 4(c)]. Matching the temporal envelope of the modes to the exponential decay at a rate $\Gamma_1$ [17], the mode is occupied by either $|0\rangle_j$ or $|1\rangle_j$ depending on the storage photon number, hence the notation $|\delta_{k,j}\rangle_j$. Finally, a single
<figure><img src="image_6.png" /><figcaption>FIG. 4. Gedanken multiplexing experiment. (a) An unconditional $\pi$ pulse is applied to the multiplexing qubit while the cavity is prepared in state $|\psi\rangle_s$. (b) The qubit is prepared in the excited state. (c) The qubit spontaneously emits a photon into the transmission line, where a multiplexer sorts the emitted radiation according to its frequency. Each port $k$ of the multiplexer is bandpass filtered around frequency $f_{MP} - k\chi_{s,MP}$ by a rectangular component displaying the frequency band. (d) Eventually, a single photodetector (detector $k=1$ in the figure) clicks with probability $|\langle k|\psi\rangle|^2$, allowing us to deduce the photon number. The storage mode is projected on the corresponding Fock state (here, $|1\rangle_s$) in a typical time $T_{1,MP}$ that does not depend on the average number of photons in the storage mode.</figcaption></figure>

photodetector clicks and reveals the number of photons $k$ with probability $|\langle k|\psi\rangle|^2$ [Fig. 4(d)]. In case of ideal detectors with zero false positives, the click detects the associated propagating mode in $|1\rangle_k$, and therefore, as the line is entangled with the storage mode, the measurement backaction projects the storage mode in Fock state $|k\rangle$. The total measurement time, a few $1/\Gamma_{1,MP}$, corresponds to the time it takes for one photodetector to click. The time is thus independent of the number of photons stored in the storage mode. Note that in order to avoid spectral leakage into other ports, $\Gamma_{1,MP}$ is limited by $\chi_{s,MP}$ so that the shortest measurement time is limited to a few $1/\chi_{s,MP}$.

In contrast to sequential measurements for which increasing the maximal number of photons that can be detected requires additional temporal resources [of the order of $\log_2(N)/\chi_{s,MP}$], this gedanken experiment shows that the multiplexed measurement is able to operate in a constant time at the expense of additional spectral resources.

## V. MULTIPLEXED PHOTON COUNTING

In practice, building such an array of frequency-sensitive photodetectors remains an open challenge in the microwave domain, despite encouraging recent progress toward this goal [5,17–22]. In this section, we demonstrate an actual experiment that implements a continuous version of the multiplexing measurement using heterodyne detectors instead of photodetectors. We multiplex the fluorescence photon counting of Sec. III by sending a pulse containing a comb with nine frequencies corresponding to photon numbers from 0 to 8. We then demultiplex the reflected pulse at the same nine frequencies $\{f_{MP} - k\chi_{s,MP}\}_{0\le k\le 8}$

and extract a reflection coefficient $r_k$ for each of them [Fig. 1(d)]. The measurement consists of simultaneously measuring the emission coefficients $1 - \text{Re}(r_k)$ for each peak in Figs. 3(a) and 3(b), which is much faster than measuring them one at a time. Figure 5 shows these emission coefficients as a function of the average initial photon number for a drive strength $\Omega = \chi_{s,MP}/2$ and a measurement duration of $2\ \mu\text{s}$. For a given $\bar{n}$, every measurement channel $k$ gives an average signal that is proportional to the probability of having $k$ photons in the storage mode. As $\bar{n}$ is varied, the shape of the average signal of channel $k$ reproduces a Poisson distribution distorted by relaxation processes and channel cross talk that increases with driving strength (see Appendix G 3). This multiplexed photon counting signal can be reproduced using a master equation approach (solid lines in Fig. 5) using the photon number calibration of the standard photon counting approach. This result thus demonstrates the applicability of our approach to photon counting by simultaneously probing information about the presence of nine possible photon numbers in the resonator.

So how does this proof-of-principle experiment compare with standard photocounting? Each method has its own advantages and drawbacks. The multiplexed photon counting scheme trades off the temporal constraint and complexity of optimal control of the standard approach for the need of an efficient quantum measurement on a large frequency bandwidth. Indeed, the efficient measurement of the reflected pulse requires the use of a near-quantum-limited amplifier with a dynamical bandwidth of at least a dozen of $\chi_{s,MP}$, which is now possible using a traveling-wave parametric amplifier (TWPA) [11]. It is a comparable technical requirement to the recently demonstrated high-efficiency
<figure><img src="image_7.png" /><figcaption>FIG. 5. Multiplexed photon counting. Dots: simultaneously measured average emission coefficients corresponding to every photon number *k* from 0 to 8 as a function of the initial mean photon number $\bar{n}$ in the storage mode. Here, $r_k$ is the reflection coefficient at $f_{MP} - k\chi_{s,MP}$. Solid lines: prediction based on a master equation without free parameters (see Appendix G 3).</figcaption></figure>

multiplexed readout of as many as six qubits coupled to a
single feedline [23–26]. The number-independent meas-
urement time in the case of ideal photodetectors (Sec. IV)
relies on the absence of noise when measuring a mode in
the vacuum state. Instead, heterodyne detectors produce at
least vacuum fluctuations in each frequency band, as
Heisenberg uncertainty relations command. For this reason,
the measurement time is expected to scale as the logarithm
of the photon number similar to state-of-the-art sequential
measurement schemes. The two main advantages of our
multiplexing technique over the standard approaches is that
the prefactor is not limited by feedback latency or optimal
control duration (see Appendix F 1) and that it is a

continuous measurement that does not require any subtle
temporal control.

VI. MEASUREMENT BACKACTION

The measurement strength of the multiplexing measure-
ment can be characterized using the yes-no qubit to observe
the dynamics of the cavity state under the action of the
continuous multiplexed measurement. The advantage of
this method is that it does not require a single-shot
measurement of the photon number, which we cannot
reach owing to the limited efficiency of our amplifier and
the too short lifetime of the storage mode. In the reciprocal
case of measuring a qubit using a cavity as a probe, the
measurement rate is bounded by the dephasing rate of the
qubit, which grows as the square of the cavity-driving
strength [27,28]. Thus, characterizing the measurement rate
of our multiplexed photon counting can be done by
observing how the storage mode dephases for a given
driving strength Ω. Indeed, owing to the inherent quantum
backaction of the photon number measurement, the meas-
urement rate is bounded by how fast the conjugated
operator, here the mode phase, diffuses. As the probe is
based on a qubit driven by a frequency comb, one expects a
different dependence of the measurement rate on Ω than for
standard dispersive qubit readout using a single tone
driving a probe cavity.

In order to measure this dephasing rate, we use the yes-
no qubit to perform a direct Wigner tomography [29–31] of
the storage mode at various times t. It provides a repre-
sentation of the state ρ in the phase space of the mode and
can be expressed as

<figure><img src="image_8.png" /><figcaption>FIG. 6. Measurement backaction. Direct Wigner tomography of the storage mode at time $t = 0.5$ µs after initialization in a coherent state with amplitude $\beta = -1.55$. Left: free evolution without driving. Middle: evolution with a single tone at $f_{MP} - \chi_{s,MP}$ that probes whether there is one photon in the storage mode. Right: evolution in the presence of the frequency comb $\{f_{MP} - k\chi_{s,MP}\}_{0\le k\le 8}$. Each bottom panel represents the density matrix $\{\rho_{nm}\}_{n,m}$ that is calculated from the Wigner function of its top panel. The cylinder heights represent $|\rho_{nm}|$ and the color encodes $\arg(\rho_{nm})$.</figcaption></figure>
$$W(\alpha) = \frac{2}{\pi} \mathrm{Tr}[D^{\dagger}(\alpha)\rho D(\alpha)\mathcal{P}]. \qquad (1)$$

Here, $D(\alpha) = e^{\alpha a_s^\dagger - \alpha^* a_s}$ is the storage displacement operator, $\mathcal{P} = e^{i\pi a_s^\dagger a_s}$ is the photon number parity operator, and $a_s$ is the canonical annihilation operator of the storage mode. Preparing the storage mode in a coherent state $|\beta = -1.55\rangle$, the Wigner function starts as a Gaussian distribution centered at $\alpha = \beta$. On the left of Fig. 6, one can see how the bare dephasing rate and the self-Kerr effect of the storage mode (frequency shift of 0.02 MHz per photon) distort the Gaussian distribution toward a torus with no phase when time increases even without any photon counting drives. Using a single drive with $\Omega = \chi_{s,\mathrm{MP}}/2$ to measure whether there is one photon, the phase diffuses faster and the Wigner function exhibits negativities in the middle of Fig. 6. As seen in the corresponding density matrix, a tone at $f_{\mathrm{probe}} = f_{\mathrm{MP}} - \chi_{s,\mathrm{MP}}$ notably induces dephasing between states $|1\rangle$ and all other states $|m \neq 1\rangle$ (see density matrix as a function of the drive frequency in Appendix H 4). The phase diffusion is more intense when all the tones of the multiplexed readout are turned on than for a single tone with the same drive strength $\Omega$ (right of Fig. 6). Likewise, all off-diagonal elements of the density matrix are then reduced.

To be more quantitative, the dephasing rate $\Gamma_{d,s}$ of the cavity is accessed through the decay of the mean quadrature of the storage mode [32] (see Appendix B 3)

$$\langle a_s + a_s^\dagger \rangle = \int 2xW(x + iy)dxdy. \qquad (2)$$

In Fig. 7(a), we show $\langle a_s + a_s^\dagger \rangle$ as a function of time under the multiplexed drive with a strength $\Omega = 0.23\chi_{s,\mathrm{MP}}$. Repeating this experiment for various values of the multiplexed driving strength $\Omega$ allows us to determine how it affects the dephasing rate $\Gamma_{d,s}$, and thus the measurement rate. The dephasing rate is nonmonotonic in the drive strength [Fig. 7(b)]. Noticeably, it reaches a maximum when $\Omega = \chi_{s,\mathrm{MP}}/2$ for which information is extracted at a rate approximately 5 times larger than the natural dephasing rate. It is possible to understand this behavior by considering a model system where the comb has an infinite number of Dirac peaks $\{f_{\mathrm{MP}} + k\chi_{s,\mathrm{MP}}\}_{k \in \mathbb{Z}}$ (see Appendix C). The Fourier transform of a comb being a comb, the drive performs sudden rotations by an angle $2\pi\Omega/\chi_{s,\mathrm{MP}}$ of the Bloch vector of the qubit every time step $1/\chi_{s,\mathrm{MP}}$. When $\Omega/\chi_{s,\mathrm{MP}}$ is an integer, the comb does not affect the qubit and thus $\Gamma_{d,s}$ vanishes. Conversely, the maximum measurement rate corresponds to half-integer $\Omega/\chi_{s,\mathrm{MP}}$ for which the effect of the comb on the qubit dynamics is maximum and leads to the strongest qubit emission. With the finite comb used in the experiment, this maximum persists and is reproduced by a model based on a

<figure><img src="image_9.png" /><figcaption>FIG. 7. Measurement-induced dephasing of the multiplexing measurement. (a) Measured exponential decay of the average storage-mode quadrature $\langle a_s + a_s^\dagger \rangle$ in the case of a driving comb of strength $\Omega = 0.23\chi_{s,\mathrm{MP}}$, and in a frame rotating at the storage-mode resonant frequency. Insets show Wigner tomography of the storage mode for two values of $t$. (b) Measured dephasing rate induced on the storage mode as a function of the drive strength $\Omega$ of the frequency comb. Insets show Wigner tomographies of the storage mode at $t = 1$ µs for four values of $\Omega$.</figcaption></figure>

master equation without any free parameter [line in Fig. 7(b)].

## VII. CONCLUSION

We experimentally demonstrate that a single qubit can be used to continuously probe a multidimensional system by encoding information about its quantum state in the frequency domain. Improving further the detection efficiency $\eta$, the dispersive shift $\chi_{s,\mathrm{MP}}$ and the coupling rate $\Gamma_{1,\mathrm{MP}}$ between the qubit and the transmission line (while protecting the storage lifetime with a notch filter at its frequency) should enable single-shot photon counting by multiplexing. To be more accurate, if the parameter $\eta\Gamma_{1,\mathrm{MP}}/\Gamma_{1,s} = 17$ increases by an order of magnitude, single-shot measurements would be possible. Interestingly, assuming perfect detectors, our gedanken experiment shows how the multiplexed measurement can determine a photon number $N$ in a time that is $\log(N)$ faster than the best standard sequential approach. Our continuous measurement opens new possibilities in terms of feedback control of the quantum state of a cavity. It can readily be applied to stabilize quantum states by feedback control [33], probe quantum trajectories of microwave modes [34],
observe quantum Zeno dynamics [35], or engineer desired decoherence channels by varying in time the amplitude of the probe tones. This measurement scheme enables the future implementation of a large class of measurement operators that would be useful to stabilize bosonic codes [36], stabilize a Fock state parity by autonomous feedback [37], or extend the reach of simultaneous probing of a single quantum system by multiple observers [38, 39] to larger systems and arbitrarily many observers. Our photo-counter for stationary modes can also be converted into a photocounter for propagating modes using a catch-and-count protocol [5]. Moving further, one could extend this frequency domain measurement to more complex probes than a single qubit and many possible physical systems beyond superconducting circuits.

## ACKNOWLEDGMENTS

This work is made possible through the support of Fondation Del Duca, IDEX Lyon (Contract No. ANR-16-IDEX-0005) and Contract No. ANR JCJC-HAMROQS. We acknowledge IARPA and Lincoln Labs for providing a Josephson traveling-wave parametric amplifier. The device was fabricated in the clean rooms of Collège de France, ENS Paris, CEA Saclay, and Observatoire de Paris. We are grateful to Rémi Azouit, Michel Devoret, and Mazyar Mirrahimi for discussions.

# APPENDIX A: DEVICE AND MEASUREMENT SETUP

## 1. Design

The circuit is composed of four electromagnetic modes whose parameters can be found in Table I. A high-Q harmonic oscillator called the storage mode is composed of a $\lambda/2$ coplanar waveguide resonator (green in Fig. 8). The storage resonator is capacitively coupled to two transmon qubits. The multiplexing qubit (orange) has a high spontaneous photon emission rate $\Gamma_{1,MP} = (44 \text{ ns})^{-1}$ into a transmission line compared to other modes. In contrast, the yes-no qubit is capacitively coupled to a low-Q readout resonator and has a long coherence time $T_{2,YN} = 27 \text{ µs}$. As required by Wigner tomography, the yes-no qubit coherence time and the lifetime of storage mode are larger than the time needed to measure the parity of storage photon number $1/2\chi_{s,YN} \ll T_{1,s}, T_{2,YN}$. As we use the multiplexing qubit to count the photon number in the storage mode, we need it to be photon number resolved [8], otherwise each record of the multiplexing measurement cannot be associated with a single specific photon number. This photon-number-resolved constraint imposes that the multiplexing qubit decoherence rate must be smaller than the cross-Kerr rate between the multiplexing qubit and the storage mode $\Gamma_{2,MP} < 2\pi\chi_{s,MP}$. This resolution constraint is not critical, as in fact, a finite amount of photon number information can be extracted as soon as $\chi_{s,MP}$ is nonzero,

<figure><img src="image_10.png" /><figcaption>FIG. 8. Optical image of the circuit. The readout resonator is colored in purple, storage mode in green, yes-no qubit in blue, and multiplexing qubit in orange. All dark gray areas are silicon, gray areas are niobium on silicon, and Josephson junctions are made of Al/AlOx/Al.</figcaption></figure>

but the decoding is much simpler if we can reason in terms of well-separated resonance peaks.

## 2. Device fabrication

The length of the readout resonator and storage mode is designed to obtain the resonant frequencies $f_{RO} \sim 7$ GHz and $f_s \sim 4.5$ GHz. The circuit consists of a sputtered 120-nm-thick niobium film deposited on a 280-µm-thick undoped silicon wafer. The resonators and feedlines are dry etched after optical lithography. After an ion milling step, the Josephson junctions are made out of e-beam-evaporated Al/AlOx/Al through a PMMA and MAA resist mask patterned in a distinct e-beam lithography step. For each transmon qubit, a single Dolan bridge is used to make the junctions.

## 3. Measurement setup

The readout resonator, the yes-no qubit, and the multiplexing qubit are driven by pulses that are generated using a Tektronix® arbitrary waveform generator (AWG) AWG5014C with a sample rate of 1 GS/s. Storage-mode pulses are generated using a Zurich Instruments® ultra-high frequency lock-in amplifier (UHFLI) with a sample rate of 1.8 GS/s. The UHFLI allows us to change the pulse amplitude and phase without recompiling the sequence. This feature decreases the time needed for Wigner tomography compared to a standard AWG and makes the pulse sequence simple with the drawback of having to synchronize the two AWGs. AWG pulses are modulated at a frequency 25 MHz for readout, 100 MHz for yes-no qubit, and 75 MHz for storage and multiplexing qubit. They are
up-converted using single sideband mixers for the readout resonator and multiplexing qubit and regular mixers for the storage resonator and yes-no qubit, with continuous micro-wave tones produced, respectively, by AnaPico® APSIN12G, Agilent® E8257D, WindFreak® SynthHD, and AnaPico® APSIN20G sources that are set at the frequencies $f_{RO} + 25$ MHz, $f_{MP} + 75$ MHz, $f_s + 75$ MHz, and $f_{YN} + 100$ MHz.

The two reflected signals from the readout and multiplexing qubit are combined with a diplexer and then amplified with a TWPA provided by Lincoln Labs. We tune the pump frequency ($f_{TWPA} = 5.998$ GHz) and power in order to reach a gain of 20.7 dB at 7.138 GHz and 18.2 dB at 4.238 GHz. The quantum efficiency of the yes-no readout signal is measured to be $18.7\% \pm 0.4\%$ and should be close to the efficiency $\eta$ of the multiplexing detection. We estimate that this efficiency is the product of the efficiency of the microwave components before the TWPA (25% to 60%), the efficiency of the TWPA itself (33% to 83%), and the (90% to 95%) efficiency coming from what is above the high-electron-mobility-transistor (HEMT) amplifier. The follow-up amplification is performed by a HEMT amplifier from Low Noise Factory® at 4 K and by two room-temperature amplifiers. The two signals are down-converted using image reject mixers before digitization by an Alazar® acquisition board and numerical demodulation. Actually, for the multiplexed

<figure><img src="image_11.png" /><figcaption>FIG. 9. Schematic of the setup. Each electromagnetic mode of the experiment is driven by a rf generator detuned by the modulation frequency and whose color matches the color of the corresponding mode in Fig. 8. Room-temperature isolators are not represented for the sake of clarity.</figcaption></figure>

signal, nine demodulation operations are performed at each of the down-converted frequencies $75$ MHz + $k\chi_{s,MP}$ for $0 \le k \le 8$. The full setup is shown in Fig. 9. The Tektronix® AWG is used as the master that triggers the UHFLI and the Alazar® board.

The frequency comb that is used for the multiplexing measurement is generated and demodulated using the following method. Nine cosine functions at frequencies $\{75$ MHz + $k\chi_{s,MP}\}_{0 \le k \le 8}$ are summed and multiplied by a Gaussian envelope numerically with a sampling rate of 1 GS/s over the duration of the pulse. A waveform is then generated by the AWG following this list of values. This method ensures a good phase coherence between all the comb frequencies. The AWG output is up-converted using a single sideband mixer whose local-oscillator (LO) port is driven at frequency $f_{MP} + 75$ MHz.

### 4. Circuit parameters and master equation

All parameters of the four modes can be measured using standard-circuit-QED measurement (see Table I). The frequencies of the readout mode and multiplexing qubit are measured by spectroscopy. The frequencies of the storage mode and yes-no qubit are measured using two-tone spectroscopy with the readout mode. The yes-no qubit decay and decoherence rate are measured with the time evolution of the probability to find the qubit excited after a $\pi$ pulse and using Ramsey oscillations. The readout-mode decay rate and cross-Kerr rate between the readout mode and yes-no qubit are measured using the measurement-induced dephasing rate by the readout mode on the yes-no qubit. The cross-Kerr rate between the storage mode and the two qubits is measured using qubit spectroscopy with the storage state initialized in various coherent states. Anharmonicities are measured using spectroscopy of the qubit excited state. The decay and decoherence rates of the storage mode are measured with the time evolution of the probability to have zero photons in the storage mode after a displacement and storage Ramsey interferometry experiment. The multiplexing qubit decay and decoherence rates are measured by fitting the qubit spectroscopy for various drive amplitudes. All those parameters enable us to write a master equation model based on the Lindblad equation with the Hamiltonian

$$ \hat{H} = h f_{RO} \hat{n}_{RO} + h f_s \hat{n}_s + h f_{YN} \hat{n}_{YN} + h f_{MP} \hat{n}_{MP} \\ - h \chi_{RO,YN} \hat{n}_{RO} \hat{n}_{YN} - h \chi_{s,YN} \hat{n}_s \hat{n}_{YN} - h \chi_{s,MP} \hat{n}_s \hat{n}_{MP} \\ - h \chi_{YN,YN} \hat{n}_{YN} (\hat{n}_{YN} - 1) - h \chi_{MP,MP} \hat{n}_{MP} (\hat{n}_{MP} - 1), \quad (A1) $$

where $\hat{n}_{RO}$, $\hat{n}_s$, $\hat{n}_{YN}$, and $\hat{n}_{MP}$ are the photon number operators, respectively, for the readout, storage, yes-no qubit, and multiplexing qubit. $\chi_{a,b}$ is the cross-Kerr rate between modes a and b. $\chi_{a,a}$ is the anharmonicity of the
<table><caption>TABLE I. Table of circuit parameters.</caption><thead><tr><th>Circuit parameters</th><th>Symbol</th><th>Hamiltonian term</th><th>Value</th></tr></thead><tbody><tr><td>Readout resonator frequency</td><td>f<sub>RO</sub></td><td>hf<sub>RO</sub>hat<sub>RO</sub></td><td>7.138 GHz</td></tr><tr><td>Storage-mode frequency</td><td>f<sub>s</sub></td><td>hf<sub>s</sub>hat<sub>s</sub></td><td>4.558 GHz</td></tr><tr><td>Yes-no qubit frequency</td><td>f<sub>YN</sub></td><td>hf<sub>YN</sub>hat<sub>YN</sub></td><td>3.848 GHz</td></tr><tr><td>Multiplexing qubit frequency</td><td>f<sub>MP</sub></td><td>hf<sub>MP</sub>hat<sub>MP</sub></td><td>4.238 GHz</td></tr><tr><td>Readout or yes-no qubit cross-Kerr rate</td><td>&chi;<sub>RO,YN</sub></td><td>-h&chi;<sub>RO,YN</sub>hat<sub>RO</sub>hat<sub>YN</sub></td><td>0.4 MHz</td></tr><tr><td>Storage or yes-no qubit cross-Kerr rate</td><td>&chi;<sub>s,YN</sub></td><td>-h&chi;<sub>s,YN</sub>hat<sub>s</sub>hat<sub>YN</sub></td><td>1.4 MHz</td></tr><tr><td>Storage or multiplexing qubit cross-Kerr rate</td><td>&chi;<sub>s,MP</sub></td><td>-h&chi;<sub>s,MP</sub>hat<sub>s</sub>hat<sub>MP</sub></td><td>4.9 MHz</td></tr><tr><td>Yes-no qubit anharmonicity</td><td>&chi;<sub>YN,YN</sub></td><td>-h&chi;<sub>YN,YN</sub>hat<sub>YN</sub>(hat<sub>YN</sub> - 1)</td><td>160 MHz</td></tr><tr><td>Multiplexing qubit anharmonicity</td><td>&chi;<sub>MP,MP</sub></td><td>-h&chi;<sub>MP,MP</sub>hat<sub>MP</sub>(hat<sub>MP</sub> - 1)</td><td>116 MHz</td></tr></tbody></table><table><thead><tr><th>Circuit parameters</th><th>Symbol</th><th>Jump operator</th><th>Value</th></tr></thead><tbody><tr><td>Readout decay rate</td><td>&Gamma;<sub>RO</sub></td><td>&Gamma;<sub>RO</sub>L(hat<sub>RO</sub>)&rho;</td><td>(40 ns)<sup>-1</sup></td></tr><tr><td>Storage decay rate</td><td>&Gamma;<sub>1,s</sub></td><td>&Gamma;<sub>1,s</sub>L(hat<sub>s</sub>)&rho;</td><td>(3.8 ns)<sup>-1</sup></td></tr><tr><td>Storage decoherence rate</td><td>&Gamma;<sub>2,s</sub></td><td>2&Gamma;<sub>&phi;,s</sub>L(hat<sub>RO</sub>)&rho;</td><td>(2 &mu;s)<sup>-1</sup></td></tr><tr><td>Yes-no decay rate</td><td>&Gamma;<sub>1,YN</sub></td><td>&Gamma;<sub>1,YN</sub>L(hat<sub>YN</sub>)&rho;</td><td>(20 &mu;s)<sup>-1</sup></td></tr><tr><td>Yes-no decoherence rate</td><td>&Gamma;<sub>2,YN</sub></td><td>2&Gamma;<sub>&phi;,YN</sub>L(hat<sub>YN</sub>)&rho;</td><td>(27 &mu;s)<sup>-1</sup></td></tr><tr><td>Multiplexing decay rate</td><td>&Gamma;<sub>1,MP</sub></td><td>&Gamma;<sub>1,MP</sub>L(hat<sub>MP</sub>)&rho;</td><td>(42 ns)<sup>-1</sup></td></tr><tr><td>Multiplexing decoherence rate</td><td>&Gamma;<sub>2,MP</sub></td><td>2&Gamma;<sub>&phi;,MP</sub>L(hat<sub>MP</sub>)&rho;</td><td>(84 ns)<sup>-1</sup></td></tr></tbody></table>

mode *a*. The master equation on the system density matrix $\rho$ reads

$$ \dot{\rho} = -\frac{i}{\hbar}[\hat{H}, \rho] + \Gamma_{RO} \mathcal{L}(\hat{a}_{RO})\rho + 2\Gamma_{\phi,s} \mathcal{L}(\hat{n}_s)\rho + \Gamma_{1,s} \mathcal{L}(\hat{a}_s)\rho \\ + 2\Gamma_{\phi,YN} \mathcal{L}(\hat{n}_{YN})\rho + \Gamma_{1,YN} \mathcal{L}(\hat{a}_{YN})\rho + 2\Gamma_{\phi,MP} \mathcal{L}(\hat{n}_{MP})\rho \\ + \Gamma_{1,MP} \mathcal{L}(\hat{a}_{MP})\rho, \quad (A2) $$

where $\mathcal{L}$ is the Lindblad superoperator defined as $\mathcal{L}(\hat{L})\rho = \hat{L}\rho\hat{L}^\dagger - \{\hat{L}^\dagger\hat{L}, \rho\}/2$ and $\hat{a}_b$ is the annihilation operator of mode *b*. For a qubit mode *b*, the dephasing rate $\Gamma_{\phi,b}$ is linked to the decoherence rate by $\Gamma_{2,b} = \Gamma_{1,b}/2 + \Gamma_{\phi,b}$.

# APPENDIX B: CALIBRATIONS

## 1. Calibration of the storage-mode-displacement amplitude

The storage mode can be displaced by driving it on resonance with a voltage $V_s(t) \cos(2\pi f_s t + \phi_s)$, where $V_s(t)$ is the pulse envelope. The driving Hamiltonian of the storage mode reads $\hbar[\epsilon_s(t)\hat{a}_s^\dagger + \epsilon_s^*(t)\hat{a}_s]$, where $\epsilon_s(t) = \mu V_s(t)e^{i\phi_s}$. The scaling factor $\mu = 1.45$ (mV $\mu$s)$^{-1}$ is calibrated by fitting the photocounting measurement results obtained using the yesno qubit with the master equation simulation (see Appendix G 1). Figure 10(a) shows the evolution of $\epsilon_s$ with $V_s$. For every experiment, the storage-mode displacements are realized using a Gaussian pulse shape $\epsilon_s(t) = \lambda(t)\epsilon_{max}$ with a maximum amplitude $\epsilon_{max}$, a width 25 ns, and a duration 100 ns. We simulate the dynamics of the storage mode under this Gaussian displacement taking into account the couplings, relaxation,

and decoherence rates (see Appendix G 2) for various amplitudes $\epsilon_{max}$. We then compute the expectation value of the number of photon operator $\langle\hat{n}_s\rangle$ at the end of the pulse. Figure 10(b) shows the square root of $\langle\hat{n}_s\rangle$ as a function of $\epsilon_{max}$. Fitting with a linear function, we find that $\sqrt{\langle\hat{n}_s\rangle} = 59.1\epsilon_{max}$. As $\epsilon_s$ increases linearly with $V_s$, $\epsilon_{max}$ increases linearly with the maximum voltage amplitude $V_{max,s}$ of the Gaussian pulse $V_s(t) = \lambda(t)V_{max,s}$. Using the two linear regressions, we can express the photon number of the storage mode as $\sqrt{\langle n_s \rangle} = (85.9 \text{ V}^{-1})V_{max,s}$.

## 2. Wigner tomography calibration

The Wigner function of a harmonic oscillator with density matrix $\rho$ is defined as $W(\alpha) = 2\text{Tr}[D^\dagger(\alpha)\rho D(\alpha)\mathcal{P}] / \pi$ where $D(\alpha) = e^{i\hat{a}_s^\dagger\alpha^*\hat{a}_s}$ is the displacement operator of the storage mode by a coherent field $\alpha$ and $\mathcal{P} = e^{i\pi\hat{a}_s^\dagger\alpha^*\hat{a}_s}$ is the photon number parity operator. A Wigner function is the expectation value of the parity after a displacement by an amplitude $\alpha$. The Wigner tomography sequence is represented on Fig. 11(a). It starts by realizing a displacement on the storage mode with a 100-ns-long Gaussian pulse at frequency $f_s$ (or detuned for Ramsey interferometry of the storage mode; see Appendix B 3) with a width of 25 ns. Then, two successive $\pi/2$ Gaussian pulses of 18 ns with a width of 4.5 ns are sent to the yes-no qubit at $f_{YN}$ and are separated by a waiting time $\Delta\tau = 337$ ns $\approx 1/2\chi_{s,YN}$. It implements a parity measurement and maps the parity of the storage mode onto the z axis on the yes-no qubit [29–31]. The sequence terminates by a 2-$\mu$s-long square pulse on the readout resonator to read out the state of the yes-no qubit. For high amplitude $\alpha$, higher-order Kerr terms
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 10. Calibration of the average number of photons $\langle n_s \rangle$ in the storage mode as a function of the displacement amplitude. (a) Evolution of the displacement amplitude $\epsilon_s$ with the pulse envelope $V_s$. The calibration is obtained by comparing the results of a photon counting experiment using the yes-no qubit with a master equation simulation (see Appendix G 1). (b) Square root of the average photon number $\langle n_s \rangle$ in the storage mode as a function of the drive amplitude. The storage is displaced by a 100-ns-long Gaussian pulse with a width of 25 ns. The same pulse shape is used in the simulation. From the two linear fits, we extract the evolution of the mean number of photons with the amplitude of the pulse $\sqrt{\langle n_s \rangle} = (85.9 \text{ V}^{-1})V_{\text{max},s}$.</figcaption></figure>

distort the Wigner function. To mitigate this effect, we interleave two sequences with a final pulse of phase either $+\pi/2$ or $-\pi/2$ [see Fig. 11(a)]. The difference between the two signals gives us the Wigner function without the distortion due to the storage-mode anharmonicity and enables us to remove low-frequency noise. The z axis of Fig. 11(b) is calibrated using the yes-no qubit Rabi oscillation amplitude to express the signal using Pauli operators. Multiplying the result by $2/\pi$ yields the Wigner function $W(\alpha)$ in Fig. 11(c).

The axes of the phase space $x, p$ are calibrated using the same pulse sequence. The photon number calibration realized before (see Appendix B 1) cannot be used here for two reasons. First, Ramsey oscillations of the storage mode impose to play the Wigner sequence with displacement pulses detuned from the storage-mode frequency, while the photon number calibration is only valid for resonant pulses. Second, high-order Kerr interaction affects

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>FIG. 11. Direct Wigner tomography of the storage mode. (a) Circuit diagram of a Wigner tomography using a parity measurement based on dispersive interaction. After a 100-ns displacement pulse on the storage mode, an unconditional $\pi/2$ pulse is applied to the yes-no qubit. The qubit evolves freely during a time $\Delta\tau = 337$ ns $\approx 1/2\chi_{s,YN}$ before a new $\pm\pi/2$ pulse is sent and the state of the yes-no qubit is measured using the readout resonator. (b) Calibration of the quadrature axes for Wigner tomography. Blue dots represent the standard deviation of the quadratures of the displaced thermal equilibrium state of the storage mode as a function of the drive amplitude for various detuning using only the photon number calibration (see Appendix B 1). In contrast, yellow dots show the same standard deviation with the noise-based quadrature calibration. (c) Wigner tomography of the storage mode. Here, the mode is prepared in two steps. First, the storage mode is displaced by a pulse with an amplitude 1.7 and then the multiplexing qubit is driven at a single tone at $f_{\text{MP}} - 1.4\chi_{s,\text{MP}}$ during 750 ns with an amplitude $\Omega = \chi_{s,\text{MP}}/2$. The appearance of negative values in the Wigner function demonstrates that one can prepare nonclassical states in the storage mode using the multiplexing qubit backaction alone.</figcaption></figure>

the calibration when the storage mode hosts a large number of photons. We decide to use the width of the Wigner function when the storage mode is in the thermal equilibrium state to calibrate the phase-space axes. For a thermal state with a thermal photon number $n_{\text{th}}$, the Wigner
function is a 2D Gaussian function with a width $\sqrt{n_{\text{th}} + 1/2}$ [34],

$$W_{\rho}(n_{\text{th}})(\alpha = x + ip) = \frac{2}{\pi} \frac{1}{2n_{\text{th}} + 1} e^{-2|\alpha|^2/(2n_{\text{th}}+1)}. \quad (B1)$$

For a thermal state displaced by an amplitude $\beta$, the Wigner function is still a 2D Gaussian function with a width $\sqrt{n_{\text{th}} + 1/2}$ but centered on $\beta$. In thermal equilibrium, the storage mode has an average photon number $n_{\text{th}} = 0.03$, which is measured using the photon counting experiment. We calibrate the quadrature axes in order to get the expected geometrical mean $\sqrt{\overline{\sigma_x \sigma_p}} = 0.53$ of the spread along the quadratures $x$ and $p$ when the storage mode is at thermal equilibrium. To take into account high-order Kerr effects, we displace the storage-mode equilibrium state and measure its Wigner function. We adjust the calibration to still find a spread of $\sqrt{\overline{\sigma_x \sigma_p}} = 0.53$. The function used for the calibration is a third-order polynomial function which gives $|\alpha|$ as a function of the pulse amplitude $V_{\text{max},s}$. We repeat this protocol for three detuning values $\delta f_s$ between the displacement pulse and storage-mode frequencies. Figure 11(b) shows the mean quadrature spread of the displaced storage-mode thermal-state Wigner function as a function of the drive amplitude $V_{\text{max},s}$ for the photon number calibration and the Wigner phase-space calibration. For example, the polynomial function for a detuning of 4 MHz reads $\alpha = x + ip = e^{i\phi_s}(77.3V_{\text{max},s} + 86.7V_{\text{max},s}^2 - 1343V_{\text{max},s}^3)$ where $V_{\text{max},s}$ is expressed in volts and $\phi_s$ is the phase of the pulse. For a typical value $V_{\text{max},s} = 20$ mV, the second-order term is a correction of about 2% and the third one is a correction of about 0.07%.

The duration $\Delta\tau$ is calibrated using qubit-state revival during Ramsey interferometry (see the Supplemental Material of Ref. [35]). We use a Ramsey interferometry sequence [Fig. 12(a)] for the yes-no qubit at its resonance frequency for various coherent states in the storage mode.

Revivals happen every $1/\chi_{s,\text{YN}}$, which allows us to set $\Delta\tau$ as half the revival time in Fig. 12(b). The signal difference between the final $-\pi/2$ and $+\pi/2$ pulses can be expressed as

$$\begin{aligned} S(t) &= \frac{\overline{\sigma_{z,\text{YN}}^+} - \overline{\sigma_{z,\text{YN}}^-}}{2} \\ &= e^{|\alpha|^2[\cos(2\pi\chi_{s,\text{YN}}t) - 1]} \\ &\quad \times \cos[|\alpha|^2 \sin(2\pi\chi_{s,\text{YN}}t)] e^{-t\Gamma_{2,\text{YN}} - \gamma|\alpha|^2 t}. \end{aligned} \quad (B2)$$

This expression is derived in the Supplemental Material of Ref. [35]. The last exponential decay factor is added to take into account the intrinsic decoherence of the yes-no qubit and the measurement-induced dephasing rate of the storage mode on the yes-no qubit during the waiting time. We also take into account a second-order Kerr correction that shifts the revival time with the amplitude of the coherent state [35]. At first order, this shift is given by

$$t_{\text{revival}} = 2\Delta\tau(1 + 2|\alpha|^2\chi_{s,s,\text{YN}}\Delta\tau). \quad (B3)$$

By adjusting the above parameters to allow the model to match the measured signal shown in Fig. 12(b), we find $\Delta\tau = 337$ ns, $\gamma = 0.23$ µs$^{-1}$, and $\chi_{s,s,\text{YN}} = 14$ kHz. However, this simple expression does not take into account the finite lifetime of the storage mode, and we prefer not to consider these values as accurate enough compared to what we obtain with the other methods presented in this work.

### 3. Measuring the mean quadratures of the resonator field

The two mean quadratures of the storage mode are computed from the Wigner tomography as follows. For any operator $\hat{O}$, one can apply the Wigner transform to obtain the operator Wigner map $W_{\hat{O}}$ [34] as

<figure><img src="image_17.png" /><figcaption>FIG. 12. Revival of the Ramsey interferometry on the yes-no qubit. (a) Circuit diagram for Ramsey interferometry in the presence of storage photons. After a 100-ns displacement pulse at the storage frequency, an unconditional $\pi/2$ pulse is applied to the yes-no qubit. We then let the qubit evolve freely during a time $t$ before doing a new $\pm\pi/2$ pulse and measure the state of the yes-no qubit. The signal $S(t)$ is half the difference between the average outcomes of the two sequences. (b) Measured (dots) and predicted (lines) signal $S$ as a function of waiting time $t$. Predicted signal is computed from Eq. (B2). Yes-no qubit revivals occur every $1/\chi_{s,\text{YN}} \approx 0.7$ µs.</figcaption></figure>
$$ 
\begin{aligned}
W_{\hat{\partial}}(\alpha = x + ip) &= W_{\hat{\partial}}(x, p) \\
&= \frac{1}{\pi} \int dy e^{-2ipy} \langle x + y/2 | \hat{\partial} | x - y/2 \rangle \\
&= \frac{2}{\pi} \text{Tr}[D^{\dagger}(\alpha) \hat{\partial} D(\alpha) \mathcal{P}],
\end{aligned} 
\quad (B4) $$

where {|x⟩} is the eigenbasis of the quadrature operator $\hat{X}$. With this tool, the Wigner function of a state |Ψ⟩ (respectively, a density matrix $\rho$) is simply given by $W_{|\Psi\rangle\langle\Psi|}(\alpha)$ [respectively, $W_{\rho}(\alpha)$]. The mean value of an operator $\hat{O}$ can be derived from the integral over the phase space of the product of the two Wigner distributions multiplied by $\pi$,

$$ 
\begin{aligned}
\pi \int dx \int dp W_{\rho}(x, p) W_{\hat{\partial}}(x, p) &= \frac{1}{\pi} \int dx \int dp \int dy \int dy' e^{-2ip(y+y')} \langle x+y/2|\rho|x-y/2\rangle \langle x+y'/2|\hat{O}|x-y'/2\rangle \\
&= \int dx \int dy \int dy' \delta(y+y') \langle x+y/2|\rho|x-y/2\rangle \langle x+y'/2|\hat{O}|x-y/2\rangle \\
&= \int dx \int dy \langle x+y/2|\rho|x-y/2\rangle \langle x-y/2|\hat{O}|x+y/2\rangle \\
&= \int du \int dv \langle u|\rho|v\rangle \langle v|\hat{O}|u\rangle \\
&= \text{Tr}(\rho \hat{O}) = \langle \hat{O} \rangle_{\rho}.
\end{aligned} 
\quad (B5) $$

In the case of $\hat{X}$ and $\hat{P}$ operators, Wigner maps take a simple expression

$$ 
\begin{aligned}
W_{\hat{X}}(\alpha = x + ip) &= x/\pi, \\
W_{\hat{P}}(\alpha = x + ip) &= p/\pi.
\end{aligned} 
\quad (B6) $$

For any density matrix $\rho$, one can extract $\langle \hat{X} \rangle = \text{Tr}(\hat{X}\rho)$ and $\langle \hat{P} \rangle = \text{Tr}(\hat{P}\rho)$ from the Wigner function $W \equiv W_{\rho}$ as

$$ 
\begin{aligned}
\langle \hat{X} \rangle &= \int dx \int dp W(x, p)x, \\
\langle \hat{P} \rangle &= \int dx \int dp W(x, p)p.
\end{aligned} 
\quad (B7) $$

### 4. Measuring the decoherence of the storage mode

For a qubit, Ramsey oscillations correspond to the evolution of the real part of the coherence between the $|g\rangle$ and $|e\rangle$ states. A typical sequence starts by a $\pi/2$ pulse detuned from resonance by $\delta f$ to create a coherent superposition of $|g\rangle$ and $|e\rangle$ states. Then, the qubit evolves freely before its state tomography. Both $\sigma_x$ and $\sigma_y$ oscillate at $\delta f$ while decaying at the decoherence rate $\Gamma_2$.

We decide to realize an analogous sequence based on the same idea for a harmonic oscillator (a similar sequence was recently performed in Ref. [32]). The first $\pi/2$ pulse is replaced by a detuned displacement pulse $D(\beta)$ on the storage mode. The field then evolves during a time $t$ (during which the multiplexing measurement can be applied) before a Wigner tomography is realized [see Fig. 13(a)].

The expectation value of $\hat{X} = (\hat{a}_s + \hat{a}_s^\dagger)/2$ and $\hat{P} = (\hat{a}_s - \hat{a}_s^\dagger)/2i$ quadratures are computed from the Wigner tomography (see Appendix B 3). The time trace of $\langle \hat{X} \rangle$ and $\langle \hat{P} \rangle$ is what we call the Ramsey oscillations for the storage mode. As in the qubit case, the frequency of the oscillations is set by the detuning $\delta f_s$ between the drive and the resonant frequency of the mode, which allows us to extract the frequency of the storage mode. At this point, a distinction has to be made between the detuning $\delta f_s^0 = f_{\text{drive}} - f_s$ between the drive and the bare storage-mode frequency (the resonant frequency when the multiplexed qubit and the storage are undriven) and the detuning $\delta f_s$ between the drive and the resonant frequency of the storage mode, which depends on the multiplexed measurement strength in perfect analogy with the ac-Stark effect for a qubit readout. Note that the Wigner tomography sequence uses the same detuned frequency $\delta f_s$ for its displacement pulse $D^\dagger(\alpha)$ in order to keep the same phase reference. The measurement of Ramsey oscillations of a harmonic oscillator takes longer than the ones of a qubit because we fully determine the quantum state of an oscillator at each time instead of a simple Bloch vector. From Eq. (A2), one finds that $\langle \hat{X} \rangle$ and $\langle \hat{P} \rangle$ evolve as

$$ 
\begin{aligned}
\langle \hat{X} \rangle &= |\beta| \cos(2\pi\delta f_s t + \phi) e^{-i\Gamma_{d,s}}, \\
\langle \hat{P} \rangle &= |\beta| \sin(2\pi\delta f_s t + \phi) e^{-i\Gamma_{d,s}},
\end{aligned} 
\quad (B8) $$

where $\beta = |\beta|e^{i\phi} = \langle \hat{a}_s \rangle (t = 0)$. For each time $t$, we compute $\langle \hat{X} \rangle$ and $\langle \hat{P} \rangle$ and define the storage-mode dephasing rate as $\Gamma_{d,s}$ which contains the intrinsic decoherence rate $\Gamma_{2,s}$. Data-extended Fig. 13(b) shows an example of
<figure><img src="image_18.png" /><figcaption>FIG. 13. Ramsey oscillations of the storage mode. (a) Circuit diagram for Ramsey oscillations of a harmonic oscillator. All storage displacement pulses are performed in 100 ns with a Gaussian envelope of 25 ns width. In this experiment, the amplitude of the prepared coherent state $\beta$ is set to -1.55. The detuning between displacement pulse and bare storage frequencies is $\delta f_s^0 = 3.96$ MHz. (b) Measured (dots) and expected (lines) signals for $\langle \hat{X} \rangle$ (blue) and $\langle \hat{P} \rangle$ (orange). The expected signals are matched to the experiment using Eq. (B8) with a frequency detuning of $\delta f_s = 3.96$ MHz and a decay rate $\Gamma_{2,s} = (2 \mu s)^{-1}$.</figcaption></figure>

measured Ramsey oscillations. In the main text, Fig. 3(b) does not exhibit oscillations because it is the mean value $\langle \hat{a}_s + \hat{a}_s^\dagger \rangle$ in the frame rotating at the resonant frequency of the storage mode. In practice, we measure them with a detuning and numerically compute the nonoscillating quantity $2\text{Re}[\langle \hat{X} \rangle + i\langle \hat{P} \rangle] \exp(-2i\pi\delta f_s t)$.

## 5. Storage-mode frequency shift and induced dephasing rate by driving the multiplexing qubit with a comb

In analogy with the ac-Stark shift of the frequency of a qubit coupled to a driven resonator, we also call ac-Stark shift the frequency shift of the storage mode induced by driving the multiplexing qubit. In order to measure this frequency shift and the dephasing rate that is induced by the multiplexing qubit on the storage mode, we realize the reciprocal protocol for a qubit measured by a cavity. We use a Ramsey interferometry sequence on the storage mode during which the multiplexing qubit is driven with a frequency comb [see Fig. 14(a)]. The drive amplitude is given by the sum of nine sine functions at the frequencies $[f_{\text{MP}}, f_{\text{MP}} - \chi_{\text{s,MP}}, \dots, f_{\text{MP}} - 8\chi_{\text{s,MP}}]$ multiplied by a Gaussian envelope of duration $t$ and width $t/4$.

For small measurement strength $\Omega/\chi_{\text{s,MP}} < 0.9$, we generate the Ramsey sequence with a displacement pulse of amplitude $\beta = -1.55$ detuned from the base storage mode by $\delta f_s^0 = 3.96$ MHz. We fit the time evolution of $\langle \hat{X} \rangle$ and $\langle \hat{P} \rangle$ using the damped sine function

$$ \begin{aligned} \langle \hat{X} \rangle &= A \cos(2\pi\delta f_s t + \phi) e^{-t\Gamma_{d,s}}, \\ \langle \hat{P} \rangle &= A \sin(2\pi\delta f_s t + \phi) e^{-t\Gamma_{d,s}}. \end{aligned} \quad (B9) $$

The parameters $A$, $\delta f_s$, $\phi$, and $\Gamma_{d,s}$ are determined all together by fitting the model to the measured oscillations. $\Gamma_{d,s}$ is the sum of the intrinsic storage dephasing rate $\Gamma_{2,s}$ and of the measurement-induced dephasing rate. $\delta f_s$ is the sum of the detuning from the bare-storage-mode frequency $\delta f_s^0$ and of the ac-Stark shift of the storage mode. Both parameters are shown in Figs. 14(c) and 14(d) (blue dots) as a function of $\Omega/\chi_{\text{s,MP}}$.

For larger measurement strength $\Omega/\chi_{\text{s,MP}} > 0.9$, we generate the Ramsey sequence with a displacement pulse detuning of $\delta f_s^0 = 5.96$ MHz, an amplitude of $\beta = -1.27$, and we model the time evolution of $\langle \hat{X} \rangle$ and $\langle \hat{P} \rangle$ by the sum of two sine functions with an exponential decay

$$ \begin{aligned} \langle \hat{X} \rangle &= A(\cos(2\pi\delta f_s t + \phi) + \zeta \cos(2\pi\nu t + \psi_X))e^{-t\Gamma_{d,s}}, \\ \langle \hat{P} \rangle &= A(\sin(2\pi\delta f_s t + \phi) + \zeta \sin(2\pi\nu t + \psi_P))e^{-t\Gamma_{d,s}}. \end{aligned} \quad (B10) $$

This empirical model originates from three ideas. The first term is identical to the simple model in Eq. (B9). Second, the measured Ramsey oscillations seem to show a small modulation in amplitude, which we try to capture with a second sine function. Third, we try to keep the model as simple as possible.

Figure 14(b) shows an example of Ramsey oscillations of the storage mode with a large amplitude of measurement. The two signals are fitted simultaneously to extract the parameters $A$, $\delta f_s$, $\nu$, $\phi$, $\psi_X$, $\psi_P$, and $\Gamma_{d,s}$. The frequency $\nu$ varies from 2.15 to 2.5 MHz. The parameter $\zeta$ is roughly constant; it varies between 0.2 to 0.27. Figure 14(c) shows measurement-induced detuning as a function of the measurement drive amplitude.

## 6. Rabi frequency calibration

We observe Rabi oscillations of the multiplexing qubit by applying a 1-μs-long square pulse at $f_{\text{MP}}$ with a varying amplitude $V_{\text{MP}}$. The reflected signal demodulated by time steps of 10 ns displays damped oscillations given by [40]

$$ \begin{aligned} \text{Re}(r(t)) - \text{Re}(r_{\text{SS}}) &= A \cos \left[ \sqrt{(2\pi\xi V_{\text{MP}})^2 - \left(\frac{\Gamma_{1,\text{MP}} - 2\Gamma_{\phi,\text{MP}}}{16}\right)^2} (t - t_0) + \phi \right] \\ & \times e^{-(t-t_0)/T}, \end{aligned} \quad (B11) $$
<figure><img src="image_19.png" /></figure>

FIG. 14. Frequency shift and dephasing rate of the storage mode induced by the multiplexed photocounting measurement. (a) Circuit diagram of the protocol used to determine the dephasing rate and frequency shift of the storage mode induced by the multiplexed photocounting measurement. The amplitude $\beta$ of the initial displacement is set at -1.55 for small measurement amplitude $\Omega/\chi_{s,MP} < 0.9$ and at -1.27 for large measurement amplitude $\Omega/\chi_{s,MP} > 0.9$. The blocks linking the multiplexing qubit and the storage mode represent the multiplexed measurement during a time $t$ made by the qubit on the storage mode. This measurement is realized by driving the qubit with a frequency comb [$f_{MP}, f_{MP}-\chi_{s,MP}, \dots, f_{MP}-8\chi_{s,MP}$] within a Gaussian envelope. (b) Ramsey oscillations of the storage mode for “large” measurement amplitude $\Omega/\chi_{s,MP} = 1$. One can observe that the dynamics $\langle\hat{X}\rangle$ and $\langle\hat{P}\rangle$ are not governed by a simple decaying sine function. The theory does not reproduce quantitatively the measurement when using the naive version of the model Eq. (B9). We use the simple model Eq. (B10) to capture this modulation. (c), (d) ac-Stark shift and measurement-induced dephasing rate measured (dots) and simulated (line) as a function of the multiplexing qubit drive amplitude $\Omega$ in units of $\chi_{s,MP}$. The evolution of the detuning and dephasing rate is strongly nonlinear with the drive amplitude.

where $r_{SS}$ is the value of the reflection coefficient in the steady state. We obtain $\xi = 0.543$ GHz $V^{-1}$ so that the Rabi frequency is calibrated as $\Omega = \xi V_{MP} = (0.543 \text{ GHz V}^{-1})V_{MP}$ (see Fig. 15).

<figure><img src="image_20.png" /></figure>

FIG. 15. Multiplexing qubit Rabi oscillations for various driving amplitudes. The measured Rabi oscillations observed in the reflection coefficient (dots) are reproduced by theory [solid line from Eq. (B11)]. The vertical axis represents the deviation of the real part of the reflection coefficient to its steady-state value. This calibration allows us to extract the scaling parameter $\xi$ such that $\Omega = \xi V_{MP} = (0.543 \text{ GHz V}^{-1})V_{MP}$.

## APPENDIX C: MAXIMAL MEASUREMENT STRENGTH

The dephasing rate $\Gamma_{d,s}$ depends on the driving amplitude $\Omega$, the dispersive shift $\chi_{s,MP}$, and the relaxation rate $\Gamma_{1,MP}$. To gain insight into its dependence and learn how to maximize the measurement strength, we start by considering the case of an isolated single qubit before extending the model to a bipartite qubit-resonator system.

### 1. Dynamics of a qubit driven by a comb

#### a. Hamiltonian evolution

We consider a single qubit driven by a frequency comb with $2p + 1$ frequency peaks at every $f_{MP} + k\chi$ for $-p \le k \le p$. In the frame rotating at the qubit frequency, the Hamiltonian reads
$$H(t) = \pi\Omega\left(\sum_{k=-p}^{p} \cos(2\pi k\chi t)\right)\sigma_x. \quad (C1)$$

After a time $t$, the qubit state will thus be rotated around the $x$ axis of the Bloch sphere by an angle $f(t)$ with

$$f(t) = 2\pi\Omega t \sum_{k=-p}^{p} \mathrm{sinc}(2\pi k\chi t), \quad (C2)$$

where $\mathrm{sinc}(x) = \sin(x)/(x)$. For large integers $p$, we can approximate the sum as $\sum_{k=-\infty}^{+\infty} \mathrm{sinc}(\pi kT) = (1/T)$, which is valid for $0 < T < 2$. Note that for a small number of peaks $2p + 1$, this approximation is invalid close to $T = 0$ or $2$. The expression allows us to approximate $f(t)$ for $0 < \chi t < 1$. It is then simple to derive $f(t)$ at any time $t$ since it is periodic up to the term in $k = 0$. We obtain

$$f(t) \approx \pi \frac{\Omega}{\chi} + 2\pi \frac{\Omega}{\chi} \lfloor t\chi \rfloor := \bar{f}(t), \quad (C3)$$

where $\lfloor x \rfloor$ is the integer part of $x$. Therefore, the rotation angle $f(t)$ evolves by steps. A comparison of the actual $f(t)$ and of the staircase approximation for a comb with 21 frequencies ($p = 10$) is shown in Fig. 16. To put it simply, the action of the comb consists of performing a Rabi rotation on the qubit by discrete steps instead of a continuous evolution as is the case for a single driving frequency. At each period $1/\chi$, the qubit rotates almost instantaneously by an angle $2\pi(\Omega/\chi)$.

Without decoherence, if the qubit starts in state $|g\rangle$ at time $t_0$, the qubit state after a time $t$ reads

$$\begin{aligned} |\psi(t)\rangle &= \cos\left(\frac{f(t) - f(t_0)}{2}\right)|g\rangle + i\sin\left(\frac{f(t) - f(t_0)}{2}\right)|e\rangle \\ &\approx \cos\left(\frac{\Omega\pi}{\chi}\lfloor(t-t_0)\chi\rfloor\right)|g\rangle + i\sin\left(\frac{\Omega\pi}{\chi}\lfloor(t-t_0)\chi\rfloor\right)|e\rangle. \end{aligned} \quad (C4)$$

Let us focus on some particular values of $\Omega/\chi$.

(i) If $\Omega/\chi$ is an integer, the staircase approximation with $\bar{f}(t)$ keeps the qubit in $|g\rangle$ at all times, just performing a full rotation on the Bloch sphere at each Rabi pulse. In the presence of relaxation, a photon loss can happen only during the short duration of the Rabi pulse, which decreases as $1/(p+1)$. The qubit

<figure><img src="image_21.png" /><figcaption>FIG. 16. Exact rotation angle $f(t)$ in blue and its staircase approximation $\bar{f}(t)$ in red. The duration of each step is equal to $1/\chi$ and its height to $2\pi\Omega/\chi$. The quality of the approximation improves as the number of peaks $2p + 1$ in the comb gets larger. The fact that the trajectory starts with a half-jump is a particularity of having assumed that all comb components have the same phase at $t = 0$. Random initial phases of the signals (e.g., due to initializing the qubit after its photon emission into the measuring transmission line at a random time) would most often position $t = 0$ on a flat portion of the stairs.</figcaption></figure>

remaining in the ground state cannot encode information on the resonator. This explains why there are minima of the measurement-induced dephasing rate when $\Omega/\chi$ is an integer [see Fig. 14(d)].

(ii) If $\Omega/\chi$ is a half-integer, the staircase approximation with $f = \bar{f}$ makes the qubit state jump periodically between $|g\rangle$ and $|e\rangle$. This maximal extent of the evolution on the Bloch sphere intuitively explains the maximal measurement strength at the driving amplitude. In the following, we explain why this situation corresponds to a larger extraction of information than any steps with a different rotation angle.

b. Integrated qubit dynamics in the presence of relaxation

In the following, we use the infinite-comb approximation $f(t) \approx \bar{f}(t)$. This allows us to integrate the qubit dynamics exactly. The continuous photon decay at rate $\Gamma_1$ is interrupted by discrete Rabi rotations at discrete times. The qubit state is confined in the $y-z$ plane of the Bloch sphere. Under this approximation, after the Rabi pulse number $k + 1$, the qubit state is given by

$$\begin{pmatrix} y(kT + T) \\ z(kT + T) \end{pmatrix} = \begin{pmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix} \begin{pmatrix} e^{-\Gamma_1 T/2} y(kT) \\ e^{-\Gamma_1 T} z(kT) + (e^{-\Gamma_1 T} - 1) \end{pmatrix}, \quad (C5)$$
where $T = 1/\chi$ is the period, and $\theta = (2\pi\Omega/\chi)$ is the angle spanned in the Bloch sphere during a discrete jump. The origin of time $t = 0$ is chosen to start just after a Rabi jump.

The permanent solution of this discrete-time map (Poincaré map) right after $k$ steps is

$$ \begin{aligned} \bar{y}(kT) &= \frac{e^{\Gamma_1 T/4} \sinh(\Gamma_1 T/2) \sin \theta}{\cosh(3\Gamma_1 T/4) - \cos \theta \cosh(\Gamma_1 T/4)}, \\ \bar{z}(kT) &= \frac{\sinh(\Gamma_1 T/2)(e^{-\Gamma_1 T/4} - e^{\Gamma_1 T/4} \cos \theta)}{\cosh(3\Gamma_1 T/4) - \cos \theta \cosh(\Gamma_1 T/4)}. \end{aligned} \quad (C6) $$

The average qubit excited-state population over one period in the permanent regime reads

$$ P(\Omega/\chi) := \frac{1}{T} \int_{0}^{T} \langle e | \rho(t) | e \rangle = \frac{1 - e^{-\Gamma_{1}T} [1 + \bar{z}(kT)]}{\Gamma_{1}T} \quad (C7) $$

When using the parametrization $c \equiv \cos(\theta)$, one easily checks that $P$ is a strictly decreasing function on $c \in [-1, 1]$. As a function of $\theta$, it has maxima for $\theta = (2k+1)\pi$ and minima for $\theta = 2k\pi$. The positions of the minima are no surprise and give $P=0$ as the Rabi pulse takes the state from $|g\rangle$ back to $|g\rangle$; for a finite number of peaks in the comb $2p+1$, the Rabi pulse is not instantaneous

and $P > 0$ at these minima. The value of the maxima would be $P = [\tanh(\Gamma_1 T/2)/\Gamma_1 T]$ with the infinite-comb approximation.

According to this approximation, the average rate of photon emission, which is linked to the measurement strength (each photon reveals information about the qubit frequency and hence the photon number), is thus $P\Gamma_1 = \tanh(\Gamma_1 T/2)/T$ with the optimal choice of $\Omega = \chi/2 + k\chi$, where $k$ is an integer.

(i) At fixed $T = 1/\chi$, the emission rate increases with $\Gamma_1$ and converges toward $\chi$. For $\Gamma_1 \gg \chi$, the qubit has the time to fully relax during one period. Therefore, in simple terms, at each period in the stepwise evolution, the qubit is excited and then releases deterministically a single photon into the output transmission line.

(ii) Likewise, for a fixed $\Gamma_1/\chi$ (thus, fixed probability $P$ to emit a photon during a period), the average emission rate increases when $T$ decreases. Therefore, the average emission rate increases as $\chi$.

(iii) For a fixed $\Gamma_1$, the largest average emission rate is obtained for $\chi = 1/T$ as large as possible, but it saturates at $P\Gamma_1 = \Gamma_1/2$. This behavior is consistent with the fact that $\Gamma_1$ is a hard limit on the photon emission rate.

## 2. Dynamics of the qubit-resonator bipartite system driven by a comb

For the sake of simplicity, we consider a lossless storage mode in this section. In this case, each $2 \times 2$ submatrix $\rho_{n_1,n_2} = \langle n_1 | \rho | n_2 \rangle$, with $|n_i\rangle$ the $n_i$ Fock state, evolves independently of the others $\rho_{m_1,m_2}$ similar to a collection of qubitlike systems. The submatrix $\rho_{n_1,n_2}$ is non-normalized because it is an off-diagonal submatrix of the storage-qubit system. We have

$$ \begin{aligned} \frac{d}{dt} \rho_{n_1, n_2} = & -i2\pi\chi \frac{n_1 + n_2}{2} \left[ \frac{\sigma_z}{2}, \rho_{n_1, n_2} \right] \\ & - i2\pi\chi \frac{n_1 - n_2}{2} \left\{ \frac{\sigma_z}{2}, \rho_{n_1, n_2} \right\} \\ & - \frac{i}{\hbar} [H(t), \rho_{n_1, n_2}] + \Gamma_1 \mathcal{L}(\sigma_-) \rho_{n_1, n_2}, \end{aligned} \quad (C8) $$

where $H(t)$ is the frequency comb drive from the previous section.

### a. Computing the decoherence rate of the $n_1, n_2$ component

The infinite-comb approximation again helps. We view $H(t)$ as applying a Rabi pulse of angle $\theta = (2\pi\Omega/\chi)$ at each period $T = 1/\chi$, without any effect for the rest of the time. Over one period, we thus have

$$ \rho_{n_1, n_2}(kT + T) = K_0 \circ K_1 \rho_{n_1, n_2}(kT), \quad (C9) $$

where $K_0$ applies the Rabi pulse, while $K_1$ contains dynamics associated with the dispersive coupling and with qubit decay. During each period between Rabi jumps, denoting $\rho_{n_1,n_2} = (x\sigma_x + y\sigma_y + z\sigma_z + \eta I/2)$, the dynamics $K_1$ corresponds to the integration of the set of equations

$$ \begin{aligned} \frac{d}{dt}x &= -\frac{\Gamma_1}{2}x - \frac{2\pi\chi(n_1 + n_2)}{2}y, \\ \frac{d}{dt}y &= -\frac{\Gamma_1}{2}y + \frac{2\pi\chi(n_1 + n_2)}{2}x, \\ \frac{d}{dt}z &= -\Gamma_1(z + \eta) - i\frac{2\pi\chi(n_1 - n_2)}{2}\eta, \\ \frac{d}{dt}\eta &= -i\frac{2\pi\chi(n_1 - n_2)}{2}z. \end{aligned} $$

After one period $T$, since the peaks in the comb are exactly separated by the dispersive shift $\chi$, the effect of the precession at a frequency $2\pi\chi(n_1 + n_2)/2$ is canceled out (modulo a possible phase flip every period when $n_1 + n_2$ is odd). Note that the infinite-comb approximation differs from the usual rotating-wave approximation that would lead to a similar disabling of the precession for $2\pi\chi \gg \Gamma_1$. We then obtain, in the above coordinates, the matrix expression
$$ K_1 = (-1)^{n_1+n_2} \times \begin{pmatrix} e^{-\Gamma_1 T/2} & 0 & 0 & 0 \\ 0 & e^{-\Gamma_1 T/2} & 0 & 0 \\ 0 & 0 & (e^{-\Gamma_1 T} - G) & (e^{-\Gamma_1 T} - 1 - G) \\ 0 & 0 & G & (G+1) \end{pmatrix} \quad (C10) $$

with $G=[-i\pi\chi(n_1-n_2)(1-e^{-\Gamma_1 T})]/[\Gamma_1+i2\pi\chi(n_1-n_2)]$. In addition, the Rabi rotation corresponds to

$$ K_0 = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & \cos\theta & -\sin\theta & 0 \\ 0 & \sin\theta & \cos\theta & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}. \quad (C11) $$

This expression further simplifies for the two drive strengths $\Omega = \chi/2 + k\chi$ or $\Omega = k\chi$, which, respectively, lead to $\theta = 0$ (no qubit emission) and $\theta = \pi$ (maximal qubit emission), since the $(z, \eta)$ variables of interest decouple from $(x, y)$.

(i) For $\theta = 0$, we tend to a stationary regime $(z, \eta)_{kT+T} = (z, \eta)_{kT}$. Note that these are the values associated with a coherence between Fock states $n_1$ and $n_2$ of the storage mode. This steady value thus confirms the intuition developed in the single-qubit case: There is no change in the coherences between Fock states in the resonator, which means that no measurement is performed [minima in $\Gamma_{d,s}$ in Fig. 14(d)].

(ii) For $\theta = \pi$, we can compute an analytical expression for the factor $R$ by which the trace $\eta$ decreases every period $T$ in the permanent regime [41]. Thus, the average decay rate is given by $-\log(|R|)/T$ with

$$ R = \frac{1}{2[\Gamma_1 + i2\pi\chi(n_1 - n_2)]} \left[ \Gamma_1(1 - e^{-\Gamma_1 T}) + 2\sqrt{B} \left( \sqrt{1 + A^2} + A + i\sqrt{1 + A^2} - A \right) \right], $$ 
 where $A = \frac{\Gamma_1^2(1 + e^{-\Gamma_1 T})^2 - 4(2\pi\chi)^2(n_1 - n_2)^2 e^{-\Gamma_1 T}}{8B}$, 
 $B = 2\pi\chi(n_1 - n_2)\Gamma_1 e^{-\Gamma_1 T}$,

and we recall $T = 1/\chi$.

(C12)

We now analyze this last expression.

### b. Optimal decoherence rate of the $n_1, n_2$ component

The decoherence rate $-\log(|R|)\chi$ can be calculated for any value of the dispersive shift $\chi$ and of the photon numbers $n_1$ and $n_2$. In Fig. 17 are shown the rates corresponding to several values of $n_1 - n_2$ for a driving amplitude $\Omega = \chi/2$ (maximal measurement strength). As expected, the decoherence is stronger when the photon numbers are farther apart. In addition, for $2\pi\chi|n_1 - n_2| \gg \Gamma_1$ the rate saturates to $\Gamma_1/2$ similar to the emission rate of the single qubit.

### 3. Upper bound on the total measurement rate as a function of the photon number

Beyond the determination of the decoherence rate between two Fock states, we are interested in the maximal information extraction rate of the storage state in the

multiplexed measurement scheme. In particular, we discuss how this maximal total measurement rate depends on the maximum number of photons $N_{max}$ that are probed by the

<figure><img src="image_22.png" /><figcaption>FIG. 17. Decoherence rate between Fock states $|n_1\rangle$ and $|n_2\rangle$ in units of $\Gamma_1$ for a driving amplitude $\Omega = \chi/2$ as a function of $2\pi\chi/\Gamma_1$ and for various values of $n_1 - n_2$.</figcaption></figure>
multiplexed scheme. In the following, we assume a perfect measurement apparatus giving us access to all the information extracted by the measurement process, which is not necessarily the case of the heterodyne measurement on each peak of the comb. In Sec. IV, we propose such a measurement apparatus.

We assume the number-resolved regime $2\pi\chi \gg \Gamma_2, \Gamma_1$. Thus, the decoherence rate between two Fock states is independent of the Fock state numbers and is equal to $\Gamma_1/2$ (see Fig. 17). In the following, we show that under these assumptions, the total measurement rate does not depend on $N_{max}$.

Since the multiplexed measurement operates by entangling the storage mode with $N_{max} + 1$ frequency modes of the transmission line, we can describe the system and the extraction of information without the multiplexing qubit and consider only its effect, which is the entanglement operation. Each Fock state $|n\rangle$ of the storage mode ($0 \le n \le N_{max}$) is associated to one of the $N_{max} + 1$ modes of the transmission line at $f_{mp} - n\chi$. Every mode is driven so that at the input it is in a coherent state (it can even be the vacuum as in the gedanken experiment). At the output, if we change the reference frame by displacing the outgoing modes $\hat{a}_{out,n}$ by the opposite of the input coherent state, a single mode will be excited, and all the nonresonant modes will be in the vacuum state. Therefore, any quantum state of the outgoing modes can be expressed as a superposition of $N_{max} + 2$ states only. States $|n\rangle_m$ correspond to all modes in the vacuum except the one at frequency $f_{MP} - n\chi$ and $|\perp\rangle_m$ is the vacuum state.

Thus, we can describe the system using two modes only: the storage mode and a simplified measurement mode. The storage mode is described using the Fock state basis $\{|n\rangle_s\}_{0\le n \le N_{max}}$. The measurement mode has the $N_{max} + 2$ states discussed above. The bipartite system starts in the state

$$|\Psi(0)\rangle = |\Psi_{\text{storage}}\rangle_s \otimes |\Psi_{\text{meas}}\rangle_m = \left(\sum_{n=0}^{N_{\text{max}}} \psi_n |n\rangle_s\right) \otimes |\perp\rangle_m. \qquad (C13)$$

After a measurement time $t$, the storage mode and the measurement mode become entangled, and the readout of the measurement mode extracts information about the storage photon number. As the decoherence rate between every storage Fock state pair is $\Gamma_1/2$, one can write the state of the bipartite system as

$$|\Psi(t)\rangle = \sqrt{e^{-\Gamma_1 t/2}} |\Psi(0)\rangle + \sqrt{1 - e^{-\Gamma_1 t/2}} \sum_{n=0}^{N_{\text{max}}} \psi_n |n\rangle_s \otimes |n\rangle_m. \qquad (C14)$$

As expected, if we trace over the measurement mode, the diagonal of the density matrix of the storage mode remains unchanged, while the off-diagonal terms decrease at a rate $\Gamma_1/2$.

One can calculate the mutual information [42] of the bipartite system in order to obtain the information extraction rate and study how it scales with $N_{max}$. The mutual information is defined as

$$I(s, m) = S(\rho_s) + S(\rho_m) - S(\rho_{s,m}), \qquad (C15)$$

where $S(\rho) = -\text{Tr}[\rho \log_2(\rho)]$ is the von Neumann entropy, $\rho_{s,m}$ is the bipartite density matrix, and $\rho_s$ (respectively, $\rho_m$) is the density matrix of the storage mode (respectively, measurement mode) obtained by tracing out on the other mode. If the state of the bipartite system is pure, then its entropy is zero, and the entropy of the storage mode is equal to the entropy of the measurement mode. Thus, the mutual information can be written as

$$I(s, m) = 2S(\rho_s). \qquad (C16)$$

The mutual information is twice the amount of information that we can acquire about the storage mode using the measurement mode. Thus, it is twice the amount of information extracted by the measurement process. The density matrix $\rho_s$ of the storage mode depends on the initial photon number distribution. As our goal is to measure the probability to have $n$ photons for all $n$ simultaneously, we consider an initially uniform photon number distribution, such as when $\psi_n = (N_{max} + 1)^{-1/2}$ for all $n$. For this photon number distribution, the storage density matrix becomes

$$\rho_s = \frac{e^{-\Gamma_1 t/2}}{N_{\text{max}} + 1} \sum_{n,l=0}^{N_{\text{max}}} |n\rangle_s \langle l|_s + \frac{1 - e^{-\Gamma_1 t/2}}{N_{\text{max}} + 1} \sum_{n=0}^{N_{\text{max}}} |n\rangle_s \langle n|_s. \qquad (C17)$$

The eigenvalues of $\rho_s$ are $e^{-\Gamma_1 t/2} + (1 - e^{-\Gamma_1 t/2})/(N_{max} + 1)$ with degeneracy 1 and $(1 - e^{-\Gamma_1 t/2})/(N_{max} + 1)$ with degeneracy $N_{max}$. Thus, one can derive the mutual information of the bipartite system

$$I(s, m) = -2 \left( r + \frac{1-r}{N_{\text{max}}+1} \right) \log_2 \left( r + \frac{1-r}{N_{\text{max}}+1} \right) \\ - 2 \frac{N_{\text{max}}(1-r)}{N_{\text{max}}+1} \log_2 \left( \frac{1-r}{N_{\text{max}}+1} \right), \qquad (C18)$$

with $r = e^{-\Gamma_1 t/2}$. The time evolution of the mutual information $I(s, m)$ is shown in Fig. 18. At short times, the mutual information increases with time $t$ at a speed which depends on $N_{max} + 1$. As the mutual information is the information we extract out of the system, the time
<figure><img src="image_23.png" /><figcaption>FIG. 18. Mutual information between the storage and measurement modes for various maximum photon number $N_{max} = 2^{n_b} - 1$ as a function of the measurement time. (a) Mutual information in bits. (b) Mutual information divided by the number of bits $n_b$ as a function of time. The dashed black curve is the lower bound.</figcaption></figure>

derivative of the mutual information at short times gives the rate at which information is extracted (e.g., the total measurement rate). In order to determine how the total measurement rate depends on the maximum photon number $N_{max}$, we look at the mutual information per bit of information $n_b = \log_2(N_{max} + 1)$. The mutual information per bit decreases with $n_b$ for small photon numbers but converges to a lower bound when $n_b$ goes to infinity

$$ \lim_{n_b \to +\infty} \frac{I(s, m)}{n_b} = 2(1 - e^{-\Gamma_1 t/2}). \quad (C19) $$

The total measurement time is of the order of the number of bits $n_b$ we have to measure divided by the total measurement rate in bits per second. As the mutual information per bit is always bigger than $2(1 - r)$, the total measurement time will be at least about $1/\Gamma_1$. Because the mutual information per bit is always bigger than the lower bound $2(1 - r)$, which does not depend on $n_b$, the total measurement time for the multiplexing protocol is independent of the maximum photon number $N_{max}$. Figure 18(b) shows the mutual information per bit as a function of the time and its lower bound.

## APPENDIX D: MEASUREMENT BANDWIDTH

The asymptotic former results are possible as long as the total measurement bandwidth increases with the photon

number. We consider two main limitations to the maximal number of photons that can be measured. First, above a certain number of photons, the qubit frequency overlaps the transition between the first and second excited states of the transmon. Indeed, the $|e\rangle$-to-$|f\rangle$ transition at zero photon in the resonator becomes resonant with the qubit $|g\rangle$-to-$|e\rangle$ transition for $\chi_{s,MP}/\chi_{s,MP}$ photons, which complicates the analysis. In our experiment, it would set a limit of about 20 photons. However, this limit can be bypassed using a qubit with a much larger anharmonicity, such as a fluxonium qubit.

Second, the higher-order terms in the Hamiltonian tend to reduce the cross-Kerr rate $\chi_{s,MP}$ when the photon number increases. Beyond some critical photon number $n_{crit}$, the dispersive shift $\chi_{s,MP}$ will be smaller than the decoherence rate of the qubit, which will escape from the number-resolved regime. This limit occurs around 1000 photons in our device, which is therefore not the dominant limitation. Indeed, one can show that

$$ \chi_{s,MP} = p_{MP} p_s \frac{h f_{MP} f_s}{4 E_J}, \quad (D1) $$

$$ \chi_{s,s,MP} = p_{MP} p_s^2 \frac{h^2 f_{MP} f_s^2}{32 E_J^2}, \quad (D2) $$

where $\chi_{s,s,MP}$ is the decrease of the cross-Kerr rate $\chi_{s,MP}$ when a photon is added to the storage mode, $E_J$ is the junction energy of the multiplexing qubit, and $p_{MP/s}$ is the fraction of the energy of the multiplexing qubit mode (resp., the storage mode) stored in the junction of the multiplexing qubit [43]. Thus, the critical photon number $n_{crit}$ in the storage mode is given by

$$ \chi_{s,MP} - n_{crit}\chi_{s,s,MP} = \Gamma_{2,MP}, \quad (D3) $$

which leads to

$$ n_{crit} = \left(1 - \frac{\Gamma_{2,MP}}{\chi_{s,MP}}\right) \frac{8E_J}{p_s h f_s}, \quad (D4) $$

and a upper bound of $8E_J/(p_s h f_s) \approx 10^3$ photons for $E_J/h \approx 15$ GHz, $p_s \approx 10^{-2}$, and $f_s \approx 5$ GHz.

Therefore, we do not believe that these limitations are hard enough to prevent the multiplexing scheme from counting very large photon numbers.

## APPENDIX E: MODELING OF THE MEASUREMENT OPERATORS

We introduce here a simple model to characterize the measurement and its backaction on the resonator. The measurement uses a phase-preserving amplifier in order to amplify the signal at each frequency $f_{MP} - k\chi_{s,MP}$ in the comb and record a complex amplitude $I_m^{(k)} + iQ_m^{(k)}$. We
assume that each of these measurement records extracts only the information on the occupation of the Fock state $|k\rangle$, which is experimentally valid in the limit $2\pi\chi_{s,MP} \gg \Gamma_{2,MP}$. Without decoherence and in the limit of long measurement time, its backaction on the storage mode would project the storage state on Fock state $|k\rangle$ or on the complementary subspace $\Pi_{\perp}^{(k)}\mathcal{H}_s$, where $\Pi_{\perp}^{(k)} = \mathbf{1} - |k\rangle\langle k|$ and $\mathcal{H}_s$ is the Hilbert space of the storage resonator.

In practice, the measurement proceeds by first entangling the resonator, which is in a state $|\psi\rangle$, and the signal mode of the phase-preserving amplifier. For the measurement channel $k$, the entangled state reads $|\alpha, 0\rangle \otimes \Pi^{(k)}|\psi\rangle + |\alpha_{\perp}, 0\rangle \otimes \Pi_{\perp}^{(k)}|\psi\rangle$, where $\Pi^{(k)} = |k\rangle\langle k|$, and states denoted as $|\alpha, \beta\rangle$ are the coherent states of the signal and the idler modes at the input of the amplifier. We distinguish two cases: the case where the probe is resonant with the multiplexing qubit, leading to a reflected amplitude $\alpha$, and the case where it is off resonant leading to a reflected amplitude $\alpha_{\perp}$. The resonance frequency of the qubit depends on the number of photons in the resonator so that the reflected amplitude $\alpha$ indicates $k$ photons, while $\alpha_{\perp}$ indicates that there are not $k$ photons. For an incoming amplitude $\alpha_{\text{in}}$ onto the multiplexing qubit, we get

$$ \alpha_{\perp} = \alpha_{\text{in}}, \\ \alpha = \alpha_{\text{in}} \left( 1 - \frac{\Gamma_{1,MP}}{\pi\Omega} \frac{\langle \sigma_{-\text{MP}} \rangle_{\text{SS}}}{i} \right), \quad (\text{E1}) $$

where $\langle \sigma_{-\text{MP}} \rangle_{\text{SS}}$ is the steady-state mean value of the multiplexing qubit-lowering operator.

If the qubit is driven by a single tone, the maximum of $|\langle \sigma_{-\text{MP}} \rangle_{\text{SS}}|$ is reached for $2\pi\Omega = \Gamma_{1,MP}/\sqrt{2}$. However, in the case of a qubit driven by an infinite-frequency comb, the time average of the lowering operator is

$$ \langle \sigma_{-\text{MP}} \rangle_{\text{comb}} = i\bar{y} \frac{1 - e^{-\Gamma_{1,MP}/2\chi_{s,MP}}}{\Gamma_{1,MP}/\chi_{s,MP}} \quad (\text{E2}) $$

with $\bar{y}$ the steady-state solution defined in Eq. (C6). The fraction on the right corresponds to the average emission between two jumps in the time-domain version of the comb.

The measurement operator (Kraus operator) corresponding to the heterodyne detection of a propagating field encoding the information on the $|k\rangle$ state thus reads

$$ M^{(k)}(I_m^{(k)}, Q_m^{(k)}) = \langle \Psi_{I_m^{(k)}, Q_m^{(k)}} | \alpha, 0 \rangle \otimes \Pi^{(k)} \\ + \langle \Psi_{I_m^{(k)}, Q_m^{(k)}} | \alpha_{\perp}, 0 \rangle \otimes \Pi_{\perp}^{(k)}, \quad (\text{E3}) $$

where $|\Psi_{I_m^{(k)}, Q_m^{(k)}}\rangle$ is the state on which the propagating field is projected after the heterodyne measurement performed

by the phase-preserving amplifier followed by a heterodyne detection setup.

Following the Supplemental Material of Ref. [44], in the case of a phase-preserving amplifier, the inner product $\xi(\beta, I_m, Q_m) = \langle \Psi_{I_m, Q_m} | \beta, 0 \rangle$ is given up to a global phase factor (independent of $\beta, I_m$, and $Q_m$) by

$$ \xi(\beta, I_m, Q_m) = \frac{1}{\sqrt{2\pi}\sigma_0} e^{-\frac{|\beta|^2}{2(2\sigma_0)^2}} \times e^{-\frac{(I_m - \beta)^2}{2(2\sigma_0)^2}} e^{-\frac{(Q_m + i\beta)^2}{2(2\sigma_0)^2}}, \quad (\text{E4}) $$

where $\sigma_0$ is the amplitude of the zero-point fluctuations (the variance of the measured $I_m$ is $2\sigma_0^2$ in the quantum limit of phase-preserving amplification). Therefore, we finally get the following analytical expression of the measurement operators for each channel $k$, in the case of $\Gamma_{2,MP} \ll 2\pi\chi_{s,MP}$:

$$ M^{(k)}(I_m^{(k)}, Q_m^{(k)}) = \frac{1}{\sqrt{2\pi}\sigma_0} e^{-\frac{|\alpha_{\perp}|^2}{2(2\sigma_0)^2}} e^{-\frac{(I_m - \alpha_{\perp})^2}{2(2\sigma_0)^2}} e^{-\frac{(Q_m + i\alpha_{\perp})^2}{2(2\sigma_0)^2}} \Pi_{\perp}^{(k)} \\ + \frac{1}{\sqrt{2\pi}\sigma_0} e^{-\frac{|\alpha|^2}{2(2\sigma_0)^2}} e^{-\frac{(I_m - \alpha)^2}{2(2\sigma_0)^2}} e^{-\frac{(Q_m + i\alpha)^2}{2(2\sigma_0)^2}} \Pi^{(k)}. \quad (\text{E5}) $$

## APPENDIX F: HOW FAST CAN THE MULTIPLEXED MEASUREMENT BE?

In principle, multiplexing the measurement enables us to determine the photon number in a constant time, no matter the maximal number of photons one wants to resolve. Note that this measurement rate is not in contradiction to the amount of information that would be contained in a qubit. Indeed, in the present scheme, we are using the two-level system not as a memory, but rather as a device whose frequency is changed by the Fock state to be measured. We are thus somehow replacing a communication channel faithfully sending a qubit state by a communication channel faithfully sending a two-level system in one out of many propagating microwave modes. In Sec. IV, we show how to obtain the photon number in a time set by the lifetime of the qubit. Here, we investigate the measurement time in the case of our experiment.

### 1. Multiplexing with quadrature measurements

The gedanken experiment requires components that are out of reach with current technologies. Yet, the experiment we perform and present in the main text implements a similar experiment that replaces the multiplexer and photo-detector array by signal processing of propagating microwave modes.

It is out of the scope of this work to derive an exact expression of the photon number measurement time in our experiment. However, we can determine whether the measurement time depends on the maximum photon
number $N_{max}$ one wants to measure. In contrast with the photodetectors, quadrature measurements are inherently noisy. Identifying the photon number in the storage mode consists of determining which channel contains an amplitude $\alpha$ while all the others contain an amplitude $\alpha_{\perp}$ (see Appendix E). The measurement records $\{I_m^{(k)}, Q_m^{(k)}\}_k$ are stochastic processes centered on $\alpha_{\perp}$ (except for one value of $k$, where it is centered on $\alpha$). Determining the photon number $k$ comes down to discriminating which record is centered on $\alpha$ using only the ensemble of noisy records $\{I_m^{(k)}, Q_m^{(k)}\}_k$.

After a measurement time $t$, the measurement records $\{I_m^{(k)}, Q_m^{(k)}\}_k$ are averaged along a duration $t$. Thus, the time-averaged intrinsic noise contained in the measurement records scales as $1/\sqrt{t}$ and the time average value is independent of $t$. The problem can be mapped onto the following game. $N_{max}$ stochastic variables $\{u_i\}_{1\le i \le N_{max}}$ are each randomly chosen using a Gaussian distribution centered on 0 with a width $1/\sqrt{t}$. Another stochastic variable $u_0$ is randomly chosen using a Gaussian distribution centered on 1 with a width $1/\sqrt{t}$. The list $\{u_i\}_{0\le i \le N_{max}}$ is scrambled randomly into a list $l$, and the goal consists of identifying the variable $u_0$ using only the list $l$. The optimal strategy is to pick the highest element of the list $l$. The probability to make an error and lose the game is then given by the probability that the maximum of the $\{u_i\}_{1\le i \le N_{max}}$ is higher than $u_0$,

$$P_{\text{error}} = P(\max_{N_{\text{max}} \ge i \ge 1}(u_i) > u_0). \quad (\text{F1})$$

We can rescale all the distribution by $\sqrt{t}$; thus, $u_0$ are chosen randomly using a Gaussian distribution centered on $\sqrt{t}$ with a width of 1 and each of the $\{u_i\}_{1\le i \le N_{max}}$ using a Gaussian distribution centered on 0 with a width of 1. One can show that the mean of the maximum of $\{u_i\}_{1\le i \le N_{max}}$ tends toward $\sqrt{2 \ln(N_{max})}$ [45],

$$\text{mean}[\max_{N_{\text{max}} \ge i \ge 1}(u_i)] \sim \sqrt{2 \ln(N_{\text{max}})}. \quad (\text{F2})$$

In addition, the median of the max of $\{u_i\}_{1\le i \le N_{max}}$ is equal to the mean value within an error scaling as $1/\sqrt{N_{max}}$:

$$\text{median}[\max_{N_{\text{max}} \ge i \ge 1}(u_i)] = \text{mean}[\max_{N_{\text{max}} \ge i \ge 1}(u_i)] + O(1/\sqrt{N_{\text{max}}}). \quad (\text{F3})$$

Since the error probability is between 1/4 and 3/4 if the median of $u_0$ is equal to the median of the maximum of $\{u_i\}_{1\le i \le N_{max}}$, it leads to

$$\begin{aligned} 1/4 < P_{\text{error}} < 3/4 &\Rightarrow \sqrt{t} \sim \sqrt{2 \log N_{\text{max}}} \\ &\Rightarrow t \sim \log(N_{\text{max}}). \end{aligned} \quad (\text{F4})$$

From this expression, we understand that the measurement time for a fixed error probability scales as $\log(N_{max})$.

## 2. Comparison among measurement schemes

In this section, we compare the following various photon number measurement schemes using a qubit of frequency $f_q$ that is dispersively coupled to a storage mode of frequency $f_s$. We assume the cross-Kerr rate $\chi$ between the storage mode and the qubit to be greater than the decoherence rate of the qubit $\Gamma_2$. The goal is to measure the photon number $N$ assuming it is smaller than $N_{max}$.

### (i) Sequential brute force

The brute-force approach consists of measuring whether or not there are $k$ photons in the storage mode for all possible values of $k$ from 0 to $N_{max}$ [2]. For each $k = 0, 1, 2, 3, \dots$, we apply a photon number conditional $\pi$ pulse to the qubit at frequency $f_q - k\chi$ so that the qubit is excited only if there are $k$ photons in the storage mode. Reading out the qubit state gives the answer to the question, “Are there $k$ photons?” The full measurement stops as soon as this binary answer is positive so that it takes a time given by $(T_{\pi} + T_{\text{RO}})(N+1)$. The time $T_{\pi}$ is the time of an conditional $\pi$ pulse; hence, it is at least $1/\chi$, while the qubit readout time $T_{\text{RO}}$ is limited by other parameters in order to get a single-shot readout.

### (ii) Passive photon number decimation using weak measurement

This approach, which was implemented with Rydberg atoms in a cavity [1], consists of encoding the photon number in the phase of the qubit by waiting a time $1/(N_{max}\chi)$ after the qubit has been prepared in state $(|g\rangle + |e\rangle)/\sqrt{2}$. The protocol is composed of a series of $p$ sequences, where each sequence encodes the photon number into the phase of the qubit and realizes a $\pi/2$ pulse on the qubit with a phase $2\pi p/N_{max}$ followed by a qubit readout. Using generalized measurement theory, one infers the probability that the cavity is in a given Fock state.

After $p$ qubit readouts, the variance on the photon number is $\sigma = N_{max}/(\sqrt{p}\pi)$ (see Appendix A in Ref. [46]). Therefore, the required number of repetitions to get a fixed error probability of the photon number scales as $p \propto N_{max}^2$. Since each measurement takes at least $1/(N_{max}\chi)$, the total measurement time scales at least as $N_{max}/\chi$.

### (iii) Active photon number decimation

The previous protocol can be improved by optimizing the phase of the final $\pi/2$ pulse to maximize
the amount of information extracted on the cavity photon number. It was realized in Ref. [3] using Rydberg atoms in a cavity. Because of the use of feedback on a weak measurement, we cannot find a closed form for the measurement time in this case [46]. However, it is shown that the total time is larger than the total time taken by a binary decimation with feedback (see below).

#### (iv) Binary decimation with feedback

This method was shown to provide the least number of steps for sequential photocounting [4]. Each step consists of applying an unconditional $\pi/2$ pulse to the qubit, wait a time $1/2^{k+1}\chi$, and apply a new unconditional $\pi/2$ pulse with a phase $\phi_k$ that encodes the least significant [47] $k$th bit $b_k$ of the photon number $N = \sum_k b_k 2^k$ into the qubit state. Importantly, the phase $\phi_k$ depends on the results of the $k-1$ former measurements. The sequence needs to be repeated $p = \log_2(N_{max} + 1)$ times with $k$ going from 0 to $p-1$. This procedure was recently implemented in Ref. [5].

The measurement time is at least given by the sum of the total interaction time between the qubit and cavity and of the total feedback latency. The total interaction time is bounded by $\sum_p 1/(2^{p+1}\chi) = 1/\chi$. However, the feedback latency scales as $p$ and can be written as $T_{FB} \log_2(N_{max} + 1)$.

#### (v) Binary decimation with optimal pulse control

An optimal binary decimation can also be implemented without using a feedback loop by measuring a series of generalized parity operators which yields the bit values of the binary decomposition of the photon number in the storage mode. The $k$th generalized parity measurement consists of an optimal pulse that excites the qubit conditioned on the value of the $k$th bit. The $p = \log_2(N_{max} + 1)$ parity measurements are performed in a time sequence. A

subsequent measurement and dynamic reset of the qubit state completes the sequence [6,7]. Such an optimal pulse can be performed only in a time of the order of the dispersive interaction time $1/\chi$. It leads to a total measurement time scaling as $\log_2(N_{max} + 1)(1/\chi + T_{reset})$ where $T_{reset}$ is the duration of the active reset protocol.

In Table II, we provide a summary of the various advantages and drawbacks of the photocounting methods. No time-sequence measurement is able to provide a measurement time that does not depend on the photon number. Using a multiplexed quantum measurement thus enables a qualitative improvement on the measurement time. In practice, this drastic improvement in the scaling with $N_{max}$ requires a detection setup that is out of reach currently. Our experiment demonstrates that multiplexing is possible using a heterodyne detection setup instead. The scaling of the measurement time is then in $\ln(N_{max})$ as in state-of-the-art sequential measurements. In addition, we depart the complexity of optimal control or feedback into the challenge to reach large measurement efficiencies on a large frequency band (many $\chi$'s).

## APPENDIX G: MASTER EQUATION SIMULATIONS

In this section, we briefly describe the master equation simulations used to understand our experimental results. We simulate the main photon counting experiments with both qubits as well as the photon number calibration of the storage mode and the dephasing rate induced by the multiplexed measurement of the storage mode.

All simulations are performed using the `PYTHON` package `OUTIP` [48]. We simulate the complete system composed of the storage mode, the yes-no qubit, and the multiplexing qubit with all couplings, except in the simulation of the measurement-induced dephasing rate for which we take into account only the storage mode and the multiplexing qubit. The storage mode is modeled as a harmonic

<table><caption>TABLE II. Protocols for photocounting using a qubit.</caption><thead><tr><th>Protocol</th><th>t<sub>meas</sub> ∝</th><th>Complexity</th><th>Check N > N<sub>max</sub></th><th>Qubit role</th><th>Error propagation</th></tr></thead><tbody><tr><td>Sequential brute force</td><td>(N + 1)/χ</td><td>N + 1 gates</td><td>Yes</td><td>Pointer</td><td>Yes</td></tr><tr><td>Passive decimation</td><td>N max /χ</td><td>N<sub>max</sub><sup>2</sup> gates and<br>complex analysis</td><td>No</td><td>Pointer</td><td>No</td></tr><tr><td>Binary code feedback</td><td>T<sub>FB</sub> log<sub>2</sub>(N<sub>max</sub> + 1)</td><td>Feedback</td><td>No</td><td>Pointer</td><td>Yes</td></tr><tr><td>Binary code<br>optimal control</td><td>(T<sub>reset</sub> + 1/χ) log<sub>2</sub>(N<sub>max</sub> + 1)</td><td>Optimal control</td><td>No</td><td>Pointer</td><td>No</td></tr><tr><td>Multiplexed<br>quadrature<br>measurement</td><td>ln(N<sub>max</sub>)/χ</td><td>Relies on near-<br>quantum-limited<br>broadband amplifier</td><td>Yes</td><td>Encoder</td><td>No</td></tr><tr><td>Gedanken<br>multiplexed<br>measurement</td><td>1/χ</td><td>Hardware to develop</td><td>Yes</td><td>Encoder</td><td>Gedanken</td></tr></tbody></table>
oscillator while the transmons are replaced by two-level systems. The Hilbert space of the storage mode is truncated at a photon number ranging from ten to 25 photons depending on the simulation. In this section, we use Pauli matrices to describe operators acting on qubits.

# 1. Photocounting simulations

## a. Photocounting with the yes-no qubit

Both photon counting approaches are simulated in a very similar manner. The first simulation (yes-no simulation) describes the use of conditional operations on the yes-no qubit. This experiment serves as a calibration of the number of photons in the storage mode and of all relevant parameters. This experiment starts with a displacement of the storage mode followed by a conditional $\pi$ pulse on the yes-no qubit at frequency $f_{YN} - \delta f_{YN}$ before detecting the expectation value of the Pauli operator $\sigma_{z,YN}$.

We write the Hamiltonian of the system in a frame rotating at $f_s - \chi_{s,MP}/2 - \chi_{s,YN}/2$ for the storage mode, $f_{YN} - \delta f_{YN}$ for the yes-no qubit mode, and $f_{MP}$ for multiplexing qubit mode as follows:

$$ \hat{H}_1/h = \delta f_{YN} \frac{\hat{\sigma}_{z,YN}}{2} - \chi_{s,YN} \hat{n}_s \frac{\hat{\sigma}_{z,YN}}{2} - \chi_{s,MP} \hat{n}_s \frac{\hat{\sigma}_{z,MP}}{2} \\ - \chi_{s,s} \hat{n}_s (\hat{n}_s - 1) - \chi_{s,s,YN} \hat{n}_s (\hat{n}_s - 1) \frac{\hat{\sigma}_{z,YN}}{2} \\ - \chi_{s,s,MP} \hat{n}_s (\hat{n}_s - 1) \frac{\hat{\sigma}_{z,MP}}{2} + \frac{\epsilon_{YN}(t)}{h} \hat{\sigma}_{x,YN} \\ + \frac{\lambda(t)}{2\pi} (\epsilon_{max} e^{i\pi(\chi_{s,MP} + \chi_{s,YN})t} \hat{a}_s \\ + \epsilon_{max}^* e^{-i\pi(\chi_{s,MP} + \chi_{s,YN})t} \hat{a}_s^\dagger), \quad (G1) $$

where $\lambda(t)$ is a Gaussian function with duration 100 ns, width 25 ns, and a maximum of 1 so that the storage-mode displacement pulse reads $\epsilon_s(t) = \lambda(t)\epsilon_{max}$, and $\epsilon_{YN}(t)$ is the time envelope of a Gaussian pulse with duration 1.9 $\mu$s and width 475 ns. The amplitude of the pulse is chosen to obtain a $\pi$ rotation on the yes-no qubit. The term $-\delta f_{YN}(\hat{\sigma}_{z,YN}/2)$ takes into account the detuning between the $\pi$ pulse and the yes-no qubit frequency. $\epsilon_{YN}(t)$ is delayed with respect to $\lambda(t)$ to match the experimental pulse sequence. In comparison with Hamiltonian [Eq. (1)], this simulation adds cross-Kerr interactions between each qubit and the storage mode, a self-Kerr term on the storage mode, but it does not take into account the readout resonator.

In addition to the Hamiltonian (G1), we supply the solver with eight collapse operators to simulate the dynamics of the following master equation:

<table><caption>TABLE III. Parameters extracted from the photocounting simulations using the multiplexing or yes-no qubit. All parameters except those related to the multiplexing qubit are determined using a fit of the yes-no qubit simulation to Figs. 2(a) and 2(c). Parameters related to multiplexing qubit are obtained using a fit of the simulation to Figs. 2(b) and 2(d).</caption><thead><tr><th>Parameter</th><th>Fitted values</th></tr></thead><tbody><tr><td>μ</td><td>1.45 (mV μs)<sup>-1</sup></td></tr><tr><td>χ<sub>s,YN</sub></td><td>1.42 MHz</td></tr><tr><td>χ<sub>s,MP</sub></td><td>4.9 MHz</td></tr><tr><td>χ<sub>s,s</sub></td><td>-0.02 MHz</td></tr><tr><td>χ<sub>s,s,YN</sub></td><td>-0.003 MHz</td></tr><tr><td>χ<sub>s,s,MP</sub></td><td>-0.08 MHz</td></tr><tr><td>n<sub>th,s</sub></td><td>0.03</td></tr></tbody></table>

$$ \dot{\rho} = -\frac{i}{\hbar} [\hat{H}_1, \rho] + 2\Gamma_{\phi,s} \mathcal{L}(\hat{n}_s)\rho \\ + (1 + n_{\text{th},s}) \Gamma_{1,s} \mathcal{L}(\hat{a}_s)\rho + n_{\text{th},s} \Gamma_{1,s} \mathcal{L}(\hat{a}_s^\dagger)\rho \\ + \frac{1}{2} \Gamma_{\phi,YN} \mathcal{L}(\hat{\sigma}_{z,YN})\rho + \Gamma_{1,YN} \mathcal{L}(\hat{\sigma}_{YN}^-)\rho \\ + \frac{1}{2} \Gamma_{\phi,MP} \mathcal{L}(\hat{\sigma}_{z,MP})\rho + \Gamma_{1,MP} \mathcal{L}(\hat{\sigma}_{MP}^-)\rho \quad (G2) $$

with $n_{th,s}$ the expectation values of $\hat{n}_s$ when the system is at rest due to thermal occupation. All decoherence and relaxation rates are measured using previously explained calibration.

The master equation is solved using the function “mesolve” of QUTIP starting from a thermal state with the $n_{th,s}$ photon in the storage mode, the yes-no qubit in the ground state $|g\rangle$, and the multiplexing qubit also in the ground state $|g\rangle$. The solver iteratively computes the density matrix with a 10-ns time step during the displacement pulse and the $\pi$ pulse. We compute the expectation value $\langle\hat{\sigma}_{z,YN}\rangle$ at the end of the sequence and convert it into a probability $P_e$ of finding the yes-no qubit in the $|e\rangle$ state.

This simulation can be used to reproduce the experiment in Figs. 2(a) and 2(b) of the main text by adjusting the following parameters: $\{\mu = \epsilon_{max}/V_{max,s}, \chi_{s,YN}, \chi_{s,s}, \chi_{s,s,YN}, n_{th,s}\}$. Note that we need to run the simulation for every couple of parameters ($V_{max,s}, \delta f_{YN}$). Table III compiles the values of the fitted parameters.

## b. Photocounting with the multiplexing qubit

A second simulation (fluorescence simulation) is carried out to compare the photon counting experiment in Figs. 3(a) and 3(b) using a single drive on the multiplexing qubit with theory. This experiment also starts with a storage-mode displacement, but it is followed by a 2-(µs) Gaussian pulse on the multiplexing qubit at the frequency $f_{MP} - \delta f_{MP}$ with an amplitude expressed as a Rabi frequency $\Omega = \chi_{s,MP}/4$. The measured reflection coefficient of the multiplexing
qubit $r(\delta f_{MP})$ can be expressed using input-output theory as [49]

$$r(\delta f_{MP}) = \frac{\langle a_{\text{out}} \rangle}{\langle a_{\text{in}} \rangle} = \frac{\langle a_{\text{in}} \rangle - \sqrt{\Gamma_{1,\text{MP}} \langle \sigma_{-,\text{MP}} \rangle}}{\langle a_{\text{in}} \rangle} \\ = 1 - \frac{\sqrt{\Gamma_{1,\text{MP}}}}{\langle a_{\text{in}} \rangle} \langle \hat{\sigma}_{-,\text{MP}} \rangle. \qquad (G3)$$

And since the Rabi frequency is given by $\Omega = \sqrt{\Gamma_{1,MP}}|\langle a_{in} \rangle|/\pi$, we get an emission coefficient

$$1 - \text{Re}[r(\delta f_{\text{MP}})] = \frac{\Gamma_{1,\text{MP}}}{\pi \Omega} \text{Re}(e^{-i \arg(\langle a_{\text{in}} \rangle)} \langle \hat{\sigma}_{-,\text{MP}} \rangle)$$

in the frame rotating at $f_{MP} - \delta f_{MP}$ for the multiplexing qubit. If we set the phase of the drive so that $i\langle a_{in} \rangle \ge 0$, meaning we drive the qubit along $\sigma_{x,MP}$, the emission coefficient becomes

$$1 - \text{Re}[r(\delta f_{\text{MP}})] = \frac{\Gamma_{1,\text{MP}}}{2\pi\Omega} \langle \hat{\sigma}_{y,\text{MP}} \rangle.$$

The Hamiltonian of the problem in the frame rotating at $f_s - \chi_{s,MP}/2 - \chi_{s,YN}/2$ for the storage mode, $f_{YN}$ for the yes-no qubit, and $f_{MP} - \delta f_{MP}$ for the multiplexing qubit reads

$$\hat{H}_2/h = \delta f_{MP} \frac{\hat{\sigma}_{z,MP}}{2} - \chi_{s,YN} \hat{n}_s \frac{\hat{\sigma}_{z,YN}}{2} - \chi_{s,MP} \hat{n}_s \frac{\hat{\sigma}_{z,MP}}{2} \\ - \chi_{s,s} \hat{n}_s (\hat{n}_s - 1) - \chi_{s,s,YN} \hat{n}_s (\hat{n}_s - 1) \frac{\hat{\sigma}_{z,YN}}{2} \\ - \chi_{s,s,MP} \hat{n}_s (\hat{n}_s - 1) \frac{\hat{\sigma}_{z,MP}}{2} + \frac{\Omega}{2} \epsilon_{MP}(t) \hat{\sigma}_{x,MP} \\ + \frac{\lambda(t)}{2\pi} (\epsilon_{\text{max}} e^{i\pi(\chi_{s,MP}+\chi_{s,YN})t} \hat{a}_s \\ + \epsilon_{\text{max}}^* e^{-i\pi(\chi_{s,MP}+\chi_{s,YN})t} \hat{a}_s^\dagger), \qquad (G4)$$

where $\epsilon_{MP}(t) \ge 0$ is a Gaussian function of duration $2 \mu s$, width 250 ns, and amplitude 1. $\epsilon_{MP}(t)$ is delayed compared to $\lambda(t)$ to reproduce the experimental pulse sequence. We add to this Hamiltonian the same relaxation and decoherence channels as for the yes-no simulation [see Eq. (G2)] for which the decoherence and relaxation rates are measured independently. The resulting master equation differs only from the yes-no simulation by the Rabi drive that addresses the multiplexing qubit instead of the yes-no qubit. The master equation is solved using the mersolve function of QUTIP with a time step of 5.25 ns starting from a thermal state with $n_{th,s}$ photons for storage and the yes-no qubit and the multiplexing qubit in the ground state $|g\rangle$. Finally, the expectation value $\langle \hat{\sigma}_{y,MP} \rangle$ is computed and integrated during the $2 \mu s$ of the pulse.

We compare the measured emission coefficient in Figs. 2(b) and 2(d) to the simulated signal $A\langle \hat{\sigma}_{y,MP} \rangle$, where

A is left as a free parameter due to a small parasitic reflection in the measurement setup and thermal population. The parameters $\{\mu, \chi_{s,t}, \chi_{s,s}, \chi_{s,s,t}, n_{th,s}\}$ are already set by the calibration above using the simulation of the yes-no qubit. From the fluorescence simulation, we thus extract the parameters $\{\chi_{s,MP}, \chi_{s,s,MP}, A\}$ by comparing the experimental observations in Figs. 2(b) and 2(d) with the simulation for various $V_{max,s}$ and $\delta f_{MP}$. Fitted values are given in Table III. Finally, we run the yes-no simulation again taking into account the updated multiplexing qubit parameters. As expected, only small changes in the results of the yes-no qubit simulation are observed.

## 2. Evolution of the average photon number in the storage mode

We simulate the filling of the storage mode by a displacement pulse on the resonator. We simulate the same master equation used for the photocounting simulations with parameters obtained from the photocounting simulations (see Table III) but without applying any drive on the qubits. Only the displacement pulse on the storage mode is modeled, i.e., $\epsilon_{MP}(t) = 0$, $\delta f_{MP} = 0$, $\epsilon_{YN}(t) = 0$, and $\delta f_{YN} = 0$.

The mersolve function of QUTIP computes the density matrix with a time step of 10 ns and returns the mean number of photons in the storage mode at the end of the displacement pulse for various drive amplitudes. Figure 10 shows the square root of the expected mean photon number as a function of the amplitude $\epsilon_{max}$. We obtain a scaling factor $\sqrt{\langle n_s \rangle} = 85.9 \text{ V}^{-1} \text{ V}_{max,s}$ used in the photon number calibration of the storage mode.

## 3. Simulation of multiplexed readout

In this subsection, we simulate how a frequency comb reflects off the multiplexing qubit. We write the Hamiltonian in the frame rotating at $f_s - \chi_{s,MP}/2 - \chi_{s,YN}/2$ for the storage mode and at the qubit frequencies for the qubits as

$$\hat{H}_3/h = -\chi_{s,YN} \hat{n}_s \frac{\hat{\sigma}_{z,YN}}{2} - \chi_{s,MP} \hat{n}_s \frac{\hat{\sigma}_{z,MP}}{2} - \chi_{s,s} \hat{n}_s (\hat{n}_s - 1) \\ - \chi_{s,s,YN} \hat{n}_s (\hat{n}_s - 1) \frac{\hat{\sigma}_{z,YN}}{2} - \chi_{s,s,MP} \hat{n}_s (\hat{n}_s - 1) \frac{\hat{\sigma}_{z,MP}}{2} \\ + \frac{\Omega}{2} [\epsilon_{\text{comb}}(t) \hat{\sigma}_{\text{MP}}^{+} + \epsilon_{\text{comb}}^{*}(t) \hat{\sigma}_{\text{MP}}] \\ + \frac{\lambda(t)}{2\pi} (\epsilon_{\text{max}} e^{i\pi(\chi_{s,MP}+\chi_{s,YN})t} \hat{a}_s \\ + \epsilon_{\text{max}}^{*} e^{-i\pi(\chi_{s,MP}+\chi_{s,YN})t} \hat{a}_s^{\dagger}), \qquad (G5)$$

where $\Omega = \chi_{s,MP}/2$, and $\epsilon_{comb}(t)$ is the product of a Gaussian function with the sum of nine complex exponential $\sum_{k=0}^{8} \exp(2i\pi\chi_{s,MP}kt)$. The Gaussian envelope of $\epsilon_{comb}(t)$ has a duration of $2 \mu s$, a width of 250 ns, and a
maximum amplitude of 1, and the delay between $\epsilon_{comb}(t)$ and $\lambda(t)$ reproduces the experimental sequence. The master equation (G2) is used with a time step of 1 ns for various amplitude $\epsilon_{max}$. We obtain the time evolution of $\langle\sigma_{y,MP}\rangle$ enabling us to compare the experimental measurements of Fig. 2(e) to the model. To do so, we integrate the simulated function $\langle\sigma_{y,MP}\rangle \times \cos(2\pi\chi_{s,MP}kt)$ for each integer k, similar to the demultiplexing processing we perform on the multiplexed experimental signal. Note that, in the case k = 0, we need to divide the integral by 2 in order to perform a proper demultiplexing. By combining this simulation with the photon number calibration, we get the expected values of the nine multiplexing readout signals as a function of the mean number of photons in the storage mode used in Fig. 2(e).

### 4. Simulation of measurement-induced dephasing on the storage mode

In this part, we simulate only the multiplexing qubit and the storage mode to decrease the computational cost of the simulation. The Hamiltonian of the simulation in the frame rotating at the multiplexing qubit-resonant frequency and at $f_s + \delta f_s^0$ for the storage mode is

$$ \hat{H}_4/h = -\chi_{s,MP} \frac{\hat{\sigma}_{z,MP} + 1}{2} \hat{n}_s - \delta f_s^0 \hat{n}_s \\ - \chi_{s,s,MP} \hat{n}_s (\hat{n}_s - 1) \frac{\hat{\sigma}_{z,MP} + 1}{2} \\ + \frac{\Omega}{2} [\epsilon_{comb}(t) \hat{\sigma}_{MP}^+ + \epsilon_{comb}^*(t) \hat{\sigma}_{MP}^-], \quad (G6) $$

where $\epsilon_{comb}(t)$ is the product of a Gaussian function with the sum of nine complex exponential $\sum_{k=0}^{8} \exp(2i\pi\chi_{s,MP}k\tau)$. The width of the Gaussian function is equal to one quarter of the duration $t$ of the pulse. We add four dephasing and relaxation channels to this Hamiltonian to obtain the master equation

$$ \dot{\rho} = -\frac{i}{\hbar} [\hat{H}_4, \rho] + 2\Gamma_{\phi,s}\mathcal{L}(\hat{n}_s)\rho + \Gamma_{1,s}\mathcal{L}(\hat{a}_s)\rho \\ + \frac{1}{2}\Gamma_{\phi,MP}\mathcal{L}(\hat{\sigma}_{z,MP})\rho + \Gamma_{1,MP}\mathcal{L}(\hat{\sigma}_{MP}^-)\rho. \quad (G7) $$

The storage is initialized in a coherent state of amplitude $\beta = 1.55$ and the multiplexing qubit is initialized in state $|g\rangle$. We simulate the dynamics of the system for a pulse duration $t$ going from 100 ns to 5 µs and for $\Omega$ ranging from 0 to $2\chi_{s,MP}$. We compute the expectation value of $\hat{X} = (\hat{a}_s + \hat{a}_s^\dagger)/2$ at the end of each simulation. For a given $\Omega$, we extract the time evolution of $\langle\hat{X}\rangle$ under the influence of the multiplexed measurement as shown on Fig. 19(a). This decaying sinusoid is fitted using Eq. (B9) to obtain the oscillation frequency $\delta f_s$ and the decay rate $\Gamma_{d,s}$.

<figure><img src="image_24.png" /><figcaption>FIG. 19. Simulations of the measurement-induced dephasing rate and of the ac-Stark shift induced by a frequency comb. (a) Ramsey-like oscillations of the storage mode for $\Omega = \chi_{s,MP}/2$ and an initial coherent field amplitude $\beta = -1.55$. Blue dots are the simulated expectation values of $\hat{X}$ and the red line is the theory given by Eq. (B9). (b) Simulated measurement-induced dephasing rate $\Gamma_{d,s}$ and ac-Stark shift as a function of $\Omega/\chi_{s,MP}$ for various values of $\chi_{s,MP}$. Simulations show the same pattern with maxima and minima for some specific values of $\Omega/\chi_{s,MP}$ as in the experiment in Fig. 3(c). (c) Simulated measurement-induced dephasing rate and ac-Stark shift as a function of $\Omega/\chi_{s,MP}$ for various initial coherent-state amplitudes $\beta$ in the storage mode. For $\Omega/\chi_{s,MP} > 1$, we see a difference of about 20% between $\beta = 1.6$ and $\beta = 1.2$. (d) Simulations for $\chi_{s,MP} \gg \Gamma_{1,MP}$. The evolution of the measurement-induced dephasing and ac-Stark shift with $\Omega$ is different compared to the case of (b). The evolution of the measurement-induced dephasing rate and the ac-Stark shift seems to be given by the ratio $2\pi\Omega/\Gamma_{1,MP}$. The red line is a guide for the eyes representing a quadratic function. It shows that the measurement-induced dephasing rate increases linearly with $\Omega^2$ for small drive amplitudes. For the ac-Stark shift, on the contrary to (b), the detuning is constant at the small drive amplitudes, and then two frequencies appear with comparable contributions to the Ramsey oscillations. The two frequencies evolve linearly with $\Omega$. (e) Example of simulated Ramsey oscillations exhibiting two frequencies.</figcaption></figure>

Figures 14(c) and 14(d) show the measurement-induced dephasing and ac-Stark shift as a function of the amplitude of the comb $\Omega$ for two sets of coherent-state amplitudes $\beta$ and detuning $\delta f_s^0$.
We identify three interesting features. The first one is the evolution of the shape of the curves $\delta f_s(\Omega)$ and $\Gamma_{d,s}(\Omega)$ with $\chi_{s,MP}$. We repeat the simulation using a square pulse envelope instead of Gaussian pulse for $\epsilon_{comb}$ to make the simulation faster for several values of $\chi_{s,MP}$ from 1.5 to 13.2 MHz in steps of 1.4 MHz. We observe that $\delta f_s(\Omega)$ and $\Gamma_{d,s}(\Omega)$ increase as $\chi_{s,MP}$ becomes larger but that the maxima and minima of the curve are always found for the same $\Omega/\chi_{s,MP}$ ratio [Fig. 19(b)] as predicted (see Appendix C).

The second observation is that $\delta f_s(\Omega)$ and $\Gamma_{d,s}(\Omega)$ vary with the initial coherent-state amplitude $\beta$ [Fig. 19(c)]. As the frequency comb contains only a finite number of frequencies $f_{MP} - n\chi_{s,MP}$ with $0 \le k \le 8$, the decoherence rate between two Fock states $|i\rangle$ and $|j\rangle$ depends on $i$ and $j$. As the probability distribution of photon numbers depends on the amplitude of the coherent state, the mean decoherence rate also depends on the coherent-state amplitude.

The third observation is that in the regime $\chi_{s,MP} \gg \Gamma_{1,MP}/2\pi$, the dephasing rate and ac-Stark shift are a function of the ratio $2\pi\Omega/\Gamma_{1,MP}$ as shown in Fig. 19 (d). The dephasing rate increases as $\Omega^2$ until a plateau is reached for $2\pi\Omega/\Gamma_{1,MP} = 0.7$; this plateau is predicted (see Appendix C). In contrast, the Stark shift is constant for $2\pi\Omega/\Gamma_{1,MP} < 0.3$ and splits into two frequencies (two oscillations on top of each other in Ramsey interferometry) with a splitting proportional to $2\pi\Omega/\Gamma_{1,MP}$. Since there are two frequencies, we use Eq. (B10) to fit the simulated Ramsey oscillations for $\chi_{s,MP} \gg \Gamma_{1,MP}/2\pi$. In practice, Eq. (B10) is a good fit function because a Fourier analysis shows that the signal is composed of two frequencies with the same amplitude. Figure 19(e) shows an example of simulated Ramsey oscillations for $\chi_{s,MP} \gg \Gamma_{1,MP}/2\pi$.

## APPENDIX H: DENSITY MATRIX ELEMENTS

In this part, we explain how one can calculate the density matrix of the storage mode from the measured Wigner function. It is the recipe we use to produce the bottom part of Fig. 3(a) in the main text. We further present original results on the decay of density matrix elements when the multiplexing qubit is driven by a single tone or by the comb of frequencies. We characterize the quantum nondemolition nature of our photocounter. Finally, we present an experiment in which we show revivals of density matrix elements as a function of time and show simulations that reproduce them qualitatively. We discuss a new quantity called the mean coherence and show its measured evolution in various measurement configurations.

### 1. Density matrix reconstruction

The Wigner tomography contains all the information about the state of the storage mode. We explain below how we reconstruct the density matrix from the measured

Wigner function. We compute the Wigner map for every operator $|n\rangle\langle m|$ with $|n\rangle$ and $|m\rangle$ two Fock states with $n$ and $m$ photons. The mean value of those operators is equal to the $(n, m)$ element $\rho_{nm}$ of the density matrix. Using the mathematical expression of $\langle x|n\rangle$,

$$ \langle x|n\rangle = \psi_n(x) = \left(\frac{2}{\pi}\right)^{1/4} \frac{1}{\sqrt{2^n n!}} \mathcal{H}_n(\sqrt{2x})e^{-x^2} \quad (H1) $$

with $\mathcal{H}_n(x) = (-1)^n e^{x^2} (d^n/dx^n) e^{-x^2}$ the Hermite polynomial function of order $n$, and $|x\rangle$ the eigenvector of the quadrature $(\hat{a}_s + \hat{a}_s^\dagger)/2$ associated with the eigenvalue $x$. The Wigner map of the operators $|n\rangle$ and $|m\rangle$ becomes

$$ W_{|n\rangle\langle m|}(x, p) = \frac{1}{\pi} \int dy e^{-2ipy} \psi_n(x + y/2) \overline{\psi}_m(x - y/2), \quad (H2) $$

and the matrix element $\rho_{nm}$ of the storage mode is given by Eq. (B5)

$$ \rho_{nm} = \pi \iint dx dp W_{|n\rangle\langle m|}(x, p) W_{\rho}(x, p). \quad (H3) $$

### 2. Accessing the measurement-induced dephasing rate

In order to characterize the decoherence due to the multiplexed measurement, we use a renormalization of the density matrix elements in order to remove most of the effects of the storage-mode relaxation. Let us now show that in the absence of Hamiltonian evolution and measurement backaction, the quantity $|\rho_{nm}|/\sqrt{\rho_{nn}\rho_{mm}}$ evolves only because of dephasing and that its dynamics is not affected by relaxation. We consider the storage mode alone under the influence of its relaxation and dephasing channels in a frame rotating at $f_s$,

$$ \dot{\rho} = \Gamma_{1,s} \mathcal{L}(\hat{a}_s) \rho + 2\Gamma_{\phi,s} \mathcal{L}(\hat{a}_s^\dagger \hat{a}_s) \rho. \quad (H4) $$

From this equation, we can compute the time derivative of the density matrix element

$$ \dot{\rho}_{nm} = \Gamma_{1,s} \left( \rho_{n+1,m+1} \sqrt{(n+1)(m+1)} - \frac{n+m}{2} \rho_{nm} \right) - \Gamma_{\phi,s} \rho_{nm} (n-m)^2. \quad (H5) $$

If the storage mode is initialized in a coherent state $|\alpha_o\rangle$, the solution of the equation is

$$ \rho_{nm}(t) = e^{-|\alpha_o|^2 e^{-\Gamma_{1,s}t}} \frac{\alpha_o^m e^{-m\Gamma_{1,s}t/2} (\alpha_o^*)^n e^{-n\Gamma_{1,s}t/2}}{\sqrt{n!m!}} e^{-\Gamma_{\phi,s}(n-m)^2t}, \quad (H6) $$
and we get

$$ \frac{|\rho_{nm}|}{\sqrt{\rho_{nn}\rho_{mm}}} (t) = e^{-\Gamma_{\phi,s}(n-m)^2 t}. \qquad (H7) $$

Thus, indeed, the renormalization removes the effect of the relaxation rate $\Gamma_{1,s}$ and characterizes only the dephasing rate. Under the action of measurement, the dephasing rate $\Gamma_{\phi,s}$ is increased by the measurement-induced dephasing rate. We use this property in the following sections and study the evolution of $\rho_{nm}/\sqrt{\rho_{nn}\rho_{mm}}$.

### 3. Decoherence of the storage mode induced by a single measurement drive

Measuring whether there are $n$ photons spoils the coherence of the superposition between Fock state $|n\rangle$ and Fock state $|m \neq n\rangle$. We evidence this dephasing by observing the evolution of $\rho_{nm}$.

$$ \begin{pmatrix} -\Gamma_{1,\text{MP}}/2 & -2\pi\Delta + \frac{n+m}{2}2\pi\chi_{s,\text{MP}} \\ 2\pi\Delta - \frac{n+m}{2}2\pi\chi_{s,\text{MP}} & -\Gamma_{1,\text{MP}}/2 \\ 0 & 2\pi\Omega \\ 0 & 0 \end{pmatrix} \begin{pmatrix} 0 & 0 \\ -2\pi\Omega & 0 \\ -\Gamma_{1,\text{MP}} & -\Gamma_{1,\text{MP}} - i \frac{n+m}{2}2\pi\chi_{s,\text{MP}} \\ -i \frac{n+m}{2}2\pi\chi_{s,\text{MP}} & 0 \end{pmatrix}. \qquad (H8) $$

Figure 20 shows the measured density matrix decoherence rates $\Gamma_{d,s}^{nm}$ and the above theory for $n$ and $m$ going from 0 to 4 [with an offset corresponding to the natural dephasing rate in Eq. (H7)]. As expected in a regime with resolved resonance peaks $(2\pi\chi_{s,\text{MP}}|m-n| > \Gamma_{2,\text{MP}})$, the decoherence rate $\Gamma_{d,s}^{nm}$ between Fock states $|n\rangle$ and $|m\rangle$ is larger when the single drive probes whether there are $n$ photons or $m$ photons with a moderate drive amplitude $\Omega$ (dependence on $\Omega$ not shown here). For much larger drive amplitude $\Omega$, one can increase the decoherence rate $\Gamma_{d,s}^{nm}$ by driving with a detuning $\Delta_{\text{MP}} = (n+m)\chi_{s,\text{MP}}/2$, similar to dispersive qubit readout, which is optimal for information extraction at large drive power and for a drive frequency detuned by $\chi_{s,\text{MP}}/2$. In fact, this regime would become particularly attractive for poorly resolved resonances as a function of the photon number $(2\pi\chi_{s,\text{MP}}|m-n| < \Gamma_{2,\text{MP}})$. Premises of this effect are visible in Fig. 20, as the maximal decoherence rate occurs at a detuning slightly closer to $(n+m)\chi_{s,\text{MP}}/2$, with a stronger effect for small $|m-n|$, both in theory and in the experimental observations. The small discrepancy between theory and experiment, in particular, the asymmetry as a function of $n$ and $m$, may be explained by the photon loss rate of the storage mode, which is not captured in the simplified theoretical model.

We prepare the storage-mode state in a coherent state with an amplitude $\beta = -1.7$ and probe the multiplexing qubit during a time $t$ with a drive at the frequency $f_{\text{MP}} - \Delta_{\text{MP}}$ before doing Wigner tomography. For various times $t$ and detunings $\Delta_{\text{MP}}$, we compute the density matrix of the storage mode using Eq. (H3). One can fit the time evolution of $|\rho_{nm}|/\sqrt{\rho_{nn}\rho_{mm}}$ with a decreasing exponential function. The extracted decoherence rate $\Gamma_{d,s}^{nm}(\Delta_{\text{MP}})$ is then compared to the theoretical value.

In Ref. [50], we show that an exact, infinite-order adiabatic elimination of the multiplexing qubit probed with a single frequency drive is possible under the assumption that there is no photon loss in the storage mode. It shows that the decoherence rate between the Fock states $|n\rangle$ and $|m\rangle$ due to the measurement is given by the highest eigenvalue, which are all negatives, of the following matrix:

$$ \begin{pmatrix} 0 & 0 \\ -2\pi\Omega & 0 \\ -\Gamma_{1,\text{MP}} & -\Gamma_{1,\text{MP}} - i \frac{n+m}{2}2\pi\chi_{s,\text{MP}} \\ -i \frac{n+m}{2}2\pi\chi_{s,\text{MP}} & 0 \end{pmatrix}. \qquad (H8) $$

### 4. Multiplexed measurement vs single tone measurement

In Fig. 3(a) of the main text, one sees that the dephasing of the storage mode induced by the measurement is stronger for multiplexed measurement than for single tone measurement. This conclusion is based on the Wigner tomography of the storage mode in three distinct cases. The storage mode is initialized in a coherent state of amplitude $\beta = -1.55$. Then, before performing the tomography of the storage mode, we either (i) wait for a time $t$, (ii) probe the multiplexing qubit for a time $t$ at a single frequency $f_{\text{MP}} - \chi_{s,\text{MP}}$ corresponding to one photon, or (iii) with a frequency comb.

From the measured Wigner functions, we compute the density matrix of the storage mode $\rho(t)$ for various times $t$ for the three cases and compare the evolution of the normalized elements $\rho_{nm}(t)$ (see Fig. 21). Without any drive on the multiplexing qubit [circles and case (i)], the density matrix elements decay due to natural dephasing only. Clearly, the drive on the multiplexing qubit induces a decay of the coherences, with a stronger effect when the comb is turned on than when a single drive is turned on. We conclude that a multiplexing measurement extracts more information than a single measurement.

The effect on $\rho_{02}$ when probing with a resonant drive for $n=1$ is consistent with the significant measurement-induced detuning that can be read off the top right plot
<figure><img src="image_25.png" /></figure>

$$ (n, m) = (0, 2), (1, 3), (2, 4) $$

<figure><img src="image_26.png" /></figure>

$$ (n, m) = (0, 3), (1, 4), (0, 4) $$

<figure><img src="image_27.png" /></figure>

FIG. 20. Decoherence rate of superpositions between Fock states induced by a single drive on the multiplexing qubit. In each panel, dots are obtained using Eq. (H3) on the measured Wigner function of the storage mode when driven by a single drive at $f_{MP} - \Delta_{MP}$ with an amplitude $\Omega = \chi_{s,MP}/2$. Lines represent the highest eigenvalue of (H8) without any free parameters. An offset equal to $\Gamma_{\phi,s}(n-m)^2$, which is the intrinsic dephasing of the storage mode, is added to obtain the total dephasing rate.

of Fig. 20 (blue, value 1 on the horizontal axis). Apparently, when driving with a comb, such an effect combines with the ones on the $n=0$ and $n=2$ resonances, and other components to induce a stronger overall measurement rate. We investigate this comb effect more precisely in Appendix H 6.

<figure><img src="image_28.png" /></figure>

FIG. 21. Dynamics of the storage-mode coherences under various measurement schemes. Normalized off-diagonal elements of the density matrix extracted from the measured Wigner function of the storage mode as a function of time. The figure focuses on three elements $\rho_{01}$ (blue), $\rho_{12}$ (orange), and $\rho_{02}$ (green). Circles: case (i) without driving the multiplexing qubit. Squares: case (ii) where a single tone at $f_{MP} - \chi_{s,MP}$ drives the multiplexing qubit with a strength $\Omega = \chi_{s,MP}/2$. Triangles: case (iii) where the multiplexing qubit is driven by a comb of nine peaks with the same strength $\Omega$ each.

### 5. Quantum nondemolition nature of the multiplexed measurement

The goal of this subsection is to quantify the quantum nondemolition (QND) nature of our multiplexed measurement. A measurement is said to be QND if

(i) the measurement time is very short compared to the timescale of evolution of the system under study,

<figure><img src="image_29.png" /></figure>

FIG. 22. Impact of multiplexed measurement on the occupation of the Fock states. (a) Measured probability to find the storage mode in Fock state $|0\rangle$ as a function of time $t$ for various comb drive amplitudes $\Omega/\chi_{s,MP}$. (b) Measured diagonal elements of the density matrix integrated during 5 µs as a function of the drive amplitude $\Omega/\chi_{s,MP}$. (c) Decay of the average photon number in time for various drive amplitudes. (d) Average photon number evolution as a function of $\Omega/\chi_{s,MP}$.
(ii) the interaction with the probe does not disturb the quantum state of the system if it belongs to the measurement basis.

If a photocounter is QND, the diagonal elements of the density matrix, i.e., the average Fock-state populations, of the resonator are unchanged (on average on all measurements) by the measurement process. In our experiment, we observe that the diagonal elements of the density matrix in the energy basis evolve owing to the decay of the storage mode but do not strongly depend on the measurement strength. To be more accurate, we notice that for large probe amplitude [red and purple points in Fig. 22(a)], the probability of finding the storage mode with 0 photons is slightly lower. This dependence on the amplitude of the drive Ω is best characterized by extracting the populations [Fig. 22(b)] and photon number [Fig. 22(d)] integrated during T = 5 µs as a function of Ω/χs,MP. For small drive amplitude Ω/χs,MP < 0.1, the probability to find a given number of photons does not change with Ω/χs,MP, but for larger drive amplitude the resonator gets populated probably because of Zeno effect due to the non-Markovian environment originating from the multiplexing qubit.

In practice, for small drive amplitude and a measurement time of 5 µs, the relaxation dynamics of the system during the measurement process increases the probability of having 0 photons at the end of the measurement by approximately 10%. We find that the mean photon number is decreased by the same percentage.

**6. Off-diagonal density matrix elements and revivals of the coherences**

In the main text, we use Wigner tomography in order to observe Ramsey-like oscillations of the storage mode. In fact, using Eq. (H3), the Wigner function allows us to visualize the dynamics of every off-diagonal element of the storage density matrix to gain insight into the physics of the dephasing process.

Figures 23(a) and 23(b) show the decay of off-diagonal elements $\rho_{12}$ and $\rho_{13}$ as a function of time. For small drive amplitudes $\Omega < 0.5\chi_{s,MP}$, off-diagonal elements decay faster when $\Omega$ is increased since more and more information is extracted per unit time by the drive. As larger drive amplitudes are reached, off-diagonal elements start oscillating. The contrast of these coherence revivals becomes more pronounced as the drive amplitude becomes larger and they exhibit a quasiperiodicity (see Fig. 24). Note that the 10% deviation to exact periodicity may originate from the Gaussian pulse shaping of the comb. This behavior is qualitatively reproduced by our simulations [Fig. 23(c)].

**7. Mean coherence between Fock states**

Since the driving frequency comb holds the promise to probe how many photons are in the storage mode, it should affect all coherences $\rho_{nm}$. In this section, we introduce two

<figure><img src="image_30.png" /><figcaption>FIG. 23. Normalized off-diagonal elements of the storage matrix density. (a) Measured normalized coherence $|\rho_{12}|/\sqrt{\rho_{11}\rho_{22}}$ between Fock states $|1\rangle$ and $|2\rangle$ of the storage mode as a function of time and for various amplitudes $\Omega$ of the driving frequency comb. (b) Similar plot for $|\rho_{13}|/\sqrt{\rho_{11}\rho_{33}}$. (c) Results of the simulation of the master equation (G7).</figcaption></figure>

ways of characterizing the impact of the multiplexed photocounting on the global coherence of the storage mode.

The first one shown in the main text in Fig. 3(b) is the quadrature of the storage mode in the frame rotating at the

<figure><img src="image_31.png" /><figcaption>FIG. 24. Normalized off-diagonal elements of the storage matrix density for the largest measurement amplitudes. Measured normalized coherence $|\rho_{12}|/\sqrt{\rho_{11}\rho_{22}}$ between Fock states $|1\rangle$ and $|2\rangle$ of the storage mode as a function of time and for various amplitudes $\Omega > 0.9\chi_{s,MP}$ of the driving frequency comb.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>FIG. 25. Impact of the multiplexed drive on storage-mode quadrature and mean coherence between Fock states. Left column: Ramsey oscillations measured (dots) and fits (line) using Eqs. (B9) and (B10) for five multiplexing drive amplitudes. Right column: quadrature of storage mode in the frame rotating at its resonance frequency (blue dots) and mean coherence between Fock states $C_ρ(t)$ (orange squares) defined in Eq. (H9) for five different drive amplitudes. The standard deviations of each data point are about 0.07 for the quadrature and 0.015 for the mean coherence. The red solid lines are exponential fits of the quadrature decay.</figcaption></figure>

frequency of this mode when the qubit is probed by a comb. It can be expressed as Re[(⟨$\hat{X}$⟩ + i⟨$\hat{P}$⟩)e⁻²iπδf_s t] with ⟨$\hat{X}$⟩ and ⟨$\hat{P}$⟩ the expectation values of the quadratures in the frame rotating at the frequency of the storage drive. This quantity is related to the first off-diagonal of the density matrix.

We introduce a second quantity: the mean coherence $C_ρ$ between Fock states 0 to 4. It is defined as

$$C_{\rho} = \text{Mean} \left[ \frac{|\rho_{ij}|}{\sqrt{\rho_{ii}\rho_{jj}}} \right]. \qquad (H9)$$

$C_ρ(t)$ contains the information about the dephasing between every different Fock state. The left part of Fig. 25 shows oscillations of the storage-mode quadratures in the frame of the drive on the storage mode (for state preparation and Wigner tomography) for various multiplexing qubit drive amplitudes. On the right part of Fig. 25, we display the mean quadrature ⟨$\hat{a}_s$ + $\hat{a}_s^†$⟩ = Re[(⟨$\hat{X}$⟩ + i⟨$\hat{P}$⟩)e⁻²iπδf_s t] in the frame rotating at the storage-mode frequency and the mean coherence $C_ρ$. Those two quantities show the same dynamics leading to the same dephasing rate and both quantities can be used to characterize it. The revivals that can be seen on each of the density matrix off-diagonal elements (see Fig. 23) also appear in the evolution of the quadrature and of the mean coherence between Fock states.

[1] C. Guerlin, J. Bernu, S. Deleglise, C. Sayrin, S. Gleyzes, S. Kuhr, M. Brune, J.-M. Raimond, and S. Haroche, *Progressive Field-State Collapse and Quantum Non-Demolition Photon Counting*, Nature (London) **448**, 889 (2007).

[2] B. R. Johnson, M. D. Reed, A. A. Houck, D. I. Schuster, L. S. Bishop, E. Ginossar, J. M. Gambetta, L. DiCarlo, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, *Quantum Non-Demolition Detection of Single Microwave Photons in a Circuit*, Nat. Phys. **6**, 663 (2010).

[3] B. Peaudecerf, T. Rybarczyk, S. Gerlich, S. Gleyzes, J. M. Raimond, S. Haroche, I. Dotsenko, and M. Brune, *Adaptive Quantum Nondemolition Measurement of a Photon Number*, Phys. Rev. Lett. **112**, 080401 (2014).

[4] S. Haroche, M. Brune, and J. M. Raimond, *Measuring Photon Numbers in a Cavity by Atomic Interferometry: Optimizing the Convergence Procedure*, J. Phys. II (France) **2**, 659 (1992).

[5] R. Dassonneville, R. Assouly, T. Peronnin, P. Rouchon, and B. Huard, *Number-Resolved Photocounter for Propagating Microwave Mode*, Phys. Rev. Applied **14**, 044022 (2020).

[6] J. C. Curtis, C. T. Hann, S. S. Elder, C. S. Wang, L. Frunzio, L. Jiang, and R. J. Schoelkopf, *Single-Shot Number-Resolved Detection of Microwave Photons with Error Mitigation*, Phys. Rev. A **103**, 023705 (2021).

[7] C. S. Wang, J. C. Curtis, B. J. Lester, Y. Zhang, Y. Y. Gao, J. Freeze, V. S. Batista, P. H. Vaccaro, I. L. Chuang, L. Frunzio, L. Jiang, S. M. Girvin, and R. J. Schoelkopf, *Efficient Multiphoton Sampling of Molecular Vibronic Spectra on a Superconducting Bosonic Processor*, Phys. Rev. X **10**, 021060 (2020).

[8] D. I. Schuster, A. A. Houck, J. A. Schreier, A. Wallraff, J. M. Gambetta, A. Blais, L. Frunzio, J. Majer, B. Johnson, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, *Resolving Photon Number States in a Superconducting Circuit*, Nature (London) **445**, 515 (2007).
[9] M. F. Gely, M. Kounalakis, C. Dickel, J. Dalle, R. Vatré, B. Baker, M. D. Jenkins, and G. A. Steele, *Observation and Stabilization of Photonic Fock States in a Hot Radio-Frequency Resonator*, Science **363**, 1072 (2019).

[10] J. G. Proakis and M. Salehi, *Communication Systems Engineering*, 2nd ed. (Prentice-Hall, Englewood Cliffs, NJ, 2002).

[11] C. Macklin, D. Hover, M. E. Schwartz, X. Zhang, W. D. Oliver, and I. Siddiqi, *A Near Quantum-Limited Josephson Traveling-Wave Parametric Amplifier*, Science **350**, 307 (2015).

[12] A. A. Houck, D. I. Schuster, J. M. Gambetta, J. A. Schreier, B. R. Johnson, J. M. Chow, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, *Generating Single Microwave Photons in a Circuit*, Nature (London) **449**, 328 (2007).

[13] O. Astafiev, A. M. Zagoskin, A. A. Abdumalikov, Y. A. Pashkin, T. Yamamoto, K. Inomata, Y. Nakamura, and J. S. Tsai, *Resonance Fluorescence of a Single Artificial Atom*, Science **327**, 840 (2010).

[14] A. A. Abdumalikov, O. V. Astafiev, Y. A. Pashkin, Y. Nakamura, and J. S. Tsai, *Dynamics of Coherent and Incoherent Emission from an Artificial Atom in a 1D Space*, Phys. Rev. Lett. **107**, 043604 (2011).

[15] P. Campagne-Ibarcq, L. Bretheau, E. Flurin, A. Auffèves, F. Mallet, and B. Huard, *Observing Interferences between Past and Future Quantum States in Resonance Fluorescence*, Phys. Rev. Lett. **112**, 180402 (2014).

[16] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, *Atom-Photon Interactions: Basic Processes and Applications* (Wiley-VCH, 2001).

[17] J.-C. Besse, S. Gasparinetti, M. C. Collodo, T. Walter, P. Kurpiers, M. Pechal, C. Eichler, and A. Wallraff, *Single-Shot Quantum Nondemolition Detection of Individual Itinerant Microwave Photons*, Phys. Rev. X **8**, 021003 (2018).

[18] Y.-F. Chen, D. Hover, S. Sendelbach, L. Maurer, S. T. Merkel, E. J. Pritchett, F. K. Wilhelm, and R. McDermott, *Microwave Photon Counter Based on Josephson Junctions*, Phys. Rev. Lett. **107**, 217401 (2011).

[19] K. Inomata, Z. Lin, K. Koshino, W. D. Oliver, J.-S. Tsai, T. Yamamoto, and Y. Nakamura, *Single Microwave-Photon Detector Using an Artificial Λ-Type Three-Level System*, Nat. Commun. **7**, 12303 (2016).

[20] A. Narla, S. Shankar, M. Hatridge, Z. Leghtas, K. M. Sliwa, E. Zalys-Geller, S. O. Mundhada, W. Pfaff, L. Frunzio, R. J. Schoelkopf, and M. H. Devoret, *Robust Concurrent Remote Entanglement between Two Superconducting Qubits*, Phys. Rev. X **6**, 031036 (2016).

[21] S. Kono, K. Koshino, Y. Tabuchi, A. Noguchi, and Y. Nakamura, *Quantum Non-Demolition Detection of an Itinerant Microwave Photon*, Nat. Phys. **14**, 546 (2018).

[22] R. Lescanne, S. Deléglise, E. Albertinale, U. Réglade, T. Capelle, E. Ivanov, T. Jacqmin, Z. Leghtas, and E. Flurin, *Irreversible Qubit-Photon Coupling for the Detection of Itinerant Microwave Photons*, Phys. Rev. X **10**, 021038 (2020).

[23] V. Schmitt, X. Zhou, K. Juliusson, B. Royer, A. Blais, P. Bertet, D. Vion, and D. Esteve, *Multiplexed Readout of Transmon Qubits with Josephson Bifurcation Amplifiers*, Phys. Rev. A **90**, 062333 (2014).

[24] 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, *Rapid High-Fidelity Multiplexed Readout of Superconducting Qubits*, Phys. Rev. Applied **10**, 034040 (2018).

[25] S. Kundu, N. Gheeraert, S. Hazra, T. Roy, K. V. Salunkhe, M. P. Patankar, and R. Vijay, *Multiplexed Readout of Four Qubits in 3D Circuit QED Architecture Using a Broadband Josephson Parametric Amplifier*, Appl. Phys. Lett. **114**, 172601 (2019).

[26] F. Arute et al., *Quantum Supremacy Using a Programmable Superconducting Processor*, Nature (London) **574**, 505 (2019).

[27] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, *Introduction to Quantum Noise, Measurement, and Amplification*, Rev. Mod. Phys. **82**, 1155 (2010).

[28] M. Boissonneault, J. M. Gambetta, and A. Blais, *Dispersive Regime of Circuit QED: Photon-Dependent Qubit Dephasing and Relaxation Rates*, Phys. Rev. A **79**, 013819 (2009).

[29] L. G. Lutterbach and L. Davidovich, *Method for Direct Measurement of the Wigner Function in Cavity QED and Ion Traps*, Phys. Rev. Lett. **78**, 2547 (1997).

[30] P. Bertet, A. Auffèves, P. Maioli, S. Osnaghi, T. Meunier, M. Brune, J. M. Raimond, and S. Haroche, *Direct Measurement of the Wigner Function of a One-Photon Fock State in a Cavity*, Phys. Rev. Lett. **89**, 200402 (2002).

[31] B. Vlastakis, G. Kirchmair, Z. Leghtas, S. E. Nigg, L. Frunzio, S. M. Girvin, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, *Deterministically Encoding Quantum Information Using 100-Photon Schrodinger Cat States*, Science **342**, 607 (2013).

[32] P. Campagne-Ibarcq, A. Eickbusch, S. Touzard, E. Zalys-Geller, N. E. Frattini, V. V. Sivak, P. Reinhold, S. Puri, S. Shankar, R. J. Schoelkopf, L. Frunzio, M. Mirrahimi, and M. H. Devoret, *Quantum Error Correction of a Qubit Encoded in Grid States of an Oscillator*, Nature (London) **584**, 368 (2020).

[33] H. M. Wiseman and G. J. Milburn, *Quantum Measurement and Control* (Cambridge University Press, Cambridge, England, 2009).

[34] S. Haroche and J. Raimond, *Exploring the Quantum: Atoms, Cavities, and Photons*, Oxford Graduated Texts, 1st ed. (Oxford University Press, New York, 2006).

[35] L. Bretheau, P. Campagne-Ibarcq, E. Flurin, F. Mallet, and B. Huard, *Quantum Dynamics of an Electromagnetic Mode that Cannot Have N Photons*, Science **348**, 776 (2015).

[36] W. Cai, Y. Ma, W. Wang, C.-L. Zou, and L. Sun, *Bosonic Quantum Error Correction Codes in Superconducting Quantum Circuits*, Fundam. Res. **1**, 50 (2021).

[37] J. M. Gertler, B. Baker, J. Li, S. Shirol, J. Koch, and C. Wang, *Protecting a Bosonic Qubit with Autonomous Quantum Error Correction*, Nature (London) **590**, 243 (2021).

[38] S. Hacohen-Gourgy, L. S. Martin, E. Flurin, V. V. Ramasesh, K. B. Whaley, and I. Siddiqi, *Quantum Dynamics of Simultaneously Measured Non-Commuting Observables*, Nature (London) **538**, 491 (2016).
[39] Q. Ficheux, S. Jezouin, Z. Leghtas, and B. Huard, *Dynamics of a Qubit While Simultaneously Monitoring Its Relaxation and Dephasing*, Nat. Commun. **9**, 1926 (2018).

[40] Q. Ficheux, *Quantum Trajectories with Incompatible Decoherence Channels*, Ph.D. thesis, Universite Paris Sciences et Lettres, 2018.

[41] We find two solutions where $z/\eta$ is stationary, one stable and one unstable. Considering the stable solution only, we compute $R$ such that $\eta(kT + T) = R\eta(kT)$. The computation boils down to solving order-two polynomials.

[42] I. L. Chuang and M. A. Nielsen, *Contemporary Physics* (Cambridge University Press, Cambridge, England, 2010), Vol. 52.

[43] Z. K. Minev, Z. Leghtas, S. O. Mundhada, L. Christakis, I. M. Pop, and M. H. Devoret, *Energy-Participation Quantization of Josephson Circuits*, arXiv:2010.00620.

[44] M. Hatridge, S. Shankar, Mirrahimi, F. Schackert, K. Geerlings, T. Brecht, K. Sliwa, B. Abdo, L. Frunzio, S. M. Girvin, R. J. Schoelkopf, and M. H. Devoret, *Quantum Back-Action of an Individual Variable-Strength Measurement*, Science **339**, 178 (2013).

[45] A. DasGupta, S. N. Lahiri, and J. Stoyanov, *Sharp Fixed n Bounds and Asymptotic Expansions for the Mean and the Median of a Gaussian Sample Maximum*, and Applications to the Donoho–Jin Model, Stat. Methodol. **20**, 40 (2014).

[46] B. Peaudecerf, *Mesure adaptative non destructive du nombre de photons dans une cavité*, Ph.D. thesis, Pierre and Marie Curie University, 2013.

[47] Here, least significant is to be understood as the last bit in the binary decomposition, and not in terms of the amount of information.

[48] J. Johansson, P. Nation, and F. Nori, *QUTIP2: A PYTHON Framework for the Dynamics of Open Quantum Systems*, Comput. Phys. Commun. **184**, 1234 (2013).

[49] N. Cottet, *Energy and Information in Fluorescence with Superconducting Circuits.*, Ph. D. thesis, PSL Research University, 2018.

[50] A. Sarlette, P. Rouchon, A. Essig, Q. Ficheux, and B. Huard, *Quantum Adiabatic Elimination at Arbitrary Order for Photon Number Measurement*, IFAC-PapersOnLine **53**, 250 (2020).