## Abstract

Nonlinear systems, whose outputs are not directly proportional to their inputs, are well known to exhibit many interesting and important phenomena that have profoundly changed our technological landscape over the last 50 years. Recently, the ability to engineer quantum metamaterials through hybridization has allowed us to explore these nonlinear effects in systems with no natural analog. 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 centers. One of the exciting properties of this spin system is its long spin lifetime, which is many orders of magnitude longer than other relevant time scales of the hybrid system. This allows us to dynamically explore this nonlinear regime of cavity quantum electrodynamics and demonstrate a critical slowing down of the cavity population on the order of several tens of thousands of seconds—a time scale much longer than observed so far for this effect. Our results provide a foundation for future quantum technologies based on nonlinear phenomena.

## INTRODUCTION

In nature, most physical systems are inherently nonlinear, giving rise to effects, such as bistability (*1*), chaos (*2*), solitons (*3*), and superradiance (*4*), and often appear counterintuitive when contrasted with much simpler linear systems. Amplitude bistability, one of the basic nonlinear phenomena [commonly used in optical switches nowadays (*5*)], has been extensively investigated both theoretically (*6*–*10*) and experimentally (*11*, *12*). It occurs in any medium where strong nonlinearities in the interaction between a radiation field and a polarizable medium, such as spins, exist. The nonlinearity in these 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 de-excited, ordered ground state occurs (*13*). The coupled system switches between these two branches and shows a hysteresis depending on the history of the system.

The usual cavity quantum electrodynamics (cQED) demonstrations use atoms or trapped ions coupled to optical light fields to investigate these nonlinear effects, but the short atomic lifetimes have made it difficult to truly observe the temporal dynamics of amplitude bistability (*14*, *15*) and restricted previous studies to the steady-state behavior. In contrast, quantum engineering by hybridizing different physical systems and by exploiting their advantages (*16*–*21*) allows us to create systems with extremely long-lived emitters, making it possible to observe the bistable system during evolution. Here, 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 [negatively charged nitrogen-vacancy (NV^{−}) centers in diamond (*22*, *23*)]. This type of system has been studied extensively (*17*, *18*), but experiments to date have been mainly carried out in the linear regime. The present work goes beyond this linear regime and thus provides the foundation for further understanding of nonlinear physics in this type of solid-state hybrid quantum system. Moreover, the long lifetime of the spin system (*24*) allows us to study the temporal behavior of the presented effect, a regime experimentally just recently accessed (*25*, *26*).

## RESULTS

### Equations of motion from the Tavis-Cummings Hamiltonian

An ensemble of spins in a cavity is characterized by the three quantities: polarization, inversion, and the cavity amplitude. Their dynamics can be derived from the driven Tavis-Cummings Hamiltonian (*27*) for *N* spins under the rotating wave approximation as(1)with *a*^{†} and ** a** being the creation and annihilation operators for the cavity mode of frequency ω

_{c}, respectively, and σ

_{j}

^{z}, σ

_{j}

^{+}, and σ

_{j}

^{−}as the spin inversion, raising, and lowering operators, respectively, for the

*j*-th spin of frequency ω

_{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 η and frequency ω

_{p}.

Using a mean-field approximation, valid in the limit of large spin ensembles, 〈*a*^{†}σ^{−}〉 ≈ *a*^{†}σ^{−} (in the following, 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 (*28*) as(2)with cavity dissipation rate κ, transversal spin relaxation rate γ_{⊥} = 1/*T*_{2}, and longitudinal spin relaxation rate γ_{∥} = 1/*T*_{1}. The relaxation rates are ordered as κ *>* γ_{⊥} ≫ γ_{∥} such that the longitudinal decay of the spin inversion is by far the slowest process. Θ_{j} are the frequency detunings with respect to the ensemble central frequency to account for inhomogeneous broadening. Setting the time derivatives to zero, the steady state of this system can be written as(3)where the dimensionless parameter *C*_{j} = *g*_{j}^{2}/[κ γ_{⊥}(1 + Θ_{j}^{2}/γ_{⊥}^{2})] is the single-spin cooperativity. The collective system cooperativity is given accordingly by *C*_{coll} = ∑_{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 σ^{j−} variables (*29*), which results in a first-order differential equation for the intracavity intensity. For the giant *S*_{z} = ∑_{j} σ_{j}^{z} spin in resonance (Θ_{j} = 0) with the cavity mode, it can be written as(4)

For our typical system parameters and a strong enough coupling, 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 bifurcations (*30*) at which the transition between both branches occurs.

In contrast to amplitude bistability using a small number of emitters (*25*, *26*), in our case, the sizeable number of spins allows us to neglect quantum fluctuations, which means that a blinking between the two stable solutions does not occur on an experimentally accessible time scale.

### The hybrid quantum system

The hybrid system is shown in Fig. 1 and is composed of an electron spin ensemble formed by NV centers in diamond, loaded onto a superconducting λ/2 resonator. To thermally polarize the *N* ≈ 10^{12} electron spins to their ground state (≥99%), we put the system in a dilution refrigerator at 25 mK. Each electron spin has a zero-field splitting of *D*/2π ≈ 2.878 GHz and an average coupling rate of *g*_{0}/2π ≈ 12 Hz to the cavity mode. We estimate the transversal relaxation rate as γ_{⊥}/2π ≤ 33 kHz and the longitudinal relaxation rate as γ_{∥}/2π ≤ 3.6 mHz (see Materials and Methods) (*18*, *24*). The superconducting resonator has a cavity linewidth of κ/2π = (440 ± 10) kHz (half width at half maximum) with a fundamental resonance frequency at ω_{c}/2π = 2.691 GHz and a loaded quality factor of *Q* = 3300. An external microwave field with frequency ω_{p} is used to probe the hybrid system.

### Steady-state amplitude bistability

First, we search for the bistable behavior in the steady state by measuring the transmitted intensities through the cavity, defined by |T|^{2} = P_{out}/P_{in} as a function of the input drive intensity P_{in} ≈ η^{2}/κ and outgoing intensity P_{out} ≈ |*a*|^{2} κ. The drive power is raised in a stepwise manner, which is slow enough to allow the system to reach a steady state for each stimulus P_{in}. For small excitations, the intracavity intensity is not sufficient to saturate the spin ensemble (σ_{j}^{z} ≈ −1) and is thus given by |*a*|^{2} = η^{2}/[κ(1 + *C*_{coll})^{2}]. As the power level increases, the cavity field bleaches the spins (σ_{j}^{z} ≈ σ_{j}^{−} ≈ 0) such that the Rabi splitting vanishes and the spin system decouples from the cavity (Fig. 2). The intracavity intensity |*a*|^{2} = η^{2}/κ^{2} is that of an empty cavity from which spins are completely decoupled.

This nonlinear saturation behavior is a necessary precursor to the observation of amplitude bistability. However, whether 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 separation.

In contrast to a homogeneously broadened spin ensemble where analytic expressions for the critical cooperativity to observe bistability exist (*C*_{coll} ≥ 8) (*10*), the inhomogeneous broadening requires numerical solutions to determine the bistability threshold (see Materials and Methods). The finite width of the spin distribution markedly increases the required collective cooperativity for which bistability can be observed (*31*). In the present case of an inhomogenously broadened line with Γ/2π = 9.5 MHz (full width at half maximum), we predict a critical cooperativity of *C*_{coll} = 42.4. For lower values of the collective cooperativity, the intracavity intensity as a function of the input drive is a continuous function, whereas at the critical cooperativity, the system response becomes a step function and, at two critical drive values P^{u}_{crit} and P^{d}_{crit} (see Fig. 3, B and C), the system switches between these two branches undergoing a first-order phase transition.

In Fig. 3, we show steady-state bistability measurements for three cooperativity values *C*_{coll} = 18, 49, and 78. The lowest value *C*_{coll} = 18 does not show bistability (Fig. 3A), but increasing the cooperativity to *C*_{coll} = 49 (see Materials and Methods) allows us to observe the first signs of bistable behavior (Fig. 3B). This value is close to the expected value for the critical cooperativity of *C*_{coll} = 42.2 for our system parameters. Increasing the cooperativity further to *C*_{coll} ≈ 78, we observe amplitude bistability (Fig. 3C) within a 2-dB range. This steady-state bistability behavior is well reproduced by a full numerical simulation, with inhomogeneous broadening taken into account (dashed lines in Fig. 3, A to C).

### Quench dynamic measurements

Given this evidence of amplitude bistability, we focus next on the temporal behavior of the hybrid system using quench dynamic measurements. We start by preparing the spin ensemble in one of the two extremal states, either polarized in the ground state or completely saturated and decoupled from the cavity. These initial states are prepared by setting the cavity input power to P_{in} = 0 or P_{in} ≫ P^{d}_{crit} for several minutes, respectively. The drive power is then nonadiabatically switched to a different 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 target drive powers. When the system is driven close to the bifurcation point (P_{in} ≈ P^{d}_{crit}), the time scales needed to settle in a stationary state become as long as 4 × 10^{4} s, as depicted in Fig. 4.

The behavior can be linked to our model given in Eq. 4, which predicts that our system features two critical drive values P^{u}_{crit} and P^{d}_{crit} at which a saddle-node bifurcation occurs (*30*). For input powers between P^{u}_{crit} and P^{d}_{crit}, two attractors coexist, and the system evolves to one of these attractors, depending on whether they are approached from below or above the bistable region. These two attractors, one with polarized and ordered spins and the other one with unordered and saturated spins, are connected by a first-order phase transition if the system cooperativity is large enough (see Fig. 3C). When driven far away from the critical drive values either in the strongly or in the weakly driven branch, the system approaches a steady state on a characteristic time scale determined by the slowest decay rate in the system—given by the longitudinal decay γ_{∥} for our implementation. The system settles in a stationary state at which the external drive and dissipation are equal and opposite in effect. Close to the critical points P^{d}_{crit} and P^{u}_{crit}, the system becomes scale-invariant and is characterized by an infinite correlation time (*32*)—an effect referred to as “critical slowing down” (*15*, *33*). In the presented experiment, we deal with a saddle-node bifurcation, where the dynamics exhibits power law divergence close to the critical drive value.

This behavior is shown in Fig. 4 (A to C) where, close to the critical drive, the system evolves toward 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 the evolution until the system relaxes to a real steady state. The time it takes to go from the upper to the lower branch diverges close to the critical drive according to *t*_{switch} ≈ |P_{in} − P^{d}_{crit}|^{−α}, as shown in Fig. 4C, with α ≈ 1.20 ± 0.04. For the simplest case of a saddle-node bifurcation after a cubic function, the phase transition shows an algebraic divergence with a critical exponent α = 1. The more complicated set of equations in the present case changes this critical exponent to the larger value α ≈ 1.20 ± 0.04. The exact value depends on the precise structure of the so-called normal form (*34*), as given by Eq. 4 for the homogeneously broadened case. Comparing these experimental results with the full numerical solutions of Eq. 2, including inhomogeneous broadening, we observe excellent agreement (see Fig. 4A). For the quench from low-power levels to the high-power levels, we observe a similar behavior with the same critical exponent.

## DISCUSSION

In summary, 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 lifetime is much longer than other decay constants in the system. This regime allows us to study the temporal evolution of the phase transition explicitly, something experimentally difficult to achieve using the standard cQED implementations. We observe a critical slowing down of the cavity population on the order of 11 hours, a time scale several orders of magnitude longer than observed so far for this effect and many orders of magnitude longer than other time scales associated with the system. Our experiment provides a foundation for the exploration of additional nonlinear phenomena in quantum metamaterials and future quantum technologies that may arise from it. One of the possible applications is microwave isolators and diodes that make use of the fact that the transmission intensity is different depending on the history of the system. For a large enough value of the collective cooperativity, our system provides an isolation of more than four orders of magnitude for a given input power in the bistable regime. By reducing the number of emitters, the nonlinearity in the system provides a way to create nonclassical states, such as spin-squeezed states, which are impossible to realize in a purely linear system. This paves the way toward possible applications in high-sensitivity magnetic field sensing and quantum metrology.

## MATERIALS AND METHODS

### Sample

The spin system was realized by enhancing a type Ib high-pressure high-temperature diamond crystal containing an initial concentration of 200 parts per million (ppm) of nitrogen, with a natural abundance of ^{13}C nuclear isotopes. We achieved NV^{−} centers with a total density of ≈6 ppm by 50 hours of neutron irradiation with a fluence of 5 × 10^{17} cm^{−2} and by annealing the crystal for 3 hours at 900°C. Excess nitrogen P1 centers (*S* = 1/2), uncharged NV^{0} centers, and additional lattice stress are the main source of inhomogeneous broadening, which exceeds decoherence because of the naturally abundant 1.1% ^{13}C spin bath. The characteristics of the diamond crystal and NV ensemble were initially determined at room temperature using an optical confocal microscope.

### Spin system

The NV^{−} center is a paramagnetic point-defect center in the 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* = *ħDS*_{z}^{2} + *ħ*μ*B*_{z}*S*_{z}, with μ = 28 MHz/mT and a large zero-field splitting of *D*/2π = 2.877 GHz. The splitting corresponds to a temperature of *D*/*hk*_{b} = 138 mK, which allows to thermally polarize the spins in the ground state at the fridge base temperature of 25 mK with up to 99% fidelity. Because of its crystallographic diamond structure, four different NV subensembles, with equal abundance (pointing in the [1, 1, 1] direction), exist. By applying *B* ≈ 30 mT with either 0° or 45° relative to the [1, 0, 0] direction in the NV resonator plane, we could Zeeman-tune four or two NV subensembles into resonance with the cavity mode.

### Superconducting resonator

The microwave cavity was loaded by placing the diamond sample on top of a λ/2 transmission line resonator. The superconducting microwave cavity was fabricated by optical lithography and reactive ion etching of a 200-nm-thick niobium film sputtered on a 330-nm-thick sapphire substrate. The loaded chip was hosted and bonded to a printed circuit board enclosed in a copper sarcophagus and connected to microwave transmission lines. The cavity exhibited a linewidth of κ/2π = 440 kHz, which we could increase to κ/2π = 1.2 MHz by applying weak magnetic fields perpendicular to the resonator plane that partially quenched the superconducting material. The cavity was coupled to the environment such that the internal losses of the cavity were much smaller than the coupling losses (κ_{int} ≪ κ_{ext}). This allowed us to approximate the total losses as κ ≈ κ_{ext}.

### Transmission measurements

Transmission measurements were 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 probed the system first with increasing power levels. For each power level, we monitored the transition until a steady state was reached. We increased the drive power until the steady state lay far in the high driving branch, after which the power was lowered again in a stepwise manner. We identified bistability if the system showed different steady states when driven with increasing and decreasing power. From Eq. 3, we immediately found that the asymptotic behavior of the transmission intensity |T|^{2} = |*a*|^{2}κ^{2}/η^{2} is described by |T_{low}|^{2} = (1 + *C*_{coll})^{−2} in the single excitation regime and by |T_{high}|^{2} = 1 in the large excitation regime. The difference in the transmission between these two regimes is therefore only determined by the collective cooperativity *C*_{coll} = ∑_{j} *g*_{j}^{2}/[κ γ_{⊥}(1 + Θ_{j}^{2}/γ_{⊥}^{2})].

### Quench dynamic measurements

Quench dynamics measurements were used to measure the temporal behavior of the observed effect. For this, we initialized 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 was completely decoupled from the cavity, the drive power was nonadiabatically switched to a lower drive, with transmission monitored for different target drive levels. We monitored the transmission until the time derivative of the transmission amplitude became smaller than an arbitrarily chosen threshold, which we then identified as our steady state.

### Transversal decay rate

We used Car-Purcell-Meiboom-Gill–like sequences to get an estimate for the spin-spin relaxation time (*T*_{2} = 1/γ_{⊥}). The best achievable echo times in our experiment were *T*_{2} = (4.8 ± 1.6) μs, which we identified as a lower bound for our relaxation times. The real spin-spin relaxation times were potentially longer, but misalignment of the external dc magnetic field with respect to the NV^{−} axis and a bath of excess electron and nuclear spins in the host material limited the echo time to times shorter than the real relaxation times.

### Longitudinal decay rate

To get a value for γ_{∥}, we used the dispersive shift of the cavity mode coupled to a detuned spin system. To enter the dispersive regime, we detuned the spin system such that the detuning was much larger than the collective coupling strength Ω. The spin system acted as a refractive medium that shifted the resonance frequency if the spin system was polarized in the ground state. By applying a strong microwave tone, we excited a fraction of the NV^{−} ensemble, which led to a shift of the resonator frequency. We monitored this frequency shift over time, while the spin system relaxed back into its thermal equilibrium state with the characteristic rate γ_{∥}. This gave a lower bound for the longitudinal relaxation time of *T*_{1} = 44 s (*18*), justifying the adiabatic elimination technique. The real longitudinal relaxation times were potentially longer, but spin diffusion processes by spin-spin interaction between neighboring NV^{−} center spins and additional electron spins limited the measured relaxation times to times shorter than the intrinsic longitudinal relaxation times.

### Theoretical modeling

To calculate the quench dynamics displayed in Fig. 4, we numerically solved the Maxwell-Bloch equations (Eq. 2) 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 took the steady state given by Eq. 3, which lay on the upper branch depicted in Fig. 3C and corresponded to the limit of strong driving with |σ_{j}^{z}|, |σ_{j}^{−}| ≪ 1. To accurately describe the dynamics and to achieve a good correspondence with experimental data, we took into account the effect of an inhomogeneous broadening by modeling the spin density with a *q*-Gaussian shape for the spin density, ρ(ω) = ζ [1 − (1 − *q*) (ω − ω_{s})^{2}/Δ^{2}]^{1/(1 − q)}, distributed around the mean frequency ω_{s}/2π = 2.6915 GHz, with the parameter *q* = 1.39, the width Δ/2π = 5.3 MHz, and a normalization constant ζ. Such a shape for ρ(ω) 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 (*35*). We then straightforwardly discretized our problem by performing the transformation *g*_{j} = Ω[ρ(ω_{j})/∑_{l} ρ(ω_{l})]^{1/2}, where Ω^{2} = ∑_{j} *g*^{2}_{j} stands for the collective coupling strength [Ω/2π = 9.6 MHz (Fig. 3, A and B) and 12.6 MHz (Fig. 3C)]. Because, in total, we dealt with a sizable number of spins (*N* ≈ 10^{12}), we made 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. 2 represent a coupling strength within each subgroup rather than an individual coupling strength.

### Critical cooperativity

From Eq. 3, it was straightforward to derive a condition for the threshold of bistability, which was accompanied by a negative slope of the driving strength η as a function of the transmission amplitude |a|. It is given by the following inequality

For the simple case of a homogeneous spin ensemble, this condition could be solved analytically, giving the well-known threshold for bistability of *C*_{coll} > 8. In the case of inhomogeneous broadening, the required collective cooperativity to observe bistability increased markedly as compared to the homogeneous case (*31*). Numerical simulations show that, for the *q*-Gaussian spin distribution used in our manuscript, the effect of bistability can be observed for *C*_{coll} > 42.2, which corresponds to a collective coupling of Ω/2π = 8.86 MHz. The threshold for bistability depends not only on the width of the distribution but also on its specific shape. Changing from a *q*-Gaussian to a Gaussian or a Lorentzian spin distribution changes the threshold to *C*_{coll} > 40.8 or *C*_{coll} > 45.2, respectively (by changing only the coupling strength while keeping all other parameters constant).

### Adiabatic elimination

Using the fact that γ_{∥}≪ κ, γ_{⊥}, and Ω, the dynamics at large times (when *t* ≫ 1/κ, 1/γ_{⊥}, and 2π/Ω) could be considerably simplified because the cavity amplitude *a* and the spin lowering expectation values σ_{j}^{−} adiabatically follow the evolution of the *z* component of the spin operator expectation value σ_{j}^{z}. By introducing the small parameter ϵ = γ_{∥}/γ_{⊥} and the slow dimensionless time τ = γ_{∥}*t*, we finally derived a first-order differential Eq. 4 for the intracavity intensity |*a*|^{2}.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank H. Ritsch, M. Trupke, and A. Amo for the helpful discussions.

**Funding:**The experimental effort led by J.M. was supported by the Top-/Anschubfinanzierung grant of Technische Universität Wien. S.P., A.A., and T.A. acknowledge support from 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 from the Austrian Science Fund (FWF) through the Spezialforschungsbereich NextLite Project no. F49-P10. K.N. acknowledges support from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) through the 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. performed 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 modeling the experiment. K.N. and W.J.M. provided support for the implementation of the theoretical methods. A.A. wrote the manuscript to which all authors suggested improvements.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the authors.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).