# Protecting collective qubits from non-Markovian dephasing

Antoine Covolo, Valentin Magro, Mathieu Girard, Sébastien Garcia, Alexei

Ourjoumtsev

► To cite this version:

Antoine Covolo, Valentin Magro, Mathieu Girard, Sébastien Garcia, Alexei Ourjoumtsev. Protecting collective qubits from non-Markovian dephasing. Optica, 2025, 12 (9), pp.1427. 10.1364/OP-TICA.563676. hal-05357497

HAL Id: hal-05357497

https://hal.science/hal-05357497v1

Submitted on 10 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Protecting collective qubits from non-Markovian dephasing

Antoine Covolo,¹ Valentin Magro,¹ Mathieu Girard,¹ Sébastien Garcia,¹ and Alexei Ourjoumtsev¹, *

¹JEIP, UAR 3573 CNRS, Collège de France, PSL University,
11, place Marcelin Berthelot, 75231 Paris Cedex 05, France

In quantum technologies, qubits encoded as collective atomic or solid-state excitations present important practical advantages, such as strong directional coupling to light. Unfortunately, they are affected by inhomogeneities between the emitters, which make them decay into dark states. In most cases, this process is non-Markovian. Through a simple formalism, we unveil a regime where this decay is suppressed by a combination of driving and non-Markovianity. We experimentally demonstrate this “driving protection” using a Rydberg superatom, making its coherent Rabi oscillations last about 14 times longer than the free lifetime of the collective qubit.

# I. INTRODUCTION

Quantum systems are never fully closed. Describing their interactions with their classical environment in an accurate and tractable way is essential in various contexts, in particular in quantum technologies. The textbook scenario assumes that these interactions are constant, and that the resulting correlations between the system and the environment instantly vanish without affecting the latter. For a qubit carried by a two-level atom, these interactions lead to three effects: the decay of the excited state, its excitation by thermal photons, and pure dephasing. In this Born-Markov approximation, the coherence of the system always decays exponentially.

Qubits can also be encoded as collective excitations of many emitters, acting together like a quantum phased array antenna. The excited state of the qubit is then an exchange-symmetric Dicke state [1], with a single excitation delocalized over the whole ensemble. This encoding is the cornerstone of heralded quantum information processing with atomic ensembles [2], of quantum memories using rare-earth ion-doped crystals or atomic gases [3], and of recent quantum optics experiments using Rydberg superatoms [4–6]. Among other practical advantages, it offers a directional coupling to light enhanced with respect to that of a single emitter [7]. However, collective encoding suffers from an inherent drawback: crystalline defects, inhomogeneous external fields, or thermal atomic motion randomly shift the transition frequencies of the emitters. With time, the phases between the elements of the “antenna” become random, scrambling its emission. In other words, inhomogeneous dephasing makes the qubit’s excited state decay in a quasi-continuum of dark, asymmetric Dicke states.

The Born-Markov approximation generally fails to describe this effect. To see this, let us consider $N \gg 1$ two-level emitters with angular frequencies $\omega_n$, carrying a qubit encoded in the states $|G\rangle = \prod_{n=1}^N |g^{(n)}\rangle$ where all the emitters are in the ground state $g$, and $|E_0\rangle = \sum_{n=1}^N \hat{\sigma}_{eg}^{(n)} |G\rangle / \sqrt{N}$ where $\hat{\sigma}_{eg}^{(n)} = |e^{(n)}\rangle\langle g^{(n)}|$ promotes the emitter $n$ into the excited state $e$, such that the

excitation is fully delocalized. After a time $t$, each term in this sum accumulates a phase factor $e^{-i\omega_n t}$. If the frequencies $\omega_n$ follow a probability distribution $P(\omega)$, an ergodic approximation shows that the overlap between the initial state $|E_0\rangle$ and the final state $|E(t)\rangle$ obeys

$$s(t) = \overline{\langle E_0 | E(t) \rangle} = \frac{1}{N} \sum_{n=1}^{N} \overline{e^{-i\omega_n t}} = \int_{-\infty}^{+\infty} d\omega P(\omega) e^{-i\omega t}. \quad (1)$$

Unless $P(\omega)$ is a Lorentzian, the decay of $s(t)$ is not exponential: the quasi-continuum of dark states has a memory, allowing some information to flow back into the system.

Equation 1 is sometimes used as-is [8, 9]. However, it assumes that the system is excited at a well-known moment, quasi-instantaneously with respect to the dephasing. While reasonable for heralded collective atomic spin qubits [2], this assumption fails for temporally multimode quantum memories [10], where the preparation time is undefined, and for Rydberg superatoms, where the dephasing is currently not negligible at the timescales of single-qubit gates. Moreover, as shown below, this decay can dramatically change in presence of driving. One is then left with two unappealing choices: keeping $N-1$ irrelevant states in the density matrix, or using a generalized master equation [11] to describe the evolution. Both are resource-intensive, and obscure much of the relevant physics.

Here, we show that the non-Markovianity of this decay process can not only be efficiently modeled, but also used to suppress it with remarkably little experimental effort. By mapping the inhomogeneous dephasing onto a displacement in time-frequency phase space, we construct a ladder of sequentially coupled Dicke states (Fig. 1(a),(b),(c)) and show that they can be determined analytically in many physically-relevant cases. Including the first few of them in the system’s Hilbert space and treating the others as a broad continuum provides a precise and resource-efficient description of the non-Markovian dynamics. In an analytically-tractable limit, we identify a process akin to cavity protection [12–15] where the qubit’s decay is inhibited by a combination of non-Markovianity and strong driving (Fig. 1(d-g)). We experimentally demonstrate this effect by using a Ryd-

<footnote>* Corresponding author: alexei.ourjoumtsev@college-de-france.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. Modeling inhomogeneous dephasing and driving protection. The effect of an inhomogeneous dephasing spectrum $P(\omega)$ (a) is mapped to a displacement in a phase space where frequency and time play the roles of the usual position and momentum quadratures $x$ and $p$ (b). This allows one to derive the corresponding ladder of $n+1$ dephasing-coupled Dicke states $|E_{k \le n}\rangle$ and their coupling coefficients (c), included in the Hamiltonian to describe the non-Markovian features, and to determine the irreversible decay of $|E_n\rangle$ into the “continuum” $C = \{|E_{k>n}\rangle\}$. In the Markovian case, $|E_0\rangle$ decays directly in the continuum (d). When driving the qubit transition with a Rabi frequency $\Omega$ (e), the dressed eigenstates $|\psi_\pm\rangle = (|G\rangle \pm |E_0\rangle)/\sqrt{2}$ also decay in the continuum, owing to its broad spectral extent. In a non-Markovian case, the decay out of the qubit’s subspace decreases drastically between the free dephasing from $|E_0\rangle$ (f) and the driven case (g), due to the energy difference between $|\psi_\pm\rangle$ and the $|E_{k \ge 1}\rangle$ asymmetric Dicke states.</figcaption></figure>

berg superatom undergoing two spectrally different non-Markovian processes. Consistently with the predictions of our model, we observe that its coherent Rabi oscillations can last up to ~14 times longer than the characteristic inhomogeneous dephasing time.

## II. THEORETICAL MODEL

Using the definitions above, let us consider the effect of the inhomogeneous dephasing Hamiltonian $\hat{H} = \sum_{n=1}^{N} \omega_n \hat{\sigma}_{ee}^{(n)}$ on a single excitation in the ensemble ($\hbar = 1$ in the following). We rewrite $\omega_n = \sqrt{2\omega_S} x_n$ where $x_n$ is a dimensionless random variable governed by a probability distribution $\rho(x) = \sqrt{2\omega_S}P(\sqrt{2\omega_S}x)$. The scaling rate $\omega_S$ leaves some freedom for optimizing the latter. Assuming ergodicity with $N \gg 1$, we approximate the sum over $n$ by $\int_{-\infty}^{+\infty} dx$ written as $\int dx$ to simplify notations, and replace the discrete operators $\hat{\sigma}_{ij}^{(n)}$ with $\hat{\sigma}_{ij}(x)$ satisfying $[\hat{\sigma}_{ij}(x), \hat{\sigma}_{kl}(x')] = \delta(x - x')(\delta_{jk}\hat{\sigma}_{il}(x) - \delta_{il}\hat{\sigma}_{kj}(x))$. We then rewrite $\hat{H} = \sqrt{2\omega_S} \int dx x \hat{\sigma}_{ee}(x)$ and $|E_0\rangle = \int dx \psi_0(x) \hat{\sigma}_{eg}(x) |G\rangle$ where the wavefunction $\psi_0(x) = \sqrt{\rho(x)}e^{i\phi(x)}$ satisfies $|\psi_0(x)|^2 = \rho(x)$, the arbitrary phase $\phi(x)$ leaving the choice of the most convenient gauge.

In the single excitation regime considered here, the di-

mensionless frequency $x$ and its conjugate dimensionless time $p$ can be treated as operators equivalent to position and momentum, with $[\hat{x}, \hat{p}] = i$ [16]. Then $\hat{\sigma}_{eg}(x) |G\rangle = |x\rangle$ is formally equivalent to a usual continuous-variable position eigenstate, giving $|E_0\rangle = \int dx \psi_0(x) |x\rangle$ and $\hat{H} = \sqrt{2\omega_S} \int dx x |x\rangle \langle x| = \sqrt{2\omega_S} \hat{x}$, which allows us to interpret the evolution $\hat{U}(t) = \exp(-i\hat{H}t)$ as a phase-space displacement in the direction $-p$ by a distance $\sqrt{2\omega_S} t$. Eq. 1 is recovered via

$$s(t) = \langle E_0 | \hat{U}(t) | E_0 \rangle = \int dx \rho(x) e^{-i\sqrt{2\omega_S} x t}. \quad (2)$$

We consider three physical examples:

* A Markovian process, with a Lorentzian spectrum $\rho_M(x) = 2/(\pi(1+4x^2))$ and a population decay rate $\sqrt{2}\omega_M$, where we can choose $\psi_{M,0}(x) = \sqrt{2/\pi}/(1+i2x)$. As expected, Eq. 2 leads to an exponential decay of the amplitude, $s_M(t) = e^{-\omega_M t}/\sqrt{2}$.

* Gaussian inhomogeneous dephasing, with $\rho_G(x) = e^{-x^2}/\sqrt{\pi}$ and $\psi_{G,0}(x) = \sqrt{\rho_G(x)}$, encountered for example in atomic ensembles where the velocities $v_n$ obey Maxwell-Boltzmann statistics. The Doppler shifts $\omega_n = k v_n$ then follow a Gaussian distribution where the width $\omega_G = k\sqrt{k_B T}/m$ depends on the effective spin-wave vector $k$, the Boltzmann constant $k_B$, the temperature $T$ and the atomic mass $m$. The state $|E_0\rangle$, which corresponds to a vacuum state $|0\rangle$ in this phase space, is displaced by $\hat{U}(t)$ into a coherent state $|-i\omega_G t\rangle$, with a Gaussian decay of the amplitude $s_G(t) = \langle 0|-i\omega_G t\rangle = e^{-\omega_G^2 t^2/2}$.

* Differential light shifts between ground and excited atomic states in a three-dimensional harmonic trap, where the Gaussian thermal distribution of atomic positions and the quadratic position-dependence of the light-shift potential lead to $\rho_L(x) = \sqrt{x}e^{-x}\theta(x)2/\sqrt{\pi}$ and $\psi_{L,0}(x) = \sqrt{\rho_L(x)}$, where $\theta$ is the Heaviside step function. The scaling rate is $\omega_L = \delta\alpha_r k_B T/\sqrt{2\hbar}$ where $\delta\alpha_r = (\alpha_e - \alpha_g)/\alpha_g$ is the relative difference of polarizabilities of the excited and ground states. The amplitude $s$ decays here as $s_L(t) = (1 + i\sqrt{2\omega_L t})^{-3/2}$.

Unless stated otherwise, we set $\omega_S = 1$ as a unit for angular frequencies and rates, and $1/\omega_S$ as a unit for time.

The interest of this mapping becomes apparent when the qubit experiences classical or quantum driving and Eqs. 1 or 2 no longer hold. As shown below, it allows one to explicitly construct the ladder of Dicke states sequentially coupled by the inhomogeneous dephasing. By truncating this ladder at a given level and treating its coupling with higher-order states as an irreversible decay, one can use a standard Gorini - Kossakowski - Sudarshan
- Lindblad (GKSL) master equation [17, 18] to describe the non-Markovian dynamics via a small extension of the Hilbert space.

We start by constructing the ladder of Dicke states $\{|E_k\rangle\}$ from $|E_0\rangle$ by successive applications of $\hat{H}$. As $\hat{H} \propto \hat{\lambda}$, the corresponding wavefunctions have the form $\psi_k(x) = Q_k(x)\psi_0(x)$ where $Q_k$ is a polynomial of degree $k$. The problem then reduces to finding a family of orthogonal polynomials for the weight function $\rho: \langle E_k | E_l \rangle = \int \rho(x) Q_k^*(x) Q_l(x) = \delta_{k,l}$. This problem is generically solved by Gram-Schmidt orthonormalization but, given its importance in mathematics, engineering or data science, explicit solutions are known in many cases [19]. By construction, $\hat{H}$ is tridiagonal in this basis, with nonzero elements $H_k = \langle E_k | \hat{H} | E_k \rangle$ and $V_k = \langle E_{k-1} | \hat{H} | E_k \rangle$. For the Gaussian $\rho_G$, $Q_k$ is a Hermite polynomial up to a normalization factor and $|E_{G,k}\rangle$ is the Fock state $|k\rangle$, which directly gives $H_{G,k} = 0$ and $V_{G,k} = \sqrt{k}$. For the light shift $\rho_L$, $Q_k(x) = \sqrt{2^k k!/(2k+1)!!} L_k^{(1/2)}(x)$ where $L_k^{(1/2)}$ is a generalized Laguerre polynomial. Their recurrence relation [19] gives $H_{L,k} = \sqrt{2}(2k+3/2)$ and $V_{L,k} = -\sqrt{k}(2k+1)$. For a generic $\rho(x)$ we initialize the recurrence with $\psi_{-1}(x) = 0$ and $\psi_0(x)$, calculate $a_k \psi_{k+1}(x) = (\sqrt{2}x - H_k)\psi_k(x) - V_k \psi_{k-1}(x)$, and determine $a_k$ by normalization.

If $\rho(x) \sim x^{-d}$ at infinity, the ladder of Dicke states terminates at step $n = \lfloor(d-1)/2\rfloor$ as $x \psi_n(x)$ becomes non-normalizable: the effect of $\hat{H}$ on $|E_n\rangle$ is then treated in a Born-Markov approximation as an irreversible decay into a spectrally-broad continuum. For the Lorentzian $\rho_M$ this happens at $n = 0$. If this “quantum Markov order” [20] is impractically large or infinite, as it is the case for $\rho_G$ and $\rho_L$, one can choose $n$ as a compromise between the precision and the size of the simulations and approximate the coupling of $|E_n\rangle$ to higher-lying states by a Markovian process.

To terminate this ladder, we thus need to calculate the amplitude decay rate $\gamma_n$ and the Lamb-like energy shift $\delta_n$ due to the coupling of $|E_n\rangle$ with the “broad continuum” of states $\{|E_{k>n}\rangle\}$. A straightforward use of the resolvent $\hat{G}(z) = (z-\hat{H})^{-1}$ [21], detailed in Appendix A, shows that $\Delta_n = \delta_n - i\gamma_n$ obeys the recurrence relation

$$ \Delta_n = H_0 - H_n - \frac{|V_n|^2}{\Delta_{n-1}}, \quad (3) $$

the initial value $\Delta_0$ being given by $\Delta_0 = -1/G_0(H_0)$ where, in the limit $\text{Im}(z) \to 0^+$,

$$ G_0(z) = \langle E_0 | \hat{G}(z) | E_0 \rangle = \int \frac{\rho(x)dx}{z - \sqrt{2}x} = -i \int_0^{+\infty} e^{itz} s(t)dt \quad (4) $$

For the Markovian case $\rho_M$, $G_{M,0}(z) = 1/(z + i/\sqrt{2})$ and $H_{M,0} = 0$ gives, self-consistently, $\delta_{M,0} = 0$ and $\gamma_{M,0} = 1/\sqrt{2}$. For the Gaussian $\rho_G$ with $H_{G,0} = 0$ and $G_{G,0}(z) = -i\sqrt{\pi/2}w(z/\sqrt{2})$ where $w(z) =$

<figure><img src="image_3.png" /><figcaption>FIG. 2. Dicke basis truncation. (a) and (b) Energy levels $H_k$ of truncated basis states with $n = 3$ for the Gaussian and lighshift-broadened cases, respectively. Arrow linewidths are proportional to coupling coefficients $|V_k|$ or decay rates $\gamma_3$ towards the continuum $|C\rangle$. The dashed line for $E_3$ represent the energy $H_3$ before the addition of the Lamb-like energy shift $\delta_3$. (c) and (d) Related evolutions of the population $|s(t)|^2$. The black curves are the analytic solutions given by Eq. 2. The blue dotted and dashed lines are computed in a truncated basis with $n = 1$, without and with the continuum state, respectively. The green curves use the basis with $n = 3$, again without or with the continuum.</figcaption></figure>

$e^{-z^2}$ erfc($-iz$) is the Faddeeva function, $\delta_{G,n} = 0$ and $\gamma_{G,n} = (2/\pi)^{(-1)^n/2} n!! / (n-1)!!$. For the light-shift spectrum $\rho_L$, $G_{L,0}(z) = -\sqrt{2} - i\sqrt{\pi z} \sqrt{2} w(\sqrt{z}/\sqrt{2})$ and $H_{L,0} = 3/\sqrt{2}$; Eq. 3 becomes $\Delta_{L,n} = -2^{3/2}n - n(2n+1)/\Delta_{L,n-1}$ and gives $\delta_{L,n}$ and $\gamma_{L,n}$ by induction.

As shown in Fig. 2, numerical simulations in a basis of size $n+1$ rapidly converge to the analytic solutions of Eq. 2 with increasing $n$. At long times, approximating the higher-lying states by a Markovian continuum significantly improves the precision compared to a simple phase-space truncation. This approach enables resource-efficient numerical simulations in more complex situations where several collective qubits are coupled to time-dependent classical or quantum light.

### III. DRIVING PROTECTION EFFECT

This formalism predicts an interesting effect when the qubit is resonantly driven with a constant Rabi angular frequency $\Omega$. A similar resolvent-based approach, detailed in Appendix B together with the derivations of the expressions below, shows that Eq. 2, describing the excited-state amplitude dynamics, is then replaced with

$$ s_{\Omega}(t) = \int_{-\infty}^{+\infty} \frac{i dz \, z G_0(z + H_0) e^{-itz}}{2\pi \, z - \frac{\Omega^2}{4} G_0(z + H_0)}. \quad (5) $$

In the following we focus on the strong-driving regime where $\Omega \gg 1$.
For the Lorentzian spectrum $\rho_M$, $H_{M,0} = 0$ and the simple form of $G_{M,0}(z)$ gives the exact solution

$$s_{M,\Omega}(t) = e^{-\gamma_M t/2} \left( \cos(\Omega_M t/2) - \frac{\gamma_M}{\Omega_M} \sin(\Omega_M t/2) \right) \quad (6)$$

with $\Omega_M = (\Omega^2 - \gamma_M^2)^{1/2}$. In this Markovian case, a conventional master equation approach leads to the same result. For $\Omega \ge \gamma_M = 1/\sqrt{2}$ the decay rate $\gamma_M/2$, reduced by half by the driving, becomes independent on $\Omega$: no matter the driving strength, the decay persists.

The situation radically changes for the Gaussian spectrum $\rho_G$. The asymptotic expansion $\text{Re}(G_{G,0}(z)) = z^{-1} + z^{-3} + O(z^{-5})$ gives

$$s_{G,\Omega}(t) \approx e^{-\gamma_{G,\Omega}t} \cos(\Omega_G t/2), \\ \Omega_G \approx \Omega + \frac{2}{\Omega}, \qquad \gamma_{G,\Omega} \approx \frac{\Omega^2}{8} \sqrt{\frac{\pi}{2e}} e^{-\Omega^2/8}. \quad (7)$$

The decay profile changes from Gaussian to exponential, its rate $\gamma_{G,\Omega}$ vanishing exponentially with $\Omega^2$. Here, strong driving combined with non-Markovianity allows high-visibility Rabi oscillations to occur over times strongly exceeding the natural dephasing time, as shown in Fig. 3(a). This “driving protection” of the collective qubit is closely related to cavity protection [12–15], the collective Rabi frequency $\Omega$ replacing the collective cavity coupling.

For the light-shift spectrum $\rho_L$, the expansion $\text{Re}(G_{L,0}(z + H_{L,0})) = z^{-1} + 3z^{-3} + O(z^{-4})$ leads to

$$s_{L,\Omega}(t) \approx \frac{1}{2} \left( e^{-\gamma_{L,\Omega}t} e^{-i\Omega_L t/2} + e^{i\Omega_L t/2} \right), \\ \Omega_L \approx \Omega + \frac{6}{\Omega}, \qquad \gamma_{L,\Omega} \approx \frac{\Omega^2}{8} \sqrt{\frac{\pi\Omega}{2e^3}} e^{-\Omega/\sqrt{8}}. \quad (8)$$

As visible in Fig. 3(b), the driving protection persists but, since $\gamma_{L,\Omega}$ decreases slower than $\gamma_{G,\Omega}$ with $\Omega$, it is weaker than for a Gaussian spectrum. While in the Gaussian case the system fully decays in the continuum, here the population $|s_{L,\Omega}(t)|^2$ tends to 1/4 at long times. In fact, the two terms in the expression of $s_{L,\Omega}(t)$ respectively correspond to the propagators $\langle \psi_{\pm} | \hat{U}_{\Omega}(t) | \psi_{\pm} \rangle$ with $|\psi_{\pm}\rangle = (|G\rangle \pm |E_0\rangle)/\sqrt{2}$: the population in $|\psi_{\pm}\rangle$ decays with a rate $2\gamma_{L,\Omega}$ while $|\psi_{-}\rangle$ is very efficiently protected by the drive.

Figure 3 shows that the approximate equations 7 and 8 correctly reproduce the dynamics computed in sufficiently large Dicke bases. For the light-shift spectrum, the numerical integration confirms that $|\psi_{-}\rangle$ is efficiently protected. Poles of the integrand in Eq. 5 neglected in the approximations above account for transient dynamics at $0 < t \lesssim 1$.

## IV. PHYSICAL INTERPRETATION

While this “driving protection” qualifies as dynamical decoupling, it is a fundamentally collective effect, dif-

<figure><img src="image_4.png" /><figcaption>FIG. 3. Driving protection. Evolution of the population $|s_{\Omega}(t)|^2$ (a) for the Gaussian spectrum $\rho_G$ and (b) for the lightshift-broadened spectrum $\rho_L$, with and without driving. The black curves are the analytic solutions of Eq. 2 for an undriven qubit ($\Omega = 0$). The blue lines are computed in a Dicke basis with a Rabi frequency $\Omega = 10$. The dotted red curves are the corresponding approximations given by Eqs. 7 and 8. In (b), the dotted and dashed green curves are the populations of $|\psi_{+}\rangle$ and $|\psi_{-}\rangle$ states, respectively.</figcaption></figure>

ferent from the techniques generally used to counteract dephasing in two-level systems. The simplest of them is spin echo [22], where the dephasing due to an unknown shift $\omega$ of the resonance frequency can be compensated by periodically swapping the ground and excited amplitudes via $\pi$ pulses. This mitigation persists when the qubit is driven continuously with a Rabi frequency $\Omega_0 \gg \omega$ [23].

The process considered here is not a mere sum of such effects acting on every emitter in the ensemble. Indeed, the Rabi frequency $\Omega = \Omega_0\sqrt{N}$ is collectively enhanced with respect to the individual Rabi frequency $\Omega_0$. The latter can be much smaller than $\omega_S$, and thus insufficient to mitigate dephasing in a single emitter. Moreover, discrete or continuous spin echo fundamentally relies on swapping the roles of the ground and the excited states. Here the excited population of each emitter remains below $1/N$. For Rydberg superatoms, single-atom spin echo requires bringing $\sim N-1$ atoms in a Rydberg state which is prevented by dipole blockade. Finally, unlike discrete or continuous spin echo, driving protection crucially relies on non-Markovianity.

The mechanism behind driving protection can be sketched in the following way. In the rotating frame, the eigenstates of the driven collective qubit are $|\psi_{\pm}\rangle = (|G\rangle \pm |E_0\rangle)/\sqrt{2}$, with eigenenergies $\pm\Omega/2$. The dephasing $\sqrt{2}\hat{x}$ couples them to the continuum of states $\{|x\rangle\}$ via $g(x) = \langle x|\sqrt{2}\hat{x}|\psi_{\pm}\rangle = \pm x \psi_0(x + H_0/\sqrt{2})$, the offset $H_0$ accounting for the resonant driving. For a Lorentzian spectrum, $g_M(x) \propto x/(1+i2x)$ is spectrally broad [Fig. 4]: for any $\Omega$, the states $|\psi_{\pm}\rangle$ are resonantly coupled to the continuum, which results in irreversible decay. In contrast, for a Gaussian spectrum, $g_G(x) \propto x e^{-x^2/2}$ is spectrally narrow and localized around 0. Strong driving makes the dephasing-induced coupling off-resonant and thus inhibits the decay of the qubit. For the light-shift spectrum, $g_L(x) \propto x \psi_{L,0}(x + 3/2)$ is asymmetric, and $|\psi_{-}\rangle$ is more strongly protected than $|\psi_{+}\rangle$.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Explanation of the driving protection mechanism. The eigenstates $|\psi_{\pm}\rangle$ of the driven qubit, shifted by $\pm\Omega/2$, remain resonantly coupled to a “continuum” $\{|x\rangle\}$ for a Lorentzian spectrum (blue, left). This coupling becomes off-resonant for a Gaussian spectrum (red, center), which decreases the decay rate of the Rabi oscillations. For a light-shift spectrum, $|\psi_+\rangle$ decays faster than $|\psi_-\rangle$.</figcaption></figure>

V. EXPERIMENTAL METHODS

To test our predictions experimentally, we use a Rydberg superatom: a small cold atomic ensemble where interactions between Rydberg-state atoms block the number of excitations to one at most, thus defining a collectively-encoded qubit [4–6]. As shown in Fig. 5, the superatom is coupled to a running-wave optical cavity enabling an efficient mapping of its excited state $|E_0\rangle$ to a single-photon Fock state $|1\rangle$. Therefore, by driving the superatom to its excited state, waiting for a storage duration $t_s$, mapping $|E_0\rangle$ to a photon, and measuring the probability to detect the latter with a photon counter, we obtain a signal proportional to $|s(t_s)|^2$.

The superatom is produced in 0.1 s through the sequential use of two-dimensional (2D) magneto-optical trapping (MOT), 3D MOT, optical molasses cooling, transporting in a one-dimensional optical lattice, and degenerate Raman sideband cooling, as described in detail in Refs. [24, 25]. It is made of $N \simeq 500$ atoms with a Gaussian spatial distribution, loaded in a crossed dipole trap. In the last initialization stage, the atoms are optically pumped in the state $|g\rangle : 5S_{1/2}, F=1, m_F=1$. By adjusting the depth of the trapping potential during either the preparation or the experimental run, we can control the light-shifts in our system, as well as increase the compression and modify the temperature of the atomic cloud. This allows us to prepare inhomogeneous broadening spectra with similar characteristic widths but with different shapes. To obtain a purely Gaussian Doppler broadening, we switch off the trap during the experimental protocol and increase the confinement potential during the preparation, creating a cloud with a root-mean-square radius $\sigma_r = 3.7(2)$ µm and a temperature $T = 6.2(3)$ µK. To probe a combination of light-shift and Doppler broadening mechanisms, we do not compress the cloud, which makes it colder ($\sigma_r = 4.4(2)$ µm and $T = 2.5(1)$ µK), but apply an average light-shift of $4.8(6)10^2$ kHz during the experimental protocol. We adjust the temperature and the light shift to achieve similar

<figure><img src="image_6.png" /><figcaption>FIG. 5. Experimental protocol. (a) Schematic of the experimental setup featuring a Rydberg superatom in an 4-mirror optical cavity, with the laser pulses and the transition scheme used in the experiment, described in Section V. (b) Measured pulse shapes (dots) and best-fit functions (lines) used in the simulations. The lower panel is a zoom on the times between the dotted vertical lines in the upper panel. The blue and dark red curves represent the blue and red drive pulses, respectively. In the lower panel, we represent the short red drive pulse used to prepare the cloud in the singly-excited symmetric Dicke state $|E_0\rangle$ with a fast $\pi$ pulse. The cyan curve represents the blue read beam mapping $|E_0\rangle$ to a photon, while the red line indicates the photon's temporal intensity profile seen by the single-photon detector. The gray dashed lines represent the average times for the $\pi$ pulse and the detected photon.</figcaption></figure>

coherence times for both configurations. Then, we repeat 10 times the excitation and measurement protocol described below. To remove any residual Rydberg-excited atoms, we apply an electric field pulse using in-vacuum electrodes before each repetition.

In the first “drive” step of our protocol, we use a “red drive” laser at 795 nm and “blue drive” laser at 475 nm to excite a resonant two-photon transition to the Rydberg state $|e\rangle : 109S, J=1/2, m_J=1/2$, with a detuning $\Delta = 2\pi \times -500$ MHz from the intermediate state $|f\rangle : 5P_{1/2}, F=2, m_F=2$ [see Fig. 5]. The $\sim 20$ µm Rydberg blockade radius [26] largely exceeds the 4.4 µm Gaussian radius of the cloud, strongly inhibiting multiple Rydberg excitations. Thus, the cloud behaves as an effective two-level system with a ground superatom state $|G\rangle$ and a collective singly-excited Rydberg state
$|E_0\rangle$. We adjust the collective two-photon Rabi frequency $\Omega = \sqrt{N}\Omega_r\Omega_b/(2\Delta)$ via the Rabi frequency $\Omega_r$ of the “red drive” beam, keeping constant the Rabi frequency $\Omega_b = 2\pi \times 7.8(4)$ MHz of the “blue drive”. The drive duration $t_d$ is set by the “red drive” pulse duration, encompassed by the longer “blue drive” pulse, as shown in Fig. 5(b).

During the subsequent “read” step, we quasi-adiabatically turn on a “blue read” laser resonant with the $|e\rangle-|f\rangle$ transition at 475 nm. As the collective coupling strength $g = 2\pi \times 8.5(3)$ MHz between the $|f\rangle-|g\rangle$ transition and the resonant cavity mode exceeds the decay rates of both the cavity field ($2\pi \times 2.9$ MHz) and the atomic dipole ($2\pi \times 3$ MHz), this efficiently transfers the excitation from $|E_0\rangle$ to an intra-cavity photon [27]. This photon then leaks out of the cavity mostly through the output-coupler mirror and is directed to a single-photon detector.

To measure the free inhomogeneous dephasing of $|E_0\rangle$, we prepare it by a fast $\pi$ pulse with a strong drive $\Omega/(2\pi) = 1.0(1)$ MHz. Then, we map $|E_0\rangle$ to a photon after a variable storage duration $t_s$, defined as the difference between the average intensity-weighted times of the $\pi$ pulse and of the extracted photon. To investigate the drive protection, we vary the duration $t_d$ of the drive pulse between 0.1 and 11.9 µs and map $|E_0\rangle$ to a photon directly afterwards, revealing point by point the Rabi oscillations.

## VI. EXPERIMENTAL RESULTS

The results of the experiment are presented in Fig. 6. Panels (a) and (b) show the decay of $|E_0\rangle$ without driving, fitted with the solutions of Eq.2. The Gaussian decay in panel (a), with a mean lifetime $\sqrt{\pi}/(2\omega_G) = 2.33(7)$ µs, results from pure Doppler thermal broadening. Panel (b) shows a decay with a similar mean lifetime of 2.03(9) µs, but with a different shape, arising from a convolution between a Gaussian, due to residual thermal motion, and a light-shift spectrum slightly modified by the anharmonicity of our optical dipole trap [details in Appendix C].

Panels (c) and (d) correspond to the same configurations in presence of Rabi flopping. They confirm the predicted driving protection: the visibility of Rabi oscillations increases with their frequency and reaches coherence times strongly exceeding the free decay lifetimes measured above. The observed dynamics are reproduced by numerical simulations on a desktop computer, using a truncated basis of Dicke states [details in Appendix C]. By fitting an exponential envelope $e^{-\Gamma t_d}$ on the measured oscillations, with a fit function accounting for shot-to-shot fluctuations of the number of atoms in the prepared superatom [details in Appendix C 1], we obtain the decay rate $\Gamma$ presented in panel (e). It decreases rapidly with the Rabi frequency, down to a factor 14(5) below the characteristic free decay rate of $|E_0\rangle$. As expected

<figure><img src="image_7.png" /><figcaption>FIG. 6. Experimental observation of driving protection in a Rydberg superatom. The detected photon population is proportional to the population $|s(t)|^2$ of the qubit's excited state $|E_0\rangle$. We realize two dephasing spectra: purely Gaussian, due to thermal Doppler broadening (left (a)(c) and red to yellow colored circles) and a convolution of a light-shift in a slightly anharmonic trap with a Gaussian due to a residual Doppler broadening (right (b)(d) and green to blue colored squares). Error bars represent the standard error. (a) and (b): Free evolution of the excited state $|E_0\rangle$, fitted with the corresponding rescaled $|s(t)|^2$ functions, yielding mean lifetimes of 2.33(7) µs and 2.03(9) µs, respectively. (c) and (d): Evolution with Rabi frequencies $\Omega = 2\pi \times (239, 406, 244, 414)$ kHz for (red, orange, green, blue), respectively. Points connected by lines show the results of corresponding numerical simulations using our dephasing model, parametrized with independently measured values. (e): Fitted exponential decay rate $\Gamma$ of the observed oscillation amplitude [fit function of Eq. C2 in Appendix C 1] as a function of the Rabi frequency $\Omega/(2\pi)$. The horizontal dashed red and green lines show inverses of the free-decay lifetimes given by fits in (a) and (b), respectively. The dashed orange line is the expression of $\gamma_{G,\Omega}$ given by Eq. 7. The dotted purple line is a best-fitted $\Gamma_{\Omega} \Omega^2 + \Gamma_0$ function for points where $\Omega > 2\pi \times 0.7$ MHz.</figcaption></figure>

from Sec. III, at intermediate values $5 < \Omega/\omega_G < 10$ (i.e $0.3$ MHz $< \Omega/(2\pi) < 0.6$ MHz), the driving protection is less efficient in the light-shifted case than in the Gaussian one. For the latter, the data points are close to the analytic result of Eq. 7, valid in the strong driving limit and plotted for $\Omega > 3\omega_G$. At very strong driving, $\Gamma$ is ultimately limited by experimental effects compatible with a
form $\Gamma_{\Omega} \Omega^2 + \Gamma_0$ [28]. $\Gamma_{\Omega}$ accounts for decoherence mechanisms proportional to the power of the driving beam, including spontaneous emission of the intermediate level of the two-photon transition between $|G\rangle$ and $|E_0\rangle$ and imperfect Rydberg blockade; while $\Gamma_0$ stems from defects such as the finite lifetime of the Rydberg state, imperfect optical pumping, impure laser polarizations, or residual electric fields.

## CONCLUSION

We unveiled, theoretically and experimentally, a dynamical decoupling effect for collective qubits, which enhances coherence times by a factor 14 in our experimental system. The demonstration of this driving protection effect was enabled by the simplicity of the underlying theoretical framework. Among existing descriptions of non-Markovian quantum dynamics [29, 30], it has some resemblance with the formalism of pseudomodes [31]. Applicable to a large class of physical systems, it gives access to an easy implementation, to efficient numerics, and to intuitive physical interpretations exemplified here.

It is worth comparing this protection mechanism with other ways of counteracting inhomogeneous dephasing, all of which come at some practical cost. Firstly, one can act on the dephasing's causes. For the Doppler and the light-shift broadenings, it means cooling the atoms further, which significantly degrades experimental duty cycles, and switching off the trapping light, which eventually leads to atomic loss. Alternatively, the system can be made less sensitive to dephasing, by using Doppler-free excitation or state-insensitive dipole traps. For Rydberg superatoms, such traps operate close to atomic transitions [32], while Doppler-free excitation requires three laser beams with specific angles [33], putting significant constrains on optical access. Other mitigation schemes for Doppler broadening in Rydberg gases involve a continuous off-resonant coupling [28] or a transient transfer [34] to an auxiliary state. Compared to these schemes,

the strengths and weaknesses of driving protection are similar to those of spin echo, applicable to collectively-encoded qubits [35] but, as mentioned earlier, incompatible with strong Rydberg blockade. It leads to unwanted phase accumulations, sensitive to fluctuations of the Rabi frequency, but it is not restricted to specific ranges of atomic levels and it requires no resources other than those already available to manipulate the qubit.

Besides its potential applications in quantum information processing with collective qubits, the driving protection process is interesting in its own right. Compared to cavity protection, where the collective coupling between the cavity and the emitters is generally fixed [14], the temporal control of the collective Rabi frequency gives access to dynamics yet to be explored, such as the active stabilization of the eigenstates of the coupled system.

## FUNDING

This work was funded by the ERC Starting Grant 677470 SEAQUEL, the Plan France 2030 through the project ANR-23-PETQ-0013, the Ile-de-France DIM network SIRTEQ through the project CRIMP, and the CI-FAR Azrieli Global Scholars program.

## ACKNOWLEDGEMENTS

We thank J. Vaneeccloo for his help during early experimental developments and characterizations.

## DISCLOSURES

The authors declare no conflicts of interest.

## DATA AVAILABILITY

Data underlying the results presented in this paper are available upon reasonable request from the authors.

[1] R. H. Dicke, Coherence in spontaneous radiation processes, Phys. Rev. **93**, 99 (1954).

[2] L.-M. Duan, M. Lukin, J. I. Cirac, and P. Zoller, Long-distance quantum communication with atomic ensembles and linear optics, Nature **414**, 413 (2001).

[3] A. I. Lvovsky, B. C. Sanders, and W. Tittel, Optical quantum memory, Nature Photonics **3**, 706 (2009).

[4] L. H. Pedersen and K. Mølmer, Few qubit atom-light interfaces with collective encoding, Phys. Rev. A **79**, 012320 (2009).

[5] J. Kumlin, C. Braun, C. Tresp, N. Stiesdal, S. Hofferberth, and A. Paris-Mandoki, Quantum optics with Rydberg superatoms, Journal of Physics Communications **7**, 052001 (2023).

[6] X.-Q. Shao, S.-L. Su, L. Li, R. Nath, J.-H. Wu, and W. Li, Rydberg superatoms: An artificial quantum system for quantum information processing and quantum optics, Applied Physics Reviews **11**, 031320 (2024).

[7] Y. O. Dudin, L. Li, F. Bariani, and A. Kuzmich, Observation of coherent many-body Rabi oscillations, Nature Physics **8**, 790 (2012).

[8] B. Zhao, Y.-A. Chen, X.-H. Bao, T. Strassel, C.-S. Chuu, X.-M. Jin, J. Schmiedmayer, Z.-S. Yuan, S. Chen, and J.-W. Pan, A millisecond quantum memory for scalable quantum networks, Nature Physics **5**, 95 (2009).

[9] S. Schmidt-Eberle, T. Stolz, G. Rempe, and S. Dürr, Dark-time decay of the retrieval efficiency of light stored as a Rydberg excitation in a noninteracting ultracold gas, Phys. Rev. A **101**, 013421 (2020).
[10] J. Nunn, K. Reim, K. C. Lee, V. O. Lorenz, B. J. Sussman, I. A. Walmsley, and D. Jaksch, Multimode memories in atomic ensembles, Phys. Rev. Lett. **101**, 260502 (2008).

[11] H.-P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems* (Oxford University Press, 2007).

[12] I. Diniz, S. Portolan, R. Ferreira, J. M. Gérard, P. Bertet, and A. Auffèves, Strongly coupling a cavity to inhomogeneous ensembles of emitters: Potential for long-lived solid-state quantum memories, Phys. Rev. A **84**, 063810 (2011).

[13] Z. Kurucz, J. H. Wesenberg, and K. Mølmer, Spectroscopic properties of inhomogeneously broadened spin ensembles in a cavity, Phys. Rev. A **83**, 053852 (2011).

[14] M. Baghdad, P.-A. Bourdel, S. Schwartz, F. Ferri, J. Reichel, and R. Long, Spectral engineering of cavity-protected polaritons in an atomic ensemble, Nature Physics **19**, 1104 (2023).

[15] N. Sauerwein, F. Orsi, P. Uhrich, S. Bandyopadhyay, F. Mattiotti, T. Cantat-Moltrecht, G. Pupillo, P. Hauke, and J.-P. Brantut, Engineering random spin models with atoms in a high-finesse cavity, Nature Physics **19**, 1128 (2023).

[16] N. Fabre, A. Keller, and P. Milman, Time and frequency as quantum continuous variables, Phys. Rev. A **105**, 052429 (2022).

[17] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of N-level systems, Journal of Mathematical Physics **17**, 821 (1976).

[18] G. Lindblad, On the generators of quantum dynamical semigroups, Communications in Mathematical Physics **48**, 119 (1976).

[19] G. Szegö, *Orthogonal Polynomials* (American Mathematical Society, Providence, RI, 1975).

[20] P. Taranto, F. A. Pollock, S. Milz, M. Tomamichel, and K. Modi, Quantum Markov order, Phys. Rev. Lett. **122**, 140401 (2019).

[21] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, Nonperturbative calculation of transition amplitudes, in *Atom-Photon Interactions* (John Wiley & Sons, Ltd, 1998) Chap. 3, pp. 165–255.

[22] E. L. Hahn, Spin echoes, Phys. Rev. **80**, 580 (1950).

[23] L. Viola and S. Lloyd, Dynamical suppression of decoherence in two-state quantum systems, Phys. Rev. A **58**, 2733 (1998).

[24] J. Vaneeecloo, S. Garcia, and A. Ourjoumtsev, Intracavity Rydberg superatom for optical quantum engineering: Coherent control, single-shot detection, and optical π phase shift, Phys. Rev. X **12**, 021034 (2022).

[25] V. Magro, J. Vaneeecloo, S. Garcia, and A. Ourjoumtsev, Deterministic freely propagating photonic qubits with negative Wigner functions, Nature Photonics **17**, 688 (2023).

[26] D. Tong, S. M. Farooqi, J. Stanojevic, S. Krishnan, Y. P. Zhang, R. Côté, E. E. Eyler, and P. L. Gould, Local blockade of Rydberg excitation in an ultracold gas, Phys. Rev. Lett. **93**, 063001 (2004).

[27] A. V. Gorshkov, A. André, M. D. Lukin, and A. S. Sørensen, Photon storage in Λ-type optically dense atomic media. i. cavity model, Phys. Rev. A **76**, 033804 (2007).

[28] R. Finkelstein, O. Lahad, I. Cohen, O. Davidson, S. Kiriati, E. Poem, and O. Firstenberg, Continuous protection of a collective state from inhomogeneous dephasing, Phys. Rev. X **11**, 011008 (2021).

[29] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Colloquium: Non-Markovian dynamics in open quantum systems, Rev. Mod. Phys. **88**, 021002 (2016).

[30] I. de Vega and D. Alonso, Dynamics of non-Markovian open quantum systems, Rev. Mod. Phys. **89**, 015001 (2017).

[31] L. Mazzola, S. Maniscalco, J. Piilo, K.-A. Suominen, and B. M. Garraway, Pseudomodes as an effective description of memory: Non-Markovian dynamics of two-state systems in structured reservoirs, Phys. Rev. A **80**, 012104 (2009).

[32] J. Lampen, H. Nguyen, L. Li, P. R. Berman, and A. Kuzmich, Long-lived coherence between ground and Rydberg levels in a magic-wavelength lattice, Phys. Rev. A **98**, 033411 (2018).

[33] I. I. Ryabtsev, I. I. Beterov, D. B. Tretyakov, V. M. Entin, and E. A. Yakshina, Doppler- and recoil-free laser excitation of Rydberg states via three-photon transitions, Phys. Rev. A **84**, 053409 (2011).

[34] Y. Jiao, C. Li, X.-F. Shi, J. Fan, J. Bai, S. Jia, J. Zhao, and C. S. Adams, Suppression of motional dephasing using state mapping, Phys. Rev. Lett. **134**, 053604 (2025).

[35] J. Rui, Y. Jiang, S.-J. Yang, B. Zhao, X.-H. Bao, and J.-W. Pan, Operating spin echo in the quantum regime for an atomic-ensemble quantum memory, Phys. Rev. Lett. **115**, 133002 (2015).

[36] J. Johansson, P. Nation, and F. Nori, Qutip 2: A Python framework for the dynamics of open quantum systems, Computer Physics Communications **184**, 1234 (2013).
Appendix A: Lamb-like shift and decay rate

In the model described in the manuscript, we identify a ladder of Dicke states {|E<sub>k</sub>⟩} sequentially coupled by the inhomogeneous-dephasing Hamiltonian ˆH = √2 ˆx, and calculate the diagonal terms H<sub>k</sub> = 〈E<sub>k</sub>| ˆH |E<sub>k</sub>⟩ and the coupling terms V<sub>k</sub> = 〈E<sub>k-1</sub>| ˆH |E<sub>k</sub>⟩. We then restrict the system’s Hilbert space to the finite subspace {|E<sub>k≤n</sub>⟩} formed by the lowest n + 1 states, and treat the orthogonal subspace spanned by {|E<sub>k>n</sub>⟩} as a broad continuum. Here we calculate the amplitude decay rate γ<sub>n</sub> and the Lamb-like energy shift δ<sub>n</sub> of |E<sub>n</sub>⟩ induced by the coupling to this continuum. For this, we define the projectors ˆP = ∑<sub>k=0</sub><sup>n</sup> |E<sub>k</sub>⟩⟨E<sub>k</sub>| and ˆQ = 1 − ˆP into the two subspaces of interest, and derive an effective non-Hermitian Hamiltonian ˆH<sub>e</sub> from the resolvent ˆG(z) = (z − ˆH)<sup>-1</sup> which obeys [21]

$$
\hat{P}\hat{G}(z)\hat{P} = \frac{1}{z\hat{P} - \hat{H}_e(z)}, \qquad \hat{H}_e(z) = \hat{P}\hat{H}\hat{P} + \hat{P}\hat{H}\hat{Q}\frac{1}{z\hat{Q} - \hat{Q}\hat{H}\hat{Q}}\hat{Q}\hat{H}\hat{P}. \qquad (A1)
$$

To obtain a parameter-free ˆH<sub>e</sub>, the energy z is fixed to a physically-relevant value. As |E<sub>0</sub>⟩ is the only experimentally-
accessible state, we set z to H<sub>0</sub>. Since ˆH is tridiagonal in this basis,

$$
\hat{H}_e(z) = \hat{P}\hat{H}\hat{P} + \Delta_n |E_n\rangle\langle E_n|, \qquad \Delta_n = \delta_n - i\gamma_n = |V_{n+1}|^2 \langle E_{n+1} | \hat{G}_{n+1}(z) | E_{n+1} \rangle \quad (\text{A2})
$$

where $\hat{G}_{n+1}(z) = (z\hat{Q} - \hat{Q}\hat{H}\hat{Q})^{-1}$ is the resolvent of $\hat{Q}\hat{H}\hat{Q}$. $\Delta_n$ can be calculated inductively by replacing $G$ with $\hat{G}_n$ and $\hat{P}$ with $\hat{P}_n = |E_n\rangle\langle E_n|$ in Eq. A1:

$$
\hat{P}_n \hat{G}_n(z) \hat{P}_n = \hat{P}_n \frac{\Delta_{n-1}}{|V_n|^2} = \frac{\hat{P}_n}{z - H_n - \Delta_n} \qquad \Rightarrow \Delta_n = z - H_n - \frac{|V_n|^2}{\Delta_{n-1}}. \tag{A3}
$$

The initial value $\Delta_0$ is given by the limit $\mathrm{Im}(z) \to 0^+$ of

$$
G_0(z) = \frac{1}{z - H_0 - \Delta_0} = \langle E_0 | \hat{G}(z) | E_0 \rangle = \int \frac{\rho(x)dx}{z - \sqrt{2}x} = -i \int_0^{+\infty} e^{itz} s(t)dt, \quad (A4)
$$

the noise spectral density ρ(x) and the amplitude s(t) = 〈E₀|e⁻iHT|E₀〉 = ∫ dxρ(x)e⁻i√2xt being defined in the main text.

We can now apply this procedure to the three specific cases discussed in the manuscript.

• For the Markovian case ρ<sub>M</sub> where H<sub>M,0</sub> = 0,

$$
s_M(t) = \int_{-\infty}^{+\infty} dx \frac{2e^{-i\sqrt{2}xt}}{\pi(1+4x^2)} = \frac{1}{2i\pi} \int_{-\infty}^{+\infty} dx e^{-ixt} \left( \frac{1}{x - i/\sqrt{2}} - \frac{1}{x + i/\sqrt{2}} \right) = e^{-t/\sqrt{2}} \quad (A5)
$$

$$
\Rightarrow G_{M,0}(z) = -i \int_{0}^{+\infty} dt e^{-t/\sqrt{2} + itz} = \frac{1}{z + i/\sqrt{2}}. \qquad (A6)
$$

Self-consistently, this gives $\delta_{M,0} = 0$ and $\gamma_{M,0} = 1/\sqrt{2}$ independently on $z$.

• For the Gaussian spectrum ρ<sub>G</sub>,

$$
s_G(t) = \frac{1}{\sqrt{\pi}} \int_{-\infty}^{+\infty} dx e^{-x^2 - i\sqrt{2}xt} = e^{-t^2/2}, \qquad (A7)
$$

$$
\Rightarrow G_{G,0}(z) = -i \int_{0}^{+\infty} dt \, e^{-t^2/2+itz} = -i\sqrt{2}e^{-z^2/2} \int_{-iz/\sqrt{2}}^{+\infty} dt \, e^{-t^2} = -i\sqrt{\frac{\pi}{2}} w\left(\frac{z}{\sqrt{2}}\right), \quad (A8)
$$

where $w(z) = e^{-z^2} \operatorname{erfc}(-iz)$ is the Faddeeva function. With $z = H_{G,0} = 0$, we get $\Delta_{G,0} = -1/G_{G,0}(0) = -i\sqrt{2}/\pi$ and, from Eq. A3,

$$
\Delta_{G,n} = \frac{n}{\Delta_{G,n-1}} = \frac{n}{n-1} \Delta_{G,n-2} = \dots \quad \Rightarrow \quad \delta_{G,n} = 0, \quad \gamma_{G,n} = \frac{(n)!!}{(n-1)!!} \left(\frac{2}{\pi}\right)^{(-1)^n/2}. \quad (A9)
$$
• For the light-shift spectrum $\rho_L$, assuming $\mathrm{Im}(z) \to 0^+$,

$$
\begin{align*}
G_{L,0}(z) &= \frac{2}{\sqrt{\pi}} \int_{0}^{+\infty} dx \frac{\sqrt{x}e^{-x}}{z - \sqrt{2}x} = -\sqrt{\frac{2}{\pi}} \int_{-\infty}^{+\infty} du \frac{u^2 e^{-u^2}}{u^2 - z/\sqrt{2}} \\
&= -\sqrt{2} - \sqrt{\frac{2}{\pi}} \int_{-\infty}^{+\infty} du e^{-u^2} \frac{1}{2} \sqrt{\frac{z}{\sqrt{2}}} \left( \frac{1}{u - \sqrt{z/\sqrt{2}}} - \frac{1}{u + \sqrt{z/\sqrt{2}}} \right) \\
&= -\sqrt{2} - i \sqrt{z\sqrt{2}} \frac{1}{\sqrt{\pi}} \int_{-\infty}^{+\infty} du \int_{0}^{+\infty} dt e^{-u^2 - it(u - \sqrt{z/\sqrt{2}})} = -\sqrt{2} - i \sqrt{z\sqrt{2}} \int_{0}^{+\infty} dt e^{-t^2/4 + it\sqrt{z/\sqrt{2}}} \\
&= -\sqrt{2} - i \sqrt{\pi} \sqrt{z\sqrt{2}} w \left( \sqrt{\frac{z}{\sqrt{2}}} \right). \tag{A10}
\end{align*}
$$

With $z = H_{L,0} = 3/\sqrt{2}$, $\Delta_{L,0} = 1/(\sqrt{2} + i\sqrt{3\pi}w(\sqrt{3}/2))$. Eq. A3 becomes $\Delta_n = -2^{3/2}n - n(2n+1)/\Delta_{n-1}$ and gives $\delta_n$ and $\gamma_n$ by induction.

Appendix B: Driving protection mechanism

Let us consider here that the collective qubit is resonantly driven between $|G\rangle$ and $|E_0\rangle$ with a constant Rabi frequency $\Omega$, assumed to be real. In the rotating frame, writing $\hat{\sigma}_{G0} = |G\rangle\langle E_0|$, the Hamiltonian becomes

$$
\hat{H}_{\Omega} = \hat{H}_{r} + \hat{D} \qquad \hat{H}_{r} = \sqrt{2}\hat{x} - H_{0}(1 - \hat{\sigma}_{GG}), \qquad \hat{D} = \frac{\Omega}{2}(\hat{\sigma}_{G0} + \hat{\sigma}_{0G}) \quad (B1)
$$

where the dephasing term $\hat{H}_r$ is offset by $H_0$ to make the driving $\hat{D}$ resonant. We define the resolvents $\hat{G}_r$ and $\hat{G}_\Omega$ of $\hat{H}_r$ and $\hat{H}_\Omega$, and use the general relation $\hat{B}^{-1} = \hat{A}^{-1} + \hat{A}^{-1}(\hat{A} - \hat{B})\hat{B}^{-1}$ to rewrite [21]

$$
\hat{G}_{\Omega}(z) = \hat{G}_{r}(z) + \hat{G}_{r}(z) \hat{D} \hat{G}_{\Omega}(z) = \hat{G}_{r}(z) + \hat{G}_{r}(z) \hat{D} \hat{G}'(z) + \hat{G}_{r}(z) \hat{D} \hat{G}'(z) \hat{D} \hat{G}_{\Omega}(z) = \dots = \sum_{n=0}^{+\infty} (\hat{G}_{r}(z) \hat{D})^{n} \hat{G}_{r}(z). \quad (B2)
$$

As $\hat{x}|G\rangle = 0$, we have

$$
\hat{G}_r(z) \hat{D} = \frac{\Omega}{2} \left( \frac{1}{z} \hat{\sigma}_{G0} + \frac{1}{z - \sqrt{2}\hat{x} + H_0(1 - \hat{\sigma}_{GG})} \hat{\sigma}_{0G} \right), \quad (B3)
$$

$$
\left(\hat{G}_{r}(z) \hat{D}\right)^{2} = \left(\frac{\Omega}{2}\right)^{2} \frac{1}{z} \left(G_{0}(z+H_{0})\hat{\sigma}_{GG} + \frac{1}{z-\sqrt{2}\hat{x}+H_{0}(1-\hat{\sigma}_{GG})}\hat{\sigma}_{00}\right), \quad (B4)
$$

$$
\left(\hat{G}_{r}(z) \hat{D}\right)^{3} = \left(\frac{\Omega}{2}\right)^{3} \frac{G_{0}(z+H_{0})}{z} \left(\frac{1}{z}\hat{\sigma}_{G0} + \frac{1}{z-\sqrt{2}\hat{x}+H_{0}(1-\hat{\sigma}_{GG})}\hat{\sigma}_{0G}\right) = \frac{\Omega^{2}G_{0}(z+H_{0})}{4z}\hat{G}_{0}(z)\hat{D} \quad (B5)
$$

$$
\Rightarrow \left(\hat{G}_{r}(z) \hat{D}\right)^{2n+k} = \left(\frac{\Omega^2 G_0(z + H_0)}{4z}\right)^n \left(\hat{G}_{r}(z) \hat{D}\right)^k, \qquad (B6)
$$

$$
\Rightarrow \hat{G}_{\Omega}(z) = \hat{G}_{r}(z) + \frac{\hat{G}_{r}(z) \hat{D} \hat{G}_{r}(z) + (\hat{G}'(z) \hat{D})^2 \hat{G}_{r}(z)}{1 - \Omega^2 G_0(z + H_0)/(4z)} \quad (B7)
$$

where $G_0$ is given by Eq. A4. The forward propagator [21]

$$
\hat{U}_{\Omega}(t > 0) = e^{-i\hat{H}_{\Omega}t} = \lim_{\eta \to 0^+} \int_{-\infty}^{+\infty} \frac{i dz}{2\pi} \hat{G}_{\Omega}(z + i\eta)e^{-itz}, \quad (B8)
$$

gives access to the dynamics of the qubit through the restriction of $\hat{G}_{\Omega}$ to the subspace $\{|G\rangle, |E_0\rangle\}$,

$$
\hat{G}_{\Omega,0}(z) = \frac{\hat{\sigma}_{GG} + G_0(z + H_0)(z \hat{\sigma}_{00} + \hat{D})}{z - \frac{\Omega^2}{4} G_0(z + H_0)}. \quad (B9)
$$
In particular, the amplitude of the qubit's excited state now evolves as

$$s_{\Omega}(t) = \langle E_0 | \hat{U}_{\Omega}(t) | E_0 \rangle = \int_{-\infty}^{+\infty} \frac{i dz}{2\pi} \frac{z G_0(z+H_0) e^{-itz}}{z - \frac{\Omega^2}{4} G_0(z+H_0)}. \quad (B10)$$

For $\Omega \to 0$ we recover the undriven case $s(t)$ discussed in the theoretical section of the main manuscript, up to a phase factor $e^{iH_0t}$ corresponding to the resonance condition. In the following we focus on the opposite strong-driving regime where $\Omega \gg 1$.

## 1. Markovian process

For a Lorentzian spectrum $\rho_M$, the simple form of $G_{M,0}$ in Eq. A5 gives

$$s_{M,\Omega}(t) = \int_{-\infty}^{+\infty} \frac{i dz}{2\pi} \frac{ze^{-itz}}{z(z + i\gamma_M) - \frac{\Omega^2}{4}}. \quad (B11)$$

where $\gamma_M = 1/\sqrt{2}$ is the driving-less decay rate derived in the theoretical section of the manuscript. Writing $\Omega_M = \sqrt{\Omega^2 - \gamma_M^2}$,

$$s_{M,\Omega}(t) = \int_{-\infty}^{+\infty} \frac{i dz}{2\pi} \frac{e^{-itz}}{2} \left( \frac{1 - i\gamma_M/\Omega_M}{z - \Omega_M/2 + i\gamma_M/2} + \frac{1 + i\gamma_M/\Omega_M}{z + \Omega_M/2 + i\gamma_M/2} \right) \\ = e^{-\gamma_M t/2} \left( \cos\left(\frac{\Omega_M t}{2}\right) - \frac{\gamma_M}{\Omega_M} \sin\left(\frac{\Omega_M t}{2}\right) \right). \quad (B12)$$

When $\Omega > \gamma_M$, the decay rate $\gamma_M/2$, two times smaller than without driving, becomes independent on the driving strength.

We can verify these results using standard Bloch equations. For consistency with the model above, let us assume that the qubit oscillates between the states $|g\rangle$ and $|e\rangle$ and decays into a different state $|d\rangle$. The master equation

$$\frac{d}{dt}\hat{\rho} = -i\frac{\Omega}{2}[\hat{\sigma}_{ge} + \hat{\sigma}_{eg}, \hat{\rho}] + \gamma(2\hat{\sigma}_{de}\hat{\rho}\hat{\sigma}_{ed} - \hat{\sigma}_{ee}\hat{\rho} - \hat{\rho}\hat{\sigma}_{ee}) \quad (B13)$$

leads to

$$\frac{d}{dt} \begin{pmatrix} \hat{\rho}_{ee} \\ \mathrm{Im}\hat{\rho}_{ge} \\ \hat{\rho}_{dd} \end{pmatrix} = - \begin{pmatrix} 2\gamma & -\Omega & 0 \\ \Omega & \gamma & \Omega/2 \\ -2\gamma & 0 & 0 \end{pmatrix} \begin{pmatrix} \hat{\rho}_{ee} \\ \mathrm{Im}\hat{\rho}_{ge} \\ \hat{\rho}_{dd} \end{pmatrix} + \begin{pmatrix} 0 \\ \Omega/2 \\ 0 \end{pmatrix}. \quad (B14)$$

For $\Omega > \gamma$, the real parts of the eigenvalues $\{\gamma, \gamma \pm i\sqrt{\Omega^2 - \gamma^2}\}$ of the matrix are all equal to $\gamma$, independently on $\Omega$. Setting $\gamma = \gamma_M$ allows one to recover the results above.

## 2. Gaussian process

For a Gaussian spectrum, $H_{G,0} = 0$ and $G_{G,0}(z)$ is given by Eq. A8, which can be rewritten as

$$G_{G,0}(z) = \sqrt{2D} \left( \frac{z}{\sqrt{2}} \right) - i \sqrt{\frac{\pi}{2}} e^{-z^2/2}, \qquad D(z) = e^{-z^2} \int_0^z e^{t^2} dt, \quad (B15)$$

where $D(z)$ is the Dawson integral. The even function $D(z)/z$ monotonically decreases from 1 to 0 on $\mathbb{R}_+$, which means that for $\Omega \gg 1$ the real part of the denominator in the integrand of Eq. B10 vanishes in $z_{\pm} = \pm \Omega_G/2$ with $\Omega_G \gg 1$, where its imaginary part becomes very small. We can recursively use the asymptotic expansion

$$\mathrm{Re}(G_{G,0}(z)) = \sum_{k=0}^{N} \frac{(2k-1)!!}{z^{2k+1}} + O\left(\frac{1}{z^{2N+3}}\right), \quad (B16)$$
to solve $\Omega_G/2 = (\Omega^2/4)\text{Re}(G_{G,0}(\Omega_G/2))$ and find $\Omega_G$ to an arbitrary order in $\Omega$. To the lowest non-trivial order, $\Omega_G = \Omega + 2/\Omega + O(\Omega^{-3})$. Near $z_{\pm}$, to the first order in $z - z_{\pm}$ and in $1/\Omega$,

$$z - \frac{\Omega^2}{4} \text{Re}((G_{G,0}(z)) \approx (z - z_{\pm}) \left(1 - \frac{\Omega^2}{4} \text{Re}((G'_{G,0}(z_{\pm}))\right) \approx (z - z_{\pm}) \left(1 + \frac{\Omega^2}{4} \frac{1}{z_{\pm}^2}\right) \approx 2(z - z_{\pm}). \quad (B17)$$

As the imaginary part of the integrand's denominator is suppressed exponentially in $\Omega^2$ by $\text{Im}(G_{G,0}(z_{\pm})) = -\sqrt{\pi/(2e)}e^{-\Omega^2/8}(1 + O(\Omega^{-2}))$, we find

$$\frac{zG_{G,0}(z)}{z - \frac{\Omega^2}{4}G_{G,0}(z)} \approx \frac{1}{2} \left( \frac{1}{z - \frac{\Omega_G}{2} + i\gamma_{G,\Omega}} + \frac{1}{z + \frac{\Omega_G}{2} + i\gamma_{G,\Omega}} \right), \quad \Omega_G \approx \Omega + \frac{2}{\Omega}, \quad \gamma_{G,\Omega} \approx \frac{\Omega^2}{8}\sqrt{\frac{\pi}{2e}}e^{-\Omega^2/8}. \quad (B18)$$

Eq. B10 then integrates into $s_G(t) \approx e^{-\gamma_{G,\Omega}t} \cos(\Omega_G t/2)$, showing that strong driving turns the Gaussian decay into an exponential one, and that the rate $\gamma_G$ of this exponential decay decreases extremely fast when $\Omega$ increases. Even for a moderate $\Omega = 10$, $\gamma_G$ drops below $10^{-4}$, allowing Rabi oscillations to persist for much longer than the decay time of the undriven qubit.

### 3. Light shift in a harmonic trap

The case of the harmonic light shift spectrum $\rho_L$ resembles the Gaussian one. Here $H_{L,0} = 3/\sqrt{2}$ and the asymptotic expansion of the real part of Eq. A10 gives

$$\text{Re}(G_{L,0}(z)) = \frac{1}{z} \sum_{k=0}^{N} \frac{(2k+1)!!}{(\sqrt{2}z)^k} + O\left(\frac{1}{z^{N+2}}\right) \quad \Rightarrow \quad \text{Re}(G_{L,0}(z+H_{L,0})) = \frac{1}{z} + \frac{3}{z^3} + O\left(\frac{1}{z^4}\right). \quad (B19)$$

This shows that real part of the denominator of the integrand in Eq. B10 vanishes in two nearly opposite points $z_{\pm} = \pm(\Omega/2 + 3/\Omega) + 12\sqrt{2}/\Omega^2 + O(1/\Omega^3)$. Using $w(iz) = e^{iz^2}\text{erfc}(z)$ to rewrite

$$G_{L,0}(-z) = -\sqrt{2} + \sqrt{\pi z \sqrt{2} e^{z/\sqrt{2}}} \text{erfc}\left(\sqrt{\frac{z}{\sqrt{2}}}\right), \quad (B20)$$

we see that, for $z > 0$, $\text{Im}(G_{L,0}(z)) = -\sqrt{2}\sqrt{\pi z}e^{-z/\sqrt{2}}$ but $\text{Im}(G_{L,0}(-z)) = 0$. With $G'_{L,0}(z+H_{L,0}) = -z^{-2}+O(z^{-4})$, the same approximations as in the Gaussian case give again

$$z - \frac{\Omega^2}{4} \text{Re}((G_{L,0}(z + H_{L,0})) \approx (z - z_{\pm}) \left(1 - \frac{\Omega^2}{4} \text{Re}((G'_{L,0}(z_{\pm}))\right) \approx (z - z_{\pm}) \left(1 + \frac{\Omega^2}{4} \frac{1}{z_{\pm}^2}\right) \approx 2(z - z_{\pm}). \quad (B21)$$

Near $z_+$, the imaginary part of the integrand's denominator now decreases exponentially with $\Omega$, while near $z_-$ it vanishes:

$$\frac{\Omega^2}{4} \text{Im}((G_{L,0}(z_+ + H_{L,0})) = \frac{\Omega^2}{4} \sqrt{\frac{\pi \Omega}{\sqrt{2} e^3}} e^{-\Omega/\sqrt{8}} (1 + O(\Omega^{-2})), \quad \frac{\Omega^2}{4} \text{Im}((G_{L,0}(z_- + H_{L,0})) = 0 \quad (B22)$$

This gives

$$\frac{zG_{L,0}(z)}{z - \frac{\Omega^2}{4}G_{L,0}(z)} \approx \frac{1}{2} \left( \frac{1}{z - \frac{\Omega_L}{2} + i\gamma_{L,\Omega}} + \frac{1}{z + \frac{\Omega_L}{2}} \right), \quad \Omega_L \approx \Omega + \frac{6}{\Omega}, \quad \gamma_{L,\Omega} \approx \frac{\Omega^2}{8}\sqrt{\frac{\pi\Omega}{\sqrt{2}e^3}}e^{-\Omega/\sqrt{8}}, \quad (B23)$$

and integrates into

$$s_L(t) = \langle E_0 | \hat{U}_{L,\Omega}(t) | E_0 \rangle \approx (e^{-\gamma_{L,\Omega}t}e^{-i\Omega_L t/2} + e^{i\Omega_L t/2})/2. \quad (B24)$$

Using Eqs. B8 and B9, we can apply the same procedure to the three other elements of the restricted propagator $\hat{U}_{L,\Omega}(t > 0)$:

$$\langle G | \hat{U}_{L,\Omega}(t) | G \rangle = \int_{-\infty}^{+\infty} \frac{i dz}{2\pi} \frac{e^{-itz}}{z - \frac{\Omega^2}{4} G_{L,0}(z + H_{L,0})} \approx \frac{e^{-\gamma_{L,\Omega}t} e^{-i\Omega_L t/2} + e^{i\Omega_L t/2}}{2} = \langle E_0 | \hat{U}_{L,\Omega}(t) | E_0 \rangle, \quad (B25)$$

$$\langle E_0 | \hat{U}_{L,\Omega}(t) | G \rangle = \int_{-\infty}^{+\infty} \frac{i dz}{2\pi} \frac{\Omega^2 G_{L,0}(z + H_{L,0}) e^{-itz}}{z - \frac{\Omega^2}{4} G_{L,0}(z + H_{L,0})} \approx \frac{e^{-\gamma_{L,\Omega}t} e^{-i\Omega_L t/2} - e^{i\Omega_L t/2}}{2} = \langle G | \hat{U}_{L,\Omega}(t) | E_0 \rangle. \quad (B26)$$
Thus, in the strong driving limit, the propagator restricted to the qubit’s subspace can be written as

$$ \hat{U}_{L,\Omega}(t) \approx e^{-\gamma_{L,\Omega}t} e^{-i\Omega_L t/2} |\psi_+\rangle\langle\psi_+| + e^{i\Omega_L t/2} |\psi_-\rangle\langle\psi_-| \quad (B27) $$

where $\psi_{\pm} = (|G\rangle \pm |E_0\rangle)/\sqrt{2}$ are the eigenstates of the Rabi flopping Hamiltonian $\hat{D}$. The “driving protection” of $|\psi_+\rangle$ is weaker than for the Gaussian spectrum, as $\gamma_{L,\Omega}$ decreases exponentially with $\Omega$ instead of $\Omega^2$, while $|\psi_-\rangle$ does not decay at all.

## Appendix C: Simulations

The simulations presented together with our experimental results in the main manuscript use the Python package QuTiP [36] to integrate the GKSL master equation

$$ \frac{d\hat{\rho}}{dt} = -i \left[ \hat{H}_{\text{Drive}} + \hat{H}_{\text{Deph}}, \hat{\rho} \right] + \mathcal{L}_{\text{Deph}}[\hat{\rho}] + \mathcal{L}_{\text{Loss}}[\hat{\rho}], \quad (C1) $$

incorporating the following terms.

### 1. Driving term $\hat{H}_{\text{Drive}}$ and fluctuations of $\Omega$.

In $\hat{H}_{\text{Drive}} = (\Omega/2)(|E_0\rangle\langle G| + |G\rangle\langle E_0|)$, describing the Rabi flopping of the qubit, the temporal profile of the collective Rabi frequency $\Omega = \sqrt{N\Omega_r\Omega_b/(2\Delta)}$ is obtained from “red drive” and “blue drive” laser pulse shapes measured with photodiodes and plotted in Fig. 5(b) in the main manuscript. The detuning $\Delta$ and the “blue drive” Rabi frequency $\Omega_b$ being fixed (see “Experimental Methods” in the manuscript), the proportionality factor between $|\Omega|^2$ and the controlled power of the “red drive” beam is extracted from experimental Rabi oscillations in the regime where $\Omega \gg \omega_G, \omega_L \sim 0.1$ MHz. The Poissonian shot-to-shot fluctuations of the atomic number $N$, which can be approximated by a Gaussian as $N$ is large, result in fluctuations of $\Omega$ with a relative standard deviation $\sigma_{\Omega} \approx 0.024\Omega$, increased to the experimentally-measured value $\sigma_{\Omega} \approx 0.028\Omega$ by residual fluctuations of laser power. We account for these fluctuations by averaging the results of each simulation over 500 runs, drawing $\Omega$ from a Gaussian distribution.

We characterize these fluctuations by fitting the experimental Rabi oscillations in the regime $\Omega \gg \omega_G, \omega_L$, where the oscillation envelope is primarily governed by Rabi frequency fluctuations. Gaussian-distributed fluctuations in $\Omega$ are incorporated through the following fit function:

$$ f(t) = \frac{A}{2} \left( 1 - \cos(\Omega t) e^{-\frac{1}{2} \sigma_{\Omega}^2 t^2} \right) e^{-\Gamma t}, \quad (C2) $$

where $A$ is an amplitude scaling factor and $\Gamma$ the exponential decay rate represented in Fig. 6(e). Given the linear dependence of $\sigma_{\Omega}$ on $\Omega$, we fix the fluctuation parameter in the fit to extract $\Gamma$ for any Rabi frequencies. However, when the decay rate $\Gamma$ becomes small, due to the protection mechanism, the envelope becomes dominated by fluctuations of $\Omega$. In this case, $\Gamma$ is difficult to estimate, leading to increased uncertainties for low values of $\Gamma$ visible in Fig. 6(e) of the main text.

### 2. Inhomogeneous dephasing terms $\hat{H}_{\text{Deph}}$ and $\mathcal{L}_{\text{Deph}}$

The Hamiltonian $\hat{H}_{\text{Deph}}$ and the Liouvillian decay term $\mathcal{L}_{\text{Deph}}$, describing inhomogeneous dephasing, are determined using the procedure described in the theoretical section of the manuscript and in the section A above. We account for non-Markovianity by including in $\hat{H}_{\text{Deph}}$ $n = 15$ Dicke states for the purely Gaussian Doppler broadening, and $n = 20$ for the light-shifted one.

For Gaussian dephasing, we extract $\omega_G = 2\pi \times 61(8)$ kHz from the fit of the measured population decay $|s_G(t)|^2$ presented in “Experimental results” section of the main manuscript. This value is consistent with the temperature obtained by time-of-flight measurement with absorption imaging, yielding $\omega_G = 2\pi \times 66(1)$ kHz.

For the light-shifted case, the scaling angular frequency $\omega_S$ is obtained by independently measuring the temperature, the average light shift, and the width of the Gaussian atomic density profile in the trap. The anharmonicity of the trap, formed by crossed Gaussian beams, leads to a noticeable deviation between $\rho_L$ and the actual light-shift distribution $\rho_{L'}$, shown in Fig. C.1(b). We reproduce $\rho_{L'}$ via Monte Carlo sampling, by randomly drawing a large number of atoms
<figure><img src="image_8.png" /><figcaption>FIG. C.1. **Light-shift and Gaussian inhomogeneous dephasing spectra.** (a) Light shift induced by a Gaussian dipole trap. (b) Inhomogeneous dephasing spectra. Red: Gaussian experimental Doppler shift distribution. Green: light-shift distribution $P_L(\omega)$ for a thermal cloud, assuming that the trap is harmonic. Blue: light-shift distribution $P_{L'}(\omega)$ accounting for the trap's anharmonicity, obtained by Monte Carlo sampling using experimental parameters. Grey and black: full inhomogeneous dephasing spectrum combining a Gaussian Doppler broadening with a light shift in an anharmonic trap. Grey circles correspond to the Monte-Carlo-sampled spectrum, the full black line is a fit by $P_{GL'}(\omega) = |\psi_{GL',0}(\omega/(\sqrt{2}\omega_{GL'}))^2/(\sqrt{2}\omega_{GL'})$ where $\psi_{GL',0}(x)$ is given by Eq. C3, the frequency scale $\omega_{GL'}$ and the Fock state coefficients $\{\beta_k\}$ being used as free parameters.</figcaption></figure>

from the spatial distribution obtained by absorption imaging, calculating the light shift distribution knowing the trap's geometry, and rescaling it using the measured mean value of the light shift. The velocity-dependent Doppler shifts being statistically independent on the position-dependent light shifts in a thermal cloud, the total inhomogeneous broadening distribution $\rho_{GL'}$ is obtained either by a convolution of $\rho_{L'}$ with a Gaussian $\rho_G$ or by an additional Monte Carlo sampling of Gaussian-distributed atomic velocities. We center $\rho_{GL'}$ and approximate the wavefunction of the symmetric Dicke state $\psi_0(x) = \sqrt{\rho_{GL'}(x)}$ by a superposition of the first five Fock states,

$$ \psi_{GL',0}(x) = \sum_{k=0}^{4} \beta_k \frac{1}{\pi^{1/4} \sqrt{2^k k!}} \mathcal{H}_k(x) e^{-x^2/2}, \quad (C3) $$

where $\beta_k$ are the coefficients and $\mathcal{H}_k$ are the Hermite polynomials. The procedure described in the theoretical section of the main manuscript is then applied algebraically, involving only the list of coefficients $\{\beta_k\}$.

### 3. Additional loss term $\mathcal{L}_{Loss}$

Additional losses from the state $|E_0\rangle$, described by the Liouvillian term $\mathcal{L}_{Loss}$, arise from a combination of several effects left aside in the idealized model above. Among them, imperfections in the Rydberg blockade or a residual population of the short-lived intermediate state $|f\rangle$ stemming from the finite detuning $\Delta$ are expected to induce loss rates scaling as $\Omega^2$, making the overall additional damping rate evolve as $\Gamma_\Omega\Omega^2 + \Gamma_0$. We estimate this rate in the strong driving regime with no light shifts, such as $\Omega \gg \omega_G$, where driving protection makes inhomogeneous dephasing negligible. We fit the experimental oscillations to extract $\Gamma_{\Omega,\text{fit}} \times 2\pi = 4.6(1)10^{-3} \text{ MHz}^{-1}$ and $\Gamma_{0,\text{fit}} = 2\pi \times 2.6(1.5) \text{ kHz}$.