\section{Introduction}
\label{sec:intro}

Objective audio quality measurement systems are tools for evaluating or monitoring the perceived quality of processed music or speech signals by means of signal analysis and feature extraction. These systems are expected to save time and resources by predicting results of subjective listening tests, which are considered the standard method for estimating quality.

The majority of the current objective quality systems incorporate models of human auditory perception, mostly peripheral effects like auditory masking, motivated by their successful use in the field of audio coding \cite{RixTaslp}. Quality degradation is a multidimensional phenomenon, and it can be driven by multiple factors such as perceived roughness, noisiness and linear distortions ---among others--- in the processed signals \cite{erne2001perceptual, dick2017generation}. In many popular subjective quality assessment methods, quality degradation in multiple dimensions is rated using a single mean opinion score (MOS) \cite{MUSHRA,ITUP800}. Likewise, objective perceptual quality systems map different distortion measures (DM) derived from the perceptual model into one single quality score by means of multivariate regression, using subjective MOS as training data. The combination of different DM into an objective quality score is considered to be a basic data-driven (and task-dependent) model of cognition \cite{ThiedePHD}. The relative importance of each of the factors influencing perceived quality degradation is related to distortion \textit{salience}, and is driven by cognitive effects \cite{ThiedePHD,bregman1994auditory}. Salience influences the perceived severity of an artifact in an audio signal and is usually implicitly modeled by the mapping stage, although it has also been explicitly modeled \cite{vincent2012improved}.

Cognitive models have been considered in objective perceptual quality measurement systems  to model the salience of a DM \cite{beerends1996the,Par2056_2019}. The cognitive model outputs (here, cognitive effect metrics - CEM) can interact with DM by scaling their measured intensity. Systems with cognitive-corrected DM have shown inconsistent performance \cite{ThiedePHD}, since salience rarely depends on a single cognitive effect \cite{sporer1997objective}.
Ideally, performing multivariate statistical learning with the DM and CEM \cite{barbedo2005a} as direct inputs can lead to a more complete model of salience. The learning models may find meaningful CEM/DM interactions describing overall quality given enough training data available. However, training data is usually costly and scarce \cite{DelgadoPEAQ}, and learning models may not have enough data to perform DM-to-quality-score mapping and determine meaningful CEM/DM interactions at the same time, causing the model to overfit \cite{ThiedePHD}.

Based on the assumption that particular cognitive effects influence distortion salience, we propose a novel cognitive salience model (CSM) as an extension to the multivariate statistical learning approach. While also data-driven, the CSM separates the quality mapping from the CEM/DM interaction model. The CEM/DM interactions are informed and optimized using a salience cost function, and not the cost function used by the multivariate regression mapping. The training data in this stage only informs a small number of CEM/DM interaction model parameters. In this way, the CSM model is expected to outperform multivariate learning models that use DM and CEM as direct inputs for the training dataset sizes typically available in practice.

\section{Method}
\label{sec:method}

\begin{figure}[htb]

  \makebox[0.8\textwidth][l]{
      \includegraphics[width=0.48\textwidth]{figures/data-driven}
  }

\caption{Block diagram of the proposed objective perceptual quality measurement system.}
\label{fig:block_diagram}
\end{figure}

As basis for our contribution, we propose a \textit{full-reference} \cite{RixTaslp} objective quality measurement system (\Figure{fig:block_diagram}) composed by a perceptual model (PM) which produces a series of perceptually-motivated DM and CEM, and a cognitive salience model  that maps said metrics into a single objective quality score. The objective quality score is expected to be a predictor of an overall MOS, as graded by listening test subjects.

\subsection{Perceptual Model}

The perceptual model used corresponds to our own implementation \cite{DelgadoPEAQ} of the one used by the Perceptual Evaluation Audio Quality (PEAQ) method \cite{PEAQ}, in its \textit{advanced} version. The perceptual model compares reference (REF) and signal under test (SUT) in the transformed psychophysical domain (time-pitch-loudness or modulation). Multiple comparisons using peripheral and cognitive models in the psychophysical domain produce the different DM.

\begin{table}[t]
  \centering
    \resizebox{\columnwidth}{!}{%
      \begin{tabular}{l|c|c}
        \textbf{DM} & \textbf{Associated MOV} &\textbf{Distortion Type} \\ \hline
        LinDist & AvgLinDist\_A & Linear distortions  \\
        ModDiff & RmsModDiff\_A & Modulation disturbances \\
        NoiseLoudness & RmsNoiseLoudAsym\_A & Added noise in SUT \\
        MissingComponents & RmsMissingComponents\_A & Missing components in SUT \\
        EHS  & EHS\_B & Harmonic structure of the error \\
        SegNMR  & Segmental NMR\_B & Noise-to-mask ratio \\
      \end{tabular}
    }
				\caption{Distortion metrics used and their related PEAQ MOVs (Tables 4 and 17 in \cite{PEAQ}).}
	\label{tab:distmetrics}
\end{table}

The distortion metrics produced by the perceptual model in PEAQ can be averaged over time and frequency to produce a single-number estimate per comparison, termed the Model Output Values (MOV). We perform frequency averaging and time averaging over 2-second segments on the distortion metrics, leaving the overall time averaging (AVG) over the whole length of the audio excerpt to a later stage in the CSM. The distortion metrics we used are based on corresponding MOVs, and are shown in \Table{tab:distmetrics}.

\begin{table}[t]
  \centering
    \resizebox{\columnwidth}{!}{%
      \begin{tabular}{l|c}
        \textbf{CEM} & \textbf{Description} \\ \hline
        EPN & Perceptual streaming measure from \cite{barbedo2005a}, Eq. (30)  \\
        PDEV & Informational masking from \cite{barbedo2005a}, Eq. (31) \\
        probSpeech & Probability of speech-like signal from \cite{SpeechVsMusicQuality} \\
      \end{tabular}
    }
				\caption{Cognitive effect metrics used in this work.}
	\label{tab:cognimetrics}
\end{table}

In addition to DM, CEM (\Table{tab:cognimetrics}) can be derived from the perceptual model. We consider a model for two important cognitive effects from the Auditory Stream Analysis field \cite{bregman1994auditory}: \textit{perceptual streaming} (PS) and \textit{informational masking} (IM). PS indicates whether some disturbance in the SUT is perceived as a separate auditory event from the masker signal itself (i.\,e.,\  REF). In this case, the disturbance is likely to be perceived as more annoying than when it forms one single percept with the masker signal. In contrast, IM increases the masking threshold when the signal and disturbance variations in time and frequency are large enough, due to an increased rate of information reaching the auditory cortex. PS is thought to weaken the effect of IM and they can be seen as competing effects \cite{beerends1996the}. We implement the measures of PS and IM as proposed in \cite{barbedo2005a}. In addition, speech and music quality perception are believed to make use of some slightly different processing mechanisms in the auditory nerve \cite{SpeechVsMusicQuality}. Therefore, the speech-music classifier proposed in \cite{USACSpeechClas} will be used as CEM. The classifier analyzes the signal in the early stages of the perceptual model, before the psychophysical representation transformations.

\subsection{Cognitive Salience Model}
\label{ref:CSM}
The proposed CSM combines the different DM into a single objective score via a weighted sum based on the outputs of different CEM. The weighting mechanism is expected to model distortion salience in quality perception.

\label{sec:BF}
In order to perform the weighted sum of DM into a single output score, the individual DM stemming from the perceptual model need to be mapped to a target subjective quality scale though a basis function (BF). The BFs were estimated using Multivariate Adaptive Regression Splines (MARS) \cite{Jekabsons_areslab}, mapping the associated MOVs in \Table{tab:distmetrics} to MOS (using the MUSHRA procedure \cite{MUSHRA}) in the database described in \cite{dick2017generation}, produced with isolated audio coding artifacts. The isolated audio coding artifacts minimize distortion type interaction, favoring the independence of each estimated BF.

As with many psychoacoustic results, cognitive effects modeled after controlled laboratory experiments with simple signals are not likely to completely generalize to overall quality perception in real use case signals \cite{ThiedePHD}. In real use cases, distortion interactions and experimental conditions in listening tests are most likely to mask the influence of small cognitive effect sizes in perceived quality, decreasing the \textit{probability} that subjects detect this effect. Likewise, perceived audio quality might no longer be influenced by larger measured effect sizes. For example, the IM effect has been reported to saturate with larger masker information amounts \cite{lutfiAndOh}.

We use psychometric functions \cite{kingdom2016psychophysics} ---in particular, the logistic function family--- to relate the CEM outputs and the actual weights to model the described phenomena. The CEM transformed by the psychometric functions are termed the \textit{detection probability weights} (DPW). Since a cognitive effect size can actually hinder the salience of a distortion type, an inverse operation (INV) in the detection probability domain can be defined as $1-DPW$ to address these cases.

\subsection{Optimization of Detection Probability Parameters}
\label{sec:optimization}

The interaction model parameters between CEM and MOVs will be estimated using a data-driven approach using listening test data. In our model, the interaction model parameters are specifically two parameters of the involved logistic functions: curve steepness and crossover midpoint. A third parameter, the curve's maximum value, will be therefore fixed at 1, emulating a detection probability. We carry out an exhaustive search optimization procedure for the two parameters in each of the DPW relating a CEM with a DM. The optimization cost function will be based on a salience measure.

The salience $\mathcal{S}$ for a DM $m,\mbox{ }m=1\ldots M$ in a signal $j, \mbox{ } j=1\ldots J$ is defined as the Pearson correlation coefficient between the respective DM basis function $BF_m$ and the MOS $y$ for all available audio treatments $i=1\ldots I$:

\begin{equation}
  \resizebox{0.9\columnwidth}{!}{
  $\mathcal{S}_m(j)

    = \frac{\sum_{i=1}^{I} ( y_{ij}-\overline{y}_{j} )( BF_{mij}-\overline{BF}_{mj} )}
  {\sqrt{\sum_{i=1}^{I} (y_{ij}-\overline{y}_{j})^2} \sqrt{\sum_{i=1}^{I} (BF_{mij}-\overline{BF}_{mj})^2}}$
  }.
\end{equation}
The salience metric is calculated for averaged values over time for the duration of the signals. The stronger the BF/MOS correlation is for a given DM in a signal, the more it is assumed that the DM describes quality degradation, and therefore the higher its salience.

The optimization cost function describes the covariance of the salience metric for a DM against the DPW output for said DM across all signals in the database. It is a metric of DM/CEM interaction. The optimization of the DPW parameters is expected to improve the interaction model. The correlation coefficient is defined as:

\begin{equation}
  \resizebox{0.9\columnwidth}{!}{
  $\mathcal{C}_m
    = \bigg|\frac{\sum_{j=1}^{J} ( S_m(j)-\overline{S}_{m} )( DPW_{m}(j)-\overline{DPW}_{m} )}
  {\sqrt{\sum_{j=1}^{J} (S_{m}(j)-\overline{S}_{m})^2} \sqrt{\sum_{j=1}^{J} (DPW_{m}(j)-\overline{DPW}_{m})^2}}\bigg|$
  }.
\end{equation}
The stronger the correlation, the better the DPW will predict a given DM salience over the signals in the database.

\section{Experiment Design}
\label{sec:Experiment}

We train and validate the proposed system with an extensive audio coding quality database of speech, music, and mixed signals at a sampling rate of 48 kHz under different treatments at bitrates ranging from 16kbps to 96kbps in stereo mode \cite{USACdatabase}. A total of 168 quality scores for 24 signals and 7 treatments for training, and a total of 216 pooled scores of 24 signals and 9 treatments (different than those of the training database) were used for validation. More than 25000 individual MUSHRA scores between the two used tests were used for the pooling.

We carry out the optimization procedure for the DPW on each CEM/BF interaction described in Section \ref{sec:optimization} using the Verification Test 2. We keep only the $DPW_m$ in the CSM that provided strong values of $\mathcal{C}_m$ after optimization, or CEM/BF interactions that have been established in psychoacoustic research. The system is then validated on the Verification Test 1 database along with some other variants.

\subsection{Evaluated System Variants}
\label{sec:proposed}

The main purpose of the proposed CSM system is to establish an interaction model between cognitive effects and distortion salience that will improve objective measurement performance. We propose three variants to account for these interactions. The first one will consider both the CEM and the DM as input variables of a mapping stage (emulating the approach used in \cite{barbedo2005a}) trained with the database used in Section \ref{sec:BF} and the training database of Section \ref{sec:Experiment}. The mapping stage consists of an Artificial Neural Network (ANN) with the same settings as listed for the advanced version of PEAQ \cite{RevisionBS1387}. The resulting model will be labeled $DM + CEM$. The second and the third variants correspond to the same implementation of the CSM of Section \ref{ref:CSM} but differing in the use of DPW for the DM/CEM interactions. The $PROPOSED$ system weights the BF directly by the CEM. The $PROPOSED\mbox{ }(Opt.)$ system weights the BF using the optimized DPW from Section \ref{sec:optimization}. The performance difference should indicate if there is any gain in modeling detection probability for cognitive effects. Additionally, the performance on the same validation data will be reported for two objective quality measurement systems: ViSQOL Audio \cite{Visqol_soft} and PEAQ Advanced DI (distortion index) \cite{RevisionBS1387}. PEAQ's DI output encompasses a broader quality range than its ODG output, the latter being quality-scale dependent \cite{ThiedePHD,DelgadoPEAQ}.

Overall performance is rated in terms of MOS/Objective score correlation $R$, and prediction error $RMSE^*$ as recommended in \cite{EvalObjective}, including techniques meant to reduce possible biases and quality gradients between the evaluated systems' training and validation databases.

\begin{figure}[htb]
  \makebox[\textwidth][l]{
    \resizebox{0.53\width}{!}{
      \input{figures/PAPER_1__plot_CEM_vs_MOV.tex}
    }
  }
\caption{Linear Distortion Salience, Mean Effect Size for a CEM (probSpeech) of the optimization database signals and associated DPW for the CEM/DM interaction.}
\label{fig:plot_CEM_vs_MOV}
\end{figure}

\section{RESULTS AND DISCUSSION}
\label{sec:results}

\begin{table*}[t!]
  \centering
    \resizebox{1.6\columnwidth}{!}{%
      \small
      \begin{tabular}{l|c|c|c|c}
        \textbf{Weight} & \textbf{CEM} &\textbf{Target DM} & \textbf{C (CEM/DPW)}& \textbf{Equation} \\ \hline
        DPW1 & probSpeech & LinDist            & -0.77 / -0.92 & DPW1 = 1- probSpeech\_th\_lin \\
        DPW2 & probSpeech & NoiseLoudness      & 0.67 / 0.80   & DPW2 = probSpeech\_th\_nl \\
        DPW3 & probSpeech & MissingComponents  & -0.20 / -0.37   & DPW3 = 1-DPW2 \\
        DPW4 & EPN           & LinDist       & -0.40 / -0.70 & DPW4 = 1-EPN\_th\_lin \\
        DPW5 & EPN           & SegmentalNMR       & 0.1 / 0.25 & DPW5 = (EPN\_th\_sgm)(1-PDEV\_th\_sgm) \\
        DPW5 & PDEV          & SegmentalNMR       & -0.18 / -0.21 & - \\
      \end{tabular}
    }
				\caption{Selected interactions and cost function values (CEM and DPW) for the proposed CSM.}
	\label{tab:resultsOptim}
\end{table*}

The interaction between the salience measure of the $LinDist$ DM, the $probSpeech$ CEM and its optimized DPW for the optimization database signals is shown in \Figure{fig:plot_CEM_vs_MOV}. A strong CEM/Salience covariance $(\mathcal{C}=-0.77)$ supports the assumption that CEM can predict a given DM importance in quality degradation. Linear distortions are more salient in music than in speech signals for the analyzed database. The optimized psychometric function improved the cognitive effect and distortion interaction $(\mathcal{C}=-0.92)$.

The results of the optimization procedure and DPW choice for the proposed models can be seen in Table \ref{tab:resultsOptim}. The $\_th$ suffix indicates the ``thresholded" CEM with the optimized logistic function considering the target DM salience. The remaining target DM that have not been listed are considered to contribute in equal parts to the final objective measure, with no DPW weighing.

Some conclusions consistent with previous reports can be drawn from Table \ref{tab:resultsOptim}. Considering DPW2 and DPW3, added noise is more salient for speech signals. The importance of added noise measures for speech has been confirmed by the objective speech quality literature \cite{PESQ}. Likewise, missing components are considered less annoying in this context. Coefficient DPW4 predicts that linear distortions are less important when perceptual streaming is likely taking place. In general, linear distortions are considered less annoying than added noise \cite{ThiedePHD}. More so if PS exacerbates the perception of added noise. Regarding DPW5, the implementation of PS/IM weights as competing effects (as in \cite{beerends1996the}) on an instantaneous measure of the noise-to-mask ratio improved performance of the model despite weak individual $\mathcal{C}$ values. Future work might contemplate how to better capture these higher order interactions between CEM and DM using data analysis.

\begin{table}[t!]
  \centering

      \begin{tabular}{l|c|c}
        \textbf{System} & $\mathbf{R}$ & $\mathbf{RMSE^*}$ \\ \hline
        ViSQOL NSIM & 0.82 & 5.6  \\
        PEAQ DI & 0.69 & 8.1  \\
        DM + CEM & 0.84 & 5.1 \\
        PROPOSED  & 0.86 & 4.6 \\
        PROPOSED (Opt.) & \textbf{0.90} & \textbf{3.7} \\
      \end{tabular}

				\caption{System validation performance metrics \cite{EvalObjective} for the systems proposed in Section \ref{sec:proposed}. Best performance metrics in bold.}
	\label{tab:resultsTotal}
\end{table}

\Table{tab:resultsTotal} shows the validation results for the proposed systems of Section \ref{sec:proposed}. The proposed CSM variant with the DPW outperforms all investigated systems for the data used. Note that the pure combination of CEM and DPW into an ANN ($DM+CEM$) did not perform as well as the proposed CSM. The ANN learning algorithm did not seem to find interactions as meaningful as is our perceptually-motivated and data-driven model.

\section{Conclusion and Future Work}

The overall inclusion of the presented cognitive effect metrics in objective perceptual quality measurement has proven to be promising. Systems including cognitive effects showed an improved performance with respect to baseline systems, encouraging the use of such metrics. The proposed CSM and the data-driven optimization of CEM/DM interactions seem to be a reasonable approach for incorporating cognitive metrics with the dataset sizes available in practice, outperforming the use of more traditional statistical learning techniques that consider CEM and DM as input features.

Future work should contemplate further validation in larger and more diverse databases, especially for use cases other than perceptual audio coding. Additionally, the inclusion of other cognitive effect models, interactions with other distortion metrics and using alternative salience measures may be further investigated.

\clearpage
\bibliographystyle{IEEEbib}


\section{Supporting information}
Schematic of microscope system, Images of fabricated devices, Room-temperature transmission spectra of GaP-on-diamond devices, Tuning curves on different devices, Resonant excitation measurements on both device and non-device SiV, Intensity dependence of SiV photoluminescence with 532\,nm and 700\,nm excitation, Theory and simulation of the SiV-cavity system, Fits and error approximation.

\section{Acknowledgements}
This material is based upon work supported by Department of Energy, Office of Science, National Quantum Information Science Research Centers, Co-design Center for Quantum Advantage (C2QA) under contract number DE-SC0012704 and National Science Foundation Grant No. ECCS-1807566. N.S.Y.\ was supported by the National Science Foundation Graduate Research Fellowship Program under Grant No.~DGE-2140004.  A.A.\ was supported by a Post Graduate Scholarship from the Natural Sciences and Engineering Research Council of Canada (PGSD3 - 545932 - 2020). D.H.\ was supported by the National Science Scholarship from A*STAR, Singapore. The photonic devices were fabricated at the Washington Nanofabrication Facility, a National Nanotechnology Coordinated Infrastructure (NNCI) site at the University of Washington which is supported in part by funds from the National Science Foundation (awards NNCI-2025489, 1542101, 1337840 and 0335765). The authors acknowledge the use of Princeton’s Imaging and Analysis Center, which is partially supported through the Princeton Center for Complex Materials (PCCM), a National Science Foundation (NSF)-MRSEC program (DMR-2011750) as well as the Princeton Micro-Nano Fabrication Lab.

\section{Optical measurements}

\begin{figure*}[ht]
\centering
\includegraphics[width=5in]{figure_si_1.pdf}
\caption{Schematic of the optical measurement setup. The GaP-on-diamond devices are placed inside a 10\,K closed cycle He cryostat.}
\label{fig:SI_1}
\end{figure*}

\newpage
\section{Fabricated devices}

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{figure_si_2.pdf}
\caption{(a) A 1.5\,mm$\times$1.5\,mm GaP membrane (300\,nm thick) transferred to a silicon oxide carrier substrate. The photonic devices are aligned to the membrane and patterned with e-beam lithography using HSQ negative-resist. (b) After etching the GaP using (Ar/Cl$_2$/N$_2$) ICP-RIE plasma processing. (c) Devices after stamp transfer to diamond (after pre-transfer characterization and GaP thinning for fine tuning cavity resonances). The stamping yield is low, however many high-Q devices are measured (indicated by the colored boxes). (d) A profilometer scan across the transferred devices showing the average final GaP thickness of 225\,nm. (e) SEM image of an array of PhC devices. (f) An angled SEM image of the PhC and grating coupler. (g) A close up SEM of the PhC holes. The final dimensions of the PhC structure are within 5\,nm of the optimized design.}
\label{fig:SI_2}
\end{figure*}

\newpage
\section{GaP-on-diamond PhC transmission}

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{figure_si_3.pdf}
\caption{Measured cavity resonances of ten GaP-on-diamond PhC devices at room temperature. The devices are excited at the cavity with 500\,\textmu W, 532\,nm CW laser and the coupled SiV PL is measured at the output grating on a spectrometer. The resonance center wavelength and $Q$ factor are extracted with a Lorentzian fit. Single SiV centers coupled to Device-3 are discussed in the main text. The cavity resonances blue-shift by 5.5\,nm after cooling down to 10\,K.}
\label{fig:SI_3}
\end{figure*}

\newpage
\section{Measured PhC Xenon tuning curves}
\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{figure_si_4.pdf}
\caption{(a--d) Show a sequence of spectra captured during the Xe gas deposition process (at 10\,K) for four cavity devices with resonances from 737\,nm (a) to 744\,nm (d). Coupled SiV transitions can be observed throughout the tuning range as the cavity resonances are gradually red-shifted.}
\label{fig:SI_5}
\end{figure*}

\newpage
\section{Resonant excitation of SiV centers}
Photoluminescence excitation (PLE) measurements of SiV centers (both device-coupled and non-device-coupled) were taken with a continuous-wave scanning laser (Sacher Lasertechnik).
We spectrally filter the SiV phonon-sideband using a combination of 760-nm long-pass and 785/62-nm band-pass filters and collect the photoluminescence on an avalanche photodiode (Excelitas).
The output data stream of the detector is temporally correlated with the laser scan signal and an additional analog power meter output which samples the laser.
This power meter monitors the intrinsic variation of power over the scan range and can signify mode hopping of the laser during individual scans as identified by large discontinuous jumps in laser power.
For all PLE data shown in this work, the laser scan range is chosen to avoid significant variation of the power and mode hops during the course of the scan and so this calibration data is not used explicitly.

\subsection{Device-coupled SiV centers}
As mentioned in the main text, PLE measurements taken on the devices can be done with either excitation through the waveguide or from the top of the cavity.
Exciting through the waveguide is ideal to spatially filter only SiV with some degree of cavity coupling.
PLE scans on the device-coupled SiV (P-1 and P-2) taken while the cavity was on resonance are shown in Fig. \ref{fig:device_ple_1}.

\begin{figure*}[ht]
\centering
\includegraphics[width=0.6\textwidth]{figure_si_5.pdf}
\caption{\textbf{\textit{(PLE measurements of device-coupled SiV)}} Resonant excitation scans of (a) P-1 and (b) P-2 taken with excitation through the grating and collection from above the cavity.
The integrated counts for the entire scan are shown below, however there is significant spectral diffusion between subsequent scans.
The structure of P-1 matches that of Fig. 5e,f in the main text indicating that the same SiV lines could be identified despite exciting from the top of the cavity.
In both cases, the PLE data shows that there are actually several lines in each set.
Specifically in the case of P-2 it is difficult to distinguish which lines contribute to the measured PL.
}
\label{fig:device_ple_1}
\end{figure*}

Because of the well-resolved transitions of P-1, we performed PLE while continuously tuning the cavity to observe a change in the linewidth.
The entire PLE measurement is shown in Fig. \ref{fig:device_ple_2}.
These scans were taken over two seconds each at 10-second intervals.
Over the course of the measurement there were a number of instances where SiV line disappeared which was attributed to ionization.
A green laser pulse was manually triggered with the software after scans denoted with the green marker in order to repump the SiV charge state.
The indicated intervals correspond to the integration regions of Fig.~5g,h in the main text.
At the end of the sequence the SiV line disappeared which is attributed to the degradation in imaging quality as a result of xenon gas condensation.
The tuning rate in this experiment was significantly faster than the the rate in the PL.
At the end of the tuning, resonance with P-1 was confirmed via spectra.

\begin{figure*}[ht]
\centering
\includegraphics[width=0.8\textwidth]{figure_si_6.pdf}
\caption{\textbf{\textit{(Continuous PLE measurement on ZPL P-1 while tuning the cavity)}} (a) The full series of scans taken while the cavity is tuned on to resonance with P-1 as described in the main text.
The green blocks in the bar above the plot indicates scans after which the green laser was pulsed to repump the charge state.
In scans where the bar is white there was no repump.
The indicated intervals of scans show the regions which were integrated to generate Fig. 5g,h in the main text.
Cavity tuning occurs from the indicated scan (scan \#78) with the preceding scans having larger background due to usage of a flashlight in the room while opening the Xe gas line.
The drift in frequency with additional Xe tuning is of unknown origin.
(b) Replicates the data binning as shown in Fig. 5g,h of the main text.
}
\label{fig:device_ple_2}
\end{figure*}

\subsection{Non-device SiV centers}
Additional PLE measurements were taken on implanted SiV that were located around 10 micrometers away from any transferred GaP.
Due to the relatively high density of SiV in the sample, it was generally possible to locate sharp SiV PLE lines within 10\,GHz detuning from the main ensemble.
Additionally, single lines of highly strained SiV with significant detuning could be located spectrally and subsequently scanned over in PLE.
As discussed in the main text, the excitation power which maximized the signal while also minimizing the ionization was generally variable between SiV, ranging from 10--100\,\textmu W.

In general, the SiV charge-state stability varied between different SiV and often would be unstable for individual SiV.
We utilized a blue LED ($\sim200$\,nW) and/or green laser ($\sim200$\,\textmu W) to repump the charge state during measurements.
The blue LED is uncollimated (used for wide-field imaging) and so the actual power delivered to the SiV is significantly less than $200$\,nW.
Fig.~\ref{fig:bulk_ple_1} shows a number of different scans on a blue-detuned SiV located near 736.07\,nm utilizing different repump mechanisms.
It was consistently observed that both the continuous blue LED and interleaved green laser pulse repump schemes resulted in a significant increase to the spectral diffusion between scans relative to the intervals without any repump mechanism.
Despite the increase in spectral diffusion, the linewidth in any particular scan did not appear to change significantly when a repump mechanism was used.

\begin{figure*}[t]
\centering
\includegraphics[width=0.8\textwidth]{figure_si_7.pdf}
\caption{\textbf{\textit{(Charge state repumping mechanisms)}} PLE measurements taken on a blue-detuned SiV (near 736.07\,nm) using different repump mechanisms. As before, the bar above the PLE scan figure indicates the scans in which a particular repump method was used.
Here the green color indicates a green repump pulse was performed after that particular scan whereas the blue color indicates a continuous blue LED which is on throughout the indicated scans.
(a) Demonstration of both methods being able to repump the SiV charge state which had gone dark for a few scans.
(b) PLE scans taken with the green pulse occuring after each scan.
(c) PLE scans with the blue LED on for the entirety of the measurement.
The usage of either repump method repeatedly results in a noticable increase of spectral diffusion between subsequent scans.
}
\label{fig:bulk_ple_1}
\end{figure*}

Finally, Fig.~\ref{fig:bulk_ple_2} shows different sets of PLE scans taken on three different SiV with varying degrees of detuning.
We show here two cases of SiV with linewidths near the thermally broadened limit of $250$\,MHz (at 10\,K) and one case in which the SiV had a significantly broadened linwidth of around $500$\,MHz.
A range of spectral intensity, stability, and linewidths can be observed.

\begin{figure*}[ht]
\centering
\includegraphics[width=0.9\textwidth]{figure_si_8.pdf}
\caption{\textbf{\textit{(Linewidths of implanted non-device-coupled SiV)}} PLE scans over three different SiV transitions in bulk diamond (a) without and (b,c) with a continuous blue repump.
The indicated binning ranges are integrated and a Lorentzian fit is taken to determine the linewidth of the transition.
For (b), the fit was taken over the portion of data for which the curve is drawn in order to avoid noise from the additional peaks.
The linewidths of the SiV vary from near the thermally broadened limit of 250\,MHz for scans (a) and (b) to more than 500\,MHz in (c).
}
\label{fig:bulk_ple_2}
\end{figure*}

\newpage
\section{Intensity dependence of SiV with 532 nm and 700 nm excitation}

\begin{figure*}[ht]
\centering
\includegraphics[width=\textwidth]{figure_si_9.pdf}
\caption{Variation of SiV photo-luminescence intensity under 532\,nm and 700\,nm excitation separately and simultaneously.  When excited simultaneously, the observed intensity is more than the sum of the intensities from each excitation separately. The 700\,nm excitation comes from a pulsed Tsunami laser while the 532\,nm excitation is the output of a continuous diode laser.}
\label{fig:SI_6}
\end{figure*}

Optically exciting defects at different wavelengths can reveal complex behavior resulting from that defect's excited state structure and charge exchange with nearby charge traps. Here we excite with 700\,nm excitation, which excites SiV from the ground state into the phonon side-band of the optically-active excited state. 532\,nm excitation excites the SiV from the ground state into a different, higher lying excited state \cite{haussler}. We observe that the SiV exhibit a non-linear increase in PL intensity when both lasers are on simultaneously, a clear signature of laser-induced deshelving of at least one meta-stable state. We suspect that the metastable states are likely the alternate charge states SiV$^{0}$, and SiV$^{2-}$ formed either from direct laser ionization or charge capture from nearby charge traps \cite{jayakumar_long-term_2020,dhomkar_-demand_2018}. Further understanding and engineering of the diamond surface, and SiV photo-dynamics should further stabilize and increase the SiV$^-$ PL intensity.

\newpage
\section{Modeling the PhC-SiV system}
\label{sec:model}
We follow the standard approach \cite{lukin2016notes} to model the SiV-cavity system as a two-level system with energy $\hbar\omega_a$ coupled to a cavity with resonance at $\omega_c$.
The Hamiltonian of the system is given by the Jaynes-Cummings Hamiltonian, here written in the interaction picture with constant-coefficients,
\begin{equation}
    H_{\mathrm{JC}} = \hbar\delta a^\dagger a
    -\hbar g \big( a^\dagger \dyad{g}{e}
    + a \dyad{e}{g} \big)
\end{equation}
where we define $\delta \equiv \omega_c - \omega_e$ as the detuning and $g$ as the single-photon Rabi frequency which characterizes the coupling strength.

We model dissipation from the SiV-cavity system via the open system formalism, taking the appropriate reservoir assumptions such that the evolution of the SiV-cavity system density operator $\rho$ is given by the Lindblad mater equation (at zero temperature)
\begin{equation}
    \dv{\rho}{t} = \frac{1}{i\hbar} [H_{JC},\rho]
    + \sum_j \Gamma_j \left(
    L_j \rho L_j^\dagger
    - \frac{1}{2} \{L_j^\dagger L_j, \rho\}
    \right)
\end{equation}
where $[\cdot,\cdot]$ ($\{\cdot,\cdot\}$) is the (anti-)commutator and $\Gamma_j$ is the rate of (population) dissipation via channel $j$ associated to the Lindblad operator $L_j$ defined on the SiV-cavity system Hilbert space.
For the system of interest the dissipation is described by tuples $(\Gamma_j, L_j)$ with
\begin{align*}
    (\gamma, \dyad{g}{e}) &\equiv \text{SiV \textit{decay} outside of the cavity (radiative and non-radiative),} \\
    (\gamma_d, \dyad{e}{e}) &\equiv \text{Pure dephasing of the SiV, and } \\
    (\kappa, a) &\equiv \text{Cavity loss (both radiation and coupling).}
\end{align*}
Note that these rates are given in units of angular frequency.

To model the system analytically we use the so-called stochastic wavefunction formalism.
Our aim is to describe the ensemble-average behavior and so we neglect quantum jumps and model the system using the effective Hamiltonian
\begin{equation}
    H_{\text{eff}} = H_{\text{JC}}
    - i\hbar \sum_j \frac{\Gamma_j}{2} L_j^\dagger L_j.
\end{equation}
This allows us to describe the dynamics of the system state $\ket{\psi}$ on a \textit{subspace} of the SiV-cavity system Hilbert space provided that all dissipative processes result in projections outside of the subspace under consideration (i.e. dissipation must not occur \textit{within} the subspace).
A necessary constraint, then, is that we must neglect pure dephasing $(\gamma_d, \dyad{e}{e})$ for this analysis.

\subsection{Modification of emitter lifetime}
Our aim is to describe the lifetime measurements (off-resonant, pulsed driving of the SiV).
In this case the decoherence due to pure dephasing can be ignored and the relevant subspace is spanned by the states $\{ \ket{e,0}, \ket{g,1} \}$ where $\ket{x,n}$ corresponds to the atom in state $x$ and $n$ photons in the cavity.
As required in the stochastic wavefunction formalism, the dissipation processes of interest, $(\gamma, \dyad{g}{e})$ and $(\kappa, a)$, result in projections into the state $\ket{g,0}$ which is outside of this subspace.

We consider the system in the so-called leaky-cavity regime distinguished by rates $\kappa \gg g, \gamma$, which is generally true of cavities with moderate to low quality factor $Q = \omega/\kappa$.
The effective Hamiltonian is given by
\begin{equation}
    H_{\text{eff}} = H_{\text{JC}}
    - i\hbar\frac{\gamma}{2} \dyad{e}
    - i\hbar\frac{\kappa}{2} a^\dagger a.
\end{equation}
The Schr\"odinger equation for a general state $\ket{\psi} = c_e\ket{e,0} + c_g \ket{g,1}$ yields the system of equations
\begin{align}
    \dv{c_g}{t} &= - \left(\frac{\kappa}{2}+i\delta\right) c_g
    - ig c_e , \\
    \dv{c_e}{t} &= -igc_g - \frac{\gamma}{2} c_e.
\end{align}
In the leaky-cavity regime, cavity losses dominate so the state $\ket{g,1}$ dissipates faster than the dynamics of other processes.
As a result, $c_g$ ``reaches equilibrium'' on timescales relevant to the coupling $g$ and SiV emission $\gamma$ rates so we may set $\dv*{c_g}{t} = 0$.
This process is known as adiabatic elimination \cite{lukin2016notes}.

The system of equations can then be solved with initial condition $\ket{\psi(t=0} = \ket{e,0}$ to find the population dynamics of the excited state as $|c_e(t)|^2 = \exp(-\gamma_{\text{eff}} t)$ with
\begin{equation} \label{eq:gamma_eff}
    \gamma_{\text{eff}} = \gamma
    + \frac{4g^2}{\kappa} \frac{(\kappa/2)^2}{(\kappa/2)^2 + \delta^2}.
\end{equation}
We may compare this effective decay rate to the intrinsic decay rate of the SiV dipole transition given by the Wigner-Weisskopf theory,
\begin{equation} \label{eq:wigner_weisskopf}
    \gamma_{eg} \equiv \frac{\omega_e^3 n |\mu|^2}{3\pi\epsilon_0 \hbar c^2},
\end{equation}
where $|\mu|$ is the amplitude of the dipole moment of the relative SiV transition, $\epsilon_0$ is the permittivity of free space, $c$ is the speed of light, and $n$ is the refractive index of diamond.
For an ideal emitter, when the cavity is \textit{open} such that $\gamma \approx \gamma_{eg}$ (i.e. the cavity does not suppress the emission outside of the cavity, which is the case for evanescently coupled emitters), Eq. \ref{eq:gamma_eff} shows that the cavity coupling $g$ results in an enhanced decay rate of the excited state.
This is the well-known Purcell effect \cite{purcell1946physrev}.

\subsection{Purcell enhancement}
The \textit{Purcell enhancement} is a figure of merit that quantifies the modification of the spontaneous decay rate of an emitter coupled to a cavity.
While commonly used in the literature, there exist several different definitions which obfuscate its utility as a cross-platform metric.
Specifically in the case of non-ideal emitters, the Purcell enhancement can be defined in ways that implicitly depend on the emitter properties if one does not account for additional loss mechanisms.
Here we specify a particular definition of the Purcell enhancement which quantifies the modification of the relevant transition rate and then clarify its relation to other figures of merit used in the literature.

We define the Purcell enhancement $F$ as the ratio of the cavity-enhanced and intrinsic (e.g. Eq. \ref{eq:wigner_weisskopf}) emission rates along the particular transition coupled to the cavity.
In general, $\gamma \neq \gamma_{eg}$, even in an open cavity, as $\gamma$ will include alternative decay channels.
Let $\gamma_{eg,\perp}$ be the contribution to $\gamma$ corresponding to decay through the relevant transition ($\ket{e}\to\ket{g}$) but with emission outside of the cavity.
The Purcell enhancement is then defined as
\begin{equation} \label{eq:purcell_enhance_def}
    F \equiv \frac{\gamma_{eg,\perp}}{\gamma_{eg}}
    + \frac{4g^2}{\kappa \gamma_{eg}} \frac{(\kappa/2)^2}{(\kappa/2)^2 + \delta^2}.
\end{equation}
The advantage of this definition is that it is essentially emitter independent --- the cavity parameters and emitter orientation/placement determine $F$.
We can see this by noting that
\begin{equation}
    g = \frac{|\mu|}{n}\sqrt{\frac{\omega}{4\pi\epsilon_0 \hbar V}} \left(\frac{\vb E \cdot \mu}{|\vb E_{\text{max}}||\mu|} \right)
\end{equation}
where $\vb E / |\vb E_{\text{max}}|$ is the vector electric field profile of the cavity mode normalized to a maximum of unity and $V$ is the corresponding mode volume.
Using this and Eq. \ref{eq:wigner_weisskopf} we can write the Purcell enhancement as
\begin{equation}
    F = \frac{\gamma_{eg,\perp}}{\gamma_{eg}}
    + F_P
    \left(\frac{\vb E \cdot \mu}{|\vb E_{\text{max}}||\mu|} \right)
    \frac{(\kappa/2)^2}{(\kappa/2)^2 + \delta^2}.
\end{equation}
where
\begin{equation}
    F_P \equiv \frac{3}{4\pi^2} \left(\frac{\lambda}{n}\right)^3 \frac{Q}{V}
\end{equation}
is the well-known \textit{Purcell factor} for a cavity with quality factor $Q=\omega_c / \kappa$.

From this analysis we see that the Purcell enhancement is a metric that describes the effect of the cavity on the emitter, encapsulating both the effectiveness of the coupling (spatially and spectrally) where as the Purcell factor is a figure of merit for the cavity itself, independent of any emitter.
These parameters are related with $F_P$ essentially giving the maximum Purcell enhancement that one could expect from a cavity.
From an engineering standpoint, it is clear that the Purcell enhancement $F$ is the correct metric for comparing designs given that optimization of both coupling and cavity properties is needed.

\subsection{Estimation of Purcell enhancement from experiment}
Estimation of the Purcell enhancement is challenging, especially when additional loss channels are significant.
Using time-resolved spectroscopy, one is able to resolve the temporal profile of emitted photons which follows the excited state population.
This allows one to determine the decay rate $\gamma_{\text{eff}}$ as a function of detuning $\delta$ as in Eq. \ref{eq:gamma_eff}, which is discussed in the main text as the excited state lifetime $\tau \equiv 1 / \gamma_{\text{eff}}$.

In general, $\gamma$ (and likewise $\gamma_{\text{eff}}$) contains contributions from a variety of different mechanisms including the aforementioned emission out of the cavity $\gamma_{eg,\perp}$, radiative transitions along other channels $\gamma_{\text{rad,other}}$ (e.g. via the phonon-sideband), and non-radiative transitions $\gamma_{\text{non-rad}}$.
These contributions are additive with $\gamma = \gamma_{eg,\perp}+\gamma_{\text{rad,other}}+\gamma_{\text{non-rad}}$.
From the analysis, the effect of the cavity is seen to be an additional decay channel with rate $(4g^2/\kappa)$ for $\delta=0$.
When $|\delta| \gg \kappa$ the cavity contribution is suppressed.

We consider the case of an open cavity ($\gamma_{eg,\perp} \approx \gamma_{eg}$) and look at the ratio of the experimentally measured lifetimes for the resonant ($\delta = 0$) and non-resonant ($|\delta| \to \infty$) cases which is given by
\begin{equation}
    \frac{\tau_{\text{off-res}}}{\tau_{\text{on-res}}}
    = \frac{\gamma + 4g^2 / k\kappa}{\gamma} = 1 + \frac{4g^2}{\kappa\gamma}.
\end{equation}
Comparing this to the definition of the Purcell enhancement (Eq. \ref{eq:purcell_enhance_def}) for an open cavity at zero detuning, $F=1+4g^2/\kappa\gamma_{eg}$, we find
\begin{equation}
    F = 1 + \frac{\gamma}{\gamma_{eg}} \left( \frac{\tau_{\text{off-res}}}{\tau_{\text{on-res}}} - 1 \right)
\end{equation}
reproducing the result of \cite{zhang2018strongly} with the correction of an additive constant of unity.

The factor $\gamma/\gamma_{eg}$ is the reciprocal of the fraction of decay via the dipole transition for an intrinsic SiV in bulk diamond.
For an ideal two-level system $\gamma/\gamma_{eg}=1$ and so the Purcell enhancement simply reduces to the ratio of the lifetimes $F = \tau_{\text{off-res}} / \tau_{\text{on-res}}$ which is the expected result.
In general however, $\gamma/\gamma_{eg} > 1$ corresponding to the fact that realistic emitters do not decay entirely through a single transition.
In such cases, estimation of this factor from the literature is necessary and can be accomplished given the branching ratio of the zero-phonon line (ZPL) $\xi\equiv \gamma_{eg} / \gamma_{\text{ZPL}}$, Debye-Waller factor $DW \equiv \gamma_{\text{ZPL}} / (\gamma_{eg} + \gamma_{\text{rad,other}})$, and quantum efficiency $QE\equiv (\gamma_{eg} + \gamma_{\text{rad,other}}) / \gamma$ with $\gamma_{\text{ZPL}}$ being the rate of emission through the ZPL. From these parameters one finds
\begin{equation}
    \frac{\gamma}{\gamma_{eg}} = \frac{1}{\xi \cdot DW \cdot QE}.
\end{equation}
This factor can then be used to estimate the Purcell enhancement as in the main text.

\subsection{Simulation of Purcell enhancement}
A useful feature of our definition of the Purcell enhancement, as in Eq. \ref{eq:purcell_enhance_def}, is that it can be directly simulated using finite-difference-time-domain (FDTD) simulation (Lumerical FDTD, MEEP).
Such methods typically involve approximating the emitter as a dipole antenna driven with a known power and then comparing this power to the actual power radiated in the simulation.
Since both the radiated power and physical emitter dipole transition rate depends on the local density of states, the relative enhancement of the radiated power is itself the Purcell enhancement as in Eq. \ref{eq:purcell_enhance_def}, provided the simulation remains in the leaky- and open-cavity limit.

The Purcell enhancement values as shown in Fig. 1f,g in the main text were computed using this simulation method in Lumerical FDTD.
For these simulations the parameters $(a_{\mathrm{mir}}, a_{\mathrm{cav}}, h_x, h_y, w) = (143, 133, 71, 183, 360)$ were used with the slight deviation of $a_{\mathrm{mir}}$ taken to handle meshing issues at lower resolutions.
The design under consideration utilized 6 mirror holes and 20 taper holes, matching the device considered in the main text.

The convergence of both the Purcell enhancement and cavity quality factor $Q$ as a function of simulation time and mesh size (via the ``mesh accuracy'' setting within Lumerical) is verified as shown in Fig.~\ref{fig:SI_9_1}a,b.
These simulations were performed with an ideally oriented dipole (axis orthogonal to the cavity) positioned 20\,nm below the GaP-diamond interface.
Incidentally, this also shows the potential of this device given proper orientation of the SiV --- here an increase in the Purcell enhancement by a factor of 4.
The simulations ultimately used in Fig. 1f-g were done with a mesh accuracy of 4 and simulation time of $20$\,ps in order to maintain high accuracy with practical computation time.

Typical Purcell enhancement and transmission spectra obtained from the simulation results are shown in Fig.~\ref{fig:SI_9_1}c,d.
The simulated $Q$ factors were determined by Lorentzian fits to the transmission spectra.
Significant oscillations can be observed in both the Purcell enhancement and transmission spectra for shorter simulation times which is characteristic of a sinc-lineshape resulting from the rectangular windowing of the finite simulation time.
Additionally, the wavelength position of the resonance converges with higher mesh accuracy due to meshing errors.
\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{figure_si_10.pdf}
\caption{\textbf{a--b.} Convergence of simulation parameters with increased mesh accuracy and simulation time. Simulations shown in the main text are done with a mesh accuracy of 4 and simulation time of 2 ps to optimize between accuracy and computation time. \textbf{c--d.} Example spectra corresponding to the indicated points in (\textbf{a--b}). Oscillatory behavior of the spectra at short simulation times is due to the finite bandwidth of the simulation. The resonance frequency shifts with improved mesh resolution. Purcell enhancement and Q are extracted from the spectra.}
\label{fig:SI_9_1}
\end{figure*}

\newpage
\section{Fits and error approximation}
Fit parameters were determined by a non-linear least-squares fitting algorithm as implemented in Python's SciPy module.
The fit parameters and corresponding covariance matrix are returned and used for subsequent calculations.
Error bounds are given by the standard deviations which are calculated as the square root of the respective diagonal element of the parameter covariance matrix.

\subsection{Lifetime measurements}
To fit the lifetime data we use an exponentially modified Gaussian (exGaussian) lineshape which is equivalent to a Gaussian impulse response (of the overall system) convolved with a single-sided exponential decay, i.e.
\begin{equation}
    f(t) = \exp(-\frac{(t-\mu)^2}{2\sigma^2}) * \Big(a \exp(-\Gamma t) H(t) \Big),
\end{equation}
where $\mu$ and $\sigma$ correspond to the center position (delay of system response after trigger) and width of the impulse response.
Likewise, $a$, and $\Gamma$ are the amplitude and decay rate of the exponential cavity-defect response.
The function $H(t)$ is the unit-step (Heaviside-step) function which is zero for $t<0$ and unity for $t>0$.
An explicit form for the exGaussian lineshape can be computed directly
\begin{equation}
    f(t) = a \int_0^\infty \dd{\tau} \exp(-\frac{(t - \tau-\mu)^2}{2\sigma^2}) \exp(-\Gamma \tau).
\end{equation}
Changing variables to $z = (\tau - (t - \mu - \sigma^2 \Gamma)) / \sqrt{2\sigma^2}$ gives
\begin{equation}
    f(t) = a \exp(-\Gamma\left( (t-\mu) - \frac{\sigma^2\Gamma}{2} \right)) \int_x^\infty (\sqrt{2\sigma^2}\dd{z}) e^{-z^2}
\end{equation}
where $x = -((t-\mu) - \sigma^2 \Gamma)/\sqrt{2\sigma^2}$.
The second term is identified as the complementary error function $\mathrm{erfc}(x) \equiv (2/\sqrt{\pi})\int_x^\infty \dd{z} e^{-z^2}$.
Thus, the final lineshape is given by
\begin{equation}
    f(t) = a\sigma \sqrt{\frac{\pi}{2}} \exp(-\Gamma\left( (t-\mu) - \frac{\sigma^2\Gamma}{2} \right) )
    \mathrm{erfc}\left( -\frac{(t-\mu) - \sigma^2 \Gamma}{\sqrt{2\sigma^2}} \right),
\end{equation}
which was used for fits of the lifetime data.
In fitting the actual data some additional modifications needed to be made.
Firstly, the width of the Gaussian input response ($\sigma$) is determined to be 228\,ps from a fit of the instrument impulse response which was measured using a reflection of the 2-ps pulse of the driving laser (which is thus far below the resolution of the timing board $\sim80$-ps).
Secondly, an additional free parameter $b$ is needed to account for the background (dark counts, light leakage, etc.) of the system such that the overall fitting function is given by $f^\prime(t) \equiv f(t) + b$.

As shown above, the original data fits the decay constant $\Gamma$ of the population decay which is given in units of inverse time.
The population lifetime is defined as $\tau \equiv 1/\Gamma$.
With the (standard-deviation) error of $\Gamma$, $\sigma_\Gamma$ computed from the self-variance of the fit as described above, the error in the lifetime is given from a Taylor expansion as
\begin{equation}
    \sigma_\tau = \abs{\frac{\sigma_\Gamma}{\Gamma^2}}.
\end{equation}
The lifetimes and corresponding errors are tabulated below in table \ref{tab:lifetime_fits_stddev}.

\begin{table}[ht]

    \caption{\label{tab:lifetime_fits_stddev}%
        Standard deviations of fitted lifetimes.
    }
    \begin{ruledtabular}
        \begin{tabular}{cccccc}
        \textrm{SiV}&
        \textrm{Meas. No.}&
        $\Gamma$ (ns$^{-1}$)&
        $\sigma_\Gamma$ (ns$^{-1}$)&
        $\tau$ (ns)&
        $\sigma_\tau$ (ns)
        \\

    \colrule

        \textrm{P-1}
         & 1  & 0.644 & 0.025 & 1.55 & 0.06 \\
         & 2  & 0.577 & 0.020 & 1.73 & 0.06 \\
         & 3  & 0.615 & 0.012 & 1.63 & 0.03 \\
         & 4  & 0.656 & 0.010 & 1.52 & 0.02 \\
         & 5  & 0.799 & 0.007 & 1.25 & 0.01 \\
         & 6  & 0.887 & 0.009 & 1.13 & 0.01 \\
         & 7  & 0.716 & 0.009 & 1.40 & 0.02 \\
         & 8  & 0.718 & 0.006 & 1.39 & 0.01 \\
         & 9  & 0.661 & 0.011 & 1.51 & 0.02 \\
         & 10 & 0.591 & 0.008 & 1.69 & 0.02 \\
         & 11 & 0.573 & 0.015 & 1.75 & 0.05 \\
         & 12 & 0.569 & 0.011 & 1.76 & 0.03 \\
         & 13 & 0.527 & 0.014 & 1.90 & 0.05 \\
         & 14 & 0.517 & 0.015 & 1.94 & 0.06 \\

    \colrule
        \textrm{P-2}
         & 1 & 0.675 & 0.023 & 1.48 & 0.05 \\
         & 2 & 0.846 & 0.013 & 1.18 & 0.02 \\
         & 3 & 1.034 & 0.013 & 0.97 & 0.01 \\
         & 4 & 1.384 & 0.015 & 0.72 & 0.01 \\
         & 5 & 2.010 & 0.022 & 0.50 & 0.01 \\
         & 6 & 2.030 & 0.034 & 0.49 & 0.01 \\
         & 7 & 1.703 & 0.022 & 0.59 & 0.01 \\
         & 8 & 1.753 & 0.032 & 0.57 & 0.01 \\
         & 9 & 1.816 & 0.040 & 0.55 & 0.01 \\

        \end{tabular}
    \end{ruledtabular}

\end{table}

To compute the propagated error of the lifetime ratio (such as $\tau_{\text{off-res}}/\tau_{\text{on-res}}$ as in equation (2) of the main text), we assume statistical independence of the lifetime decay constants (i.e. the covariance is $0$).
The error propagation of the lifetime ratio
\begin{equation}
    r = \frac{\tau_1}{\tau_2}
\end{equation}
can be computed from a Taylor expansion as
\begin{equation}
    \sigma_r = \abs{r}\sqrt{
    \left( \frac{\sigma_{\tau_1}}{\tau_1}\right)^2
    +
    \left( \frac{\sigma_{\tau_2}}{\tau_2}\right)^2
    }.
\end{equation}
The table \ref{tab:ratio_fit_stddev} shows the standard deviation of the lifetime parameter fits for SiV P-1 and P-2, as well as the corresponding error of the lifetime ratio.

\begin{table}[ht]

    \caption{\label{tab:ratio_fit_stddev}%
        Standard deviations of lifetime ratio parameters.
    }
    \begin{ruledtabular}
        \begin{tabular}{ccccccc}
        \textrm{SiV}&
        $\tau_{\text{off-res}}$ (ns)&
        $\sigma_{\tau,\text{off-res}}$ (ns)&
        $\tau_{\text{on-res}}$ (ns)&
        $\sigma_{\tau,\text{on-res}}$ (ns)&
        $r=\tau_{\text{off-res}}/\tau_{\text{on-res}}$&
        $\sigma_r$
        \\
        \colrule
        \textrm{P-1} &
        1.94 &
        0.06 &
        1.13 &
        0.01 &
        1.72 &
        0.05 \\
        \textrm{P-2} &
        1.48 &
        0.05 &
        0.50 &
        0.01 &
        2.98 &
        0.11 \\
        \end{tabular}
    \end{ruledtabular}

\end{table}

\subsection{Resonant excitation measurements: linewidth and cooperativity}
The linewidth of the SiV transitions are extracted from the resonant excitation measurements by fitting to a Lorentzian lineshape given by
\begin{equation}
    g(t) = \frac{(\Delta\nu/2)^2}{(\nu-\nu_0)^2 + (\Delta\nu/2)^2},
\end{equation}
where $\Delta\nu$ gives the full-width-half-max (FWHM) in dimensions of $\nu$ (here in units of Hertz; note that decay constants $\gamma$ from the Lindblad operators are given by $\gamma \equiv 2\pi\Delta\nu$).
As before, the error of the linewidth $\sigma_{\Delta\nu}$ is given by the square-root of the corresponding diagonal element from the fit covariance matrix.
The corresponding error bounds are directly referenced in Figure 5 of the main text.

In general, the measured linewidth from the resonant excitation measurements are broadened by spectral diffusion and pure dephasing.
We would like to compare this to the lifetime-limited linewidth which is determined by the Fourier transform of the excited state \textit{amplitude} decay.
Assuming a form of the excited-state amplitude decay $x(t) = \exp(-\Gamma t /2) H(t)$ (which corresponds to the population decay $|x(t)|^2$ with rate $\Gamma$), the Fourier transform in angular-frequency space is
\begin{equation}
    x(t)= \exp(-\Gamma t /2) H(t) \quad\xrightarrow{\mathrm{F.T.}}\quad \tilde x (\omega) = \frac{1}{(\Gamma/2) + i\omega}
\end{equation}
where $\omega \equiv 2\pi \nu$ is the angular frequency for an oscillation at frequency $\nu$ in Hertz.
The lineshape of a lifetime-limited transition follows the spectral density $|\tilde x (\omega)|^2$ which has Lorentzian profile with FWHM in $\omega$-space of $\Gamma$.
Correspondingly, the FWHM in $\nu$-space is then $\Gamma/2\pi$, or correspondingly $1/2\pi \tau$ in terms of the lifetime.
The associated error is then $\sigma_{\text{FWHM},\tau} = |\sigma_{\tau}/2\pi\tau^2|$.
Error bounds for lifetime-limited linewidth in the on-/off-resonant cases are shown below in table \ref{tab:linewidth_lifetime_limit_stddev}.

The cooperativity is given by $C = 4g^2/\kappa(\gamma + \gamma_d)$ where $\gamma_d$ describes the contribution to the loss via pure dephasing of the SiV linewidth (see Section \ref{sec:model} for additional details).
The measured linewidth in the resonant excitation measurements gives the broadened linewidth which we associate to the total loss $\gamma_{\text{tot}} \equiv \gamma + \gamma_d$ (in $\omega$-space) whose errors are tablulated in table \ref{tab:linewidth_lifetime_limit_stddev}.
As discussed in the main text, in the limit of negligible pure dephasing ($\gamma_d \to 0$), we have
\begin{equation}
    C_{\text{max}} = \frac{\tau_{\text{off-res}}}{\tau_{\text{on-res}}} - 1.
\end{equation}
As a result the error bounds on $C_{\text{max}}$ correspond to that of the lifetime ratios themselves.
In order to estimate the actual cooperativity given the measured linewidth we use the following relation
\begin{equation}
    C = \frac{\gamma}{\gamma_{\text{tot}}} C_{\text{max}}
\end{equation}
where $\gamma = 1/2\pi\tau_{\text{off-res}}$ is numerically given by the off-resonant lifetime-limited linewidth as tabulated in table \ref{tab:linewidth_lifetime_limit_stddev}.
We first compute the error of the linewidth ratio assuming statistical independence of each of the parameters such that the overall error can be computed
\begin{equation}
    \sigma_{\gamma/\gamma_{\text{tot}}} = \frac{\gamma}{\gamma_{\text{tot}}}
    \sqrt{\left( \frac{\sigma_{\gamma}}{\gamma} \right)^2 +
    \left( \frac{\sigma_{\gamma_{\text{tot}}}}{\gamma_{\text{tot}}} \right)^2}.
\end{equation}
The error associated with the cooperativity is thus
\begin{equation}
    \sigma_C = C \sqrt{
    \left( \frac{\sigma_{\gamma/\gamma_{\text{tot}}}}{\gamma/\gamma_{\text{tot}}} \right)^2 +
    \left( \frac{\sigma_{C_{\text{max}}}}{C_{\text{max}}} \right)^2}
\end{equation}
The calculated values are tabulated in table \ref{tab:linewidth_lifetime_limit_stddev}.

\begin{table}[ht]

    \caption{\label{tab:linewidth_lifetime_limit_stddev}%
        Standard deviations of parameters for cooperativity estimation. All frequencies (rates) are in units of MHz.
    }

        \begin{ruledtabular}
        \begin{tabular}{ccccccccccc}
        \textrm{SiV} &
        $\gamma/2\pi$ &
        $\sigma_{\gamma}/2\pi$ &
        $\gamma_{\text{tot}}/2\pi$ &
        $\sigma_{\gamma_\text{tot}}/2\pi$ &
        $\gamma/\gamma_{\text{tot}}$ &
        $\sigma_{\gamma/\gamma_{\text{tot}}}$ &
        $C_{\text{max}}$ &
        $\sigma_{C_{\text{max}}}$ &
        $C$ &
        $\sigma_{C}$ \\
        \colrule
        P-1 &
            82.22 &
            2.39 &
            578.81 &
            31.90 &
            0.140 &
            0.009 &
            0.72 &
            0.05 &
            0.10 &
            0.01
        \\
        P-2 &
            107.42 &
            3.64 &
            - &
            - &
            - &
            - &
            1.98 &
            0.10 &
            - &
            -
        \end{tabular}
    \end{ruledtabular}

\end{table}