Spin State Dynamics in a Bichromatic Microwave Field: Role of Bright and Dark States in coupling with Reservoir

TL;DR

This paper investigates how a bichromatic microwave field influences spin state dynamics, revealing complex resonance behaviors linked to bright and dark states and their interactions with the environment.

Contribution

It introduces a detailed analysis of spin state interactions under bichromatic excitation, highlighting the roles of bright and dark states in resonance behavior and decoherence.

Findings

Nested magnetic resonances with power-dependent widths

Field-induced stabilization of certain resonances

Light-induced narrowing of composite resonances

Abstract

Driving an open spin system by two strong, nearly degenerate fields enables addressing populations of individual spin states, characterisation of their interaction with thermal bath, and measurements of their relaxation/decoherence rates. With such addressing we observe nested magnetic resonances having nontrivial dependence on microwave field intensity: while the width of one of the resonances undergoes a strong power broadening, the other one exhibits a peculiar field-induced stabilization. We also observe light-induced narrowing of such composite resonances. The observations are explained by the dynamics of bright and dark superposition states and their interaction with reservoir.

Equations

$$\frac{d}{dt}\rho =-\frac{i}{\hslash }[H,\rho ]-(\Gamma \rho ),$$

$$\rho =\left(\begin{array}{cc}{\rho }_{11}& {\rho }_{10}\\ {\rho }_{01}& {\rho }_{00}\end{array}\right)$$

$$H=\frac{\hslash }{2}\left(\begin{array}{cc}-\omega & {\Omega }_1+{\Omega }_2{e}^{i\delta t}\\ {\Omega }_1+{\Omega }_2{e}^{-i\delta t}& \omega \end{array}\right),$$

$$(\Gamma \rho )=\left(\begin{array}{cc}{\gamma }_1({\rho }_{11}-{\rho }_{11}^0)& \Gamma {\rho }_{10}\\ \Gamma {\rho }_{01}& {\gamma }_0({\rho }_{00}-{\rho }_{00}^0)\end{array}\right),$$

$$\begin{array}{rl}{\dot{\rho }}_{11}& =-\gamma (1+ϵ)({\rho }_{11}-{\rho }_{11}^0)+\frac{i}{2}({\Omega }_1+{\Omega }_2{e}^{-i\delta t}){\rho }_{10}-\frac{i}{2}({\Omega }_1+{\Omega }_2{e}^{i\delta t}){\rho }_{01},\\ {\dot{\rho }}_{00}& =-\gamma (1-ϵ)({\rho }_{00}-{\rho }_{00}^0)-\frac{i}{2}({\Omega }_1+{\Omega }_2{e}^{-i\delta t}){\rho }_{10}+\frac{i}{2}({\Omega }_1+{\Omega }_2{e}^{i\delta t}){\rho }_{01},\\ {\dot{\rho }}_{10}& =-(\Gamma -i\omega ){\rho }_{10}+\frac{i}{2}({\Omega }_1+{\Omega }_2{e}^{i\delta t})({\rho }_{11}-{\rho }_{00}),\\ {\dot{\rho }}_{01}& =-(\Gamma +i\omega ){\rho }_{01}-\frac{i}{2}({\Omega }_1+{\Omega }_2{e}^{-i\delta t})({\rho }_{11}-{\rho }_{00}).\end{array}$$

$$\dot{\rho}_{11}+\dot{\rho}_{00}=-\gamma\Big{(}(1+\epsilon)(\rho_{11}-\rho_{11}^{0})+(1-\epsilon)(\rho_{00}-\rho_{00}^{0})\Big{)}.$$

$${\rho }_{11}+{\rho }_{00}+ϵ({\rho }_{11}-{\rho }_{00})={\rho }_{11}^0+{\rho }_{00}^0+ϵ({\rho }_{11}^0-{\rho }_{00}^0).$$

$$\left\{\begin{array}[]{rcl}\rho_{11}+\rho_{00}&=&1\\ \epsilon(\rho_{11}-\rho_{00})&=&\epsilon(\rho_{11}^{0}-\rho_{00}^{0})\end{array}\right.$$

$${\rho }_{\alpha \beta }(t)=\sum _{k\in \mathbb{Z}}{\tilde{\rho }}_{\alpha \beta }(k)e^{ik\delta t}\forall \alpha ,\beta =0,1.$$

$$\begin{split}\tilde{\rho}_{10}(k)&=\frac{i}{2}\cdot\frac{\Omega_{1}\Big{(}\tilde{\rho}_{11}(k)-\tilde{\rho}_{00}(k)\Big{)}+\Omega_{2}\Big{(}\tilde{\rho}_{11}(k-1)-\tilde{\rho}_{00}(k-1)\Big{)}}{\Gamma-i\omega+ik\delta},\\ \tilde{\rho}_{01}(k)&=-\frac{i}{2}\cdot\frac{\Omega_{1}\Big{(}\tilde{\rho}_{11}(k)-\tilde{\rho}_{00}(k)\Big{)}+\Omega_{2}\Big{(}\tilde{\rho}_{11}(k+1)-\tilde{\rho}_{00}(k+1)\Big{)}}{\Gamma+i\omega+ik\delta}.\end{split}$$

$$\begin{array}{r}\frac{ik\delta }{\gamma }{\tilde{n}}_1(k)=-(1+ϵ)({\tilde{n}}_1(k)-n_1^0{\delta }_{k0})-\frac{1}{2}\mathcal{F}(k),\\ \frac{ik\delta }{\gamma }{\tilde{n}}_0(k)=-(1-ϵ)({\tilde{n}}_0(k)-n_0^0{\delta }_{k0})+\frac{1}{2}\mathcal{F}(k),\end{array}$$

$$\begin{split}\mathcal{F}(k)=&\mathcal{K}(k)\Big{(}\tilde{n}_{1}(k)-\tilde{n}_{0}(k)\Big{)}\\ &+\frac{\mathcal{L}(k)}{2}\Big{(}\tilde{n}_{1}(k+1)+\tilde{n}_{1}(k-1)-\tilde{n}_{0}(k+1)-\tilde{n}_{0}(k-1)\Big{)}\\ &+\frac{\mathcal{M}(k)}{2}\Big{(}\tilde{n}_{1}(k+1)-\tilde{n}_{1}(k-1)-\tilde{n}_{0}(k+1)+\tilde{n}_{0}(k-1)\Big{)},\end{split}$$

$$\begin{array}{rl}\mathcal{K}(k)& =\frac{{\Omega }_1^2}{\gamma }\cdot \frac{\Gamma +ik\delta }{(\Gamma +ik\delta {)}^2+({\omega }_1-{\omega }_0{)}^2}+\frac{{\Omega }_2^2}{\gamma }\cdot \frac{\Gamma +ik\delta }{(\Gamma +ik\delta {)}^2+({\omega }_2-{\omega }_0{)}^2},\\ \mathcal{L}(k)& =\frac{{\Omega }_1{\Omega }_2}{\gamma }\cdot \left(\frac{\Gamma +ik\delta }{(\Gamma +ik\delta {)}^2+({\omega }_1-{\omega }_0{)}^2}+\frac{\Gamma +ik\delta }{(\Gamma +ik\delta {)}^2+({\omega }_2-{\omega }_0{)}^2}\right),\\ \mathcal{M}(k)& =\frac{{\Omega }_1{\Omega }_2}{\gamma }\cdot \left(\frac{-({\omega }_1-{\omega }_0)}{(\Gamma +ik\delta {)}^2+({\omega }_1-{\omega }_0{)}^2}+\frac{{\omega }_2-{\omega }_0}{(\Gamma +ik\delta {)}^2+({\omega }_2-{\omega }_0{)}^2}\right).\end{array}$$

$$\frac{\delta }{\Gamma }\simeq 0,$$

$${\tilde{n}}_i(n)\simeq 0,$$

$$\frac{{\omega }_i-{\omega }_0}{\Gamma }\simeq 0,$$

$${\Omega }_1={\Omega }_2=:\Omega .$$

$$\begin{split}\frac{d}{d(\gamma t)}n_{1}=-(1+\epsilon)(n_{1}-n_{1}^{0})-S\big{(}1+\cos\delta t\big{)}(n_{1}-n_{0}),\\ \frac{d}{d(\gamma t)}n_{0}=-(1-\epsilon)(n_{0}-n_{0}^{0})+S\big{(}1+\cos\delta t\big{)}(n_{1}-n_{0}),\end{split}$$

$$S=\frac{\Omega^{2}}{\gamma\Gamma}\Big{(}L_{\Gamma}(\omega_{1}-\omega_{0})+L_{\Gamma}(\omega_{2}-\omega_{0})\Big{)}$$

$$\dot{n}=-({\mathcal{L}}_0+{\mathcal{L}}_1(t))\cdot n+n^0,$$

$$\begin{array}{rl}n& =\left(\begin{array}{c}{n}_1\\ n_0\end{array}\right),{\mathcal{L}}_0=\left(\begin{array}{cc}1+S+ϵ& -S\\ -S& 1+S-ϵ\end{array}\right),\\ n^0& =\left(\begin{array}{c}(1+ϵ)n_1^0\\ (1-ϵ)n_0^0\end{array}\right),{\mathcal{L}}_1(t)=S\cos \delta t\cdot \left(\begin{array}{cc}1& -1\\ -1& 1\end{array}\right),\end{array}$$

$$n(t)=e^{-{\mathcal{L}}_0\cdot \gamma t}(n(0)-\bar{n})+\bar{n},$$

$$\bar{n}=\left(\begin{array}{c}{n}_1^0-\frac{S(1-ϵ)}{1+2S-{ϵ}^2}(n_1^0-n_0^0)\\ n_0^0+\frac{S(1+ϵ)}{1+2S-{ϵ}^2}(n_1^0-n_0^0)\end{array}\right).$$

$$U_R=\left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right),$$

$$\tan \theta =\sqrt{1+{\left(\frac{ϵ}{S}\right)}^2}-\frac{ϵ}{S}.$$

$${\lambda }_k=1+S-(-1{)}^k\sqrt{S^2+{ϵ}^2}\forall k=0,1.$$

$$\begin{array}{rl}{\eta }_1& =\cos \theta (n_1-{\bar{n}}_1)-\sin \theta (n_0-{\bar{n}}_0)\\ {\eta }_0& =\sin \theta (n_1-{\bar{n}}_1)+\cos \theta (n_0-{\bar{n}}_0),\end{array}$$

$$\eta (t):=U_R^{-1}(n(t)-\bar{n})=\left(\begin{array}{c}{e}^{-{\lambda }_1\cdot \gamma t}{\eta }_1(0)\\ e^{-{\lambda }_0\cdot \gamma t}{\eta }_0(0)\end{array}\right),$$

$${\mathcal{L}}_1^{\prime }(t)=U_R^{-1}{\mathcal{L}}_1(t)U_R=S\cos (\delta t)\cdot \left(1+\left(\begin{array}{cc}\sin 2\theta & -\cos 2\theta \\ -\cos 2\theta & -\sin 2\theta \end{array}\right)\right)$$

$${\mathcal{L}}_1^{\prime }(t)=\frac{S}{\sqrt{S^2+{ϵ}^2}}\cos (\delta t)\cdot \left(\begin{array}{cc}\sqrt{S^2+{ϵ}^2}+S& -ϵ\\ -ϵ& \sqrt{S^2+{ϵ}^2}-S\end{array}\right).$$

Full text

Spin State Dynamics in a Bichromatic Microwave Field: Role of Bright and Dark States in coupling with Reservoir – Supplemental Materials

Wojciech Gawlik

  

Piotr Olczykowski

  

Mariusz Mrózek

  

Adam M. Wojciechowski

Institute of Physics, Jagiellonian University, Łojasiewicza 11, 30-348 Kraków, Poland

I Introduction

The analysis is based on a simplified two-level system and is performed with the density matrix formalism. The quantum master equation yields:

[TABLE]

with the density matrix

[TABLE]

and the Hamiltonian defined in the frame rotating with frequency $\omega_{1}$:

[TABLE]

where $\Omega_{1},\,\Omega_{2}$ denote Rabi frequencies of two $MW$s which have their frequencies differing by $\delta=\omega_{1}-\omega_{2}$ and $\omega=\omega_{1}-\omega_{0}$ represents detuning of $\omega_{1}$ from the ODMR central frequency $\omega_{0}$.

The part of the evolution generator responsible for relaxation and decoherence is independent of the choice of the rotating frame and defined as:

[TABLE]

where $\rho_{kk}^{0}$ are equilibrium populations, $\gamma_{k}=\gamma(1-(-1)^{k}\epsilon)$ denote the rates with which populations $\rho_{kk}$ relax to equilibrium, and $\Gamma$ is the relaxation rate of coherences $\rho_{jk}$ (here and throughout the paper indices $j,\,k=0,\,1$). Moreover, we assume that both $\Gamma$ and $\gamma$ are nonzero. Master equation (Eq. S1) results then in the following set of equations:

[TABLE]

II Closed vs. open systems

The trace of the density matrix represents the sum of individual populations. Since the master equation describes the evolution of the density matrix, its trace represents the time derivative of total probability stored in the system. For a closed system total probability is conserved, so $\mathrm{Tr}\dot{\rho}=0$. On the contrary, the systems with $\mathrm{Tr}\dot{\rho}\neq 0$ are referred to as open ones since the probability flows in and out of the system from a reservoir of probability. To examine the conditions of openness of the system we compute the trace of Eq. S1:

[TABLE]

Conservation of the total probability requires that the left-hand side of the above equation nullifies, i.e. $\dot{\rho}_{00}+\dot{\rho}_{11}=0$, which results in:

[TABLE]

Without loss of generality, we can define a normalization condition: $\rho_{11}^{0}+\rho_{00}^{0}=1$. For $\epsilon=0$ there exists non-stationary solution for which the total probability is conserved with the normalization $\rho_{11}+\rho_{00}=1$. On the contrary, for $\epsilon\neq 0$ the set of equations:

[TABLE]

has only stationary solution with $\rho_{kk}=\rho_{kk}^{0}$. This solution, in turn, does not fulfill Equation S2 for $\Omega_{1}\neq 0$ or $\Omega_{2}\neq 0$ if only $\delta\neq 0$ and $\rho_{11}^{0}-\rho_{00}^{0}\neq 0$.

To summarize, equality of relaxation rates $\gamma_{0}=\gamma_{1}$ signifies a closed system with the normalization condition $\rho_{11}+\rho_{00}=\rho_{11}^{0}+\rho_{00}^{0}$, whereas inequality $\gamma_{0}\neq\gamma_{1}$ represented by the existence of asymmetry parameter $\epsilon\neq 0,$ requires opening of the system. In the open system probability is not conserved, flows in and out of the system from a reservoir of probability. Further, our analysis shows that the openness of the system under consideration can be detected only for two $MW$s oscillating with different frequencies $\omega_{1}\neq\omega_{0}$ and the system is prepared in such a manner that its equilibrium populations are not equal. i.e. $\rho_{11}^{0}\neq\rho_{00}^{0}$.

This can be explained in the following way, for a closed system ($\epsilon=0$) the solution is restricted to a one-dimensional space (a line) with $\rho_{11}+\rho_{00}=1$, on the contrary, for $\epsilon\neq 0$ the solution is not constrained by total probability, thus it penetrates the subset of a two-dimensional manifold parameterized by individual populations. Moreover, it is driven by $MW$ field in the environment of the stationary solution, hence it can be decomposed into two evolving states identified by different relaxation rates and coupled by $MW$ field. This very observation is sufficient to explain qualitatively the existence of two components of the resonance which reflects the dimensionality of the solution. As $\mathrm{Tr}([H,\rho])=0$, this effect does not depend on the specific form of interaction and is general.

III Adiabatic approximation

By application of the Floquet theory we obtain that the solution exponentially approaches periodic function with fundamental period $\frac{2\pi}{\delta}$. Assuming that the time of measurement is much longer that the inverse of the Floquet exponentials, we expand the solution into the Fourier series, i.e.:

[TABLE]

The last two equations of Eq. S2 in the Fourier basis yield:

[TABLE]

We use Eq. S3 to eliminate coherences from the first two equations S2 (for brevity we change the notation $n_{i}:=\rho_{ii}$):

[TABLE]

where:

[TABLE]

[TABLE]

Up to now, the calculations were exact, we will now make some simplifying assumptions:

• •

we consider effects which occur for small detunings only:

[TABLE]

• •

we assume that for very high harmonic numbers ($n\simeq\frac{\Gamma}{\delta}$) the amplitudes of the elements of the Fourier series are negligible:

[TABLE]

• •

we consider effects caused by burning a hole close to the center of the ODMR resonance:

[TABLE]

• •

both MWs are taken as equally strong:

[TABLE]

With these assumptions, the population dynamics can be cast in a simple form of first-order inhomogeneous equations (presented in the main paper as Eqs. 1):

[TABLE]

where:

[TABLE]

and we use a normalized Lorentz function $L_{a}(x)=\frac{a^{2}}{x^{2}+a^{2}}$.

IV Liouville equation, eigenvalues, dark and bright states

IV.1 Liouvillian expansion

We decompose the Liouvillian into a sum of a time independent $\mathcal{L}_{0}$ and time dependent $\mathcal{L}_{1}(t)$, i.e.:

[TABLE]

where:

[TABLE]

which become Eqs. 2 in the main paper. Operators $\mathcal{L}_{0}$ and $\mathcal{L}_{1}$ play significantly different roles in the dynamics of the system. The analysis of Eq. S6 is performed in three steps: first, to identify the states spanning the space of the solutions interpreted as dressed populations and causing the emergence of two components of the composite resonance; second, to designate their time-dependent coupling; and last, to obtain a quantitative prediction on the widths of the resonance. To this aim we interpret $\mathcal{L}_{0}$ as dressing initial states resulting in creation of two superpositions which we denote $\eta_{0}$ and $\eta_{1}$ and $\mathcal{L}_{1}$ as a time-dependent coupling between them.

IV.2 Time independent part

By zeroing $\mathcal{L}_{1}(t)$ in Eq. S6 we focus only on its time independent part. This part enables the definition of the reference system for analyzing the evolution of the full, time-dependent equation. A formal solution of Eq. S6 can be written as

[TABLE]

where point $\bar{n}\colon=\mathcal{L}_{0}^{-1}n^{0}$ denotes the stationary solution:

[TABLE]

For the sake of identification of independently evolving states and their lifetimes, we diagonalize $\mathcal{L}_{0}$ with the help of the $U_{R}$ matrix: $U_{R}^{-1}\mathcal{L}_{0}U_{R}=\mathrm{diag}(\lambda_{1},\lambda_{0})$, with $\lambda_{1},\,\lambda_{0}$ being the eigenvalues of $\mathcal{L}_{0}$. Straightforward computation yields:

[TABLE]

where we interpret $\theta$ as the mixing angle between initial populations $n_{1}$ i $n_{0}$ and $\theta$ satisfies the implicit equation:

[TABLE]

Eigenvalues of $\mathcal{L}_{0}$ are:

[TABLE]

The differences between the momentary populations and stationary ones, $n(t)-\bar{n}^{0}$, span the vector space. The eigenvectors of $\mathcal{L}_{0}$ resulting from the diagonalization can be interpreted as two independently evolving states $\eta_{1}$ i $\eta_{0}$ with coordinates $\begin{pmatrix}\eta_{1}\\ \eta_{0}\end{pmatrix}=U_{R}^{-1}(n-\bar{n})$, i.e.:

[TABLE]

which leads to the solution of the time-independent part of Eq. S6:

[TABLE]

where $\eta_{k}(0)$ are determined by initial conditions. This provides us with interpretation of $\eta_{k}$ as two dressed states with coefficients $\gamma\cdot\lambda_{k}=\gamma\cdot(1+S-(-1)^{k}\sqrt{S^{2}+\epsilon^{2}})$ as their relaxation rates. We interpret the modes as the bright and dark states, i.e. respectively short and long living.

It is important to realize that the calculated above states $\eta_{0,1}$ and corresponding rates $\lambda_{0,1}$ are very general features of mathematical properties of two coupled equations (Eq. S5). Their consequences are seen in a wide range of problems, for example in a collisionally perturbed atomic fluorescence Series (1978); Pendrill and Series (1978).

IV.3 Time dependent part of Equation S6 in the dark and bright states basis

In the previous section, the population evolution has been discussed under conditions when the evolution generator is independent on time. Such a situation is realized when there is just one driving field, like one $MW$. The analysis revealed existence of two states, i.e. linear combinations of the differences between populations and stationary solution, which tend to an equilibrium with two different rates $\gamma\cdot\lambda_{k}=\gamma\cdot(1+S-(-1)^{k}\sqrt{S^{2}+\epsilon^{2}})$ given by eigenvalues of $\gamma\cdot\mathcal{L}_{0}$. These states, derived above, are defined as the difference between initial and equilibrium values, so they can be added and multiplied and thus they span a vector space in the vicinity of stationary solution.

Here, we will present equations for a full, i.e. time-dependent evolution in the $\eta_{k}$ basis. The representation of the time-dependent Liouvillian in the this basis is:

[TABLE]

or, using the explicit form of the mixing angle $\theta$:

[TABLE]

The above representation of $\mathcal{L}^{\prime}_{1}$ yields

[TABLE]

where $\bar{\eta}_{1}$ and $\bar{\eta}_{0}$ are coordinates of $\bar{\eta}=U_{R}^{-1}\bar{n}$:

[TABLE]

Equations S9 demonstrate that states $\eta_{1}$, $\eta_{0}$ are mutually coupled by interaction with $MWs$ or, more precisely, by beating between two $MWs$, or $CPO$. Interestingly, when $S$ becomes large ($S\gg\Gamma$), the coupling affects the two states differently.

In the limit of very strong MWs ($S\rightarrow+\infty$) one obtains $\bar{\eta}_{1}\rightarrow 0$, $\bar{\eta}_{0}\rightarrow\frac{(n_{1}^{0}+n_{0}^{0})+\epsilon(n_{1}^{0}-n_{0}^{0})}{\sqrt{2}}$. Moreover, since $\tan\theta\xrightarrow{S\rightarrow+\infty}1$, the states $\eta_{1}$ and $\eta_{0}$ converges to the difference and sum of bare populations:

[TABLE]

We apply this result to the initial differential equations (Eqs. S5) which yields:

[TABLE]

The above equations induce following observations:

• •

Firstly, they show that strong bichromatic field couples directly to the population difference, and only indirectly to their sum.

• •

Secondly, in the first equation the relaxation rate in the self-coupling term ($\lambda_{1}$) grows asymptotically to $1+2S$ and dominates the cross-coupling term $\epsilon\,\eta_{0}$. Consequently, the difference of the initial populations exhibits a periodic CPO dynamics driven by oscillations proportional to $2S\cos\delta t$.

• •

Further, in the second equation the self-coupling constant converges to $1$ which is the limit of $\lambda_{0}$ for $S$ approaching infinity.

Which triggers several conclusions (considered in the three regimes of the MWs detuning):

• •

For $\delta\simeq\gamma$ both states exhibit a periodic CPO dynamics. The amplitude of $\eta_{1}$ is proportional to the source term $\frac{1-\epsilon^{2}}{\sqrt{2}}(n_{1}^{0}-n_{0}^{0})$. State $\eta_{0}$ is indirectly driven by force oscillating with the MWs detuning via the cross term $\epsilon\,\eta_{1}$. Still, its oscillations are insignificant compared to the amplitude of $\eta_{1}$ and restricted to the MWs detuning smaller than $\gamma$. This results in the composite resonance and justifies why the width of the resonance component corresponding to $\eta_{0}$ is limited from above by $\gamma$ and its amplitude decreases monotonically with $S$.

• •

For $\Gamma\gg\delta\gg\gamma$ the cross-coupling term $\epsilon\,\eta_{1}$ in the second equation oscillates with a frequency much higher than the relaxation rate determined by stabilization of the self-coupling constant $\lambda_{0}$ for strong MWs, hence it can be neglected. Since the dynamics of state $\eta_{0}$ is not driven by any source term its amplitude for $\delta\gg\gamma$ vanishes exponentially in time. The sum of populations becomes constant and the system effectively closes. Consequently, the cross term $\epsilon\,\eta_{0}$ in the first of Equations S11 becomes zero and states $\eta_{0}$ and $\eta_{1}$ can be regarded as decoupled for strong MWs.

• •

For $\delta\simeq\Gamma$ the driving term $2S\cos\delta t$ oscillates much faster than $1+2S$. Consequently, it has negligible effect on the dynamics of state $\eta_{1}$ which is governed solely by the source term $\frac{1-\epsilon^{2}}{\sqrt{2}}(n_{1}^{0}-n_{0}^{0})$. Thus, for a large detuning of the MWs the amplitude of $\eta_{1}$ becomes stationary at $\langle\eta_{1}\rangle\sim\frac{n_{1}^{0}-n_{0}^{0}}{1+2S}$ and one recovers the regular ODMR resonance shape.

Importantly, Eqs. S11 provide mathematical evidence and interpretation of the observations presented in the main paper - the effect of the power stabilization of the narrow resonance. This also justifies our interpretation of the uncoupled state as the bright one and the coupled state as the dark one.

The described effect is to some extent analogous to the effect of CPT. An important difference, though, is that here we consider combinations of populations, rather than wave functions or coherences in the familiar interpretation of the CPT effect Gray et al. (1978); Arimondo (1996).

V First harmonic approximation, resonance widths

By substitution of $n_{i}(t)=\alpha_{i}+\beta_{i}\cos(\delta t+\phi_{i})\,\,\forall i=0,1$ into Eqs. S5, ignoring the time dependence of all expressions except of constants and quantities oscillating as $e^{\pm i\delta t}$, we arrive at three pairs of equations for $\alpha_{i},\,\beta_{i},\,\phi_{i}$. As stated in the main article, we take the fluorescence intensity as proportional to the time-averaged population difference $I\sim\langle n_{1}-n_{0}\rangle=\alpha_{1}-\alpha_{0}$ and present the fluorescence as composed of several contributions, two of them exhibiting resonance dependence on $\delta$:

[TABLE]

By making substitutions

[TABLE]

the widths $w_{0}$ and $w_{1}$ and amplitudes $A_{0}$ and $A_{1}$ of the two last resonance contributions can be cast in simple forms:

[TABLE]

Eqs. S13 can be further simplified by considering a strong field limit ($S\gg$1) and making the approximation $\frac{\nu}{2\mu}\simeq 0$. Then the widths can be expressed with very simple formulae:

[TABLE]

Equations S14 represent the power dependence of linewidths of the composite resonance discussed in Figs. 3 and 5 of the main paper.