Multiplexed photon number measurement

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 back-action of the measurement as well as the optimal information extraction parameters. Our experiment unleashes the full potential of quantum meters by replacing a sequential quantum measurements with simultaneous and continuous measurements separated in the frequency domain.


I. INTRODUCTION
The most general measurement of a quantum system consists in using a quantum apparatus as a probe. The system interacts with the probe before the latter 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. 1a).
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 one 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? * Equal contribution author We propose to use a qubit as an encoder of information about the cavity state into 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 9 possible photon numbers (more than 3 bits) in a microwave resonator is simultaneously extracted by a single superconducting qubit into 9 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 9 test frequencies by multiplexing, the qubit simultaneously emits 9 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 questions used time division multiplexing while our experiment demonstrates the analogous of frequency division multiplexing, where the qubit alone acts as the frequency multiplexing transducer (Fig. 1a). This experiment directly benefits from the recent bandwidth improvements of near quantum limited amplifiers [11], which enabled us to bring the measurement process in the frequency domain. Current limitations in the cavity lifetime and detector efficiency prevented 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. Besides we manage to observe the multiplexed measurement back-action 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 2 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 fmp − kχs,mp. This multiplexing-demultiplexing process leads to reflection coefficients r k that each encodes the probability that k photons are stored.
on the resonator using direct Wigner tomography, which allowed 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 multi-frequency 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. Sec. III describes photon counting using the fluorescence of a second qubit. Sec. 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. Sec. V presents an actual experiment that implements a version of the gedanken experiment where the photo detectors 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. Sec. VI evidences the measurement back-action 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 offresonant transmon qubits are coupled to the storage mode (see Fig. 1b). 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-VI). Both qubits are dispersively coupled to the resonator so that their frequency respectively redshifts by χ s,yn = 1.4 MHz and χ 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 π-pulse at f drive = f yn − kχ s,yn . The state of the yes-no qubit is read out using a dedicated resonator (Fig. 1b). To demonstrate this photon counting ability, we use a microwave tone at f s to prepare the storage mode in a coherent state |β = e −|β| 2 /2 +∞ n=0 β n √ n! |n , which is a superposition of all Fock states with mean photon numbern = |β| 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 (Fig. 2a,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 β 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. 2b).
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 π-pulse on the qubit. Finding the photon number N , by asking successively whether the resonator is in state |k 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 one 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 adaptative 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 (GRAPE) [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 the order of log 2 (N ) measurement steps that each requires a time of the order of 1/χ s,yn at best (see Appendix F 2). The storage mode is prepared in a coherent state with an average photon number n using a microwave pulse at storage frequency and amplitude Vs. a, b, Measured probability Pe that the yes-no qubit gets excited by a π-pulse at a frequency f drive . Peaks appear at fyn − kχ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 ofn as a function of the drive amplitude Vs.

III. FLUORESCENCE PHOTON COUNTING
The intrinsic limitation of the standard approach is that measuring the qubit state can at most reveal one 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 around the qubit frequency f mp − kχ 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. 1b) [12][13][14][15]. The measured real part 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. 1c) [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 storage mode can be deduced from the relative amplitudes of the reduction of Re(r) at each frequency f mp − lχ s,mp .
In Figs. 3 a,b, we show the measured qubit emission coefficient 1 − Re(r) as a function of a single probe frequency f probe and of the initial amplitude of the storage mode coherent state √n . The measurement is performed using a drive strength Ω = χ s,mp /4 (expressed as the corresponding Rabi frequency) and pulse duration of 2 µs, which is smaller than the storage lifetime of 3.8 µs. Resolved peaks develop for every photon number up to at least 9. Using the former calibration ofn, a masterequation 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 designed the relaxation rate of the multiplexing qubit Γ 1,mp = (42 ns) −1 so that the decoherence rate Γ 2,mp = Γ 1,mp /2 is smaller than the dispersive shift 2πχ s,mp . When peaks are separated, probing the qubit at one of its resonance frequencies f mp −kχ s,mp opens a communication channel with a maximal bandwidth Γ 2,mp carrying information only about Fock state |k . We maximize the bandwidth of each channel by designing Γ 1,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 have shown 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 mp − kχ s,mp of the transmission line. It thus enables the key ingredient of our approach: the multiplexing measurement of reflection at every frequency f mp − kχ s,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 mp − lχ s,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 mp − kχ s,mp with a bandwidth χ s,mp . The protocol proceeds in three steps to count the number of photons in the storage mode starting in state |ψ s , as detailed in Fig. 4. First, the multiplexing qubit is excited with a π-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 wavepacket whose envelope decays at a rate Γ 1,mp . In the limit where Γ 1,mp χ s,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 where |· j represents the quantum state of the propagating mode going through the multiplexer on branch j corresponding to frequencies in the band [f mp − (j + 1/2)χ s,mp , f mp − (j − 1/2)χ s,mp [ (see Fig. 4c). Matching the temporal envelope of the modes to the exponential decay at a rate Γ 1 [17], the mode is occupied by either |0 j or |1 j depending on the storage photon number, hence the notation |δ k,j j . Finally, a single photodetector clicks and reveals the number of photons k with prob- FIG. 4. Gedanken multiplexing experiment. a) An unconditional π pulse is applied to the multiplexing qubit while the cavity is prepared in state |ψ 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 fmp − kχs,mp by a rectangular component displaying the frequency band. d) Eventually a single photodetector (detector k = 1 in the figure) clicks with probability | k|ψ | 2 , allowing us to deduce the photon number. The storage mode is projected on the corresponding Fock state (here |1 s) in a typical time T1,mp that does not depend on the average number of photons in the storage mode.
ability | k|ψ | 2 (Fig. 4d). In case of ideal detectors with zero false positives, the click detects the associated propagating mode in |1 k , and therefore, as the line is entangled with the storage mode, the measurement backaction projects the storage mode in Fock state |k . The total measurement time, a few 1/Γ 1,mp , corresponds to the time it takes for one photodetector to click. The time is thus independent on the number of photons stored in the storage mode. Note that in order to avoid spectral leakage into other ports, Γ 1,mp is limited by χ s,mp so that the shortest measurement time is limited to a few 1/χ 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 )/χ 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 towards this goal [5,[17][18][19][20][21][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 9 frequencies 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 numbern in the storage mode. r k is here the reflection coefficient at fmp − kχs,mp. Solid lines: prediction based on a master equation without free parameters (see Appendix G 3).
corresponding to photon numbers from 0 to 8. We then demultiplex the reflected pulse at the same 9 frequencies {f mp − kχ s,mp } 0≤k≤8 and extract a reflection coefficient r k for each of them (Fig. 1d). The measurement consists in simultaneously measuring the emission coefficients 1 − Re(r k ) for each peak in Fig. 3a,b, which is much faster than measuring them one at a time. the average initial photon number for a drive strength Ω = χ s,mp /2 and a measurement duration of 2 µs. For a givenn, every measurement channel k gives an average signal that is proportional to the probability of having k photons in the storage mode. Asn 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 9 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 χ s,mp which is now possible using a TWPA [11]. It is a comparable technical requirement to the recently demonstrated high-efficiency multiplexed readout of as many as 6 qubits coupled to a single feed line [23][24][25][26]. The number independent measurement time in the case of ideal photodetectors (Sec. IV) relies on the absence of noise when measuring a mode in the vaccum 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 similarly to state-of-theart 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 BACK-ACTION
The measurement strength of the multiplexing measurement 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 could not 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 measurement 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 yesno qubit to perform a direct Wigner tomography [29][30][31] of the storage mode at various times t. It provides a representation of the state ρ in the phase space of the mode and can be expressed as Here D(α) = e αa † s −α * as is the storage displacement operator, P = e iπa † s as 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 |β = −1.55 , the Wigner function starts as a Gaussian distribution centered at α = β. On the left of Fig. 6, one can see how the bare dephasing rate and the self-Kerr effect of the storage mode (0.02 MHz frequency shift per photon) distort the Gaussian distribution towards a torus with no phase when time increases even without any photon counting drives. Using a single drive with Ω = χ s,mp /2 to measure whether there is 1 photon, the phase diffuses faster and the Wigner function exhibits negativities in middle of Fig. 6. As seen in the corresponding density matrix, a tone at f probe = f mp − χ s,mp notably induces dephasing between states |1 and all other states |m = 1 (see density matrix as a function of 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 Ω (right of Fig. 6). Likewise all offdiagonal elements of the density matrix are then reduced.
To be more quantitative, the dephasing rate Γ d,s of the cavity is accessed through the decay of the mean quadrature of the storage mode [32] (see Appendix B 3) a s + a † s = 2xW (x + iy)dxdy.
In Fig. 7a, we show a s + a † s as a function of time under the multiplexed drive with a strength Ω = 0.23χ s,mp . Repeating this experiment for various values of the multiplexed driving strength Ω allows us to determine how the latter affects the dephasing rate Γ d,s , and thus the measurement rate. The dephasing rate is non monotonic in the drive strength (Fig. 7b). Noticeably, it reaches a maximum when Ω = χ s,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 mp + kχ s,mp } k∈Z (see Appendix C). The Fourier transform of a comb being a comb, the drive performs sudden rotations by an angle 2πΩ/χ s,mp of the Bloch vector of the qubit every time step 1/χ s,mp . When Ω/χ s,mp is integer, the comb does not affect the qubit and thus Γ d,s vanishes. Conversely, the maximum measurement rate corresponds to half-integer Ω/χ s,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 master equation without any free parameter (line in Fig. 7b).

VII. CONCLUSION
We have experimentally demonstrated 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 η, the dispersive shift χ s,mp and the coupling rate Γ 1,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 ηΓ 1,mp /Γ 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], to stabilize a Fock state parity by autonomous feedback [37], or to extend the reach of simultaneous probing of a single quantum system by multiple observers [38,39] to larger systems and arbitrarily many observers. Our photocounter for stationnary 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. The circuit is composed of 4 electromagnetic modes whose parameters can be found in Methods. A high-Q harmonic oscillator, called storage mode, is composed of a λ/2 coplanar waveguide (CPW) 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 Γ 1,mp = (44 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 µ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χ s,yn 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 could not be associated to 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 Γ 2,mp < 2πχ s,mp . This resolution constraint is not critical, as in fact a finite amount of photon number information can be extracted as soon as χ s,mp is nonzero, but the decoding is much simpler if we can reason in terms of well-separated resonance peaks.

Device fabrication.
The length of the readout resonator and storage mode was designed to obtain the resonant frequencies f ro ∼ 7 GHz and f s ∼ 4.5 GHz. The circuit consists of a sputtered 120 nm-thick Niobium film deposited on a 280 µmthick undoped silicon wafer. The resonators and feed lines 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/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.

Measurement setup
The readout resonator, the yes-no qubit and the multiplexing qubit are driven by pulses that are gener-ated using a Tektronix® Arbitrary Waveform Generator (AWG) AWG5014C with a sample rate of 1 GS/s. Storage mode pulses are generated using a Zurich In-struments® 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 yesno qubit, and 75 MHz for storage and multiplexing qubit. They are up-converted using single sideband mixers for readout resonator and multiplexing qubit and regular mixers for the storage resonator and yes-no qubit, with continuous microwave tones produced respectively by AnaPico® APSIN12G, Agilent® E8257D, Wind-Freak® 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 Travelling Wave Parametric Amplifier (TWPA) provided by Lincoln Labs. We tuned 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 was measured to be 18.7 ± 0.4%, and should be close to the efficiency η 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 HEMT amplifier. The follow-up amplification is performed by a High Electron Mobility Transistor (HEMT) amplifier from Low Noise Factory (LNF®) 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 signal, nine demodulation operations are performed at each of the down-converted frequencies 75 MHz + kχ s,mp for 0 ≤ k ≤ 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χ s,mp } 0≤k≤8 are summed and multiplied by a Gaussian envelop numerically with a sampling rate of 1 GHz/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 side band mixer whose LO port is driven at frequency f mp + 75 MHz.

Circuit parameters and master equation.
All parameters of the 4 modes can be measured using standard circuit-QED measurement (see Table I). Frequencies of the readout mode and multiplexing qubit are measured by spectroscopy. Frequencies of storage mode and yes-no qubit are measured using two-tone spectroscopy with the readout mode. Yes-no qubit decay and decoherence rate are measured with the time evolution of the probability to find the qubit excited after a π pulse and using Ramsey oscillations. Readout mode decay rate and cross-Kerr rate between readout mode and yes-no qubit are measured using the measurement in-duced dephasing rate by the readout mode on the yes-no qubit. Cross-Kerr rate between the storage mode and the two qubits are measured using qubit spectroscopy with the storage state initialized in various coherent states. Anharmonicities are measured using spectroscopy of the qubit excited state. Decay and decoherence rates of the storage mode are measured with the time evolution of the probability to have 0 photon in the storage mode after a displacement and storage Ramsey interferometry experiment. 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 wheren ro ,n s ,n yn , andn mp are the photon number operators respectively for the readout, storage, yes-no qubit and multiplexing qubit. χ a,b is the cross-Kerr rate between modes a and b. χ a,a is the anharmonicity of the mode a. The master equation on the system density matrix ρ readṡ where L is the Lindblad superoperator defined as L(L)ρ =LρL † − L †L , ρ /2 andâ b is the annihilation operator of mode b. For a qubit mode b, the dephasing rate Γ φ,b is linked to the decoherence rate by Appendix B: Calibrations

Calibration of the storage mode displacement amplitude
The storage mode can be displaced by driving it on res- is the pulse envelope. The driving Hamiltonian of the storage mode reads ( The scaling factor µ = 1.45 (mV µs) −1 is calibrated by fitting the photocounting measurement results obtained using the yes-no qubit with the master equation simulation (see section G 1). Fig. 10a shows the evolution of s with V s . For every experiment, the storage mode displacements are realized using a Gaussian pulse shape s (t) = λ(t) max with a maximum amplitude max , a width 25 ns and a duration 100 ns. We simulated the dynamics of the storage mode under this Gaussian displacement taking into account the couplings, relaxation and decoherence rates (see section G 2) for various amplitudes max . We then computed the expectation value of the number of photon operator n s at the end of the pulse. Fig. 10b shows the square root of n s as a function of max . Fitting with a linear function, we find that n s = 59.1 max . As s increases linearly with V s , max increases linearly with the maximum voltage amplitude V max,s of the Gaussian pulse V s (t) = λ(t)V max,s . Using the two linear regressions, we can express the photon number of the storage mode as n s = (85.9 V −1 ) V max,s .

Wigner tomography calibration
The Wigner function of a harmonic oscillator with density matrix ρ is defined as W(α) = 2 Tr(D † (α)ρD(α)P)/π where D(α) = e αâ † s −α * â s is the displacement operator of the storage mode by a coherent field α and P = e iπâ † sâs is the photon number parity operator. A Wigner function is the expectation value of the parity after a displacement by an amplitude α. The Wigner tomography sequence is represented on Fig. 11a . 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 section B 3) with a width of 25 ns. Then two successive π/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 ∆τ = 337 ns ≈ 1/2χ 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][30][31]. The sequence terminates by a 2 µs long square pulse on the readout resonator to read out the state of the yes-no qubit. For high amplitude α, higher order Kerr terms distort the Wigner function. To mitigate this effect, we interleave two sequences with a final pulse of phase either +π/2 or −π/2.(see Fig. 11a). 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. 11b is calibrated using the yes-no qubit Rabi oscillation amplitude to express the signal using Pauli operators. Multiplying the result by 2/π yields the Wigner function W(α) in Fig. 11c.
The axes of the phase space x, p are calibrated using the same pulse sequence. The photon number calibration realized before (see section 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 the calibration when the storage mode hosts a large number of photons. We decided 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 th the Wigner function is a 2D Gaussian function with a width n th + 1/2 [34] For a thermal state displaced by an amplitude β the Wigner function is still a 2D Gaussian function with a width n th + 1/2 but centered on β. In thermal equilibrium, the storage mode has an average photon number n th = 0.03, which is measured using the photon counting experiment. We calibrated the quadrature axes in order to get the expected geometrical mean √ σ x σ 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 √ σ x σ p = 0.53. The function used for the calibra- . 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 π/2 pulse is applied to the yes-no qubit. The qubit evolves freely during a time ∆τ = 337 ns ≈ 1/2χs,yn before a new ±π/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 drive amplitude for various detuning using only the photon number calibration (see section 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 was 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 fmp − 1.4χs,mp during 750 ns with an amplitude Ω = χs,mp/2. The appearance of negative values in the Wigner function demonstrate that one can prepare non-classical states in the storage mode using the multiplexing qubit backaction alone.
tion is a third order polynomial function which gives |α| as a function of the pulse amplitude V max,s . We repeat this protocol for 3 detuning values δf s between the displacement pulse and storage mode frequencies. Fig. 11b shows the mean quadrature spread of the displaced stor- 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 π/2 pulse is applied to the yes-no qubit. We then let the qubit evolve freely during a time t before doing a new ±π/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/χs,yn ≈ 0.7µs.
age mode thermal state Wigner function as a function of drive amplitude V 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 s is expressed in Volts and φ s is the phase of the pulse. For a typical value V 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 ∆τ is calibrated using qubit state revival during Ramsey interferometry (see supplementary material of Ref. [35]). We used a Ramsey interferometry sequence (Fig. 12a) for the yes-no qubit at its resonance frequency for various coherent states in the storage mode. Revivals happen every 1/χ s,yn which allows us to set ∆τ as half the revival time in Fig. 12b. The signal difference between the final −π/2 and +π/2 pulses can be expressed as This expression is derived in the supplementary material of Ref. [35]. The last exponential decay factor was 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 By adjusting the above parameters to allow the model to match the measured signal shown in Fig. 12b, we find ∆τ = 337 ns, γ = 0.23 µs −1 and χ s,s,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.

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Ô, one can apply the Wigner transform to obtain the operator Wigner map WÔ [34] as where {|x } is the eigenbasis of the quadrature opera-torX. With this tool, the Wigner function of a state |Ψ (respectively a density matrix ρ) is simply given by W |Ψ Ψ| (α) (respectively W ρ (α)). The mean value of an operatorÔ can be derived from the integral over the phase-space of the product of the two Wigner distributions multiplied by π, π dx dp W ρ (x, p)WÔ(x, p) = 1 π dx dp dy dy e −2ip(y+y ) x + y/2|ρ|x − y/2 x + y /2|Ô|x − y /2 = dx dy dy δ(y + y ) x + y/2|ρ|x − y/2 x + y /2|Ô|x − y /2 = dx dy x + y/2|ρ|x − y/2 x − y/2|Ô|x + y/2 = du dv u|ρ|v v|Ô|u = Tr(ρÔ) = Ô ρ . (B5) In the case ofX andP operators, Wigner maps take a simple expression For any density matrix ρ, one can extract X = Tr(Xρ) and P = Tr(P ρ) from the Wigner function W ≡ W ρ as 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 and |e states. A typical sequence starts by a π/2 pulse detuned from resonance by δf to create a coherent superposition of |g and |e states. Then the qubit evolves freely before its state tomography. Both σ x and σ y oscillate at δf while decaying at the decoherence rate Γ 2 .
We decided 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 π/2 pulse is replaced by a detuned displacement pulse D(β) on the storage mode. The field then evolves during a time t (during which the multiplexing measurement could be applied) before a Wigner tomography is realized (see Fig. 13a). The expectation value ofX = (â s +â † s )/2 andP = (â s −â † s )/2i quadratures are computed from the Wigner tomography (see Sec. B 3). The time trace of X and P 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 δ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 δf 0 s = f 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 δ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 δf s for its displacement pulse D † (α) 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 X and P evolve as where β = |β|e iφ = â s (t = 0). For each time t, we computed X and P and defined the storage mode dephasing rate as Γ d,s which contains the intrinsic decoherence rate Γ 2,s . Data Extended Fig. 13b shows an example of measured Ramsey oscillations. In the main text, Fig. 3B does not exhibit oscillations because it is the mean value â s +â † s in the frame rotating at the resonant frequency of the storage mode. In practice, we measured them with a detuning and numerically computed the non-oscillating quantity 2Re(( X + i P )exp(−2iπδf s t)).

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. 14a). The drive amplitude is given by the sum of nine sine functions at the frequencies [f mp , f mp − χ s,mp , ..., f mp − 8χ s,mp ] multiplied by a Gaussian envelope of duration t and width t/4.
For small measurement strength Ω/χ s,mp < 0.9, we generated the Ramsey sequence with a displacement pulse of amplitude β = −1.55 detuned from the base storage mode by δf 0 s = 3.96 MHz. We fit the time evolution of X and P using the damped sine function The parameters A, δf s , φ, and Γ d,s are determined altogether by fitting the model to the measured oscillations. Γ d,s is the sum of the intrinsic storage dephasing rate Γ 2,s For larger measurement strength Ω/χ s,mp > 0.9, we generated the Ramsey sequence with a displacement pulse detuning of δf 0 s = 5.96 MHz, an amplitude of β = −1.27, and we model the time evolution of X and P by the sum of two sine functions with an exponential decay (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. Fig. 14b 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, δf s , ν, φ, ψ X , ψ P , and Γ d,s . The frequency ν varies from 2.15 MHz to 2.5 MHz. The parameter ζ is roughly constant, it varies between 0.2 to 0.27. Fig. 14c shows measurement induced detuning as a function of measurement drive amplitude.

Rabi frequency calibration
We observe Rabi oscillations of the multiplexing qubit by applying a 1 µs-long square pulse at f mp with a varying amplitude V mp . The reflected signal, demodulated by time steps of 10 ns, displays damped oscillations given by [40] Re (r(t)) − Re (r ss ) = where r ss is the value of the reflection coefficient in the steady state. We obtain ξ = 0.543 GHz.V −1 so that the Rabi frequency is calibrated as Ω = ξV mp = (0.543 GHz.V −1 )V mp (see Fig. 15).

Appendix C: Maximal measurement strength
The dephasing rate Γ d,s depends on the driving amplitude Ω, the dispersive shift χ s,mp and on the relaxation rate Γ 1,mp . To gain insight in 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.

a. Hamiltonian evolution
We consider a single qubit driven by a frequency comb with 2p + 1 frequency peaks at every f mp + kχ for −p ≤ k ≤ p. In the frame rotating at the qubit frequency, the Hamiltonian reads 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 where sinc(x) = sin(x)/(x). For large integers p, we can approximate the sum as +∞ k=−∞ sinc(π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 < χ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. With this we get where x 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 Figure 16. To put it simply, the action of the comb consists in 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/χ, the qubit rotates almost instantaneously by an angle 2π Ω χ . Without decoherence, if the qubit starts in state |g at time t 0 , the qubit state after a time t reads (C4) Let us focus on some particular values of Ω χ .
• If Ω χ is integer, the staircase approximation with f (t) keeps the qubit in |g at all times, just performing a full rotation on the Bloch sphere at each Rabi pulse. In presence of relaxation, a photon loss can only happen during the short duration of the Rabi pulse, which decreases as 1/(p+1). The qubit remaining in the ground state cannot encode information on the resonator. This explains why there are minima of measurement induced dephasing rate when Ω/χ is integer (see Fig. 14d). Exact rotation angle f (t) in blue and its staircase approximationf (t) in red. The duration of each step is equal to 1/χ and its height to 2πΩ χ . 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.
• If Ω χ is half-integer, the staircase approximation with f =f makes the qubit state jump periodically between |g and |e . 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) ≈f (t). This allows us to integrate the qubit dynamics exactly. The continuous photon decay at rate Γ 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 χ 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 The average qubit excited state population over one period in the permanent regime reads (1 +z(k T )) 2 .
(C7) When using the parametrization c ≡ cos(θ), one easily checks that P is a strictly decreasing function on c ∈ [−1, 1]. As a function of θ, it has maxima for θ = (2k + 1)π and minima for θ = 2kπ. The latter are no surprise and give P = 0 as the Rabi pulse takes the state from |g back to |g ; for 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(Γ1T /2)

Γ1T
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 Γ 1 = tanh(Γ 1 T /2)/T with the optimal choice of Ω = χ/2 + kχ, where k is integer.
• At fixed T = 1/χ, the emission rate increases with Γ 1 and converges towards χ. For Γ 1 χ, 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.
• Likewise, for a fixed Γ 1 /χ (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 χ.
• For a fixed Γ 1 , the largest average emission rate is obtained for χ = 1/T as large as possible, but it saturates at P Γ 1 = Γ 1 /2. This is consistent with the fact that Γ 1 is a hard limit on the photon emission rate.

Dynamics of the qubit-resonator bipartite system driven by comb
For the sake of simplicity, we consider a lossless storage mode in this section. In this case, each 2 by 2 submatrix ρ n1,n2 = n 1 |ρ|n 2 , with |n i the n i Fock state, evolves independently of the others ρ m1,m2 similarly to a collection of qubit-like system. The sub-matrix ρ n1,n2 is non-normalized because it is an off diagonal sub-matrix of the storage-qubit system. We have where H(t) is the frequency comb drive from the previous section.
a. Computing the decoherence rate of the n1, n2 component The infinite comb approximation again helps. We view H(t) as applying a Rabi pulse of angle θ = 2πΩ χ at each period T = 1/χ, without any effect for the rest of the time. Over one period, we thus have where K 0 applies the Rabi pulse, while K 1 contains dynamics associated to the dispersive coupling and to qubit decay. During each period between Rabi jumps, denoting ρ n1,n2 = xσx+yσy+zσz+ηI 2 , the dynamics K 1 corresponds to the integration of the set of equations After one period T , since the peaks in the comb are exactly separated by the dispersive shift χ, the effect of the precession at a frequency 2πχ(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πχ Γ 1 . We then obtain, in the above coordinates, the matrix expression . Besides, the Rabi rotation corresponds to This expression further simplifies for the two drive strengths Ω = χ/2 + kχ or Ω = kχ, which respectively lead to θ = 0 (no qubit emission) and θ = π (maximal qubit emission), since the (z, η) variables of interest decouple from (x, y).
• For θ = 0, we tend to a stationary regime (z, η) kT +T = (z, η) kT . Note that these are the values associated to 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 Γ d,s in Fig. 14d).
• For θ = π, we can compute an analytical expression for the factor R by which the trace η decreases every period T in the permanent regime [41]. Thus the average decay rate is given by − log(|R|)/T with and we recall T = 1/χ . We now analyze this last expression.
b. Optimal decoherence rate of the n1, n2 component The decoherence rate − log(|R|)χ can be calculated for any value of the dispersive shift χ 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 Ω = χ/2 (maximal measurement strength). As expected, the decoherence is stronger when the photon numbers are further apart. Besides, for 2πχ|n 1 − n 2 | Γ 1 the rate saturates to Γ 1 /2 similarly to the emission rate of the single qubit.

Upper bound on the total measurement rate as a function of photon number
Beyond the determination of the decoherence rate between two Fock states, we are interested in the maximal information extraction rate about the storage state in the multiplexed measurement scheme. In particular, we will discuss how this maximal total measurement rate depends on the maximum number of photon N max that is probed by the 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 heterodyne measurement on each peak of the comb. In Sec. IV, we propose such a measurement apparatus.
We assume the number resolved regime 2πχ Γ 2 , Γ 1 . Thus the decoherence rate between two Fock states is independent of the Fock state numbers and is equal to Γ 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 only consider its effect, which is the entanglement operation. To each Fock state |n in the storage mode (0 ≤ n ≤ N max ) is associated one out of the N max + 1 modes of the transmission line at f mp − nχ. 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 of reference frame by displacing the outgoing modesâ out,n by the opposite of the input coherent state, a single mode will be excited and all the non resonant 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 m correspond to all modes in the vacuum except the one at frequency f mp − nχ and | ⊥ 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 s } 0≤n≤Nmax . The measurement mode has the N max + 2 states discussed above. The bipartite system starts in the state (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 Γ 1 /2, one can write the state of the bipartite system as 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 Γ 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 where S(ρ) = −Tr(ρ log 2 (ρ)) is the Von Neumann entropy, ρ s,m is the bipartite density matrix and ρ s ( respectively ρ m ) is the density matrix of the storage mode ( resp. measurement mode) obtained by tracing out on the other mode. If the state of the bipartite system is pure then its entropy of 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 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 ρ 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 a initially uniform photon number distribution, such as when ψ n = (N max + 1) −1/2 for all n.
with r = e −Γ1t/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 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 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/Γ 1 . Because the mutual information per bit is always bigger than the 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 . Fig. 18b shows the mutual information per bit as a function of 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 to |f transition at zero photon in the resonator becomes resonant with the qubit |g to |e transition for χ mp,mp /χ s,mp photons, which complicates the analysis. In our experiment, it would set a limit to 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 χ s,mp when the photon number increases. Beyond some critical photon number n crit , the dispersive shift χ 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 where χ s,s,mp is the decrease of the cross Kerr rate χ 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 which leads to and a upper bound of 8E J /(p s hf s ) ≈ 10 3 photons for E J /h ≈ 15 GHz, p s ≈ 10 −2 and f s ≈ 5 GHz. Therefore, we do not believe that these limitations are hard enough to prevent the multiplexing scheme to count 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χ s,mp in the comb and record a complex amplitude I We assume that each of these measurement records only extracts the information on the occupation of the Fock state |k , which is experimentally valid in the limit 2πχ s,mp Γ 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 or on the complementary subspace Π In practice, the measurement proceeds by first entangling the resonator, which is in a state |ψ , and the signal mode of the phase preserving amplifier. For the measurement channel k, the entangled state reads and states denoted as |α, β 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 α, and the case where it is off resonant leading to a reflected amplitude α ⊥ . The resonance frequency of the qubit depends on the number of photons in the resonator so that the reflected amplitude α indicates k photons while α ⊥ indicates that there are not k photons. For an incoming amplitude α in onto the multiplexing qubit, we get where σ −,mp 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 | σ −,mp ss | is reached for 2πΩ = Γ 1,mp / √ 2. However, in the case of a qubit driven by an infinite frequency comb, the time average of the lowering operator is with 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 state thus reads 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 supplementary information of Ref. [44], in the case of a phase preserving amplifier the inner product ξ(β, I m , Q m ) = Ψ Im,Qm |β, 0 is given up to a global phase factor (independent on β, I m and Q m ) by where σ 0 is the amplitude of the zero-point fluctuations (the variance of the measured I m is 2σ 2 0 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 Γ 2,mp 2πχ s,mp

(E5)
Appendix F: How fast can the multiplexed measurement be?
In principle, multiplexing the measurement enables 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 with 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 showed 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.

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 photodetector array by signal processing of propagating microwave modes.
It is out of 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 in determining which channel contains an amplitude α while all the others contain an amplitude α ⊥ (see Sec. E). The measurement records {I m } k are averaged along a duration t. Thus the time averaged intrinsic noise contained in the measurement records scales as 1/ √ t and the time average value is independent on t. The problem can be mapped onto the following game. N max stochastic variables {u i } 1≤i≤Nmax are each randomly chosen using a Gaussian distribution centered on 0 with a width 1/ √ t. Another stochastic variable u 0 is randomly chosen using a Gaussian distribution centered on 1 with a width 1/ √ t. The list {u i } 0≤i≤Nmax is scrambled randomly into a list l and the goal consists in 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≤i≤Nmax are higher than u 0 We can rescale all the distribution by √ t, thus u 0 are chosen randomly using a Gaussian distribution centered on √ t with a width of 1 and each of the {u i } 1≤i≤Nmax 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≤i≤Nmax tends towards 2 ln(N max ) [45].
Besides, the median of the max of {u i } 1≤i≤Nmax is equal to the mean value within an error scaling as 1/ √ N max : (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≤i≤Nmax , it leads to From this expression, we understand that the measurement time for a fixed error probability scales as log(N max ).

Comparison between 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 χ between the storage mode and the qubit to be greater than the decoherence rate of the qubit Γ 2 . The goal is to measure the photon number N assuming it is smaller than N max .

Sequential brute force
The brute force approach consists in 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, ..., we apply a photon number conditional π pulse to the qubit at frequency f q − kχ 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 π + T ro )(N + 1). The time T π is the time of an conditional π pulse, hence it is at least 1/χ, while the qubit readout time T ro is limited by other parameters in order to get a single shot readout.

Passive photon number decimation using weak measurement
This approach, which was implemented with Rydberg atoms in cavity [1], consists in encoding the photon number in the phase of the qubit by waiting a time 1/(N max χ) after the qubit has been prepared in state (|g + |e )/ √ 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 π/2 pulse on the qubit with a phase 2π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 σ = N max /( √ pπ) (see appendix A in Ref. [46]).
Therefore, the required number of repetitions to get a fixed error probability on the photon number scales as p ∝ N 2 max . Since each measurement takes at least 1/(N max χ), the total measurement time scales at least as N max /χ.

Active photon number decimation
The previous protocol can be improved by optimizing the phase of the final π/2 pulse to maximize the amount of information extracted on the cavity photon number. It was realized in Ref. [3] using Rydberg atoms in cavity. Because of the use of feedback on a weak measurement, we could not find a closed form for the measurement time in this case [46]. However it was shown that the total time is larger than the total time taken by a binary decimation with feedback (see below).
Binary decimation with feedback This method was shown to provide the least number of steps for sequential photocounting [4]. Each step consists in applying an unconditional π/2 pulse to the qubit, wait a time 1/2 k+1 χ, apply a new unconditional π/2 pulse with a phase φ k that encodes the least significant [47] k th bit b k of the photon number N = k b k 2 k into the qubit state. Importantly, the phase φ 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 qubit and cavity and of the total feedback latency. The total interaction time is bounded by p 1/(2 p+1 χ) = 1/χ. However the feedback latency scales as p and can be written as T fb log 2 (N max + 1).
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 in 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 only be performed in a time of the order of the dispersive interaction time 1/χ. It leads to a total measurement time scaling as log 2 (N max +1)(1/χ+T reset ) where T reset is the duration of the active reset protocol.
In the table, 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 an heterodyne detection setup instead. The scaling of the measurement time is then in ln(N max ) as in state-of-the-art sequential measurements. Besides, we deport the complexity of optimal control or feedback into the challenge to reach large measurement efficiencies on a large frequency band (many χ's).

Appendix G: Master Equation simulations
In this section, we briefly describe the master equation simulations used to understand our experimental results. We simulated 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 were performed using Python package QuTiP [48]. We simulated 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 only took into account the storage mode and the multiplexing qubit. The storage mode was modeled as an harmonic oscillator while the transmons were replaced by two level systems. The Hilbert space of the storage mode was truncated at a photon number ranging from 10 to 25 photons depending on the simulation. In this section we will use Pauli matrices to describe operators acting on qubits.

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 π pulse on the yes-no qubit at frequency f yn − δf yn before detecting the expectation value of the Pauli operator σ z,yn .
We write the Hamiltonian of the system in a frame rotating at f s − χ s,mp /2 − χ s,yn /2 for storage mode, f yn − δf yn for yes-no qubit mode and f mp for multiplexing qubit mode as followŝ , (G1) where λ(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 s (t) = λ(t) max and yn (t) is the time envelope of a Gaussian pulse with duration 1.9 µs and width 475 ns. The amplitude of the pulse is chosen to obtain a π rotation on the yes-no qubit. The term −δf ynσ z,yn 2 takes into account the detuning between the π pulse and the yes-no qubit frequency. yn (t) is delayed with respect to λ(t) to match the experimental pulse sequence. In comparison with Hamiltonian (Eq. (1) in Methods), 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.
The master equation is solved using the function "mesolve" of QuTiP starting from a thermal state with n th,s photon in storage mode, the yes-no qubit in the ground state |g and the multiplexing qubit also in the ground state |g . The solver iteratively computes the density matrix with a 10 ns time step during the displacement pulse and the π pulse. We compute the expectation value σ z,yn at the end of the sequence and convert it into a probability P e of finding the yes-no qubit in the |e state.
This simulation can be used to reproduce the experiment in Fig. 2A,B of the main text by adjusting the following parameters {µ = max /V max,s , χ s,yn , χ s,s , χ s,s,yn , n th,s }. Note that we need to run the simulation for every couple of parameters (V max,s ,δf yn ). Table III compiles the values of fitted parameters.
b. Photocounting with the multiplexing qubit A second simulation (fluorescence simulation) was carried out to compare the photon counting experiment in Fig. 3A,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 − δf mp with an amplitude expressed as a Rabi frequency Ω = χ s,mp /4. The measured reflection coefficient of the multiplexing qubit r(δf mp ) can be expressed using input-output theory as [49] r(δf mp ) = aout ain = ain − √ Γ1,mp σ−,mp ain And since the Rabi frequency is given by Ω = Γ 1,mp | a in |/π we get an emission coefficient in the frame rotating at f mp − δf mp for the multiplexing qubit. If we set the phase of the drive so that i a in ≥ 0, meaning we drive the qubit along σ x,mp , the emission coefficient becomes The Hamiltonian of the problem in the frame rotating at f s − χ s,mp /2 − χ s,yn /2 for storage mode, f yn for yes-no qubit and f mp − δf mp for multiplexing qubit readŝ , (G4) where mp (t) ≥ 0 is a Gaussian function of duration 2 µs, width 250 ns and amplitude 1. mp (t) is delayed compare to λ(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 were measured independently. The resulting master equation only differs 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 "mesolve" 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 . Finally, the expectation value σ y,mp is computed and integrated during the 2 µs of the pulse.
We compare the measured emission coefficient in Fig. 2B,D to the simulated signal A σ y,mp where A is left as a free parameter due to a small parasitic reflection in the measurement setup and thermal population. The parameters {µ, χ s,t , χ s,s , χ s,s,t , n th,s } is already set by the calibration above using the simulation of the yes-no qubit. From the fluorescence simulation, we thus extract the parameters {χ s,mp , χ s,s,mp , A} by comparing the experimental observations in Fig. 2B,D with the simulation for various V max,s and δf mp . Fitted values are given in Tab. III. Finally, we ran 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. Parameters extracted from the photocounting simulations using the multiplexing or yesno qubit. All parameters except those related to the multiplexing qubit are determined using a fit of the yes-no qubit simulation to the Fig. 2A,C. Parameters related to multiplexing qubit are obtained using a fit of the simulation to the Fig. 2B,D.

Evolution of the average photon number in the storage mode
We simulated the filling of the storage mode by a displacement pulse on the resonator. We simulated the same master equation used for the photocounting simulations with parameters obtained from the photocounting simulations (see Tab. III) but without applying any drive on the qubits. Only the displacement pulse on the storage mode is modeled i.e. mp (t) = 0, δf mp = 0, yn (t) = 0, and δf yn = 0.
The "mesolve" 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. Fig. 10 shows the square root of the expected mean photon number as a function of the amplitude max . We obtain a scaling factor n s = 85.9 V −1 V max,s used in the photon number calibration of the storage mode.

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 − χ s,mp /2 − χ s,yn /2 for the storage mode and at the qubit frequencies for the qubits aŝ ( max e iπ(χs,mp+χs,yn)tâ s + * max e −iπ(χs,mp+χs,yn)tâ † s ) , (G5) where Ω = χ s,mp /2 and comb (t) is the product of a Gaussian function with the sum of nine complex exponential 8 k=0 exp(2iπχ s,mp kt). The Gaussian envelope of comb (t) has a duration of 2 µs, a width of 250 ns, and a maximum amplitude of 1 and the delay between comb (t) and λ(t) reproduces the experimental sequence. The master equation (G2) is used with a time step of 1 ns for various amplitude max . We obtain the time evolution of σ y,mp enabling us to compare the experimental measurements of Fig. 2E to the model. To do so, we integrate the simulated function σ y,mp × cos(2πχ s,mp kt) for each integer k, similarly 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 9 multiplexing readout signals as a function of the mean number of photon in the storage mode used in Fig. 2E.

Simulation of measurement induced dephasing on the storage mode
In this part, we only simulate 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 + δf 0 s for the storage mode iŝ where comb (t) is the product of a Gaussian function with the sum of nine complex exponential 8 k=0 exp(2iπχ s,mp kτ ). 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 equatioṅ The storage is initialized in a coherent state of amplitude β = 1.55 and the multiplexing qubit is initialized in state |g . We simulate the dynamics of the system for a pulse duration t going from 100 ns to 5 µs and for Ω ranging from 0 to 2χ s,mp . We compute the expectation value of X = (â s +â † s )/2 at the end of each simulation. For a given Ω, we extract the time evolution of X under the influence of the multiplexed measurement as shown on Fig. 19a. This decaying sinusoid is fitted using Eq. (B9) to obtain the oscillation frequency δf s and the decay rate Γ d,s . Fig. 14c and d show the measurement induced dephasing and ac-Stark shift as a function of amplitude of the comb Ω for two sets of coherent state amplitudes β and detuning δf 0 s . We identify three interesting features. The first one is the evolution of the shape of the curves δf s (Ω) and Γ d,s (Ω) with χ s,mp . We repeat the simulation using a square pulse envelope instead of Gaussian pulse for comb to make the simulation faster for several values of χ s,mp from 1.5 to 13.2 MHz by steps of 1.4 MHz. We observe that δf s (Ω) and Γ d,s (Ω) increase as χ s,mp becomes larger but that the maxima and minima of the curve are always found for the same Ω/χ s,mp ratio (Fig. 19b) as predicted (see Sec. C).
The second observation is that δf s (Ω) and Γ d,s (Ω) vary with the initial coherent state amplitude β (Fig. 19c). As the frequency comb contains only a finite number of frequencies f mp − nχ s,mp with 0 ≤ k ≤ 8, the decoherence rate between two Fock states, |i and |j , 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 χ s,mp Γ 1,mp /2π, the dephasing rate and ac-Stark shift are a function of the ratio 2πΩ/Γ 1,mp as shown on Fig. 19d. The dephasing rate increases as Ω 2 until a plateau is reached for 2πΩ/Γ 1,mp = 0.7, this plateau is predicted (see Sec. C). In contrast, the Stark shift is constant for 2πΩ/Γ 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πΩ/Γ 1,mp . Since there are two frequencies, we use Eq. (B10) to fit the simulated Ramsey oscillations for χ s,mp Γ 1,mp /2π. 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. Fig. 19e shows an example of simulated Ramsey oscillations for χ s,mp Γ 1,mp /2π. Simulations show the same pattern with maxima and minima for some specific values of Ω/χs,mp as in the experiment in Fig. 3C. c. Simulated measurement induced dephasing rate and ac-Stark shift as a function of Ω/χs,mp for various initial coherent state amplitudes β in the storage mode. For Ω/χs,mp > 1, we see a difference of about 20 % between β = 1.6 and β = 1.2. d. Simulations for χs,mp Γ1,mp. The evolution of the measurement induced dephasing and ac-Stark shift with Ω 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πΩ/Γ1,mp. The red line is a guide for eyes representing a quadratic function. It shows that the measurement induced dephasing rate increases linearly with Ω 2 for small drive amplitudes. For the ac-Stark shift, on the contrary with b, the detuning is constant at the small drive amplitudes, then two frequencies appear with comparable contributions to the Ramsey oscillations. The two frequencies evolve linearly with Ω. e. Example of simulated Ramsey oscillations exhibiting two frequencies.

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 used to produce the bottom part of Fig. 3A 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 non-demolition 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.

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 m| with |n and |m two fock states with n and m photons. The mean value of those operators is equal to the (n, m) element ρ nm of the density matrix. Using the mathematical expression of x|n with H n (x) = (−1) n e x 2 d n dx n e −x 2 the Hermite polynomial function of order n and |x the eigenvector of the quadrature (â s +â † s )/2 associated to the eigenvalue x. The Wigner map of the operator |n and |m becomes W |n m| (x, p) = 1 π dye −2ipy ψ n (x + y/2)ψ m (x − y/2) (H2) and the matrix element ρ nm of the storage mode is given by Eq. (B5)

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 back-action, the quantity |ρ nm |/ √ ρ nn ρ 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 ṡ From this equation, we can compute the time derivative of the density matrix elemenṫ (H5) If the storage mode is initialized in a coherent state |α o , the solution of the equation is and we get Thus indeed, the renormalization removes the effect of the relaxation rate Γ 1,s and only characterizes the dephasing rate. Under the action of measurement, the dephasing rate Γ φ,s is increased by the measurement induced dephasing rate. We will use this property in the following sections and study the evolution of ρ mn / √ ρ nn ρ mm .

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 and Fock state |m = n . We evidence this dephasing by observing the evolution of ρ nm We prepare the storage mode state in a coherent state with an amplitude β = −1.7, and probe the multiplexing qubit during a time t with a drive at the frequency f mp − ∆ mp before doing a Wigner tomography. For various times t and detunings ∆ mp , we compute the density matrix of the storage mode using Eq. (H3). One can fit the time evolution of |ρ nm |/ √ ρ nn ρ mm with a decreasing exponential function. The extracted decoherence rate Γ nm d,s (∆ 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 state |n and |m due to the measurement is given by the highest eigenvalue, which are all negatives, of the following matrix Fig. 20 shows the measured density matrix decoherence rates Γ nm d,s 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πχ s,mp |m − n| > Γ 2,mp ), the decoherence rate Γ nm d,s between Fock states |n and |m is larger when the single drive probes whether there are n photons or m photons with a moderate drive amplitude Ω (dependence on Ω not shown here). For much larger drive amplitude Ω, one can increase the decoherence rate Γ nm d,s by driving with a detuning ∆ mp = (n + m)χ s,mp /2, similarly to dispersive qubit readout which is optimal for information extraction at large drive power and for a drive frequency detuned by χ s,mp /2. In fact this regime would become particularly attractive for poorly resolved resonances as a function of photon number (2πχ s,mp |m−n| < Γ 2,mp ). Premises of this effect are visible on Fig. 20, as the maximal decoherence rate occurs at a detuning slightly closer to (n + m)χ s,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.

Multiplexed measurement vs single tone measurement
In Fig. 3A 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 β = −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 mp − χ s,mp corresponding to 1 photon or (iii) with a frequency comb.
From the measured Wigner functions, we compute the density matrix of the storage mode ρ(t) for various times t for the three cases and compare the evolution of the normalized elements ρ 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 singe drive is turned on. We conclude that a multiplexing measurement extracts more information than a single measurement.
The effect on ρ 02 when probing with a resonant drive for n = 1, is consistent with the significant measurementinduced detuning that can be read off the top right plot of Fig. 20 (blue, value 1 on the horizontal axis). Apparently, when driving with a comb, such an effect combines with the ones on n = 0 and n = 2 resonances, and other components, to induce a stronger overall measurement rate. We will investigate this comb effect more precisely in section H 6.

Quantum Non Demolition nature of the multiplexed measurement
The goal of this subsection is to quantify the Quantum Non Demolition (QND) nature of our multiplexed measurement. A measurement is said to be QND if • the measurement time is very short compared to the timescale of evolution of the system under study, • the interaction with the probe does not disturb the quantum state of the system if it belongs to the measurement basis.
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 fmp − ∆mp with an amplitude Ω = χs,mp/2. Lines represent the highest eigenvalue of (H8) without any free parameters. An offset equals to Γ φ,s (n−m) 2 , which is the intrinsic dephasing of the storage mode, is added to obtain the total dephasing rate. 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. 22a), the probability of finding the storage mode with 0 photon is slightly lower. This dependence on the amplitude of the drive Ω is best characterized by extracting the populations (Fig. 22b) and photon number (Fig. 22d) 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 photon 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 photon at the end of the measurement by approximately 10 %. We find that the mean photon number is decreased by the same percentage.

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 elements of the storage density matrix to gain insight into the physics of the dephasing process. Fig. 23a and b show the decay of off-diagonal elements ρ 12 and ρ 13 as a function of time. For small drive amplitudes Ω < 0.5χ s,mp , off-diagonal elements decay faster when Ω 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 become more pronounced as the drive amplitude becomes larger and they exhibit a quasi periodicity. 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. 23c).

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 ρ nm . In this section, we introduce two 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. 3B, is the quadrature of the storage mode in the frame rotating at the frequency of this mode when the qubit is probed by a comb. It can be expressed as Re[( X + i P )e −2iπδfst ] with X and 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 FIG. 24. Normalized off-diagonal elements of the storage matrix density for the largest measurement amplitudes. Measured normalized coherence |ρ12|/ √ ρ11ρ22 between Fock states |1 and |2 of the storage mode as a function of time and for various amplitudes Ω > 0.9χs,mp of the driving frequency comb.
C ρ (t) contains the information about the dephasing between every different Fock states. 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 â s +â † s = Re[( X +i P )e −2iπδfst ] 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.