# Multiplexed photon number measurement

Antoine Essig∗,1, Quentin Ficheux∗,1, Théau Peronnin1, Nathanaël Cottet1,
Raphaël Lescanne2,3, Alain Sarlette3, Pierre Rouchon4,3, Zaki Leghtas4,2,3,
Benjamin Huard1

$*$ The two authors contributed equally.

1Univ Lyon, ENS de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique,
F-69342 Lyon, France
2 Laboratoire de Physique de l’Ecole Normale Supérieure, ENS, Université PSL,
CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité,
Paris, France
3QUANTIC team, INRIA de Paris, 2 Rue Simone Iff, 75012 Paris, France
4Centre Automatique et Systèmes, Mines-ParisTech, PSL Research University,
60, bd Saint-Michel, 75006 Paris, France

> The evolution of quantum systems under measurement is a central aspect of
> quantum mechanics. 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 followed by its corresponding quantum collapse. Here,
> we report an experiment where a single superconducting qubit is counter-
> intuitively able to answer not a single but nine yes-no questions about the
> number of photons in a microwave resonator at the same time. The key
> ingredients are twofold. First, we exploit the fact that observing the color
> of a qubit carries additional information to the conventional readout of its
> state. The qubit-system interaction is hence designed so that the qubit
> color encodes the number of photons in the resonator. Secondly, we multiplex
> the qubit color observation by recording how the qubit reflects a frequency
> comb. Interestingly the amount of extracted information reaches a maximum at
> a finite drive amplitude of the comb. We evidence it by direct Wigner
> tomography of the quantum state of the resonator. Our experiment unleashes
> the full potential of quantum meters by bringing the measurement process in
> the frequency domain.

##

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 a series 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 (?, ?) with the possible refinement of choosing what binary question
should be optimally asked conditioned on the previous answers (?). At most,
each answer provides one bit of information about the system state (?) so that
determining an arbitrary number of photons in the cavity between $0$ and
$2^{m}-1$ takes at least $m$ consecutive probe measurements. In this work, we
experimentally show that one can ask all questions at once using a single
qubit as a probe. Daring an analogy with communication protocols (?), previous
measurement

Figure 1: Multiplexing photon counting setup and protocol. (A) Time domain
division multiplexing proceeds one binary question at a time. Frequency
division multiplexing simultaneously retrieves multiple binary answers. (B)
Scheme of the device in coplanar waveguide architecture. The storage mode
(green) is coupled to a transmon multiplexing qubit (orange), which is
directly coupled to a transmission line. 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 for standard photon counting and Wigner tomography
(?). (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 illuminated by a frequency comb of amplitude
$\Omega$. The reflected pulse is amplified and digitized before numerical
demodulation at every frequency $f_{\mathrm{mp}}-k\chi_{\mathrm{s,mp}}$. This
multiplexing-demultiplexing process leads to reflection coefficients $r_{k}$
that each contain the probability that $k$ photons are stored.

schemes with time series of binary questions used time division multiplexing
while our experiment demonstrates the analogous of frequency division
multiplexing (Fig. 1A).

Both approaches can be implemented to count the photon number of a
superconducting microwave resonator dubbed the _storage mode_ , which
resonates at $f_{\mathrm{s}}=4.558~{}\mathrm{GHz}$. We couple the latter to
two off-resonant transmon qubits (see Fig. 1B). The _yes-no qubit_ with a
frequency $f_{\mathrm{yn}}=3.848~{}\mathrm{GHz}$ can implement binary
questions about the photon number while the _multiplexing qubit_ with a
frequency $f_{\mathrm{mp}}=4.238~{}\mathrm{GHz}$ can implement frequency
multiplexed photon counting. Both qubits are dispersively coupled to the
resonator so that their frequency respectively redshifts by
$\chi_{\mathrm{s,yn}}=1.4~{}\mathrm{MHz}$ and
$\chi_{\mathrm{s,mp}}=4.9~{}\mathrm{MHz}$ per additional photon in the storage
mode.

Figure 2: Photocounting a mode in a coherent state. The storage mode is
prepared in a coherent state with $\bar{n}$ photons on average using a
microwave pulse at storage frequency and amplitude $V_{s}$. (A,C) Standard
approach to photon counting. A $\pi$-pulse at a tunable frequency
$f_{\mathrm{drive}}$ on the yes-no qubit excites the qubit with a measured
probability $\mathbb{P}_{\mathrm{e}}$ as a function of $V_{s}$. Peaks appear
at $f_{\mathrm{yn}}-k\chi_{\mathrm{s,yn}}$ and indicate the probability to
store $k$ photons. The dots in (C) are cuts along the dashed lines in (A) and
match the master equation model (solid lines) using it as a calibration of the
storage mode occupation. (B,D) Photon counting using the reflection
coefficient $r$ of the multiplexing qubit at a single frequency
$f_{\mathrm{probe}}$. The measured emission coefficient $1-\mathcal{R}e(r)$ is
shown as a function of $f_{\mathrm{probe}}$ and $V_{\mathrm{s}}$ and exhibits
resolved peaks for every photon number. The dots in (D) are cuts along the
dashed lines in (B) and are captured by a master equation model (solid lines).
(E) Multiplexed photon counting. Dots: simultaneously measured average
emission coefficients for each photon number $k$ from 0 to 8 and for various
initial state of the storage mode. $r_{k}$ is here the reflection coefficient
at $f_{\mathrm{mp}}-k\chi_{\mathrm{s,mp}}$. Solid lines: prediction based on a
master equation without free parameter (?).

Starting with the standard approach (?, ?), the probability to have $k$
photons is encoded as the probability $\mathbb{P}_{\mathrm{e}}$ to excite the
yes-no qubit by driving it with a $\pi$-pulse at
$f_{\mathrm{drive}}=f_{\mathrm{yn}}-k\chi_{\mathrm{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 initialize the storage mode in a coherent
state
$\left|\beta\right>=e^{-|\beta|^{2}/2}\sum_{n=0}^{+\infty}\frac{\beta^{n}}{\sqrt{n!}}\left|n\right>$,
which is a superposition of all Fock states with mean photon number
$\bar{n}=|\beta|^{2}$ using a microwave tone at $f_{\mathrm{s}}$. The
probability $\mathbb{P}_{\mathrm{e}}$ is then measured and shows resolved
peaks as a function of $f_{\mathrm{drive}}$ for every photon number up to
about 7 (Fig. 2A,C) and various average photon numbers $\bar{n}$. The linear
relation between $\beta$ and the amplitude $V_{\mathrm{s}}$ of the tone at
$f_{\mathrm{s}}$ was calibrated using a master-equation based model (?)
reproducing the measured $\mathbb{P}_{\mathrm{e}}$ (solid lines in Fig. 2C).

The intrinsic limitation of the standard approach is that measuring a qubit
state can at most reveal one bit of information at a time. We propose to avoid
this constraint by observing the qubit frequency (or color) directly without
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
into the mode of the transmission line that is centered around the qubit
frequency $f_{\mathrm{mp}}-k\chi_{\mathrm{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_{\mathrm{mp}}-k\chi_{\mathrm{s,mp}}\\}_{k}$, which can host much more
than one bit of information. This encoding ability can be observed by driving
the multiplexing qubit in reflection through the transmission line (Fig. 1B)
(?, ?, ?, ?). The measured real part $\mathcal{R}e(r)$ of the reflection
coefficient of a microwave pulse at frequency $f_{\mathrm{probe}}$ is reduced
when the probe resonates with the qubit, hence revealing the photon number $k$
(Fig. 1C) (?). This reduction arises from the coherent emission by the qubit
in phase opposition with the reflected drive (?). Interestingly, the roles of
qubit and resonator are here swapped compared to usual qubit readout in
circuit-QED (?).

In Fig. 2B,D we show the measured qubit _emission coefficient_
$1-\mathcal{R}e(r)$ as a function of a single probe frequency
$f_{\mathrm{probe}}$ and of the initial amplitude of the storage mode coherent
state $\sqrt{\bar{n}}$. The measurement is performed using a drive strength
$\Omega=\chi_{\mathrm{s,mp}}/4$ (expressed as the corresponding Rabi
frequency) and pulse duration of $2~{}\mu\mathrm{s}$, which is smaller than
the storage lifetime of $3.8~{}\mu\mathrm{s}$. Resolved peaks develop for
every photon number up to at least 9 which is already a clearer indicator of
the photon number than the yes-no probe. Using the former calibration of
$\bar{n}$, a master-equation based model enables us to reproduce the
measurement results (?). In order to get resolved peaks we designed the
relaxation rate of the multiplexing qubit
$\Gamma_{\mathrm{1,mp}}=(44\mathrm{~{}ns})^{-1}$ so that the decoherence rate
$\Gamma_{\mathrm{2,mp}}=\Gamma_{\mathrm{1,mp}}/2$ is smaller than the
dispersive shift $2\pi\chi_{\mathrm{s,mp}}$. Probing the qubit at one of its
resonance frequencies opens a communication channel with a maximal bandwidth
$\Gamma_{\mathrm{2,mp}}$. Hence we maximized this bandwidth by designing
$\Gamma_{\mathrm{1,mp}}$ as large as possible using the direct coupling to the
transmission line.

Therefore we have shown that both measurements in Figs. 2A,B allow us to ask
questions of the kind “are there $k$ photons?”. However the yes-no qubit is
read out using a pulse that is independent on the storage state. In contrast,
each continuous measurement of the frequency mode
$f_{\mathrm{mp}}-k\chi_{\mathrm{s,mp}}$ of the transmission line in which the
multiplexing qubit emits reveals information about a different photon number
$k$. It thus enables the key ingredient of our approach: the multiplexing
measurement of reflection at every frequency
$f_{\mathrm{mp}}-k\chi_{\mathrm{s,mp}}$.

To do so, we illuminate the multiplexing qubit with a pulse containing a comb
with $9$ frequencies corresponding to photon numbers from 0 to 8 (see Fig.
1D). The efficient measurement of the reflected pulse requires the use of a
near-quantum limited amplifier with a dynamical bandwidth of at least a dozen
of $\chi_{\mathrm{s,mp}}$ such as a TWPA (?). This technical challenge is
similar to the recently demonstrated high-efficiency multiplexed readout of as
many as 6 qubits coupled to a single feed line (?, ?, ?, ?). We then
demultiplex the reflected pulse at 9 frequencies
$\\{f_{\mathrm{mp}}-k\chi_{\mathrm{s,mp}}\\}_{0\leq k\leq 8}$ and extract a
reflection coefficient $r_{k}$ for each of them (Fig. 1D). The measurement
consists in simultaneously measuring the emission coefficients
$1-\mathcal{R}e(r_{k})$ for each peak in Fig. 2B,D, which is of course much
faster than measuring them one at a time. They are shown as a function of the
average initial photon number in Fig. 2E for a drive strength
$\Omega=\chi_{\mathrm{s,mp}}/2$ and a measurement duration of
$2~{}\mu\mathrm{s}$. For a given $\bar{n}$, every measurement channel $k$
gives an average signal that is proportional to the probability of having $k$
photons in the storage mode. As $\bar{n}$ is varied, the shape of the average
signal of channel $k$ reproduces a Poisson distribution distorted by
relaxation processes and channel cross-talk due to large driving strength (?).
This multiplexed photon counting signal can be reproduced using a master
equation approach (solid lines in Fig. 2E) without any free parameter. This
result thus demonstrates the possibility to multiplex the photon number
measurement using a single qubit.

How is the measurement affected by the strength $\Omega$ of the driving
frequency comb? 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 (?, ?). Thus,
characterizing the measurement rate of our multiplexed photon counting can
elegantly be done by observing how the storage mode dephases for a given
driving strength $\Omega$. 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 a qubit,
one expects a different dependence of the measurement rate on $\Omega$ than
for standard qubit readout using a probe cavity.

Figure 3: Storage mode dephasing due to multiplexed photon counting. (A)
Direct Wigner tomography of the storage mode at time $t=0.5~{}\mu\mathrm{s}$
after initialization in a coherent state with amplitude $\beta=-1.55$. Left:
natural evolution without driving. Middle: evolution with a single tone at
$f_{\mathrm{mp}}-\chi_{\mathrm{s,mp}}$ that probes whether there is 1 photon
in the storage mode. Right: evolution in presence of the frequency comb
$\\{f_{\mathrm{mp}}-k\chi_{\mathrm{s,mp}}\\}_{0\leq k\leq 8}$. Each bottom
panel represents the density matrix $\\{\rho_{nm}\\}_{n,m}$ that is calculated
from the Wigner function of its top panel. The cylinder heights represent
$|\rho_{nm}|$ and the color encodes $\mathrm{arg}(\rho_{nm})$. (B) Measured
exponential decay of the average storage mode quadrature $\langle
a_{\mathrm{s}}+a_{\mathrm{s}}^{\dagger}\rangle$ in the case of a driving comb
of strength $\Omega=0.23\chi_{\mathrm{s,mp}}$, and in a frame rotating at the
storage mode resonant frequency. Insets show Wigner tomography of the storage
mode for two values of $t$. (C) Measured dephasing rate induced on the storage
mode as a function of the drive strength $\Omega$ of the frequency comb.
Insets show Wigner tomographies of the storage mode at $t=1~{}\mu\mathrm{s}$
for four values of $\Omega$.

In order to measure this dephasing rate, we perform a direct Wigner tomography
of the storage mode at various times $t$ (?, ?, ?). It provides a
representation of the state $\rho$ in the phase space of the mode and can be
expressed as $W(\alpha)=\frac{2}{\pi}\mathrm{Tr}(D^{\dagger}(\alpha)\rho
D(\alpha)\mathcal{P})$. Here $D(\alpha)=e^{\alpha
a_{s}^{\dagger}-\alpha^{*}a_{s}}$ is the storage displacement operator,
$\mathcal{P}=e^{i\pi a_{s}^{\dagger}a_{s}}$ is the photon number parity
operator, and $a_{s}$ is the canonical annihilation operator of the storage
mode. Preparing the storage mode with a coherent state $|\beta=-1.55\rangle$,
the Wigner function starts as a Gaussian distribution centered at
$\alpha=\beta$. On the left of Fig. 3A, one can see how the bare dephasing
rate and the self-Kerr effect of the storage mode ($0.02~{}\mathrm{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 $\Omega=\chi_{\mathrm{s,mp}}/2$ to measure whether
there is 1 photon, the phase diffuses faster and the Wigner function exhibits
negativities in middle of Fig. 3A. As seen in the corresponding density
matrix, a tone at $f_{\mathrm{probe}}=f_{\mathrm{mp}}-\chi_{\mathrm{s,mp}}$
notably induces dephasing between states $|1\rangle$ and all other states
$|m\neq 1\rangle$ (see density matrix as a function of drive frequency in
(?)). As expected the phase diffusion is more intense when all the tones of
the multiplexed readout are turned on than for a single tone with the same
drive strength $\Omega$ (right of Fig. 3A). Likewise all off-diagonal elements
of the density matrix are then reduced.

To be more quantitative, the dephasing rate $\Gamma_{\mathrm{d,s}}$ of the
cavity is accessed through the decay of the mean quadrature of the storage
mode $\langle a_{s}+a_{s}^{\dagger}\rangle=\int
2xW(x+iy)\mathrm{d}x\mathrm{d}y$ (?). In Fig. 3B, we show $\langle
a_{s}+a_{s}^{\dagger}\rangle$ as a function of time under the multiplexed
drive with a strength $\Omega=0.23\chi_{\mathrm{s,mp}}$. Repeating this
experiment for various values of the multiplexed driving strength $\Omega$
allows us to determine how the latter affects the dephasing rate
$\Gamma_{\mathrm{d,s}}$, and thus the measurement rate. The dephasing rate is
non monotonic in the drive strength (Fig. 3C). Noticeably, it reaches a
maximum when $\Omega=\chi_{\mathrm{s,mp}}/2$ for which information is
extracted at a rate approximately 5 times larger than the natural dephasing
rate. This maximum is well captured by a measurement model based on a master
equation without any free parameter (line in Fig. 3C) (?).

In conclusion, we have demonstrated that a single qubit can be used to
continuously probe a multidimensional system by encoding the information in
the frequency domain. The measurement rate depends non-monotonically on the
drive strength. Improving further the detection efficiency, the storage
lifetime, and lowering the residual qubit thermal excitation should enable
single shot photon counting by multiplexing. Our new measurement scheme can
readily be applied to stabilize quantum states by feedback control (?), probe
quantum trajectories of microwave modes (?), Quantum Zeno dynamics (?), or
vary the driving comb tones in real time to engineer desired decoherence
channels. Furthermore, exciting the multiplexing qubit and detecting the
emitted photon by a series of frequency resolved photodetectors (?, ?, ?), one
could realize a discrete and projective photocounter. Moving further, one
could extend this frequency domain measurement to more complex probes than a
single qubit and many possible physical systems beyond superconducting
circuits.

## Acknowledgments

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

## Supplementary materials

###  A Device and Fabrication

####  A.1 Design

The circuit is composed of 4 electromagnetic modes whose parameters can be
found in Tab. 1 and 2. A high-Q harmonic oscillator, called storage mode, is
composed of a $\mathrm{\lambda/2}$ coplanar waveguide (CPW) resonator (green
in Fig. 4). The storage resonator is capacitively coupled to two transmon
qubits. The multiplexing qubit (orange) has a high spontaneous photon emission
rate $\Gamma_{\mathrm{1,mp}}=(44\mathrm{~{}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_{\mathrm{2,yn}}=27~{}\mu\mathrm{s}$. As required by Wigner tomography, the
yes-no qubit coherence time and the lifetime of storage mode are larger than
the time needed for storage photon number parity measurement
$1/2\chi_{\mathrm{s,yn}}\ll T\mathrm{{}_{1,s}},T_{\mathrm{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 (?) 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
$\mathrm{\Gamma_{2,mp}<2\pi\chi_{s,mp}}$. This resolution constraint is not
critical, as in fact a finite amount of photon number information can be
extracted as soon as $\chi_{\mathrm{s,mp}}$ is nonzero, but the decoding is
much simpler if we can reason in terms of well-separated resonance peaks.

circuit parameters | Symbol | Hamiltonian term | parameter value
---|---|---|---
readout resonator frequency | $f_{\mathrm{ro}}$ | $hf_{\mathrm{ro}}\hat{n}_{\mathrm{ro}}$ | 7.138 GHz
storage mode frequency | $f_{\mathrm{s}}$ | $hf_{\mathrm{s}}\hat{n}_{\mathrm{s}}$ | 4.558 GHz
yes-no qubit frequency | $f_{\mathrm{yn}}$ | $hf_{\mathrm{yn}}\hat{n}_{\mathrm{yn}}$ | 3.848 GHz
multiplexing qubit frequency | $f_{\mathrm{mp}}$ | $hf_{\mathrm{mp}}\hat{n}_{\mathrm{mp}}$ | 4.238 GHz
readout/yes-no qubit cross-Kerr rate | $\chi_{\mathrm{ro,yn}}$ | $-h\chi_{\mathrm{ro,yn}}\hat{n}_{\mathrm{ro}}\hat{n}_{\mathrm{yn}}$ | 0.4 MHz
storage/yes-no qubit cross-Kerr rate | $\chi_{\mathrm{s,yn}}$ | $-h\chi_{\mathrm{s,yn}}\hat{n}_{\mathrm{s}}\hat{n}_{\mathrm{yn}}$ | 1.4 MHz
storage/multiplexing qubit cross-Kerr rate | $\chi_{\mathrm{s,mp}}$ | $-h\chi_{\mathrm{s,mp}}\hat{n}_{\mathrm{s}}\hat{n}_{\mathrm{mp}}$ | 4.9 MHz
yes-no qubit anharmonicity | $\chi_{\mathrm{yn,yn}}$ | $-h\chi_{\mathrm{yn,yn}}\hat{n}_{\mathrm{yn}}(\hat{n}_{\mathrm{yn}}-1)$ | 160 MHz
multiplexing qubit anharmonicity | $\chi_{\mathrm{mp,mp}}$ | $-h\chi_{\mathrm{mp,mp}}\hat{n}_{\mathrm{mp}}(\hat{n}_{\mathrm{mp}}-1)$ | 116 MHz
Table 1: Table of circuit parameters. Cross-Kerr rates with the storage mode
are measured using qubit spectroscopy for various coherent states. The cross-
Kerr rate between yes-no qubit and readout resonator is measured by
calibrating the measurement induced dephasing by the readout resonator on the
qubit as a function of the drive amplitude and detuning of the readout
resonator. Anharmonicities are measured using spectroscopy of the qubit
excited state. Figure 4: Optical image of the circuit. The readout resonator
is colored in purple, storage mode in green, yes-no qubit in blue and
multiplexing qubit in orange. All dark grey areas are silicon, grey areas are
niobium on silicon and Josephson junctions are made of Al/AlOx/Al.

####  A.2 Measurement setup

The readout resonator, the yes-no qubit and the multiplexing qubit are driven
by pulses that are generated using a Tektronix® Arbitrary Waveform Generator
(AWG) AWG5014C with a sample rate of 1 GS/s. Storage mode pulses are generated
using a Zurich Instruments® UHFLI with a sample rate of 1.8 GS/s. The UHFLI
allows us to change the pulse amplitude and phase without recompiling the
sequence. This feature decreases the time needed for Wigner tomography
compared to a standard AWG and makes the pulse sequence simple with the
drawback of having to synchronize the two AWGs. AWG pulses are modulated at a
frequency 25 MHz for readout, 100 MHz for yes-no qubit, and 75 MHz for storage
and multiplexing qubit. They are up-converted using single sideband mixers for
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, WindFreak® SynthHD, and
AnaPico® APSIN20G sources that are set at the frequencies $f_{\mathrm{ro}}+25$
MHz, $f_{\mathrm{mp}}+75$ MHz, $f_{\mathrm{s}}+75$ MHz and
$f_{\mathrm{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_{\mathrm{TWPA}}$ = 5.998 GHz) and power in order to reach a gain of 20 dB
at 7.138 GHz and 4.238 GHz. 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. The full setup is shown in Fig. 5. The Tektronix®
AWG is used as the master that triggers the UHFLI and the Alazar® board.

Figure 5: Schematic of the setup. Each electromagnetic mode of the experiment
is driven by a RF generator detuned by the modulation frequency and whose
color matches the color of the corresponding mode in Fig. 4. Room temperature
isolators are not represented for the sake of clarity.

####  A.3 Master equation of the full system

We give the full Hamiltonian of the circuit up to the second order in photon
number, a more complex Hamiltonian is given in section C.1. We denote the
photon number operator respectively for the readout, storage, yes-no qubit and
multiplexing qubit with $\hat{n}_{\mathrm{{ro}}}$, $\hat{n}_{\mathrm{s}}$,
$\hat{n}_{\mathrm{{yn}}}$, and $\hat{n}_{\mathrm{{mp}}}$. Note that the
annihilation operator $\hat{a}_{\mathrm{s}}$ in the main text is related to
these operators by
$\hat{n}_{\mathrm{s}}=\hat{a}_{\mathrm{s}}^{\dagger}\hat{a}_{\mathrm{s}}$.
With those operators, the undriven Hamiltonian can be written as

$\begin{array}[]{rcl}\hat{H}&=&hf_{\mathrm{ro}}\hat{n}_{\mathrm{ro}}+hf_{\mathrm{s}}\hat{n}_{\mathrm{s}}+hf_{\mathrm{yn}}\hat{n}_{\mathrm{yn}}+hf_{\mathrm{mp}}\hat{n}_{\mathrm{{mp}}}-h\chi_{\mathrm{ro,yn}}\hat{n}_{\mathrm{ro}}\hat{n}_{\mathrm{yn}}-h\chi_{\mathrm{s,yn}}\hat{n}_{\mathrm{s}}\hat{n}_{\mathrm{yn}}\\\
&&-h\chi_{\mathrm{s,mp}}\hat{n}_{\mathrm{s}}\hat{n}_{\mathrm{mp}}-h\chi_{\mathrm{yn,yn}}\hat{n}_{\mathrm{yn}}(\hat{n}_{\mathrm{yn}}-1)-h\chi_{\mathrm{mp,mp}}\hat{n}_{\mathrm{mp}}(\hat{n}_{\mathrm{mp}}-1)\end{array},$
(1)

with $\chi_{a,b}$ the cross-Kerr rate between mode $a$ and $b$ and
$\chi_{a,a}$ the anharmonicity of the mode $a$. The values of all parameters
are given in Table 1.

To this Hamiltonian, we need to add seven dephasing and relaxation channels in
order to get the Lindblad equation describing the circuit dynamics.

$\begin{array}[]{rcl}\dot{\rho}&=&-\dfrac{i}{\hbar}[\hat{H},\rho]+\Gamma_{\mathrm{ro}}\mathcal{L}(\hat{a}_{\mathrm{ro}})\rho+2\Gamma_{\mathrm{\phi,s}}\mathcal{L}(\hat{n}_{\mathrm{s}})\rho+\Gamma_{\mathrm{1,s}}\mathcal{L}(\hat{a}_{\mathrm{s}})\rho\\\
&&+2\Gamma_{\mathrm{\phi,yn}}\mathcal{L}(\hat{n}_{\mathrm{yn}})\rho+\Gamma_{\mathrm{1,yn}}\mathcal{L}(\hat{a}_{\mathrm{yn}})\rho+2\Gamma_{\mathrm{\phi,mp}}\mathcal{L}(\hat{n}_{\mathrm{mp}})\rho+\Gamma_{\mathrm{1,mp}}\mathcal{L}(\hat{a}_{\mathrm{mp}})\rho,\end{array}$
(2)

where $\mathcal{L}$ is the Lindblad superoperator defined as
$\mathcal{L}(\hat{L})\rho=\hat{L}\rho\hat{L}^{\dagger}-\left\\{\hat{L}^{\dagger}\hat{L},\rho\right\\}/2$
, $\hat{a}_{b}$ is the annihilation operator of mode $b$. For a qubit mode
$b$, the dephasing rate $\Gamma_{\mathrm{\phi,b}}$ is linked to the
decoherence rate by
$\Gamma_{\mathrm{2,b}}=\Gamma_{\mathrm{1,b}}/2+\Gamma_{\mathrm{\phi,b}}$.
Table 2 gives the values of all relaxation and decoherence rates.

parameters | symbol | value
---|---|---
readout decay rate | $\Gamma_{\mathrm{ro}}$ | $(40~{}\mathrm{ns})^{-1}$
storage decay rate | $\Gamma_{\mathrm{1,s}}$ | $\mathrm{(}3.8~{}\mu s)^{-1}$
storage decoherence rate | $\Gamma_{\mathrm{2,s}}$ | $\mathrm{(}2~{}\mu s)^{-1}$
yes-no decay rate | $\Gamma_{\mathrm{1,yn}}$ | $\mathrm{(}20~{}\mu s)^{-1}$
yes-no decoherence rate | $\Gamma_{\mathrm{2,yn}}$ | $\mathrm{(}27~{}\mu s)^{-1}$
multiplexing decay rate | $\Gamma_{\mathrm{1,mp}}$ | $\mathrm{(}44~{}ns)^{-1}$
multiplexing decoherence rate | $\Gamma_{\mathrm{2,mp}}$ | $\mathrm{(}88~{}ns)^{-1}$
Table 2: Table of all relaxation and decoherence rates. The readout decay rate
is obtained from the measurement induced dephasing rate of the yes-no qubit by
readout photons. The yes-no qubit decay rate is measured using the time
evolution of probability to find the qubit excited after a $\pi$ pulse. The
yes-no qubit decoherence rate is measured using Ramsey interferometry. The
storage decay rate is measured using the time evolution of the probability to
have zero photon after populating the storage mode. Storage decoherence rate
is measured using what we define as a Ramsey interferometry experiment for an
harmonic oscillator (see section B.3). The multiplexing decay and decoherence
rates are measured by fitting the qubit resonance for various probe
amplitudes.

####  A.4 Electromagnetic simulation

The circuit was simulated using Sonnet® software. The length of the readout
resonator and storage mode was designed to obtain the resonant frequencies
$f_{\mathrm{ro}}\sim 7~{}\mathrm{GHz}$ and $f_{\mathrm{s}}\sim
4.5~{}\mathrm{GHz}$. In the simulation, the qubit Josephson junctions were
replaced by linear inductances. The value of the inductances and of the
capacitive couplings between modes were adjusted to match the resolved number
resolved condition (see section A.1) between the qubits and the storage mode.
Couplings were extracted from the electromagnetic simulation by varying the
Josephson junction inductances and by predicting the avoided crossing between
the qubits and the storage and readout resonators.

####  A.5 Device fabrication

The circuit consists of a sputtered 120 nm-thick Niobium film deposited on a
280 $\mu$m-thick undoped silicon wafer. The resonators and feed lines are dry
etched after optical lithography. 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.

The first cooldown of the device revealed that the multiplexing qubit
frequency and storage frequency were too close to be modeled by a dispersive
Hamiltonian. They were detuned by about 100 MHz and the storage mode was
significantly hybridized with the multiplexing qubit, entailing a large photon
loss rate. We estimated that a 10$\%$ increase of the junction resistance
should decrease the frequency of the multiplexing qubit by 200 MHz. We warmed
up the circuit at room temperature and placed it on a hot plate during 60 min
at 100∘C. We monitored the resistance of 9 test junctions located on the chip
every 10 minutes until we reached a 10% increase. During the second cool down
we measured a detuning between the multiplexing qubit and the storage mode of
about 320 MHz, as targeted.

###  B Calibration

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

The storage mode can be displaced by driving it on resonance with a voltage
$V_{\mathrm{s}}(t)\cos(2\pi f_{\mathrm{s}}t+\phi_{\mathrm{s}})$, where
$V_{\mathrm{s}}(t)$ is the pulse envelope. The driving Hamiltonian of the
storage mode reads
$\hbar(\epsilon_{\mathrm{s}}(t)\hat{a}_{\mathrm{s}}^{\dagger}+\epsilon_{\mathrm{s}}^{\ast}(t)\hat{a}_{\mathrm{s}})$
where $\epsilon_{\mathrm{s}}(t)=\mu V_{\mathrm{s}}(t)e^{i\phi_{\mathrm{s}}}$.
The scaling factor $\mu=1.45~{}\mathrm{(mV.\mu s)^{-1}}$ is calibrated by
fitting the photocounting measurement results obtained using the yes-no qubit
with the master equation simulation (see section C.1). Fig. 6a shows the
evolution of $\epsilon_{s}$ with $V_{\mathrm{s}}$. For every experiment, the
storage mode displacements are realized using a Gaussian pulse shape
$\epsilon_{\mathrm{s}}(t)=\lambda(t)\epsilon_{\mathrm{max}}$ with a maximum
amplitude $\epsilon_{\mathrm{max}}$, a width $25~{}\mathrm{~{}ns}$ and a
duration $100~{}\mathrm{~{}ns}$. We simulated the dynamics of the storage mode
under this Gaussian displacement taking into account the couplings, relaxation
and decoherence rates (see section C.2) for various amplitudes
$\epsilon_{\mathrm{max}}$. We then computed the expectation value of the
number of photon operator $\langle\hat{n}_{\mathrm{s}}\rangle$ at the end of
the pulse. Fig. 6b shows the square root of
$\langle\hat{n}_{\mathrm{s}}\rangle$ as a function of
$\epsilon_{\mathrm{max}}$. Fitting with a linear function, we find that
$\sqrt{\langle\hat{n}_{\mathrm{s}}\rangle}=59.1\epsilon_{\mathrm{max}}$. As
$\epsilon_{s}$ increases linearly with $V_{\mathrm{s}}$,
$\epsilon_{\mathrm{max}}$ increases linearly with the maximum voltage
amplitude $V_{\mathrm{max,s}}$ of the Gaussian pulse
$V_{\mathrm{s}}(t)=\lambda(t)V_{\mathrm{max,s}}$. Using the two linear
regressions, we can express the photon number of the storage mode as
$\sqrt{\langle
n_{\mathrm{s}}\rangle}=(85.9~{}\mathrm{V}^{-1})\,V_{\mathrm{max,s}}$.

Figure 6: Calibration of the average number of photons $\langle
n_{\mathrm{s}}\rangle$ in the storage mode as a function of the displacement
amplitude. a. Evolution of the displacement amplitude $\epsilon_{s}$ with the
pulse envelope $V_{s}$. The calibration is obtained by comparing the results
of a photon counting experiment using the yes-no qubit with a master equation
simulation (see section C.1). b. Square root of the average photon number
$\langle n_{\mathrm{s}}\rangle$ in the storage mode as a function of the drive
amplitude. The storage is displaced by 100 ns long Gaussian pulse with a width
of 25 ns. The same pulse shape is used in the simulation. From the two linear
fits we extract the evolution of the mean number of photons with the amplitude
of the pulse $\sqrt{\langle
n_{\mathrm{s}}\rangle}=\left(85.9~{}\mathrm{V}^{-1}\right)\,V_{\mathrm{max,s}}$.

####  B.2 Wigner tomography calibration

The Wigner function of a harmonic oscillator with density matrix $\rho$ is
defined as $W(\alpha)=$ $2\,\mathrm{Tr}(D^{\dagger}(\alpha)\rho
D(\alpha)\mathcal{P})/\pi$ where
$\mathrm{D}(\alpha)=e^{\alpha\hat{a}_{\mathrm{s}}^{\dagger}-\alpha^{*}\hat{a}_{\mathrm{s}}}$
is the displacement operator of the storage mode by a coherent field $\alpha$
and $\mathcal{P}=e^{i\pi\hat{a}_{\mathrm{s}}^{\dagger}\hat{a}_{\mathrm{s}}}$
is the photon number parity operator. A Wigner function is the expectation
value of the parity after a displacement by an amplitude $\alpha$. The Wigner
tomography sequence is represented on Fig. 7a . It starts by realizing a
displacement on the storage mode with a 100 ns long Gaussian pulse at
frequency $f_{\mathrm{s}}$ (or detuned for Ramsey interferometry of the
storage mode, see section B.3) with a width of 25 ns. Then two successive
$\pi/2$ Gaussian pulses of 18 ns with a width of 4.5 ns are sent to the yes-no
qubit at $f_{\mathrm{yn}}$ and are separated by a waiting time
$\Delta\tau=337~{}\mathrm{ns}\approx 1/2\chi_{\mathrm{s,yn}}$. It implements a
parity measurement and maps the parity of the storage mode onto the $z$-axis
on the yes-no qubit (?, ?, ?). The sequence terminates by a 2 $\mu$s long
square pulse on the readout resonator to read out the state of the yes-no
qubit. For high amplitude $\alpha$, higher order Kerr terms distort the Wigner
function. To mitigate this effect, we interleaves two sequences with a final
pulse of phase either $+\pi/2$ or $-\pi/2$.(see Fig. 7a). 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. 7b is calibrated using the yes-no qubit Rabi
oscillation amplitude to express the signal using Pauli operators. Multiplying
the result by $2/\pi$ yields the Wigner function W($\alpha$) in Fig. 7c.

Figure 7: Direct Wigner tomography of the storage mode a. Circuit diagram of a
Wigner tomography using a parity measurement based on dispersive interaction.
After a 100 ns displacement pulse on the storage mode, an unconditional
$\pi/2$ pulse is applied to the yes-no qubit. The qubit evolves freely during
a time $\Delta\tau=337~{}\mathrm{ns}\approx 1/2\chi_{\mathrm{s,yn}}$ before a
new $\pm\pi/2$ pulse is sent and the state of the yes-no qubit is measured
using the readout resonator. b. Calibration of the quadrature axes for Wigner
tomography. Blue dots represent the standard deviation of the quadratures of
the displaced thermal equilibrium state of the storage mode as a function of
drive amplitude for various detuning using the photon number calibration alone
(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 $f_{\mathrm{mp}}-1.4\chi_{\mathrm{s,mp}}$
during 750 ns with an amplitude $\Omega=\chi_{\mathrm{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. Figure 8: Revival of the Ramsey interferometry on the yes-no
qubit a. Circuit diagram for Ramsey interferometry in the presence of storage
photons. After a 100 ns displacement pulse at the storage frequency, an
unconditional $\pi/2$ pulse is applied to the yes-no qubit. We then let the
qubit evolve freely during a time $t$ before doing a new $\pm\pi/2$ pulse and
measure the state of the yes-no qubit. The signal S(t) is half the difference
between the average outcomes of the two sequences. b. Measured (dots) and
predicted (lines) signal S as a function of waiting time $t$. Predicted signal
is computed from Eq. (4). yes-no qubit revivals occur every
$1/\chi_{\mathrm{s,yn}}\approx 0.7\mu\mathrm{s}$.

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_{\mathrm{th}}$ the Wigner function is a
2D Gaussian function with a width $\sqrt{n_{\mathrm{th}}+1/2}$ (?)

$W_{\rho(n_{\mathrm{th}})}(\alpha=x+ip)=\dfrac{2}{\pi}\dfrac{1}{2n_{\mathrm{th}}+1}e^{-2|\alpha|^{2}/(2n_{\mathrm{th}}+1)}.$
(3)

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

The duration $\Delta\tau$ is calibrated using qubit state revival during
Ramsey interferometry (see supplementary material of Ref. (?)). We used a
Ramsey interferometry sequence (Fig. 8a) for the yes-no qubit at its resonance
frequency for various coherent states in the storage mode. Revivals happen
every $1/\chi_{\mathrm{s,yn}}$ which allows us to set $\Delta\tau$ as half the
revival time in Fig. 8b. The signal difference between the final $-\pi/2$ and
$+\pi/2$ pulses can be expressed as

$S(t)=\dfrac{\overline{\sigma_{\mathrm{z,yn}}}^{+}-\overline{\sigma_{\mathrm{z,yn}}}^{-}}{2}=e^{|\alpha|^{2}(\mathrm{cos}(2\pi\chi_{\mathrm{s,yn}}t)-1)}\mathrm{cos}(|\alpha|^{2}\mathrm{sin}(2\pi\chi_{\mathrm{s,yn}}t))e^{-t\Gamma_{\mathrm{2,yn}}-\gamma|\alpha|^{2}t}$
(4)

This expression is derived in the supplementary material of Ref. (?). 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 (?). At first order this shift is given by

$t_{\mathrm{revival}}=2\Delta\tau\left(1+2|\alpha|^{2}\chi_{\mathrm{s,s,yn}}\Delta\tau\right).$
(5)

Finding the above parameters that allow the model to match the measured signal
shown in Fig. 8b, we find $\Delta\tau=337~{}\mathrm{ns}$,
$\gamma=0.23~{}\mathrm{\mu s^{-1}}$ and $\chi_{\mathrm{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 take these values as accurate enough
compared to what we can obtain with the other methods presented in this work.

####  B.3 Ramsey oscillations of the storage mode

For a qubit, Ramsey oscillations correspond to the evolution of the real part
of the coherence between the $|g\rangle$ and $|e\rangle$ states. A typical
sequence starts by a $\pi/2$ pulse detuned from resonance by $\delta\\!f$ to
create a coherent superposition of $|g\rangle$ and $|e\rangle$ states. Then
the qubit is let to evolve freely before its state tomography is performed.
Both $\sigma_{\mathrm{x}}$ and $\sigma_{\mathrm{y}}$ oscillate at $\delta\\!f$
while decaying at the decoherence rate $\Gamma_{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. (?)).
The first $\pi/2$ pulse is replaced by a detuned displacement pulse $D(\beta)$
on the storage mode. The field then evolves freely during a time $t$ before a
Wigner tomography is realized. The expectation value of
$\hat{X}=(\hat{a}_{\mathrm{s}}+\hat{a}^{\dagger}_{\mathrm{s}})/2$ and
$\hat{P}=(\hat{a}_{\mathrm{s}}-\hat{a}^{\dagger}_{\mathrm{s}})/2i$ quadratures
are computed from the Wigner tomography.

Indeed for any operator $\hat{O}$, we can apply the Wigner transform to obtain
the operator Wigner map $W_{\hat{O}}$ (?) as

$W_{\hat{O}}(\alpha=x+ip)=W_{\hat{O}}(x,p)=\dfrac{1}{\pi}\int\mathrm{d}y\,e^{-2ipy}\langle
x+y/2|\hat{O}|x-y/2\rangle=\frac{2}{\pi}\,\mathrm{Tr}(D^{\dagger}(\alpha)\hat{O}D(\alpha)\mathcal{P})$
(6)

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

$\begin{array}[]{lllll}\pi\int\mathrm{d}x\int\mathrm{d}p\,W_{\rho}(x,p)W_{\hat{O}}(x,p)\\\
=\frac{1}{\pi}\int\mathrm{d}x\int\mathrm{d}p\int\mathrm{d}y\int\mathrm{d}y\prime\,e^{-2ip(y+y\prime)}\langle
x+y/2|\rho|x-y/2\rangle\langle x+y\prime/2|\hat{O}|x-y\prime/2\rangle\\\
=\int\mathrm{d}x\int\mathrm{d}y\int\mathrm{d}y\prime\,\delta(y+y\prime)\langle
x+y/2|\rho|x-y/2\rangle\langle x+y\prime/2|\hat{O}|x-y\prime/2\rangle\\\
=\int\mathrm{d}x\int\mathrm{d}y\,\langle x+y/2|\rho|x-y/2\rangle\langle
x-y/2|\hat{O}|x+y/2\rangle=\int\mathrm{d}u\int\mathrm{d}v\,\langle
u|\rho|v\rangle\langle v|\hat{O}|u\rangle\\\
=\mathrm{Tr}(\rho\hat{O})=\langle\hat{O}\rangle_{\rho}\end{array}.$ (7)

Figure 9: Ramsey oscillations of the storage mode a. Circuit diagram for
Ramsey oscillations of an harmonic oscillator. All storage displacement pulses
are performed in 100 ns with a Gaussian envelope of 25 ns width. For this
experiment $\beta$ is the amplitude of the prepared coherent state and is set
to -1.55. The detuning between displacement pulse and bare storage frequencies
is $\delta\\!f^{0}_{\mathrm{s}}=3.96$ MHz. b. Measured (dots) and expected
(lines) signals for $\langle\hat{X}\rangle$ (blue) and $\langle\hat{P}\rangle$
(orange). The expected signals are matched to the experiment using Eq. (11),
when we set a frequency detuning of $\delta\\!f_{\mathrm{s}}=3.96$ MHz and a
decay rate $\Gamma_{\mathrm{2,s}}=(2~{}\mu\mathrm{s})^{-1}$.

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

$\begin{array}[]{ll}W_{\hat{X}}(\alpha=x+ip)=x/\pi\\\
W_{\hat{P}}(\alpha=x+ip)=p/\pi\end{array}.$ (8)

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

$\begin{array}[]{ll}\langle\hat{X}\rangle=\int\mathrm{d}x\int\mathrm{d}p\,W(x,p)x\\\
\langle\hat{P}\rangle=\int\mathrm{d}x\int\mathrm{d}p\,W(x,p)p\end{array}.$ (9)

The time trace of $\langle\hat{X}\rangle$ and $\langle\hat{P}\rangle$ is what
we call the Ramsey oscillations for the storage mode. As in the qubit case,
the frequency of the oscillations is set by the detuning
$\delta\\!f_{\mathrm{s}}$ between the drive and the resonant frequency of the
mode, which allows us to extract the frequency of the storage mode. At this
point, a distinction has to be made between the detuning
$\delta\\!f^{0}_{\mathrm{s}}=f_{\mathrm{drive}}-f_{s}$ between the drive and
the bare storage mode frequency (the resonant frequency when the multiplexed
qubit and the storage are undriven) and the detuning $\delta\\!f_{\mathrm{s}}$
between the drive and the resonant frequency of the storage mode, which
depends on the multiplexed measurement strength in perfect analogy with the
AC-Stark effect for a qubit readout. Note that the Wigner tomography sequence
uses the same detuned frequency $\delta\\!f_{\mathrm{s}}$ for its displacement
pulse $D(\alpha)$ in order to keep the same phase reference. The measurement
of Ramsey oscillations of a harmonic oscillator takes longer than the ones of
a qubit because we fully determine the quantum state of an oscillator at each
time instead of a simple Bloch vector. From Eq. (2), one can compute the time
derivative of $\langle a\rangle$ and find

$\langle\dot{\hat{a}}_{\mathrm{s}}\rangle=(2i\pi\delta\\!f_{\mathrm{s}}-\Gamma_{\mathrm{2,s}})\langle\hat{a}_{\mathrm{s}}\rangle.$
(10)

This differential equation can be easily solved and taking the real and
imaginary parts, one finds that $\langle\hat{X}\rangle$ and
$\langle\hat{P}\rangle$ evolve as

$\begin{array}[]{ll}\langle\hat{X}\rangle=|\beta|\mathrm{cos}(2\pi\delta\\!f_{\mathrm{s}}t+\phi)e^{-t\Gamma_{\mathrm{2,s}}}\\\
\langle\hat{P}\rangle=|\beta|\mathrm{sin}(2\pi\delta\\!f_{\mathrm{s}}t+\phi)e^{-t\Gamma_{\mathrm{2,s}}}\end{array}$
(11)

where $\beta=|\beta|e^{i\phi}=\langle\hat{a}_{s}\rangle(t=0)$. For each time
$t$, we computed $\langle\hat{X}\rangle$ and $\langle\hat{P}\rangle$ and
define the storage mode decoherence rate as $\Gamma_{\mathrm{2,s}}$. Fig. 9
shows an example of measured Ramsey oscillations.

In the main text, Fig. 3B does not exhibit oscillations because it is the mean
value $\langle\hat{a}_{\mathrm{s}}+\hat{a}_{\mathrm{s}}^{\dagger}\rangle$ in
the frame rotating at the resonant frequency of the storage mode. In practice,
we measured them with a detuning shown in Fig. 10c and numerically computed
the non-oscillating quantity
$2\mathcal{R}e((\langle\hat{X}\rangle+i\langle\hat{P}\rangle)\mathrm{exp}(-2i\pi\delta\\!f_{\mathrm{s}}t))$.

####  B.4 Storage mode frequency shift and induced dephasing rate by driving a
comb on the multiplexing qubit

By analogy with the AC Stark effect that shifts a qubit frequency when its
readout cavity is driven, we will call AC Stark shift the detuning of the
storage mode when the multiplexing qubit is driven. In order to measure this
frequency shift and to measure 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 decided to use a Ramsey interferometry
sequence for the storage mode while the multiplexing qubit is continuously
driven by a frequency comb during the time of free evolution $t$ (see Fig.
10a). The frequency comb is a Gaussian envelope of duration $t$ and width
$t/4$ multiplied by the sum of nine sine functions at the frequencies
$[f_{\mathrm{mp}},f_{\mathrm{mp}}-\chi_{\mathrm{s,mp}},...,f_{\mathrm{mp}}-8\chi_{\mathrm{s,mp}}]$.

For small measurement strength ($\Omega/\chi_{\mathrm{s,mp}}<0.9$) we
generated the Ramsey sequence with a displacement pulse detuned from the bare
storage mode frequency by $\delta\\!f^{0}_{\mathrm{s}}=3.96$ MHz, and an
amplitude $\beta=-1.55$. We fit the time evolution of $\langle\hat{X}\rangle$
and $\langle\hat{P}\rangle$ using the damped sine function

$\begin{array}[]{ll}\langle\hat{X}\rangle=A\mathrm{cos}(2\pi\delta\\!f_{\mathrm{s}}t+\phi)e^{-t\Gamma_{\mathrm{d,s}}}\\\
\langle\hat{P}\rangle=A\mathrm{sin}(2\pi\delta\\!f_{\mathrm{s}}t+\phi)e^{-t\Gamma_{\mathrm{d,s}}}\end{array}.$
(12)

The parameters $A$, $\delta\\!f_{\mathrm{s}}$, $\phi$, and
$\Gamma_{\mathrm{d,s}}$ are determined altogether by fitting the model to the
measured oscillations. $\Gamma_{\mathrm{d,s}}$ is interpreted as the sum of
the intrinsic storage dephasing rate $\Gamma_{\mathrm{2,s}}$ and of the
measurement induced dephasing rate. $\delta\\!f_{\mathrm{s}}$ is the sum of
the detuning from the bare storage mode frequency
$\delta\\!f^{0}_{\mathrm{s}}$ and of the AC stark shift of the storage mode.
Both parameters are shown in Fig. 10c,d (blue dots) as a function of
$\Omega/\chi_{\mathrm{s,mp}}$.

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

$\begin{array}[]{ll}\langle\hat{X}\rangle=A(\mathrm{cos}(2\pi\delta\\!f_{\mathrm{s}}t+\phi)+\zeta\mathrm{cos}(2\pi\nu
t+\psi_{\mathrm{X}}))e^{-t\Gamma_{\mathrm{d,s}}}\\\
\langle\hat{P}\rangle=A(\mathrm{sin}(2\pi\delta\\!f_{\mathrm{s}}t+\phi)+\zeta\mathrm{sin}(2\pi\nu
t+\psi_{\mathrm{P}}))e^{-t\Gamma_{\mathrm{d,s}}}\end{array}.$ (13)

This empirical model originates from three ideas. The first term is identical
to the simple model in Eq. (12). 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. 10b shows an example of Ramsey oscillations of the storage mode with a
large amplitude of measurement. The two signals are used altogether to extract
the parameters A, $\delta\\!f_{\mathrm{s}}$, $\nu$, $\phi$,
$\psi_{\mathrm{X}}$, $\psi_{\mathrm{P}}$, and $\Gamma_{\mathrm{d,s}}$. The
frequency $\nu$ varies from 2.15 MHz to 2.5 MHz. The parameter $\zeta$ is
roughly constant, it varies between 0.2 to 0.27. We decompose
$\Gamma_{\mathrm{d,s}}$ as the intrinsic storage dephasing rate
$\Gamma_{\mathrm{2,s}}$ plus the measurement induced dephasing rate and
$\delta\\!f_{\mathrm{s}}$ as the detuning from the storage frequency
$\delta\\!f^{0}_{\mathrm{s}}$ plus the measurement induced frequency shift of
the storage mode. Fig. 10c shows measurement induced detuning as a function of
measurement drive amplitude.

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

####  B.5 Calibration of Rabi frequency $\Omega$

In order to realize the calibration of the Rabi frequency $\Omega$ for the
multiplexing qubit we recorded Rabi oscillations and measured their frequency
for various drive amplitudes. We played a square pulse of $1~{}\mu$s on the
multiplexing qubit at its bare frequency $f_{\mathrm{mp}}$ with an amplitude
$V_{\mathrm{mp}}$. We demodulated the reflected pulse by time steps of $10$ ns
at $f_{\mathrm{mp}}$. The reflection coefficient shows damped Rabi
oscillations at a frequency that depends on $V_{\mathrm{mp}}$. We matched to
the measurement the Rabi oscillations for seven different amplitudes
$V_{\mathrm{mp}}$ using the model (?)

$\mathcal{R}e\left(r(t)\right)-\mathcal{R}e\left(r_{\mathrm{ss}}\right)=A~{}\mathrm{cos}\left[\sqrt{(2\pi\xi
V_{\mathrm{mp}})^{2}-\left(\dfrac{\Gamma_{\mathrm{1,mp}}-2\Gamma_{\mathrm{2,mp}}}{16}\right)^{2}}(t-t0)+\phi\right]e^{-(t-t0)/T},$
(14)

where $r_{\mathrm{ss}}$ is the value of the reflection coefficient in the
steady state. We obtain $\xi=0.543\mathrm{~{}GHz.V^{-1}~{}}$ so that the Rabi
frequency is thus calibrated as $\Omega=\xi
V_{\mathrm{mp}}=0.543\mathrm{~{}GHz.V^{-1}~{}}V_{\mathrm{mp}}$ (see Fig. 11).

Figure 11: Multiplexing qubit Rabi oscillations for various driving amplitude.
The measured Rabi oscillations observed in the reflection coefficient (dots)
is reproduced by theory (solid line from Eq. (14). The vertical axis
represents the deviation of the real part of the reflection coefficient to its
steady state value. This calibration allows us to extract the scaling
parameter $\xi$ such that $\Omega=\xi
V_{\mathrm{mp}}=(0.543\mathrm{~{}GHz.V^{-1}~{}})V_{\mathrm{mp}}$.

###  C Qutip 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_ (?). We simulated
the complete system composed of the storage mode, the yes-no qubit and the
multiplexing qubit with all couplings, except in the case 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 qubits 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.

####  C.1 Photocounting simulations

C.1.1 Photocounting with the yes-no qubit

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

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

$\begin{array}[]{rcl}\hat{H_{1}}/h&=&\delta\\!f_{\mathrm{yn}}\dfrac{\hat{\sigma}_{\mathrm{z,yn}}}{2}-\chi_{\mathrm{s,yn}}\hat{n}_{\mathrm{s}}\dfrac{\hat{\sigma}_{\mathrm{z,yn}}}{2}-\chi_{\mathrm{s,mp}}\hat{n}_{\mathrm{s}}\dfrac{\hat{\sigma}_{\mathrm{z,mp}}}{2}-\chi_{\mathrm{s,s}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\\\
&&-\chi_{\mathrm{s,s,yn}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\dfrac{\hat{\sigma}_{\mathrm{z,yn}}}{2}-\chi_{\mathrm{s,s,mp}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\dfrac{\hat{\sigma}_{\mathrm{z,mp}}}{2}+\frac{\epsilon_{\mathrm{yn}}(t)}{h}\hat{\sigma}_{\mathrm{x,yn}}\\\
&&+\frac{\lambda(t)}{2\pi}(\epsilon_{\mathrm{max}}e^{i\pi(\chi_{\mathrm{s,mp}}+\chi_{\mathrm{s,yn}})t}\hat{a}_{\mathbf{s}}+\epsilon_{\mathrm{max}}^{\ast}e^{-i\pi(\chi_{\mathrm{s,mp}}+\chi_{\mathrm{s,yn}})t}\hat{a}^{\dagger}_{\mathbf{s}})\end{array},$
(15)

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

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

$\begin{array}[]{rcl}\dot{\rho}&=&-\dfrac{i}{\hbar}[\hat{H_{1}},\rho]+2\Gamma_{\mathrm{\phi,s}}\mathcal{L}(\hat{n}_{\mathrm{s}})\rho+(1+n_{\mathrm{th,s}})\Gamma_{\mathrm{1,s}}\mathcal{L}(\hat{a}_{\mathrm{s}})\rho+n_{\mathrm{th,s}}\Gamma_{\mathrm{1,s}}\mathcal{L}(\hat{a}^{\dagger}_{\mathrm{s}})\rho\\\
&&+\dfrac{1}{2}\Gamma_{\mathrm{\phi,yn}}\mathcal{L}(\hat{\sigma}_{\mathrm{z,yn}})\rho+\Gamma_{\mathrm{1,yn}}\mathcal{L}(\hat{\sigma}^{-}_{\mathrm{yn}})\rho+\dfrac{1}{2}\Gamma_{\mathrm{\phi,mp}}\mathcal{L}(\hat{\sigma}_{\mathrm{z,mp}})\rho+\Gamma_{\mathrm{1,mp}}\mathcal{L}(\hat{\sigma}^{-}_{\mathrm{mp}})\rho\end{array},$
(16)

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

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

This simulation can be used to reproduce the experiment in Fig. 2A,C of the
main text by adjusting the following parameters
$\left\\{\mu=\epsilon_{\mathrm{max}}/V_{\mathrm{max,s}},\chi_{\mathrm{s,yn}},\chi_{\mathrm{s,s}},\chi_{\mathrm{s,s,yn}},n_{\mathrm{th,s}}\right\\}$.
Note that we need to run the simulation for every couple of parameters
($V_{\mathrm{max,s}}$,$\delta\\!f_{\mathrm{yn}}$). The table Tab. 3 compiles
the values of fitted parameters.

C.1.2 Photocounting with the multiplexing qubit

A second simulation (fluorescence simulation) was carried out to compare the
photon counting experiment in Fig. 2B,D 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\mathrm{~{}\mu s}$ Gaussian pulse
on the multiplexing qubit at the frequency
$f_{\mathrm{mp}}-\delta\\!f_{\mathrm{mp}}$ with an amplitude expressed as a
Rabi frequency $\Omega=\chi_{\mathrm{s,mp}}/4$. The measured reflection
coefficient of the multiplexing qubit $r(\delta\\!f_{\mathrm{mp}})$ can be
expressed using input-output theory as (?)

$r(\delta\\!f_{\mathrm{mp}})=\frac{\langle a_{\mathrm{out}}\rangle}{\langle
a_{\mathrm{in}}\rangle}=\frac{\langle
a_{\mathrm{in}}\rangle-\sqrt{\Gamma_{\mathrm{1,mp}}}\langle\sigma_{-,\mathrm{mp}}\rangle}{\langle
a_{\mathrm{in}}\rangle}=1-\frac{\sqrt{\Gamma_{\mathrm{1,mp}}}}{\langle
a_{\mathrm{in}}\rangle}\langle\hat{\sigma}_{\mathrm{-,mp}}\rangle.$

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

$1-\mathcal{R}e\left(r(\delta\\!f_{\mathrm{mp}})\right)=\frac{\Gamma_{\mathrm{1,mp}}}{\pi\Omega}\mathcal{R}e\left(e^{-i\arg(\langle
a_{\mathrm{in}}\rangle)}\langle\hat{\sigma}_{\mathrm{-,mp}}\rangle\right)$

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

$1-\mathcal{R}e\left(r(\delta\\!f_{\mathrm{mp}})\right)=\frac{\Gamma_{\mathrm{1,mp}}}{2\pi\Omega}\langle\hat{\sigma}_{\mathrm{y,mp}}\rangle.$

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

$\begin{array}[]{rcl}\hat{H_{2}}/h&=&\delta\\!f_{\mathrm{mp}}\dfrac{\hat{\sigma}_{\mathrm{z,mp}}}{2}-\chi_{\mathrm{s,yn}}\hat{n}_{\mathrm{s}}\dfrac{\hat{\sigma}_{\mathrm{z,yn}}}{2}-\chi_{\mathrm{s,mp}}\hat{n}_{\mathrm{s}}\dfrac{\hat{\sigma}_{\mathrm{z,mp}}}{2}-\chi_{\mathrm{s,s}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\\\
&&-\chi_{\mathrm{s,s,yn}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\dfrac{\hat{\sigma}_{\mathrm{z,yn}}}{2}-\chi_{\mathrm{s,s,mp}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\dfrac{\hat{\sigma}_{\mathrm{z,mp}}}{2}+\dfrac{\Omega}{2}\epsilon_{\mathrm{mp}}(t)\hat{\sigma}_{\mathrm{x,mp}}\\\
&&+\frac{\lambda(t)}{2\pi}(\epsilon_{\mathrm{max}}e^{i\pi(\chi_{\mathrm{s,mp}}+\chi_{\mathrm{s,yn}})t}\hat{a}_{\mathbf{s}}+\epsilon_{\mathrm{max}}^{\ast}e^{-i\pi(\chi_{\mathrm{s,mp}}+\chi_{\mathrm{s,yn}})t}\hat{a}^{\dagger}_{\mathbf{s}}),\end{array}$
(17)

where $\epsilon_{\mathrm{mp}}(t)\geq 0$ is a Gaussian function of duration 2
$\mu$s, width 250 ns and with an amplitude 1. $\epsilon_{\mathrm{mp}}(t)$ is
delayed compare to $\lambda(t)$ to reproduce the experimental pulse sequence.
We add to this Hamiltonian the same relaxation and decoherence channels as for
the yes-no simulation (see Eq. (16)) 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\mathrm{~{}ns}$ starting
from a thermal state with $n_{\mathrm{th,s}}$ photons for storage and the yes-
no qubit and the multiplexing qubit in the ground state $|g\rangle$. Finally,
the expectation value $\langle\hat{\sigma}_{\mathrm{y,mp}}\rangle$ is computed
and integrated during the $2~{}\mathrm{\mu s}$ of the pulse.

We compare the measured emission coefficient in Fig. 2B,D to the simulated
signal $A\langle\hat{\sigma}_{\mathrm{y,mp}}\rangle$ where $A$ is left as a
free parameter due to a small parasitic reflection in the measurement setup
and thermal population. The parameters
$\left\\{\mu,\chi_{\mathrm{s,t}},\chi_{\mathrm{s,s}},\chi_{\mathrm{s,s,t}},n_{\mathrm{th,s}}\right\\}$
is already set by the calibration above using the simulation of the yes-no
qubit. From the fluorescence simulation, we thus extract the parameters
$\left\\{\chi_{\mathrm{s,mp}},\chi_{\mathrm{s,s,mp}},A\right\\}$ by comparing
the experimental observation in Fig. 2B,D with the simulation for various
$V_{\mathrm{max,s}}$ and $\delta\\!f_{\mathrm{mp}}$. Fitted values are given
in Tab. 3. 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.

parameter | fitted values
---|---
$\mu$ | 1.45 (mV.$\mu$s)-1
$\chi_{\mathrm{s,yn}}$ | 1.42 MHz
$\chi_{\mathrm{s,mp}}$ | 4.9 MHz
$\chi_{\mathrm{s,s}}$ | -0.02 MHz
$\chi_{\mathrm{s,s,yn}}$ | -0.003 MHz
$\chi_{\mathrm{s,s,mp}}$ | -0.08 MHz
$n_{\mathrm{th,s}}$ | 0.03
Table 3: Parameters extracted from the photocounting simulations using the
multiplexing or yes-no qubit. All parameters except those related to the
multiplexing qubit are determined using a fit of the yes-no qubit simulation
to the Fig. 2A,C. Parameters related to multiplexing qubit are obtained using
a fit of the simulation to the Fig. 2B,D.

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

We simulated the filling of the storage mode by a displacement pulse on the
cavity. We simulated the same master equation used for the photocounting
simulations with parameters obtained from the photocounting simulations (see
Tab. 3) but without applying any drive on the qubits. Only the displacement
pulse on the storage mode is modeled i.e. $\epsilon_{\mathrm{mp}}(t)=0$,
$\delta\\!f_{\mathrm{mp}}=0$, $\epsilon_{\mathrm{yn}}(t)=0$, and
$\delta\\!f_{\mathrm{yn}}=0$.

The ”mesolve” function of QuTiP computes the density matrix with a time step
of $10~{}\mathrm{ns}$ and returns the mean number of photons in the storage
mode at the end of the displacement pulse for various drive amplitudes. Fig. 6
shows the square root of the expected mean photon number as a function of the
amplitude $\epsilon_{\mathrm{max}}$. We obtain a scaling factor
$\sqrt{\left<n_{\mathrm{s}}\right>}=85.9~{}\mathrm{V^{-1}}V_{\mathrm{max,s}}$
used in the photon number calibration of the storage mode.

####  C.3 Simulation of multiplexed readout

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

$\begin{array}[]{rcl}\hat{H_{3}}/h&=&-\chi_{\mathrm{s,yn}}\hat{n}_{\mathrm{s}}\dfrac{\hat{\sigma}_{\mathrm{z,yn}}}{2}-\chi_{\mathrm{s,mp}}\hat{n}_{\mathrm{s}}\dfrac{\hat{\sigma}_{\mathrm{z,mp}}}{2}-\chi_{\mathrm{s,s}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)-\chi_{\mathrm{s,s,yn}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\dfrac{\hat{\sigma}_{\mathrm{z,yn}}}{2}\\\
&&-\chi_{\mathrm{s,s,mp}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\dfrac{\hat{\sigma}_{\mathrm{z,mp}}}{2}+\dfrac{\Omega}{2}(\epsilon_{\mathrm{comb}}(t)\hat{\sigma}^{+}_{\mathrm{mp}}+\epsilon_{\mathrm{comb}}^{\ast}(t)\hat{\sigma}^{-}_{\mathrm{mp}})\\\
&&+\frac{\lambda(t)}{2\pi}(\epsilon_{\mathrm{max}}e^{i\pi(\chi_{\mathrm{s,mp}}+\chi_{\mathrm{s,yn}})t}\hat{a}_{\mathbf{s}}+\epsilon_{\mathrm{max}}^{\ast}e^{-i\pi(\chi_{\mathrm{s,mp}}+\chi_{\mathrm{s,yn}})t}\hat{a}^{\dagger}_{\mathbf{s}})\end{array},$
(18)

where $\Omega=\chi_{\mathrm{s,mp}}/2$ and $\epsilon_{\mathrm{comb}}(t)$ is the
product of a Gaussian function with the sum of nine complex exponential
$\sum_{k=0}^{8}\mathrm{exp}(2i\pi\chi_{\mathrm{s,mp}}kt)$. The Gaussian
envelope of $\epsilon_{\mathrm{comb}}(t)$ has a duration of $2~{}\mu s$, a
width of $250\mathrm{~{}ns}$, and a maximum amplitude of $1$ and the delay
between $\epsilon_{\mathrm{comb}}(t)$ and $\lambda(t)$ reproduces the
experimental sequence. The master equation (16) is used with a time step of
$1\mathrm{~{}ns}$ for various amplitude $\epsilon_{\mathrm{max}}$. We obtain
the time evolution of $\langle\sigma_{y,mp}\rangle$ enabling us to compare the
experimental measurements of Fig. 2E to the model. To do so, we integrate the
simulated function
$\langle\sigma_{y,mp}\rangle\times\mathrm{cos}(2\pi\chi_{\mathrm{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 (it is
easy to understand it when demultiplexing a signal of the form $\cos(\omega
t)$ when $\omega=0$ or not). 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.

####  C.4 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_{\mathrm{s}}+\delta\\!f_{\mathrm{s}}^{0}$ for the storage mode is

$\begin{array}[]{rcl}\hat{H_{4}}/h&=&-\chi_{\mathrm{s,mp}}\dfrac{\hat{\sigma}_{\mathrm{z,mp}}+\mathbbm{1}}{2}\hat{n}_{\mathrm{s}}-\delta\\!f_{\mathrm{s}}^{0}\hat{n}_{\mathrm{s}}-\chi_{\mathrm{s,s,mp}}\hat{n}_{\mathrm{s}}(\hat{n}_{\mathrm{s}}-1)\dfrac{\hat{\sigma}_{\mathrm{z,mp}}+\mathbbm{1}}{2}\\\
&&+\dfrac{\Omega}{2}(\epsilon_{\mathrm{comb}}(t)\hat{\sigma}^{+}_{\mathrm{mp}}+\epsilon_{\mathrm{comb}}^{\ast}(t)\hat{\sigma}^{-}_{\mathrm{mp}}),\end{array}$
(19)

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

$\dot{\rho}=-\dfrac{i}{\hbar}[\hat{H_{4}},\rho]+2\Gamma_{\mathrm{\phi,s}}\mathcal{L}(\hat{n}_{\mathrm{s}})\rho+\Gamma_{\mathrm{1,s}}\mathcal{L}(\hat{a}_{\mathrm{s}})\rho+\dfrac{1}{2}\Gamma_{\mathrm{\phi,mp}}\mathcal{L}(\hat{\sigma}_{\mathrm{z,mp}})\rho+\Gamma_{\mathrm{1,mp}}\mathcal{L}(\hat{\sigma}^{-}_{\mathrm{mp}})\rho.$
(20)

The storage is initialized in a coherent state of amplitude $\beta=1.55$ and
the multiplexing qubit is initialized in state $|g\rangle$. We simulate the
dynamics of the system for a pulse duration $t$ going from $100\mathrm{~{}ns}$
to $5\mathrm{~{}\mu s}$ and for $\Omega$ ranging from $0$ to
$2\chi_{\mathrm{s,mp}}$. We compute the expectation value of
$\hat{X}=(\hat{a}_{\mathrm{s}}+\hat{a}^{\dagger}_{\mathrm{s}})/2$ at the end
of each simulation. For a given $\Omega$, we extract the time evolution of
$\langle\hat{X}\rangle$ under the influence of the multiplexed measurement as
shown on Fig. 12a. This decaying sinusoid is fitted using Eq. (12) to obtain
the oscillation frequency $\delta\\!f_{\mathrm{s}}$ and the decay rate
$\Gamma_{\mathrm{d,s}}$. Fig. 10c and d show the measurement induced dephasing
and AC Stark shift as a function of amplitude of the comb $\Omega$ for two
sets of coherent state amplitudes $\beta$ and detuning
$\delta\\!f_{\mathrm{s}}^{0}$.

We identify three interesting features. The first one is the evolution of the
shape of the curves $\delta\\!f_{\mathrm{s}}(\Omega)$ and
$\Gamma_{\mathrm{d,s}}(\Omega)$ with $\chi_{\mathrm{s,mp}}$. We repeat the
simulation using a square pulse envelope instead of Gaussian pulse for
$\epsilon_{\mathrm{comb}}$ to make the simulation faster for several values of
$\chi_{\mathrm{s,mp}}$ from $1.5$ to $13.2\mathrm{~{}MHz}$ by steps of
$1.4\mathrm{~{}MHz}$. We observe that $\delta\\!f_{\mathrm{s}}(\Omega)$ and
$\Gamma_{\mathrm{d,s}}(\Omega)$ increase as $\chi_{\mathrm{s,mp}}$ becomes
larger but that the maxima and minima of the curve are always found for the
same $\Omega/\chi_{\mathrm{s,mp}}$ ratio (Fig. 12b).

The second observation is that $\delta\\!f_{\mathrm{s}}(\Omega)$ and
$\Gamma_{\mathrm{d,s}}(\Omega)$ vary with the initial coherent state amplitude
$\beta$ (Fig. 12c).

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

Figure 12: Simulations of the measurement induced dephasing rate and of the AC
Stark shift induced by a frequency comb. a. Ramsey-like oscillations of the
storage mode for $\Omega=\chi_{s,mp}/2$ and an initial coherent field
amplitude $\beta=-1.55$. Blue dots are the simulated expectation values of
$\hat{X}$ and red line is the theory given by the Eq. (12). b. Simulated
measurement induced dephasing rate $\Gamma_{\mathrm{d,s}}$ and AC Stark shift
as a function of $\Omega/\chi_{s,mp}$ for various values of
$\chi_{\mathrm{s,mp}}$. Simulations show the same pattern with maxima and
minima for some specific values of $\Omega/\chi_{s,mp}$ as in the experiment
in Fig. 3C. c. Simulated measurement induced dephasing rate and AC Stark shift
as a function of $\Omega/\chi_{s,mp}$ for various initial coherent state
amplitudes $\beta$ in the storage mode. For $\Omega/\chi_{s,mp}>1$, we see a
difference of about $20~{}\%$ between $\beta=1.6$ and $\beta=1.2$. d.
Simulations for $\chi_{\mathrm{s,mp}}\gg\Gamma_{\mathrm{1,mp}}$. The evolution
of the measurement induced dephasing and AC Stark shift with $\Omega$ is
different compared to the case of b. The evolution of the measurement induced
dephasing rate and the AC Stark shift seems to be given by the ratio
$2\pi\Omega/\Gamma_{\mathrm{1,mp}}$. The red line is a guide for eyes
representing a square function. It shows that the measurement induced
dephasing rate increases linearly with $\Omega^{2}$ for small drive
amplitudes. For the AC Stark shift, on the contrary 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 $\Omega$. e. Example of simulated Ramsey oscillations
exhibiting two frequencies.

###  D 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.

####  D.1 Density matrix reconstruction

The Wigner tomography contains all the information about the state of the
storage mode. We explain below how we reconstruct the density matrix from the
measured Wigner function. We compute the Wigner map for every operator
$|n\rangle\langle m|$ with $|n\rangle$ and $|m\rangle$ two fock states with
$n$ and $m$ photons. The mean value of those operators is equal to the $(n,m)$
element $\rho_{\mathrm{nm}}$ of the density matrix. Using the mathematical
expression of $\langle x|n\rangle$

$\langle
x|n\rangle=\psi_{n}(x)=\left(\dfrac{2}{\pi}\right)^{1/4}\dfrac{1}{\sqrt{2^{n}n!}}\mathcal{H}_{n}(\sqrt{2}x)e^{-x^{2}}$
(21)

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

$W_{|n\rangle\langle
m|}(x,p)=\dfrac{1}{\pi}\int\mathrm{d}ye^{-2ipy}\psi_{n}(x+y/2)\overline{\psi_{m}}(x-y/2)$
(22)

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

$\rho_{\mathrm{nm}}=\pi\iint\mathrm{d}x\mathrm{d}pW_{|n\rangle\langle
m|}(x,p)W_{\rho}(x,p).$ (23)

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 the
quantity $|\rho_{\mathrm{nm}}|/\sqrt{\rho_{\mathrm{nn}}\rho_{\mathrm{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_{\mathrm{s}}$

$\dot{\rho}=\Gamma_{1,s}\mathcal{L}(\hat{a}_{\mathrm{s}})\rho+2\Gamma_{\phi,s}\mathcal{L}(\hat{a}_{\mathrm{s}}^{\dagger}\hat{a}_{\mathrm{s}})\rho.$
(24)

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

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

If the storage mode is initialized in a coherent state $|\alpha_{o}\rangle$,
the solution of the equation is

$\rho_{\mathrm{nm}}(t)=e^{-|\alpha_{o}|^{2}e^{-\Gamma_{\mathrm{1,s}}t}}\dfrac{\alpha_{o}^{m}e^{-m\Gamma_{\mathrm{1,s}}t/2}(\alpha_{o}^{\ast})^{n}e^{-n\Gamma_{\mathrm{1,s}}t/2}}{\sqrt{n!m!}}e^{-\Gamma_{\mathrm{\phi,s}}(n-m)^{2}t},$
(26)

and we get

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

Thus indeed, the renormalization removes the effect of the relaxation rate
$\Gamma_{1,s}$ and only characterizes the dephasing rate. Under the action of
measurement, we expect that the dephasing rate $\Gamma_{\mathrm{\phi,s}}$ is
increased by the measurement induced dephasing rate. In the following
sections, we will see that, indeed, the measurement process to distinguish
photon numbers $n$ and $m$ can well be captured by
$|\rho_{\mathrm{nm}}|/\sqrt{\rho_{\mathrm{nn}}\rho_{\mathrm{mm}}}=e^{-\Gamma_{\mathrm{d,s}}^{nm}t}$,
where $\Gamma_{\mathrm{d,s}}^{nm}$ is the sum of
$\Gamma_{\mathrm{\phi,s}}(n-m)^{2}$ and of a well-characterized measurement-
induced dephasing rate.

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

By essence of quantum mechanics, measuring with a single drive whether there
are $n$ photons induces a dephasing between Fock state $|n\rangle$ and Fock
state $|m\neq n\rangle$. This dephasing can be probed by observing the
dynamics of $\rho_{\mathrm{nm}}$.

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

In Ref. (?), 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. They show
that the decoherence rate between the Fock state $|n\rangle$ and $|m\rangle$
due to the measurement is given by the highest eigenvalue of the following
matrix

$\begin{pmatrix}-\Gamma_{\mathrm{1,mp}}/2&-2\pi\Delta+\tilde{\chi}&0&0\\\
2\pi\Delta-\tilde{\chi}&-\Gamma_{\mathrm{1,mp}}/2&-2\pi\Omega&0\\\
0&2\pi\Omega&-\Gamma_{\mathrm{1,mp}}&-\Gamma_{\mathrm{1,mp}}-i\tilde{\chi}\\\
0&0&-i\tilde{\chi}&0\end{pmatrix}$ (28)

with $\tilde{\chi}=\dfrac{n+m}{2}2\pi\chi_{\mathrm{s,mp}}$. Fig. 13 shows the
measured density matrix decoherence rates $\Gamma_{\mathrm{d,s}}^{nm}$ and the
above theory for $n$ and $m$ going from $0$ to $4$ (with an offset
corresponding to the natural dephasing rate in Eq. (27)). As expected in a
regime with resolved resonance peaks
$(2\pi\chi_{\mathrm{s,mp}}|m-n|>\Gamma_{\mathrm{2,mp}})$, the decoherence rate
$\Gamma_{\mathrm{d,s}}^{nm}$ between Fock states $|n\rangle$ and $|m\rangle$
is larger when the single drive probes whether there are $n$ photons or $m$
photons with a moderate drive amplitude $\Omega$ (dependence on $\Omega$ not
shown here). For much larger drive amplitude $\Omega$, one can increase the
decoherence rate $\Gamma_{\mathrm{d,s}}^{nm}$ by driving with a detuning
$\Delta_{\mathrm{mp}}=(n+m)\chi_{\mathrm{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 $\chi_{\mathrm{s,mp}}/2$. In fact this
regime would become particularly attractive for poorly resolved resonances as
a function of photon number
$(2\pi\chi_{\mathrm{s,mp}}|m-n|<\Gamma_{\mathrm{2,mp}})$. Premises of this
effect are visible on Fig. 13, as the maximal decoherence rate occurs at a
detuning slightly closer to $(n+m)\chi_{\mathrm{s,mp}}/2$, with a stronger
effect for small $|m-n|$, both in theory and in the experimental data. 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.

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

####  D.3 Multiplexed measurement vs single tone measurement

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

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 figure is based on the Wigner tomography of the
storage mode in three distinct cases. The storage mode is initialized in a
coherent state of amplitude $\beta=-1.55$, then, before performing the
tomography of the storage mode, we either (i) wait for a time $t$, (ii) probe
the multiplexing qubit for a time $t$ at a single frequency
$f_{\mathrm{mp}}-\chi_{\mathrm{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 $\rho(t)$ for various times $t$ for the three cases and compare
the evolution of the normalized elements $\rho_{\mathrm{nm}}(t)$ (see Fig.
14). 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 $\rho_{02}$ when probing with a resonant drive for $n=1$, is
consistent with the significant measurement-induced detuning that can be read
off the top right plot of Fig. 13 (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 the
section D.5.

####  D.4 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.

If a photocounter is QND, the diagonal elements of the density matrix 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 predominantly evolve owing to the decay of
the storage mode. The probability of having $0$ photon in the storage
increases in time (Fig. 15a) while the average photon number decreases with
time (Fig. 15c). Interestingly, we notice that for large probe amplitude (red
and purple points in Fig. 15a), the probability of finding the storage mode
with $0$ photon is slightly lower. This dependence on the amplitude of the
drive $\Omega$ is best characterized by extracting the populations (Fig. 15b)
and photon number (Fig. 15d) integrated during $T=5\mathrm{~{}\mu s}$ as a
function of $\Omega/\chi_{\mathrm{s,mp}}$. For small drive amplitude
$\Omega/\chi_{\mathrm{s,mp}}<0.1$, the probability to find the a given number
of photon does not change with $\Omega/\chi_{\mathrm{s,mp}}$ but for larger
drive amplitude the cavity gets populated probably because of deviations from
the dispersive approximation.

In practice, for small drive amplitude and a measurement time of
$5\mathrm{~{}\mu 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.

Figure 15: Impact of multiplexed measurement on the occupation of the Fock
states. a. Measured probability to find the storage mode in Fock state
$\left|0\right>$ as a function of time $t$ for various comb drive amplitudes
$\Omega/\chi_{\mathrm{s,mp}}$. b. Measured diagonal elements of the density
matrix integrated during $5\mathrm{~{}\mu s}$ as a function of drive amplitude
$\Omega/\chi_{\mathrm{s,mp}}$. c. Decay of the average photon number in time
for various drive amplitudes. d. Average photon number evolution as a function
of $\Omega/\chi_{\mathrm{s,mp}}$.

####  D.5 Off-diagonal density matrix elements and revivals of the coherence

In the main text we determined that the storage mode dephases faster in the
presence of measurement by measuring how much faster the Ramsey oscillations
decay. Ramsey oscillations give access to the decay of
$a_{\mathrm{s}}+a_{\mathrm{s}}^{\dagger}$, which is sensitive to the decay of
coherences of the form $\rho_{n,n+1}$. In fact it is also possible to
visualize the dynamics of individual off-diagonal elements of the storage
density matrix to gain insight into the physics of the dephasing process. For
that, we use Eq. (23) on the directly measured Wigner function of the storage
mode.

First, we study the decay of off-diagonal elements $\rho_{12}$ and $\rho_{13}$
(Fig. 16a and b). For small drive amplitude $\Omega<0.5\chi_{\mathrm{s,mp}}$,
off-diagonal elements decay faster when $\Omega$ is increased since more and
more information is extracted per unit time by the drive. As larger drive
amplitudes are reached, off-diagonal elements start oscillating. The contrast
of these coherence revivals become more pronounced as the drive amplitude
becomes larger and they exhibit a quasi periodicity. The time between the
first and second revivals is $10~{}\%$ shorter than the time between the
second and the third.

This behavior is qualitatively reproduced by our simulations (Fig. 16c). We
found that the amplitude of the revivals is strongly enhanced when the
coupling rate between the multiplexing qubit $\Gamma_{\mathrm{1,mp}}$ and the
transmission line is lowered. The revivals are periodic with a frequency equal
to $\Delta\\!f$, the frequency difference between two successive tones of the
comb. The presence of revivals in the coherences is a strong indication that
the dephasing process on the storage mode cannot be described as an open
quantum system coupled to a Markovian bath. This memory effect originates from
the finite lifetime of the multiplexing qubit. Pushing even further the
measurement amplitude $\Omega>0.9\chi_{\mathrm{s,mp}}$ in the experiment, we
observe that revivals remain but their time structure is no longer periodic
(see Fig. 17).

Figure 16: Normalized off-diagonal elements of the storage matrix density. a.
Measured normalised coherence $|\rho_{12}|$ between Fock states $|1\rangle$
and $|2\rangle$ of the storage mode as a function of time and for various
amplitudes $\Omega$ of the driving frequency comb. b. Similar plot for
$|\rho_{13}|$. c. Results of the simulation (20) for the normalized
$|\rho_{12}|$. Figure 17: Normalized off-diagonal elements of the storage
matrix density for the largest measurement amplitudes. Measured normalized
coherence $|\rho_{12}|$ between Fock states $|1\rangle$ and $|2\rangle$ of the
storage mode as a function of time and for various amplitudes
$\Omega>0.9\chi_{\mathrm{s,mp}}$ of the driving frequency comb.

####  D.6 Mean coherence between Fock states

Since the driving frequency comb holds the promise to probe how many photons
are in the storage mode, it should affect all coherences $\rho_{nm}$. In this
section, we introduce two 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
$\mathrm{Re}[(\langle\hat{X}\rangle+i\langle\hat{P}\rangle)e^{-2i\pi\delta\\!f_{\mathrm{s}}t}]$
with $\langle\hat{X}\rangle$ and $\langle\hat{P}\rangle$ 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_{\rho}$ between Fock
states $0$ to $4$. It is defined as

$C_{\rho}=\underset{4\geq i>j\geq
0}{\mathrm{Mean}}\left[\dfrac{|\rho_{ij}|}{\sqrt{\rho_{ii}\rho_{jj}}}\right].$
(29)

$C_{\rho(t)}$ contains the information about the dephasing between every
different Fock states. The left part of Fig. 18 shows oscillations of the
storage mode quadratures in the frame of the drive on the storage mode for
various drive amplitudes. On the right part of Fig. 18 we display the mean
quadrature $\langle\hat{a}_{\mathrm{s}}+\hat{a}_{\mathrm{s}}^{\dagger}\rangle$
in the frame rotating at the storage mode frequency and the mean coherence
$C_{\rho}$. 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. 16) also appear in the evolution of the quadrature and of the mean
coherence betwee Fock states.

Figure 18: Impact of the multiplexed drive on storage mode quadrature and mean
coherence between Fock states. Left column: Ramsey oscillations measured
(dots) and fits (line) using Eq. 12 and 13 for five multiplexing drive
amplitudes. Right column: quadrature of storage mode in the frame rotating at
its resonance frequency (blue dots) and mean coherence between Fock states
$C_{\rho}(t)$ (orange squares) defined in Eq. (29) for the five same drive
amplitudes. Incertiture are about $0.07$ for the quadrature and $0.015$ for
the mean coherence. The red solid lines are exponential fits of the quadrature
decay.

## References

  * 1. C. Guerlin, et al., Nature 448, 889 (2007).
  * 2. B. R. Johnson, et al., Nature Physics 6, 663 (2010).
  * 3. B. Peaudecerf, et al., Physical Review Letters 112, 080401 (2014).
  * 4. S. Haroche, M. Brune, J.M. Raimond, J. Phys. II France 2, 659 (1992).
  * 5. J. G. Proakis, M. Salehi, Communication Systems Engineering (Pentice Hall, 2002), second edn.
  * 6. Supplemental Material .
  * 7. D. I. Schuster, et al., Nature 445, 515 (2007).
  * 8. A. A. Houck, et al., Nature 449, 328 (2007).
  * 9. O. Astafiev, et al., Science (New York, N.Y.) 327, 840 (2010).
  * 10. A. A. Abdumalikov, O. V. Astafiev, Y. A. Pashkin, Y. Nakamura, J. S. Tsai, Phys. Rev. Lett. 107, 43604 (2011).
  * 11. P. Campagne-Ibarcq, et al., Phys. Rev. Lett. 112, 180402 (2014).
  * 12. M. F. Gely, et al., Science 363, 1072 (2019).
  * 13. C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, Atom-Photon Interactions: Basic Processes and Applications (2001).
  * 14. A. Wallraff, et al., Nature 431, 162 (2004).
  * 15. C. Macklin, et al., Science 350, 307 (2015).
  * 16. V. Schmitt, et al., Physical Review A 90, 062333 (2014).
  * 17. J. Heinsoo, et al., Physical Review Applied 10, 034040 (2018).
  * 18. S. Kundu, et al., Applied Physics Letters 114, 172601 (2019).
  * 19. F. Arute, et al., Nature 574, 505 (2019).
  * 20. A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, R. J. Schoelkopf, Reviews of Modern Physics 82, 1155 (2010).
  * 21. Boissonneault, Gambetta, Blais, Physical Review A 79, 13819 (2009).
  * 22. L. Lutterbach, L. Davidovich, Physical Review Letters 78, 2547 (1997).
  * 23. P. Bertet, et al., Physical Review Letters 89, 200402 (2002).
  * 24. B. Vlastakis, et al., Science 342, 607 (2013).
  * 25. P. Campagne-Ibarcq, et al., arXiv:1907.12487 (2019).
  * 26. H. M. Wiseman, G. J. Milburn, Quantum Measurement and Control (Cambridge University Press, 2009).
  * 27. S. Haroche, J. Raimond, Oxford Graduated Text p. 616 (2006).
  * 28. L. Bretheau, P. Campagne-Ibarcq, E. Flurin, F. Mallet, B. Huard, Science 348, 776 (2015).
  * 29. S. Kono, K. Koshino, Y. Tabuchi, A. Noguchi, Y. Nakamura, Nature Physics 14, 546 (2018).
  * 30. J.-C. Besse, et al., Physical Review X 8, 021003 (2018).
  * 31. R. Lescanne, et al., arXiv:1902.05102 (2019).
  * 32. Q. Ficheux, Quantum Trajectories with Incompatible Decoherence Channels, Ph.D. thesis, Universite Paris Sciences et Lettres (2018).
  * 33. J. Johansson, P. Nation, F. Nori, Computer Physics Communications 184, 1234 (2013).
  * 34. N. Cottet, Energy and Information in Fluorescence with Superconducting Circuits., Ph.D. thesis, PSL Research University; École Normale Supérieure (2018).
  * 35. A. Sarlette, P. Rouchon, A. Essig, Q. Ficheux, B. Huard, arXiv:2001.02550 (2020).