Dynamical Exploration of Amplitude Bistability in Engineered Quantum Systems

TL;DR

This paper demonstrates amplitude bistability in a hybrid quantum system with long-lived nitrogen-vacancy centers coupled to a superconducting resonator, revealing critical slowing down over unprecedented timescales, and lays groundwork for future quantum technologies.

Contribution

It introduces the experimental observation of amplitude bistability in a hybrid quantum system with long coherence times, highlighting a new regime of nonlinear quantum dynamics.

Findings

Observation of amplitude bistability in a hybrid quantum system.

Demonstration of critical slowing down over tens of thousands of seconds.

Long spin lifetime enables exploration of nonlinear quantum phenomena.

Abstract

Nonlinear systems, whose outputs are not directly proportional to their inputs, are well known to exhibit many interesting and important phenomena which have profoundly changed our technological landscape over the last 50 years. Recently the ability to engineer quantum metamaterials through hybridisation has allowed to explore these nonlinear effects in systems with no natural analogue. Here we investigate amplitude bistability, which is one of the most fundamental nonlinear phenomena, in a hybrid system composed of a superconducting resonator inductively coupled to an ensemble of nitrogen-vacancy centres. One of the exciting properties of this spin system is its extremely long spin life-time, more than ten orders of magnitude longer than other relevant timescales of the hybrid system. This allows us to dynamically explore this nonlinear regime of cavity quantum electrodynamics (cQED) and demonstrate a critical slowing down of the cavity population on the order of several tens of thousands of seconds - a timescale much longer than observed so far for this effect. Our results provide the foundation for future quantum technologies based on nonlinear phenomena.

Equations

$$\mathcal{H}=\hslash {\omega }_c{\mathbf{a}}^{†}\mathbf{a}+\frac{\hslash }{2}\sum _{j=1}^N{\omega }_j{\mathbf{\sigma }}_j^z+i\hslash \sum _{j=1}^N{g}_j\left({\mathbf{\sigma }}_j^{-}{\mathbf{a}}^{†}-{\mathbf{\sigma }}_j^{+}\mathbf{a}\right)$$

$$+i\hslash \left(\eta {\mathbf{a}}^{†}{e}^{-i{\omega }_{p}t}-h.c.\right),$$

$$\dot{a}=-\kappa a+{\sum }_j{g}_j{\sigma }_j^{-}+\eta$$

$${\dot{\sigma }}_j^{-}=-\left({\gamma }_{\perp }+i{\Theta }_j\right){\sigma }_j^{-}+g_j{\sigma }_j^{z}a$$

$${\dot{\sigma }}_j^z=-{\gamma }_{\parallel }\left(1+{\sigma }_j^z\right)-2g_j\left({\sigma }_j^{-}{a}^{†}+{\sigma }_j^{+}a\right),$$

$${\left|a\right|}^2=\frac{{\eta }^2}{{\kappa }^2}{\left(1-{\sum }_j{C}_j{\sigma }_j^z\right)}^{-2},{\sigma }_j^z=-{\left(1+\frac{4g_j^2|a{|}^2{\gamma }_{\perp }}{{\gamma }_{\parallel }({\gamma }_{\perp }^2+{\Theta }_j^2)}\right)}^{-1},$$

$$\frac{\mathrm{d}|a{|}^2}{\mathrm{d}t}$$

Full text

††footnotetext: § These authors contributed equally to this work

Dynamical Exploration of Amplitude Bistability in Engineered Quantum Systems

Andreas Angerer§

andreas.angerer@tuwien.ac.at

Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

Zentrum für Mikro- und Nanostrukturen, TU Wien, Floragasse 7, 1040 Vienna, Austria

  

Stefan Putz§

Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

Zentrum für Mikro- und Nanostrukturen, TU Wien, Floragasse 7, 1040 Vienna, Austria

Department of Physics, Princeton University, Princeton, New Jersey 08544, USA

  

Dmitry O. Krimer

Institute for Theoretical Physics, TU Wien, Wiedner Hauptstraße 8-10/136, 1040 Vienna, Austria

  

Thomas Astner

Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

Zentrum für Mikro- und Nanostrukturen, TU Wien, Floragasse 7, 1040 Vienna, Austria

  

Matthias Zens

Institute for Theoretical Physics, TU Wien, Wiedner Hauptstraße 8-10/136, 1040 Vienna, Austria

  

Ralph Glattauer

Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

  

Kirill Streltsov

Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

  

William J. Munro

NTT Basic Research Laboratories, 3-1 Morinosato-Wakamiya, Atsugi, Kanagawa 243-0198, Japan

  

Kae Nemoto

National Institute of Informatics, 2-1-2 Hitotsubashi, Chiyoda-ku, Tokyo 101-8430, Japan

  

Stefan Rotter

Institute for Theoretical Physics, TU Wien, Wiedner Hauptstraße 8-10/136, 1040 Vienna, Austria

  

Jörg Schmiedmayer

Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

  

Johannes Majer

Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, Stadionallee 2, 1020 Vienna, Austria

Zentrum für Mikro- und Nanostrukturen, TU Wien, Floragasse 7, 1040 Vienna, Austria

Nonlinear systems, whose outputs are not directly proportional to their inputs, are well known to exhibit many interesting and important phenomena which have profoundly changed our technological landscape over the last 50 years. Recently the ability to engineer quantum metamaterials through hybridisation imamoglu_cavity_2009; kubo_hybrid_2011; amsuss_cavity_2011; schuster_high-cooperativity_2010; xiang_hybrid_2013; kurizki_quantum_2015 has allowed to explore these nonlinear effects in systems with no natural analogue. Here we investigate amplitude bistability rempe_optical_1991, which is one of the most fundamental nonlinear phenomena, in a hybrid system composed of a superconducting resonator inductively coupled to an ensemble of nitrogen-vacancy centres jelezko_observation_2004; doherty_nitrogen-vacancy_2013. One of the exciting properties of this spin system is its extremely long spin life-time, more than ten orders of magnitude longer than other relevant timescales of the hybrid system jarmola_temperature-_2012. This allows us to dynamically explore this nonlinear regime of cavity quantum electrodynamics (cQED) and demonstrate a critical slowing down of the cavity population on the order of several tens of thousands of seconds - a timescale much longer than observed so far for this effect mitschke_transients_1983; garmire_transient_1979. Our results provide the foundation for future quantum technologies based on nonlinear phenomena.

In nature, most physical systems are inherently nonlinear, giving rise to effects such as bistability gibbs_optical_1979, chaos benza_symmetry_1987, solitons mollenauer_experimental_1980, superradiance dicke_coherence_1954 and often appear counter-intuitive when contrasted with much simpler linear systems. Amplitude bistability, one of the basic nonlinear phenomena (nowadays commonly used in optical switches szoke_bistable_1969) has been extensively investigated both theoretically casteels_power_2016; bonifacio_transient_1978; sawant_optical-bistability-enabled_2016; dombi_optical_2013; martin_extreme_2011 and experimentally rempe_optical_1991; kimble_strong_1998. It occurs in any medium where strong nonlinearities in the interaction between a radiation field and a polarisable medium such as spins exist. The nonlinearity in such systems arises from the two-level nature of the atoms coupled to the cavity mode but only shows up when driven beyond the single excitation regime. For a strong coupling between the spin system and the cavity mode a first order phase transition between a saturated, disordered and a de-excited, ordered ground state occurs bonifacio_cooperative_1976. The coupled system switches between these two branches and shows a hysteresis depending on the history of the system.

The usual cQED demonstrations use atoms or trapped ions coupled to optical light fields to investigate these nonlinear effects but the short atomic life-times have made it difficult to truly observe the temporal dynamics of amplitude bistability mitschke_transients_1983; garmire_transient_1979, and restricted previous studies to the steady state behaviour. In contrast emerging quantum engineering through hybridisation allows us to create systems with extremely long-lived emitters, making it possible to observe the bistable system during evolution. In this letter, we report on the observation of amplitude bistability in a cQED system composed of a superconducting resonator coupled to a long-lived electron spin ensemble formed from artificial atoms (NV- centres in diamond jelezko_observation_2004; doherty_nitrogen-vacancy_2013). Their extremely long life-times allow us to study the temporal behaviour of the presented effect, a regime experimentally not accessible so far.

An ensemble of spins in a cavity is characterised by the three quantities polarisation, inversion and the cavity amplitude. They can be derived from the driven Tavis-Cummings Hamiltonian tavis_exact_1968 for $N$ spins under a rotating wave approximation,

[TABLE]

with $\bm{a}^{\dagger}$, $\bm{a}$ as the creation and annihilation operators for the cavity mode of frequency $\omega_{c}$ and $\bm{\sigma}_{j}^{z}$, $\bm{\sigma}_{j}^{+}$, $\bm{\sigma}_{j}^{-}$ as the spin inversion, raising and lowering operators for the $j$-th spin of frequency $\omega_{j}$ coupled to the cavity with a single spin coupling strength $g_{j}$. The last term accounts for an external cavity drive with field amplitude $\eta$ and frequency $\omega_{p}$.

Using a mean-field approximation, valid in the limit of very large spin ensembles, $\braket{\bm{a}^{\dagger}\bm{\sigma}_{-}}\approx a^{\dagger}\sigma_{-}$ (from now on unbolded symbols will be used for the expectation values), we derive a set of first order differential equations, formally equivalent to the well known Maxwell-Bloch equations bonifacio_theory_1982 as

[TABLE]

with cavity dissipation rate $\kappa$, transversal spin relaxation rate $\gamma_{\perp}=1/{T_{2}}$ and longitudinal spin relaxation rate $\gamma_{\parallel}=1/{T_{1}}$. The relaxation rates are ordered as $\kappa>\gamma_{\perp}\gg\gamma_{\parallel}$, such that the spin inversion is by far the slowest process. $\Theta_{j}$ are frequency detunings with respect to the ensemble central frequency to account for inhomogeneous broadening. Setting the time derivatives to zero we obtain the steady state of this system as

[TABLE]

where the dimensionless parameter ${C_{j}=g_{j}^{2}/\left[\kappa\gamma_{\perp}\left(1+\Theta_{j}^{2}/\gamma_{\perp}^{2}\right)\right]}$ is the single spin cooperativity. The collective system cooperativity is given accordingly by $C_{\mathrm{coll}}=\sum_{j}C_{j}$.

We can classify the expected system phase transition by deriving a solution for the time dependent cavity amplitude $|a(t)|^{2}$. The difference in dissipation rates allows us to adiabatically eliminate the $a$ and ${\sigma}^{-}_{j}$ variables lugiato_nonlinear_2015, which results in a first order differential equation for the intra-cavity intensity. For the giant $S_{z}=\sum_{j}\sigma^{z}_{j}$ spin in resonance ($\Theta_{j}=0$) with the cavity mode it can be written as

[TABLE]

For our typical system parameters this equation predicts a first order phase transition that connects a strongly driven branch and a weakly driven branch, with hysteresis and two saddle-node bifurications strogatz_nonlinear_2001 at which a jump occurs between both branches.

Our hybrid system is illustrated in Fig. 1 and is composed of an electron spin ensemble formed by NV centres in diamond, loaded onto a superconducting $\lambda/2$-resonator. To thermally polarise the $N\approx 10^{12}$ electron spins to their ground state ($\geq$99 %) we put the system in a dilution refrigerator at $25\text{\,}\mathrm{m}\mathrm{K}$. Each electron spin has a zero field splitting of $D/2\pi\sim$ $2.878\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ and average coupling rate of $g_{0}/2\pi\approx$12\text{,}\mathrm{H}\mathrm{z}$$ to the cavity mode. The superconducting resonator has a cavity linewidth of $\kappa/2\pi$=$440\pm 10\text{\,}\mathrm{k}\mathrm{H}\mathrm{z}$ (HWHM) with a fundamental resonance frequency at $\omega_{c}/2\pi$=$2.691\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}$ and a loaded quality factor of $Q=3300$. An external microwave field with frequency $\omega_{p}$ is used to probe our hybrid system.

First we search for the steady state bistable behaviour by measuring the transmitted intensities through the cavity defined by $|\mathrm{T}|^{2}=\mathrm{P}_{\mathrm{out}}/\mathrm{P}_{\mathrm{in}}$ as a function of the input drive intensity $\mathrm{P}_{\mathrm{in}}\approx\eta^{2}/\kappa$ and outgoing intensity $\mathrm{P}_{\mathrm{out}}\approx\left|a\right|^{2}\kappa$. The drive power is raised in a stepwise manner, slow enough to allow the system to reach a steady state for each stimulus $\mathrm{P}_{\mathrm{in}}$. For small excitations the intra-cavity intensity is not sufficient to saturate the spin ensemble ($\sigma_{j}^{z}\sim-1$) and is thus given by $\left|a\right|^{2}\approx\frac{\eta^{2}}{\kappa^{2}}\frac{1}{(1+C_{\mathrm{coll}})^{2}}$. As the power level increases, the cavity field bleaches the spins ($\sigma_{j}^{z}\approx\sigma_{j}^{-}\sim 0$) such that the Rabi-splitting vanishes and the spin system decouples from the cavity (Fig. 2). The intra-cavity intensity $\left|a\right|^{2}\approx\frac{\eta^{2}}{\kappa^{2}}$ is that of an empty cavity from which spins are completely decoupled. We identify the power level for which the transition between both cases occur as the critical drive $\mathrm{P}_{\mathrm{crit}}$.

This nonlinear saturation behaviour is a necessary precursor to the observation of amplitude bistability. However, whether or not this is observable in the experiment is determined by the system’s collective cooperativity. This is apparent from Eq. (3), where larger cooperativity values result in stronger nonlinearity and thus a larger phase seperation. In Fig. 3 we show steady state bistability measurements for three cooperativity values $C_{\mathrm{coll}}=18,49,78$. The lowest value $C_{\mathrm{coll}}=18$ does not show bistability (Fig. 3a), but increasing the cooperativity to $C_{\mathrm{coll}}=49$ (see Methods) allows us to observe the first signs of bistable behaviour (see Fig. 3b). Increasing the cooperativity further to $C_{\mathrm{coll}}\approx 78$ shows clear amplitude bistability (Fig. 3c) within a 2dB range. This steady state bistability behaviour is well reproduced by a full numerical simulation taking inhomogeneouse broadening into account (dashed lines in Fig. 3a-c).

Given such clear evidence of amplitude bistability, we focus next on the temporal behaviour of the hybrid system using quench dynamic measurements. Here the system is prepared in an initial steady state in which the spin ensemble is completely saturated and constantly driven with $\mathrm{P}_{\mathrm{in}}\gg\mathrm{P}_{\mathrm{crit}}$. The drive power is non-adiabatically switched to a lower drive level and the system transmission is monitored. We repeat this measurement several times, always preparing the system in the same initial state, but switching to different lower drive powers. When the system is driven close to the bifurication point ($\mathrm{P}_{\mathrm{in}}\approx\mathrm{P}_{\mathrm{crit}}$) the timescales needed to settle in a stationary state become as long as $40\text{\times}{10}^{3}\text{\,}\mathrm{s}$, as depicted in Fig. 4. This is known as critical slowing down grynberg_critical_1983; garmire_transient_1979.

The behaviour can be linked to our model given in Eq. (4) which predicts that our system features two fixed points at which a saddle-node bifurication occurs. A branch with no stable solutions connects two stable branches, one weakly driven, with de-excited and ordered spins, and the other one strongly driven, with unordered and saturated spins (see Fig. 3c). Starting in the strongly driven upper branch with a large intra-cavity intensity, the saturated spin system ($\sigma_{j}^{z}\approx\sigma_{j}^{-}\sim 0$) is entirely decoupled from the cavity. This prohibits collective decay into the cavity leaving the extremely small longitudinal relaxation rate $\gamma_{\parallel}$ as the dominant decay channel. In the opposite limit of drive intensities much smaller than the critical value, the spin system starts to collectively decay through the cavity early on with a rate much larger than $\gamma_{\parallel}$ and no critical slowing down is observed.

At the critical drive intensity, the saturation due to the drive and decay of the spin system are equal and opposite in effect allowing an everlasting decay to occur without the system reaching its fixed point. If however, the drive intensity is slightly smaller than the critical drive, we observe a critical slowing down of the decay, but eventually the decay of the collectively decaying spin system becomes dominant which leads to a buildup of correlations in the system and a faster decay of the spins and thus the cavity population. This behaviour is shown in Fig. 4a-c where close to a critical drive the system evolves towards the upper unstable fixed point, with a time derivative that can approach zero arbitrarily closely (inset in Fig. 4b). Small deviations from the critical drive lead to a speed up in decay until the system relaxes to a real steady state in the end. The time it takes to go from the upper to the lower branch diverges close to the critical drive according to $t_{\mathrm{switch}}\sim|\mathrm{P}_{\mathrm{in}}-\mathrm{P}_{\mathrm{crit}}|^{-\alpha}$, as shown in Fig. 4c with $\alpha\sim 1.2$. This algebraic divergence is characteristic for nonlinear systems showing saddle-node bifurications kuehn_scaling_2009. Comparing these experimental results with the full numerical solutions of Eqs. (Dynamical Exploration of Amplitude Bistability in Engineered Quantum Systems) including inhomogeneous broadening, we observe excellent agreement (see Fig. 4a).

To summarise, we have shown how a hybrid system composed of a superconducting resonator coupled to an electron spin ensemble in diamond can be used to explore amplitude bistability in new regimes of cQED, with unusual decay rates where the spin life-time is much smaller than other decay constants in the system. We have demonstrated a critical slowing down of the cavity population on the order of eleven hours, a timescale several orders of magnitude longer than observed so far for this effect and many order of magnitude longer than other time scales associated with the system. Our experiment provides a foundation for the exploration of new nonlinear phenomena in quantum metamaterials and future quantum technologies that may arise from it.

**Acknowledgements

**We would like to thank Helmut Ritsch, Michael Trupke and Alberto Amo for helpful discussions. The experimental effort lead by J.M. has been supported by the Top-/Anschubfinanzierung grant of the TU Wien. S.P., A.A. and T.A. acknowledge support by the Austrian Science Fund (FWF) in the framework of the Doctoral School “Building Solids for Function” Project W1243. D.O.K. and S.R. acknowledge funding by the Austrian Science Fund (FWF) through the Spezialforschungsbereich (SFB) NextLite Project No. F49-P10. K.N. acknowledges support from the MEXT KAKENHI Grant-in-Aid for Scientific Research on Innovative Areas “Science of hybrid quantum systems” Grant No. 15H05870.

**Author contributions

**S.P. and D.O.K. conceived the idea and S.P., J.S, and J.M. designed and set up the experiment. A.A., S.P., R.G., T.A., and K.S. carried out the measurements under the supervision of J.M.. D.O.K, M.Z., and S.R. devised the theoretical framework and provided the theoretical support for modelling the experiment. K.N. and W.J.M provided support for implementation of the theoretical methods. A.A. wrote the manuscript to which all authors suggested improvements.

**Methods

**

**Sample

**The spin system is realised by enhancing a type Ib high-pressure high-temperature (HPHT) diamond crystal containing an initial concentration of $200\text{\,}\mathrm{p}\mathrm{p}\mathrm{m}$ nitrogen with a natural abundance of $\prescript{13}{}{\mathrm{C}}$ nuclear isotopes. We achieve a total density of $\approx$6\text{,}\mathrm{p}\mathrm{p}\mathrm{m} negatively charged NV centres by $50\text{\,}\mathrm{h}$ of neutron irradiation with a fluence of $5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-2}$ and annealing the crystal for $3\text{\,}\mathrm{h}$ at $900\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Excess nitrogen P1 centres ($S=1/2$), uncharged $\mathrm{NV}^{0}$ centres, plus additional lattice stress serves as the main source of inhomogeneous broadening, which exceeds decoherence due to the naturally abundant 1.1% $\prescript{13}{}{\mathrm{C}}$ spin bath. The characteristics of the diamond crystal and NV ensemble were initially determined at room temperature using an optical confocal microscope. In the present experiment, the broadened spin ensemble was characterised by a q-Gaussian spectral line shape $\rho(\omega)$ with a linewidth of $\gamma_{\mathrm{inh}}=$4.55\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}.

**Spin system

**The negatively charged $\mathrm{NV}^{-}$ centre is a paramagnetic point-defect centre in diamond with an electron spin $S=1$, consisting of a nitrogen atom replacing a carbon atom in the diamond lattice and an adjacent vacancy. The ground spin triplet can be described by a simplified Hamiltonian $H/h=DS_{z}^{2}+\mu_{B}S_{z}$, with $\mu=$28\text{,}\mathrm{M}\mathrm{H}\mathrm{z}\mathrm{/}\mathrm{m}\mathrm{T} and a large zero-field splitting of $D/2\pi=$2.87\text{\,}\mathrm{G}\mathrm{H}\mathrm{z}. The splitting corresponds to a temperature of $D/hk_{b}=$138\text{,}\mathrm{m}\mathrm{K}, which allows to thermally polarise the spins in the ground state at the fridge base temperature of $25\text{\,}\mathrm{m}\mathrm{K}$ with up to 99% fidelity. Due to its crystallographic diamond structure four different NV subensembles, with equal abundance, pointing in the $[1,1,1]$ direction, exist. By applying $B\approx$30\text{\,}\mathrm{m}\mathrm{T} either with 0°or 45°relative to the $[1,0,0]$ direction in the NV-resonator plane, we can Zeeman tune four or two NV subensembles into resonance with the cavity mode.

**Superconducting resonator

**The microwave cavity is loaded by placing the diamond sample on top of a $\lambda/2$-transmission line resonator. The superconducting microwave cavity is fabricated by optical lithography and reactive-ion etching of a $200\text{\,}\mathrm{n}\mathrm{m}$-thick niobium film sputtered on a $330\text{\,}\mathrm{μ}\mathrm{m}$-thick sapphire substrate. The loaded chip is hosted and bonded to a printed circuit board enclosed in a copper sarcophagus and connected to microwave transmission lines. The cavity exhibits a linewidth of $\kappa/2\pi=$440\text{,}\mathrm{k}\mathrm{H}\mathrm{z} which we can increase to $\kappa/2\pi=$1.2\text{\,}\mathrm{M}\mathrm{H}\mathrm{z} by applying weak magnetic fields perpendicular to the resonator plane that partially quenches the superconducting material. The cavity is coupled to the environment such that the internal losses of the cavity are much smaller than the coupling losses ($\kappa_{\mathrm{int}}\ll\kappa_{\mathrm{ext}}$). This allows to approximate the total losses as $\kappa\approx\kappa_{\mathrm{ext}}$.

**Transmission measurements

**Transmission measurements are performed by recording the forward scattering parameter $|S_{21}|^{2}$ of the hybrid system using a standard vector network analyzer (Agilent E5071C). To perform steady state bistability measurements we probe the system first with increasing power levels. For each power level we monitor the transition until a steady state is reached. We increase the drive power until the steady state lies far in the high driving branch, after which the power is lowered again in a stepwise manner. We identify bistability if the system shows different steady states when driven with increasing and decreasing power. From Eq. (3), we immediately see that the asymptotic behaviour of the transmission intensity $|\mathrm{T}|^{2}=|a|^{2}\kappa^{2}/\eta^{2}$ is described by $|\mathrm{T_{\mathrm{low}}}|^{2}=\frac{1}{(1+C_{\mathrm{coll}})^{2}}$ in the single excitation regime and by $|\mathrm{T_{\mathrm{high}}}|^{2}=1$ in the large excitation regime. The stroke between these two regimes is therefore only determined by the collective cooperativity $C_{\mathrm{coll}}$ as defined in the main text.

**Quench dynamic measurements

**Quench dynamics measurements are used to measure the temporal behaviour of the observed effect. For this we initialise the system in an initial state far in the large driving regime with a strong drive for several minutes. After this state preparation where the spin system is completely decoupled from the cavity, the drive power is non adiabatically switched to a lower drive, with transmission monitored for different target drive levels. We monitor the transmission until the time derivative of the transmission amplitude becomes smaller than an arbitrarily chosen threshold, which we then identify as our steady state.

**Transversal decay rate

**We employ Car-Purcell-Meiboom-Gill (CPMG)-like sequences to get an estimate for the spin-spin relaxation time ($T_{2}=1/\gamma_{\perp}$). The best achievable echo times in our experiment are $T_{2}=$4.8\pm 1.6\text{,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$$, which we identify as a lower bound for our relaxation times. The real spin-spin relaxation times are potentially longer, but misalignment of the external d.c. magnetic field with respect to the NV axis and a bath of excess electron and nuclear spins in the host material limits the echo time to times shorter than the real relaxation times.

**Longitudinal decay rate

**To get a value for $\gamma_{\parallel}$ we perform time resolved quench dynamic measurement in a regime with no bistability visible. For the limit $t\gg 1/\gamma_{\perp},1/\Theta_{j},1/\kappa$ the decay of the cavity intensity follows a single exponential decay $\mathrm{e}^{-\alpha\gamma_{\parallel}t}$ with different values $\alpha$ for different drive intensities $\mathrm{P}_{\mathrm{in}}$. By extracting these asymptotical decay constants $\Gamma=\alpha\gamma_{\parallel}$ from the experimental data, and fitting them to the numerically extracted values, we get a value for $\gamma_{\parallel}/2\pi\approx$6.25\text{\times}{10}^{-4}\text{,}\mathrm{Hz} ($T_{1}\approx$30\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}).

**Theoretical modelling

**To calculate the quench dynamics displayed in Fig. 4, we numerically solve the Maxwell-Bloch equations (Eqs. (Dynamical Exploration of Amplitude Bistability in Engineered Quantum Systems)) for the driving signals chosen nearby the first-order transition (see the main text for details) using the standard Runge-Kutta method. As an initial condition we take the steady state given by Eqs. (3), which lies on the upper branch depicted in Fig. 3c and corresponds to the limit of strong driving with $|\sigma_{j}^{z}|,\,|\sigma^{-}_{j}|\ll 1$. To accurately describe the dynamics and to achieve a good correspondence with experimental data, we take into account the effect of an inhomogeneous broadening by modelling the spin density with a $q$-Gaussian shape for the spin density, $\rho(\omega)=C\cdot\left[1-(1-q)(\omega-\omega_{s})^{2}/\Delta^{2}\right]^{1/(1-q)}$, distributed around the mean frequency $\omega_{s}/2\pi=2.6915$ GHz with the parameter $q=1.39$ and the width $\Delta/2\pi=5.3$ MHz. Such a shape for $\rho(\omega)$ was previously established by obtaining an excellent agreement between our theoretical model and the experiment, when treating the problem in the framework of the Volterra equation valid in the limit of weak driving signals krimer_non-markovian_2014. We then straightforwardly discretize our problem by performing the following transformation, $g_{j}=\Omega\cdot\left[\rho(\omega_{j})/\sum_{l}\rho(\omega_{l})\right]^{1/2}$, where $\Omega^{2}=\sum_{j}g_{j}^{2}$ stands for the collective coupling strength ($\Omega/2\pi=9.6$ MHz and $12.6$ MHz for Figs. 3ab and Fig. 3c, respectively). Since in total we deal with a sizable number of spins ($N\approx 10^{12}$), we make our problem numerically tractable by dividing spins into many subgroups with approximately the same coupling strengths, so that the numerical values for $g_{j}$ in Eq. (Dynamical Exploration of Amplitude Bistability in Engineered Quantum Systems) represent a coupling strength within each subgroup rather than an individual coupling strength.

**Adiabatic Elimination

**Using the fact that $\gamma_{\parallel}\ll\kappa,\gamma_{\perp},\Omega$, the dynamics at large times when $t\gg 1/\kappa,1/\gamma_{\perp},2\pi/\Omega$ can be considerably simplified as the cavity amplitude $a$ and the spin lowering expectation values $\sigma^{-}_{j}$ follow adiabatically the evolution of the $z$-component of the spin operator expectation value, $\sigma_{z}^{j}$. By introducing the small parameter $\epsilon=\gamma_{\parallel}/\gamma_{\perp}$ and the slow dimensionless time $\tau=\gamma_{\parallel}t$, we finally derive a first order differential Eq. (4) for the intra-cavity intensity $|a|^{2}$.