# Nonlinear Wave-Spin Interactions in Nitrogen-Vacancy Centers

Zhongqiang Hu  Department of Electrical Engineering and Computer Science,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Qiuyuan Wang
Department of Electrical Engineering and Computer Science, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  Chung-Tao Chou  Department
of Electrical Engineering and Computer Science, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Department of Physics, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  Justin T. Hou  Department
of Electrical Engineering and Computer Science, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Zhiping He  Department of Electrical
Engineering and Computer Science, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Luqiao Liu<EMAIL_ADDRESS>Department of Electrical
Engineering and Computer Science, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA

###### Abstract

Nonlinear phenomena represent one of the central topics in the study of wave-
matter interactions and constitute the key blocks for various applications in
optical communication, computing, sensing, and imaging. In this work, we show
that by employing the interactions between microwave photons and electron
spins of nitrogen-vacancy (NV) centers, one can realize a variety of nonlinear
effects, ranging from the resonance at the sum or difference frequency of two
or more waves to electromagnetically induced transparency from the
interference between spin transitions. We further verify the phase coherence
through two-photon Rabi-oscillation measurements. The highly sensitive,
optically detected NV-center dynamics not only provides a platform for
studying magnetically induced nonlinearities but also promises novel
functionalities in quantum control and quantum sensing.

##  I Introduction

Through the mixing of multiple electromagnetic waves [1, 2, 3, 4], nonlinear
processes provide useful mechanisms for frequency up- and down-conversion [5,
6, 7, 8, 9], parametric signal amplification or generation [10, 11, 12, 13],
as well as the creation of entangled photons or squeezed light [14, 15, 16,
17], the fundamental components of quantum information systems. For nonlinear
interactions between waves and matter, electric dipole transitions are
generally considered over their magnetic counterparts due to their larger
strengths [18]. However, restrained by optical selection rules, special
crystals with broken inversion symmetry are usually required for nonlinear
coefficients such as the second-order electric susceptibility to be
nonvanishing [1, 2]. On the other hand, magnetic dipole transitions can
possess nonlinearities even in centrosymmetric systems due to the inherent
breaking of time-reversal symmetry. Nonlinear magnetic dipole transitions,
particularly nonlinear spin transitions, have been touched upon in nuclear
magnetic resonance (NMR) and electron paramagnetic resonance (EPR) [19, 20,
21, 22], where more than one electromagnetic wave source is used for exciting
the resonance. However, due to the very weak wave-spin interactions, these
nonlinear signals are generally difficult to detect. To ensure measurable
resonance, very low frequencies – in the kilohertz or low-megahertz range –
have to be used for at least one of the input sources, making these
measurements effectively the same as the field-modulation scheme of magnetic
resonance. Therefore, a comprehensive study on multiphoton spin transitions
that cover a broad frequency range and that can lead to useful quantum control
and sensing protocols is highly desirable.

The nitrogen-vacancy (NV) center, an extensively studied quantum defect in
diamond, has been pursued as a magnetometer with fine spatial resolution and
high sensitivity [23, 24, 25, 26, 27, 28], and as a qubit for quantum
information processing [29, 30, 31, 32, 33, 34, 35]. To achieve quantum state
control, existing studies focus on linear processes by applying gigahertz
microwaves at or close to the intrinsic resonance frequency. Recently, quantum
frequency mixing based on sophisticated Floquet Hamiltonian engineering has
been developed for magnetic field sensing with NV centers [36]. The magnetic
field at 150 MHz has been detected using the difference frequency of two waves
through a spin-locked sensing protocol, under the assistance of a third,
control signal at the original resonance frequency. In NV-center resonance,
the detection of photons in the visible-light region rather than those in the
radio-frequency or microwave domains greatly enhances the sensitivity, and
leads to a superior platform for studying nonlinear spin transitions. In this
work, we demonstrate such opportunities by carrying out a systematic study on
nonlinear wave-spin interactions in NV centers. We show that the nonlinear
resonance condition can be reached over a broad frequency range, at the sum or
difference frequency of two waves, as well as with higher-order effects
involving three, four, or more photons. Utilizing the interference between
spin transitions, we further show that the resonance can be greatly suppressed
in the presence of a probe wave and a strong control wave, leading to
electromagnetically induced transparency (EIT). Finally, on top of continuous-
wave measurements, we also observe sum-frequency Rabi oscillations, which not
only verifies the phase coherence of these multiphoton processes but also
suggests new mechanisms for quantum control and sensing.

##  II Optical detection of multiphoton spin transitions

Figure 1: (a) The energy-level diagram of an NV center. (b)–(c) Nonlinear spin
transitions in the ground state ${}^{3}A_{2}$ when the (b) sum or (c)
difference of two applied microwave frequencies $\omega_{1}$ and $\omega_{2}$
matches the transition frequency $\omega_{\mathrm{A}}$. (d) A schematic of
nonlinear ODMR measurements. Microwaves from two signal generators are
combined through a linear power combiner, and applied onto a copper
microstrip, exciting NV-center resonance in microdiamonds on top of the strip.
The top-left inset shows an illustration of the geometrical relationship
between microwave magnetic fields $\bm{h}(t)$ and an NV spin, where an angle
of $\theta$ is formed between $\bm{h}(t)$ and the principal spin axis ($z$
axis).

Figure 1(a) illustrates the energy levels of an NV center, where both the
optical ground state ${}^{3}A_{2}$ and excited state ${}^{3}E$ are spin
triplets with spin sublevels of $\ket{m_{s}=0}$ and $\ket{m_{s}=\pm 1}$
separated by $\omega_{\mathrm{A}}/(2\pi)=2.87\ \mathrm{GHz}$ and
$\omega_{\mathrm{E}}/(2\pi)=1.42\ \mathrm{GHz}$, respectively, under zero
external static field [37, 38]. Green light can induce the transition from
${}^{3}A_{2}$ to ${}^{3}E$, and the $m_{s}$-conserving decay from ${}^{3}E$ to
${}^{3}A_{2}$ generates photoluminescence (PL) in the region of red light
[39]. The nonradiative transition path through spin singlet states
${}^{1}A_{1}$ and ${}^{1}E$ pumps the NV population into the $\ket{m_{s}=0}$
sublevel, which can be suppressed with the application of a microwave at or
close to the sublevel splitting $\omega_{\mathrm{A}}$ or
$\omega_{\mathrm{E}}$, yielding a reduction of the PL intensity [40, 41, 42].
In this work, we will delve into optically detected magnetic resonance (ODMR)
beyond the linear response regime and investigate spin transitions induced by
multiple photons, through concurrent application of two or more microwaves. In
Figs. 1(b) and 1(c), we illustrate two example scenarios where the sum or
difference of the two applied frequencies $\omega_{1}$ and $\omega_{2}$
matches the ground-state transition frequency $\omega_{\mathrm{A}}$. To excite
magnetic resonance in experiments, microwaves from two independent signal
generators are combined through a power combiner [Fig. 1(d)]. We have verified
that under our employed power levels, the external microwave circuit acts
purely linearly and is not the origin of frequency mixing (see Appendix A).
The microwaves are further applied onto a lithographically defined copper
microstrip on a silicon substrate. Diamond particles with a diameter of
approximately $1\ \rm\mu m$ and an NV-center concentration of approximately
$3.5\ \rm ppm$ (parts per million) are dispersed on top of the strip. PL
excited by a $532$-$\rm nm$ green laser is filtered and collected with a
photomultiplier tube.

Figure 2: (a) The change of the PL intensity $\Delta I$ under driving
microwaves of $\omega_{1}$ and $\omega_{2}$ with input powers $P_{1}=P_{2}=13\
\mathrm{mW}$. Nonlinear resonance signals emerge at (1)
$\omega_{2}+\omega_{1}=\omega_{\mathrm{A}}$, (2)
$\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$, (3)
$\omega_{1}-\omega_{2}=\omega_{\mathrm{A}}$, (4)
$3\omega_{2}=\omega_{\mathrm{A}}$, (5) $4\omega_{2}=\omega_{\mathrm{A}}$, (6)
$2\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$, and (7)
$2\omega_{1}-\omega_{2}=\omega_{\mathrm{A}}$. (b) $\omega_{2}$ scan in (a)
when $\omega_{1}/(2\pi)=0.6\ \mathrm{GHz}$. (c) The $\omega_{2}$ scan when
$P_{1}=0$, serving as the baseline of the measurement. The splitting of the
linear resonance dip at $\omega_{2}=\omega_{\mathrm{A}}$ (red arrow) is a
result of amplitude modulation at high applied microwave powers, which
disappears in unmodulated results (see Appendix A). The side dips denoted by
green arrows originate from the interactions between NV centers and P1 centers
[43].

In Fig. 2(a), we show the change of the PL intensity $\Delta I$ under driving
microwaves of $\omega_{1}$ and $\omega_{2}$. To enhance the signal-to-noise
ratio, we modulate the amplitude of the $\omega_{2}$ input with a frequency of
$104.42\ \mathrm{Hz}$ and detect $\Delta I$ using a lock-in amplifier. We have
compared results from this lock-in measurement with a standard unmodulated
continuous-wave measurement (see Appendix A) and we have confirmed that these
two give the consistent results and that the low-frequency amplitude
modulation is not the source of the observed nonlinear effects. Since the
amplitude modulation only acts on the $\omega_{2}$ input, it gives rise to the
asymmetry on the dependence of $\Delta I$ with respect to $\omega_{1,2}$ in
Fig. 2(a): resonance signals only show up at
$\omega_{2}=\omega_{\mathrm{A,E}}$ but not at
$\omega_{1}=\omega_{\mathrm{A,E}}$, in contrast to the unmodulated measurement
results. The large line widths associated with the
$\omega_{2}=\omega_{\mathrm{A,E}}$ resonance dips reflect the applied high
microwave power and low laser pump power (about 0.6 mW) [39, 44]. Besides the
standard, linear resonance dips, in Fig. 2(a) additional resonance signals
emerge when $\omega_{1,2}$ satisfy the relationship of
$\omega_{2}+\omega_{1}=\omega_{\mathrm{A}}$ (labeled as 1) or
$\pm(\omega_{2}-\omega_{1})=\omega_{\mathrm{A}}$ (labeled as 2 and 3). As an
example, we show the spectrum when $\omega_{1}/2\pi$ is fixed at $0.6\
\mathrm{GHz}$ and $\omega_{2}$ is swept [Fig. 2(b)] and compare it with the
baseline when the power of the $\omega_{1}$ input is set to zero [Fig. 2(c)].
In Fig. 2(b), the depths of the $\omega_{2}+\omega_{1}=\omega_{\mathrm{A}}$
and $\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$ dips reach around 50% and 20%
of that of the $\omega_{2}=\omega_{\mathrm{A}}$ dip. The two nonlinear
resonance dips can be easily detected even when $\omega_{1,2}$ are
individually far away from $\omega_{\mathrm{A}}$.

##  III Theoretical model for multiphoton spin transitions

To understand the origin of magnetic resonance occurring at the sum or
difference frequency, we next model the NV spin transitions driven by multiple
microwave photons. Here, we consider, e.g., the spin transition between
$\ket{m_{s}=0}$ and $\ket{m_{s}=+1}$ in ${}^{3}A_{2}$. Other transitions, such
as that between $\ket{m_{s}=0}$ and $\ket{m_{s}=-1}$ and those in ${}^{3}E$,
can be treated similarly. We write the Hamiltonian of the spin-photon system
as

$H(t)=\frac{\hbar\omega_{\mathrm{A}}}{2}\sigma_{z}+\gamma\mu_{0}\frac{\hbar}{2}\bm{\sigma}\cdot\bm{h}(t),$
(1)

where $\hbar$ is reduced Planck constant,
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$ are Pauli matrices, $\gamma$
is the electron’s gyromagnetic ratio, $\mu_{0}$ is the vacuum permeability,
and $\smash{\bm{h}(t)=\sum_{j=1,2}\bm{h}_{j}(t)}$ are microwave fields with
the $j$th frequency component $\bm{h}_{j}(t)$. In experiments, the two
microwave fields are launched by the same microstrip; thus $\bm{h}_{1}(t)$ and
$\bm{h}_{2}(t)$ are collinear and form an angle $\theta$ with the principal
axis of a given NV spin [see the top-left inset of Fig. 1(d)],
$\bm{h}_{j}(t)=h_{j}(\sin\theta\hat{x}+\cos\theta\hat{z})\cos(\omega_{j}t+\varphi_{j})$,
where $h_{j}$, $\omega_{j}$, and $\varphi_{j}$ are the amplitude, frequency,
and phase of the $j$th field ($j=1,2$).

Solving the quantum master equations iteratively (see the derivation details
in Appendixes B and C), we obtain the change of the PL intensity for an
ensemble of spins when the sum or difference of $\omega_{1,2}$ is close to the
resonance condition of $\omega_{2}\pm\omega_{1}=\omega_{\mathrm{A}}$:

$\Delta
I=-\eta|\chi_{xxz}^{(2)}(\omega_{2}\pm\omega_{1},\omega_{1},\omega_{2})|^{2}h_{1}^{2}h_{2}^{2}.$
(2)

Here, $\eta=|\Delta
I|_{\mathrm{max}}\Gamma_{\mathrm{2,A}}/(16\gamma^{2}\hbar^{2}\Gamma_{\mathrm{p}})$
is a factor depending on the maximum PL intensity $|\Delta I|_{\mathrm{max}}$,
the effective transverse relaxation rate $\Gamma_{\mathrm{2,A}}$ in
${}^{3}A_{2}$, as well as the laser pump rate $\Gamma_{\mathrm{p}}$.
$\smash{\chi_{xxz}^{(2)}}$, defined through
$M_{x}=\smash{\chi_{xxz}^{(2)}}h_{x}h_{z}$, is an element in the second-order
magnetic susceptibility tensor for a single spin, where $M_{x}$ is the
$x$-axis component of the magnetic moment and $\smash{h_{x(z)}}$ is the
$x(z)$-axis component of $\bm{h}(t)$. Close to the resonance, we have

$\chi_{xxz}^{(2)}(\omega_{2}\pm\omega_{1},\omega_{1},\omega_{2})=\mp\frac{\Gamma_{\mathrm{p}}\gamma^{3}\mu_{0}^{2}\hbar\omega_{\mathrm{A}}}{4\Gamma_{\mathrm{1,A}}(\Delta_{\mathrm{\pm,A}}-i\Gamma_{\mathrm{2,A}})\omega_{1}\omega_{2}},$
(3)

where the $+\ (-)$ sign is chosen when the resonance condition is satisfied by
the sum(difference) frequency, $\Gamma_{\mathrm{1,A}}$ is the effective
longitudinal relaxation rate in ${}^{3}A_{2}$, and
$\Delta_{\mathrm{\pm,A}}=\omega_{2}\pm\omega_{1}-\omega_{\mathrm{A}}$ denotes
the frequency detuning. We can infer important information on the two-photon
spin transition from $\smash{\chi_{xxz}^{(2)}}$. First,
$\smash{\chi_{xxz}^{(2)}}$ remains finite irrespective of the atom’s inversion
symmetry, in contrast to the second-order electric susceptibility which
vanishes for atoms with inversion symmetry [1, 2]. Second, the fact that the
magnetic moment resonates at the sum or difference frequency of the input
sources reflects the energy conservation. Third, angular-momentum conservation
is also respected in the two-photon spin transition described by
$\smash{\chi_{xxz}^{(2)}}$: the perpendicular microwave field $h_{x}$
comprises circularly polarized photons $\sigma^{\pm}$, each carrying an
angular momentum of $\pm\hbar$, while the parallel component $h_{z}$ comprises
$\pi$ photons that carry zero angular momentum [22]. Therefore, the transition
from $\ket{m_{s}=0}$ to $\ket{m_{s}=\pm 1}$ via $\sigma^{\pm}+\pi$ leaves the
total angular momentum of the spin-photon system unchanged.

Figure 3: (a) The sum-frequency resonance under different input microwave
powers ($P_{1}=P_{2}$ is maintained), when $\omega_{1}/(2\pi)=0.6\ \rm GHz$.
Neighboring curves and data are shifted by $-0.5$ for clear visualization. (b)
The contrast $C$ of the sum-frequency resonance as a function of $P_{1}$ when
$P_{2}$ is fixed. (c) The $C$ value of the sum-frequency resonance as a
function of $P_{2}$ when $P_{1}$ is fixed. (d) The $C$ value of the sum-
frequency resonance as a function of $P_{1}P_{2}$. (e) The $C$ value of the
sum- or difference-frequency resonance as a function of $\omega_{2}$
($P_{1}=P_{2}=9\ \rm mW$). Dashed lines in (b)–(e) represent fitting with Eq.
(4) using the same parameter $\eta$. (f) The $C$ value of the three-photon
resonance at $2\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$ as a function of
$P_{2}$ when $P_{1}$ is fixed. (g) The $C$ value of the
$2\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$ resonance as a function of
$P_{1}P_{2}^{2}$. (h) The $C$ value of the
$2\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$ resonance as a function of
$\omega_{2}$ ($P_{1}=P_{2}=11\ \rm mW$). The dashed lines in (f)-(h) represent
fitting with $\smash{\chi^{(3)}}$, the expression for which is derived in
Appendix D.

We next compare our measurement results quantitatively with the proposed
theory. Equation (2) represents a Lorentzian with a line width of
$\Gamma_{\mathrm{2,A}}$ and a contrast of

$C=\frac{|\Delta
I|_{\mathrm{max}}\Gamma_{\mathrm{p}}\omega_{\mathrm{A}}^{2}\gamma^{4}\mu_{0}^{4}h_{1}^{2}h_{2}^{2}}{256\Gamma_{\mathrm{1,A}}^{2}\Gamma_{\mathrm{2,A}}\omega_{1}^{2}\omega_{2}^{2}}.$
(4)

In Fig. 3(a), we show a series of resonance curves under varying input
microwave powers $P_{1,2}$, when $\omega_{1}/(2\pi)$ is fixed at $0.6\ \rm
GHz$ and $\omega_{2}$ is swept around $\omega_{\mathrm{A}}-\omega_{1}$. The
line shape in Eq. (2) (dashed lines) agrees well with the experimental results
(solid circles). The resonance contrast $C$ obtained under different
combinations of $P_{1,2}$ is summarized in Figs. 3(b)–3(d), in which the
theoretical curves from Eq. (4) (dashed lines) confirm the linear relationship
between $C$ and $P_{1}P_{2}$. The same value of $\eta$, the only fitting
parameter, is used across all the curves. The influence from frequencies of
the two applied microwaves is presented in Fig. 3(e), where $\omega_{1,2}$ are
varied simultaneously, and their sum or difference is maintained at
$\omega_{2}\pm\omega_{1}=\omega_{\mathrm{A}}$. $C$ has a frequency dependence
of $\omega_{2}^{-2}(\omega_{\mathrm{A}}-\omega_{2})^{-2}$, consistent with Eq.
(4). The fitting curves (dashed lines) use the same value of $\eta$ as in
Figs. 3(b)–3(d).

Besides the two-photon resonance investigated above, in Fig. 2(a) we observe
additional bright lines (labeled from 4 to 7), which can be traced to magnetic
resonance excited by even higher-order processes. Our examination shows that
signals of 4, 6, and 7 satisfy the frequency relationship of (4)
$3\omega_{2}=\omega_{\mathrm{A}}$, (6)
$2\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$, and (7)
$2\omega_{1}-\omega_{2}=\omega_{\mathrm{A}}$, corresponding to three-photon
processes. The horizontal line with a resonance frequency of $0.72\
\mathrm{GHz}$ (labeled as 5) corresponds to four-photon resonance with
$4\omega_{2}=\omega_{\mathrm{A}}$. These higher-order processes can be well
described with the theoretical framework that we have developed, by
calculating higher-order magnetic susceptibilities. For example, in Figs. 3(f)
and 3(g) we summarize the dependence on the input powers for the three-photon
resonance at $2\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$. As indicated by
$\smash{\chi^{(3)}}$ (see Appendix D), the contrast of this resonance has a
quadratic dependence on $P_{2}$ and a linear dependence on $P_{1}$. The
influence from frequencies of the two applied microwaves is presented in Fig.
3(h), which also fits well with the theory. Finally, we note that the
$2\omega_{2}=\omega_{\mathrm{A}}$ line, as a special case of the sum-frequency
resonance, falls onto the broad $\omega_{2}=\omega_{\mathrm{E}}$ resonance dip
and is difficult to distinguish due to the small difference between
$\omega_{\mathrm{A}}/2$ and $\omega_{\mathrm{E}}$.

While the results in Figs. 2 and 3 correspond to the scenario with zero
external static field, we have verified that the nonlinear processes remain in
the presence of a finite static field that lifts the degeneracy in
$\ket{m_{s}=\pm 1}$ (see Appendix A).

##  IV Electromagnetically induced transparency in NV centers

Figure 4: A demonstration of the EIT effect. (a) $\Delta I$ under $\omega_{1}$
and $\omega_{2}$ inputs when $\omega_{1,2}$ are close to
$\omega_{\mathrm{E}}$. The microwave powers are $P_{1}=13\ \mathrm{mW}$,
$P_{2}=5\ \mathrm{mW}$. (b) Destructive interference between the spin
transition with $\omega_{2}$ absorption [associated with $\smash{\chi^{(1)}}$;
see Appendix E] and that with $\omega_{2}$ absorption, $\omega_{1}$ emission,
$\omega_{1}$ absorption [associated with $\smash{\chi^{(3)}}$]. (c)–(d)
$\Delta I$ as a function of $\omega_{2}$ when $\omega_{1}/(2\pi)$ is fixed at
(c) $1.4\ \rm GHz$ and (d) $1.46\ \rm GHz$. The dashed lines represent fitting
with Eq. (5).

In Fig. 2(a), on top of the series of the bright resonance lines, we observe a
dark line satisfying $\omega_{1}=\omega_{2}$ within the broad resonance dip at
$\omega_{2}=\omega_{\mathrm{E}}$. This feature of magnetic resonance
suppression under zero detuning of two waves is very similar to the
electromagnetically induced transparency (EIT) phenomenon studied in nonlinear
optics [45], where in the presence of a strong pump wave, the interaction
between the probe wave and the matter is minimized due to interference
effects. This signal-suppression phenomenon has also been reported before
under the concept of coherent population oscillation or trapping, for nuclear
and electron spin systems [46, 47]. Treating the $\omega_{1}$ signal as the
pump and the modulated $\omega_{2}$ signal as the probe, we examine the
resonance results with a weaker probe power $P_{2}$ while maintaining a high
pump power $P_{1}$ and find that the transparency window still exists [Fig.
4(a)]. Examples of $\omega_{2}$ scans are presented in Figs. 4(c) and 4(d). On
the other hand, the transparency window disappears when a low pump power
$P_{1}$ is used (see Appendix E). The magnetic resonance suppression under
$\omega_{1}=\omega_{2}$ can be explained by considering the destructive
interference between the $\ket{m_{s}=0}\xrightarrow{\omega_{2}}\ket{m_{s}=\pm
1}$ spin transition associated with a first-order susceptibility and the
$\ket{m_{s}=0}\xrightarrow{\omega_{2}}\ket{m_{s}=\pm
1}\smash{\xrightarrow{-\omega_{1}}}\ket{m_{s}=0}\xrightarrow{\omega_{1}}\ket{m_{s}=\pm
1}$ spin transition associated with a third-order susceptibility [Fig. 4(b)].
Mathematically, when $\omega_{1}\approx\omega_{2}$ and both of them are close
to $\omega_{\mathrm{E}}$, we have (see the derivation details in Appendix E)

$\Delta
I=\frac{iA}{\omega_{2}-\omega_{\mathrm{E}}-i\Gamma_{\mathrm{2,E}}-\frac{B}{\omega_{2}-\omega_{1}-i\Gamma_{\mathrm{1,E}}}}+\mathrm{H.c.},$
(5)

where $A=|\Delta
I|_{\mathrm{max}}\Gamma_{\mathrm{p}}\gamma^{2}\mu_{0}^{2}h_{2}^{2}/(32\Gamma_{\mathrm{1,E}}^{2})$,
$B=\gamma^{2}\mu_{0}^{2}h_{1}^{2}/4$, and $\Gamma_{\mathrm{1(2),E}}$ is the
effective longitudinal(transverse) relaxation rate in ${}^{3}E$. In Eq. (5),
we see that the peak value of $|\Delta I|$ is suppressed at zero detuning
$\omega_{2}-\omega_{1}=0$ and that the EIT effect is most significant when
$\Gamma_{\mathrm{2,E}}>\gamma\mu_{0}h_{1}/2>\Gamma_{\mathrm{1,E}}$, which is
satisfied in the spin transitions in ${}^{3}E$. Comparatively, the relatively
smaller $\Gamma_{\mathrm{2,A}}$ makes the EIT feature in ${}^{3}A_{2}$ less
noticeable.

##  V Coherent control of NV spin state through multiple microwave photons

Figure 5: (a) Experimental schematic for the two-photon Rabi-oscillation
measurements. A single-crystal diamond with a diameter of 15 $\mathrm{\mu m}$
is placed at the center of a copper microstrip ring. The inner and outer
diameter of the ring is 60 $\mathrm{\mu m}$ and 100 $\mathrm{\mu m}$,
respectively. CH2 and CH3 of the AWG modulate the laser and two microwave
sources, respectively, where CH1 provides a low-frequency (200 Hz) reference
signal for the lock-in amplifier. (b) Pulse sequences for lock-in amplifier
(blue), laser (green), and two microwave sources (violet). (c) ODMR spectrum
in the presence of an external static field of 300 Oe under $\omega_{2}$
microwave source only, with $P_{2}=13\ \mathrm{mW}$. The leftmost resonance at
$\omega_{\mathrm{A1}}/(2\pi)=2.08\ \mathrm{GHz}$ is selected for Rabi-
oscillation measurements. (d) Standard Rabi oscillation results under a single
source $\omega_{2}=\omega_{\mathrm{A1}}$ with $P_{2}=200\ \mathrm{mW}$. The
single-photon Rabi frequency is 14.5 MHz. (e) Rabi oscillation results when
simultaneously applying $\omega_{1}/(2\pi)=60\ \mathrm{MHz}$, $P_{1}=13\
\mathrm{mW}$ and $\omega_{2}=\omega_{\mathrm{A1}}-\omega_{1}$, $P_{2}=200\
\mathrm{mW}$. The two-photon Rabi frequency is 2.03 MHz. In (d) and (e), the
experimental data (solid circles) are fitted with Eq. (6) (dashed lines).

Up to now, the continuous-wave measurements as described above have allowed us
to capture magnetic resonance signals at different orders for a broad
frequency range. In the following, we carry out Rabi oscillation experiments
to evaluate the phase coherence of the nonlinear multiphoton spin transitions,
which is of paramount significance in developing effective quantum control and
sensing protocols. The experimental schematic is shown in Fig. 5(a). A single-
crystal diamond with a diameter of 15 $\mathrm{\mu m}$ is placed at the center
of a copper microstrip ring. The inner and outer diameters of the ring are 60
$\mathrm{\mu m}$ and 100 $\mathrm{\mu m}$, respectively, such that the
generated microwave fields are nearly constant within the focal spot of the
objective lens, the diameter of which is $<1\ \mathrm{\mu m}$. The switch to
an individual diamond and a ring-shaped microstrip is to minimize the
inhomogeneity in the detected NV centers. An arbitrary waveform generator
(AWG) is used to program the pulse sequences for the measurements. CH2 and CH3
modulate the laser and two microwave sources, respectively, where CH1 provides
a low-frequency (200 Hz) reference signal for a lock-in amplifier. Figure 5(b)
shows the pulse sequences for the lock-in amplifier (blue), the laser (green),
and two microwave sources (violet). In our experiments,
$\tau_{\mathrm{laser}}$ is set at $5\ \mathrm{\mu s}$, $\tau_{\mathrm{mw}}$ is
varied between $0.1\ \mathrm{\mu s}$ and $2\ \mathrm{\mu s}$, and
$\tau_{\mathrm{gap}}$ is fixed at $1\ \mathrm{\mu s}$. The lock-in reference
signal (CH1) is fixed at 200 Hz with $\tau_{\mathrm{ref}}=2.5\ \mathrm{ms}$.
Within each ON half period of CH1, the laser and microwave pulses are repeated
$N=250$ times. Within each OFF half period of CH1, only the laser pulse is
repeated with the microwave pulse always OFF. This lock-in-based pulse method
avoids the necessity of high-frequency electronics for data acquisition and
possesses high sensitivity [48].

In Fig. 5(c), we show the single-frequency ODMR spectrum under an external
static field of $H_{\mathrm{ext}}=300\ \mathrm{Oe}$. We see eight resonance
dips, corresponding to four different NV-axis orientations and two spin
transitions $\ket{m_{s}=0}\rightarrow\ket{m_{s}=\pm 1}$ in the NV ground state
${}^{3}A_{2}$. The leftmost dip, with a resonance frequency of
$\omega_{\mathrm{A1}}/(2\pi)=2.08\ \mathrm{GHz}$, is selected for oscillation
measurements. By applying a single microwave source $\omega_{2}$ at the exact
resonance frequency and varying its pulse width, we summarize the oscillation
results in Fig. 5(d), which are fitted with

$\Delta I=A\sin\left(2\pi
f_{\mathrm{R}}\tau_{\mathrm{mw}}+\phi\right)\exp\left(-\frac{\tau_{\mathrm{mw}}}{T_{2}^{*}}\right)+B\tau_{\mathrm{mw}}+C,$
(6)

where $A$ is the Rabi oscillation amplitude, $B$ is a linear coefficient that
includes the heating caused shift during measurement, $C$ is the readout at
the steady state, $f_{\mathrm{R}}$ is the Rabi frequency, $\phi$ is the phase
offset, and $T_{2}^{*}$ is the transverse relaxation time. The single-photon
Rabi frequency $\smash{f_{\mathrm{R}}^{(1)}}$ is determined to be $14.5\
\mathrm{MHz}$. We next realize two-photon Rabi oscillations by applying both
microwave sources $\omega_{1}/(2\pi)=60\ \mathrm{MHz}$ and
$\omega_{2}/(2\pi)=2.02\ \mathrm{GHz}$, with synchronous pulse modulation. We
have verified that the $\omega_{1}$ or $\omega_{2}$ source at these
frequencies alone does not induce detectable oscillation signals. However, as
they satisfy the sum-frequency resonance condition, their concurrent
application leads to clear oscillations, as shown in Fig. 5(e), with a Rabi
frequency of $\smash{f_{\mathrm{R}}^{(2)}}=2.03\ \mathrm{MHz}$. The phase
coherence demonstrated by our Rabi-oscillation measurements hopefully paves
the way for the use of multiphoton spin transitions in practical quantum
control and sensing schemes.

##  VI Conclusions

In conclusion, we have conducted a systematic study on the nonlinear
interactions between electromagnetic waves and electron spins of NV centers.
With the high sensitivity of NV-center resonance, we have revealed various
nonlinear, parametric processes over a broad frequency range, ranging from
high-order magnetic resonance to EIT. We have developed a theoretical
framework based on perturbation theory to account for these nonlinear
phenomena quantitatively. In addition, we have verified the phase coherence of
the multiphoton spin transitions through Rabi-oscillation measurements, which
hopefully paves the way toward future applications in quantum control and
sensing. Furthermore, leveraging the sum- or difference-frequency resonance,
one can make a nanoscale spectrometer out of NV centers to extract the
spectrum information of oscillating magnetic fields such as those from spin-
wave excitations in magnetic materials [49, 50, 51], extending their well-
established role as a magnetometer.

###### Acknowledgements.

The work is supported by the National Science Foundation under Award No.
DMR-2104912. L.L. acknowledges support from the Sloan Research Fellowship.
Z.H. acknowledges support from the MathWorks Fellowship. The manuscript has
been accepted by Physical Review Applied (APS copyright) on Mar 29, 2024.

##  Appendix A Experimental method

Diamond samples: The samples used for continuous-wave measurements are high-
pressure and high-temperature (HPHT) microdiamonds with a diameter of 1
$\mathrm{\mu m}$ and an NV concentration of 3.5 ppm (MDNV1umHi, Adámas
Nanotechnologies). The samples used for the Rabi-oscillation measurements are
15-$\mathrm{\mu m}$ diamond particles from the same company, with the same NV
concentration and growth method (MDNV15umHi). The continuous-wave measurements
have been carried out on a cluster of diamond particles, while the Rabi
oscillations have been done on an individual diamond.

Figure 6: (a) A photograph of the device wire bonded to a home-made PCB. The
microwaves are fed into the waveguide via an SMA connector. (b) The spectrum
of the two-frequency microwaves after passing through the power combiner, when
$\omega_{1}/(2\pi)=0.6\ \mathrm{GHz}$ and $\omega_{2}/(2\pi)=2.27\
\mathrm{GHz}$. Both the output powers in the signal generators are set as $17\
\mathrm{dBm}$. Figure 7: The results of standard unmodulated continuous-wave
ODMR measurements under driving microwaves of $\omega_{1}$ and $\omega_{2}$.
The input microwave powers are $P_{1}=P_{2}=13\ \mathrm{mW}$.

Device fabrication: For continuous-wave measurements, through standard
photolithography followed by ion milling, we pattern a Cu(100 nm)/Pt(5 nm)
stack on a silicon substrate into a straight microstrip with a length and
width of 100 $\mathrm{\mu m}$ and 20 $\mathrm{\mu m}$, respectively, which is
further wire bonded onto a home-made printed circuit board (PCB) with an SMA
connector [Fig. 6(a)]. Diamonds with a 1-$\mathrm{\mu m}$ diameter are
dispersed on top of the microstrip. By calibrating the microwave signal, we
determine that when the input power is 10 mW at the input terminal of the PCB,
the microwave field is approximately 9 Oe at the microstrip surface. For the
pulsed measurements on Rabi oscillations, a Cu(500 nm)/Pt(10 nm) stack on a
sapphire substrate is patterned into a microstrip ring with an inner and outer
diameter of 60 $\mathrm{\mu m}$ and 100 $\mathrm{\mu m}$, respectively, which
is also wire bonded onto the same PCB. A single-crystal 15-$\mathrm{\mu m}$
diamond particle is placed at the center of the ring, to minimize the
inhomogeneity of oscillating fields generated by microwaves. When the input
microwave power is 200 mW, the field magnitude at the center of the ring is
calibrated to be approximately 8 Oe.

Figure 8: $\Delta I$ as a function of $H_{\mathrm{ext}}$ when $\omega_{2}$ is
swept and $\omega_{1}/(2\pi)$ is fixed at (a) $0.4\ \rm{GHz}$, (b) $0.6\
\rm{GHz}$, and (c) $0.8\ \rm{GHz}$. Both the linear resonance dips at
$\omega_{2}=\omega_{\mathrm{A,E}}$ and the nonlinear dips at
$\omega_{2}\pm\omega_{1}=\omega_{\mathrm{A}}$ manifest a cone-shaped structure
with edges separated by $2\gamma\mu_{0}H_{\mathrm{ext}}$.

Continuous-wave ODMR measurements: The setup for continuous-wave ODMR
measurements depicted in Fig. 1(d) mainly consists of a home-built confocal
microscope. To excite magnetic resonance in NV centers, continuous-wave
signals are generated from two independent microwave signal generators
(Anritsu 68369A and Anritsu 68347B) and then combined through a microwave
power combiner (CentricRF CS6072). Green light from a $532$-$\rm nm$ DPSS
laser is focused on the sample via a $1.25$-$\rm NA$ objective lens and
illuminates NV centers in an ensemble of microdiamonds. The laser power is
$0.6\ \mathrm{mW}$ measured on the sample surface. We intentionally choose a
low-power laser diode as the excitation light source, to enhance the ratio of
resonance signals from high-order effects to those from linear effects by
forcing the linear signals to approach saturation under limited optical
pumping [39, 44].The PL in the region of $600$-$800\ \mathrm{nm}$ from NV
centers is filtered and collected with a photomultiplier tube. To enhance the
signal-to-noise ratio, the amplitude of one microwave source is modulated with
a frequency of $104.42\ \mathrm{Hz}$ and the change of the PL intensity is
detected with a lock-in amplifier (EG&G 7260). We have varied the position of
the illuminated NV-center ensemble and found no qualitative difference in the
ODMR results.

Linearity verification of the circuit: With a spectrum analyzer (Anritsu
MS8609A), we verify that the external microwave circuit acts purely linearly
and is not the origin of frequency mixing. In Fig. 6(b), we show the spectrum
of microwaves after they pass through the power combiner, when
$\omega_{1}/(2\pi)=0.6\ \mathrm{GHz}$ and $\omega_{2}/(2\pi)=2.27\
\mathrm{GHz}$. Both the output powers in the signal generators are set as $17\
\mathrm{dBm}$. Peaks only appear at $\omega_{1}$ and $\omega_{2}$, verifying
the linearity of the external microwave circuit.

Comparison with unmodulated ODMR measurements: In Fig. 7, we show the PL
intensity $I$ under driving microwaves of $\omega_{1}$ and $\omega_{2}$, when
we conduct a standard unmodulated continuous-wave measurement with an Agilent
34401A multimeter. Compared to Fig. 2(a), where the amplitude of the
$\omega_{2}$ microwave is modulated with a low frequency, the unmodulated
measurement has a much lower signal-to-noise ratio but we can still clearly
identify the resonances at $\omega_{2}+\omega_{1}=\omega_{\mathrm{A}}$. The
difference frequency resonances at
$\pm(\omega_{2}-\omega_{1})=\omega_{\mathrm{A}}$ are buried in the noisy
background to some degree but we can still distinguish them with extra
attention. The feature of EIT, i.e., the dark line satisfying
$\omega_{1}=\omega_{2}$ when $\omega_{1,2}$ are both close to
$\omega_{\mathrm{E}}$, is also observed in the unmodulated measurement. In
conclusion, the unmodulated and modulated measurements show consistency in
demonstrating the resonance at the sum or difference frequency, as well as the
EIT effect. The modulated one gives a much higher signal-to-noise ratio and is
therefore applied in the main experiment.

Field-dependent ODMR measurements: The Zeeman splitting from a finite external
static field $H_{\mathrm{ext}}$ (generated by an electromagnet in our
experiments) will lift the degeneracy and result in distinct resonance
frequencies for transitions from $\ket{m_{s}=0}$ to $\ket{m_{s}=\pm 1}$. In
Figs. 8(a)–8(c), we show $\Delta I$ as a function of $H_{\mathrm{ext}}$ when
$\omega_{2}$ is swept and $\omega_{1}/(2\pi)$ is fixed at $0.4$, $0.6$, and
$0.8\ \rm{GHz}$, respectively. Due to random orientations of the NV-center
principle axes with respect to the external field, the splitting in
$\ket{m_{s}=\pm 1}$ manifests a cone-shaped structure in the resonance
spectrum, as observed near $\omega_{2}=\omega_{\mathrm{A,E}}$. In a similar
vein, the nonlinear resonance near
$\omega_{2}\pm\omega_{1}=\omega_{\mathrm{A}}$ also exhibits this conical
feature, with edges separated by $2\gamma\mu_{0}H_{\mathrm{ext}}$.

Rabi-oscillation measurements: The setup for two-photon Rabi-oscillation
measurements is depicted in Fig. 4(a). We use an AWG (Feelelec FY8300) to
modulate the laser and two microwave sources. The microwave sources and the
power combiner are the same as the ones we use for continuous-wave
measurements. The acousto-optic modulator that we use is the Isomet Model
1205C-1 with a driver of Model 532C. A permanent magnet is used to generate a
static field of approximately 300 Oe at the position of the diamond sample,
which helps us select one specific resonance frequency. The Rabi-oscillation
data are obtained by recording the lock-in readout of the PL signals for
varying pulse widths of microwaves.

##  Appendix B Solution of the density matrix for the NV spin under two
driving microwaves

When the external static field is zero, the spin sublevels $\ket{m_{s}=\pm 1}$
in the NV optical ground state ${}^{3}A_{2}$ or the excited state ${}^{3}E$
are degenerate. In what follows, we derive the microwave photon induced spin
transitions within ${}^{3}A_{2}$ or ${}^{3}E$. Throughout this section, we
discuss the general case, in which $\omega_{s}$ represents the transition
frequency between $\ket{m_{s}=\pm 1}$ and $\ket{m_{s}=0}$ either in the ground
state ${}^{3}A_{2}$ or in the excited state ${}^{3}E$. Due to the symmetry
between the $\ket{m_{s}=0}\leftrightarrow\ket{m_{s}=+1}$,
$\ket{m_{s}=0}\leftrightarrow\ket{m_{s}=-1}$ transitions and their negligible
mixing, we can focus on the $\ket{m_{s}=0}\leftrightarrow\ket{m_{s}=+1}$
transition and consider the two-level spin system on the basis of
$\ket{0}=\ket{m_{s}=0}$ and $\ket{1}=\ket{m_{s}=+1}$. The intrinsic spin
Hamiltonian is

$H_{0}=\frac{\hbar\omega_{s}}{2}\sigma_{z},$ (7)

with reduced Planck constant $\hbar$ and Pauli matrices
$\bm{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})$. We apply microwave fields

$\bm{h}(t)=\sum_{j=1,2}\bm{h}_{j}(t)=\sum_{j=1,2}h_{j}(\sin\theta\hat{x}+\cos\theta\hat{z})\cos(\omega_{j}t+\varphi_{j}),$
(8)

where $\bm{h}_{j}(t)$ is the $j$th frequency component, $h_{j}$, $\omega_{j}$,
and $\varphi_{j}$ are the amplitude, frequency, and phase of $\bm{h}_{j}(t)$,
and $\theta$ is the angle between $\bm{h}_{j}(t)$ and the principal spin axis
($z$ axis) of the NV center ($j=1,2$). The interaction Hamiltonian between the
spin and microwave photons is given by

$V(t)=\gamma\mu_{0}\frac{\hbar}{2}\bm{\sigma}\cdot\bm{h}(t),$ (9)

with electron’s gyromagnetic ratio $\gamma$ and vacuum permeability $\mu_{0}$.
The total Hamiltonian is therefore given by

$H(t)=H_{0}+V(t).$ (10)

The density matrix $\rho=\sum_{m,n=0,1}\rho_{mn}\ket{m}\bra{n}$ for the spin
can be determined by solving the quantum master equation in the Lindblad form
[39]:

$\frac{\partial\rho}{\partial
t}=\frac{i}{\hbar}[\rho,H]+\sum_{n}\left(L_{n}\rho
L_{n}^{\dagger}-\frac{1}{2}L_{n}^{\dagger}L_{n}\rho-\frac{1}{2}\rho
L_{n}^{\dagger}L_{n}\right),$ (11)

where $L_{n}$ is the operator describing a nonunitary time evolution due to
dissipative interactions between the spin and the environment. The
longitudinal spin relaxation with a rate of $\Gamma_{1}^{0}$ and the
transverse spin relaxation with a rate of $\Gamma_{2}^{0}$ can be described by

$\displaystyle L_{1}$ $\displaystyle=(\Gamma_{1}^{0}/2)^{1/2}\sigma_{x},$ (12)
$\displaystyle L_{2}$ $\displaystyle=(\Gamma_{2}^{0}/2)^{1/2}\sigma_{z}.$

Due to the photon bath of the laser pumping, $\ket{1}$ is optically pumped
into $\ket{0}$ with a pump rate of $\Gamma_{\mathrm{p}}$, which can be
described by

$\displaystyle L_{3}$ $\displaystyle=\Gamma_{\mathrm{p}}^{1/2}\ket{0}\bra{1}.$
(13)

Combining Eqs. (7)–(13), we obtain equations of motion for elements in the
density matrix

$\displaystyle\frac{\partial{\rho}_{11}}{\partial t}$
$\displaystyle=-\left(\frac{\Gamma_{1}^{0}}{2}+\Gamma_{\mathrm{p}}\right)\rho_{11}+\frac{\Gamma_{1}^{0}}{2}\rho_{00}+\frac{i}{2}\gamma\mu_{0}h_{x}(\rho_{01}-\rho_{10}),$
(14) $\displaystyle\frac{\partial{\rho}_{01}}{\partial t}$
$\displaystyle=(i\omega_{s}-\Gamma_{2})\rho_{01}-i\gamma\mu_{0}h_{z}\rho_{01}+\frac{i}{2}\gamma\mu_{0}h_{x}(\rho_{11}-\rho_{00}),$

where $h_{x}=\sum_{j=1,2}h_{j}\sin\theta\cos(\omega_{j}t+\varphi_{j})$ and
$h_{z}=\sum_{j=1,2}h_{j}\cos\theta\cos(\omega_{j}t+\varphi_{j})$ are the
$x$-axis (transverse) and $z$-axis (longitudinal) components of $\bm{h}(t)$
and $\Gamma_{2}=\Gamma_{2}^{0}+\Gamma_{\mathrm{p}}/2$ is the effective
transverse spin relaxation rate. Other two elements in $\rho$ are determined
by the constraints of $\rho_{10}=\rho_{01}^{*}$ and $\rho_{00}+\rho_{11}=1$.
We treat $V(t)$ as a perturbation to $H_{0}$ and solve for
$\rho=\sum_{n=0}^{\infty}\rho^{(n)}$, where $\rho^{(n)}$ is in the $n$th order
of $V(t)$ and is the $n$th-order correction to the zero-order solution
$\rho^{(0)}$ given by $\rho_{11}^{(0)}=\Gamma_{1}^{0}/(2\Gamma_{1})$ with the
effective longitudinal spin relaxation rate
$\Gamma_{1}=\Gamma_{1}^{0}+\Gamma_{\mathrm{p}}$ and $\rho_{01}^{(0)}=0$.
Specifically, we solve for $\rho^{(n)}$ ($n\geq 1$) using iterative equations

$\displaystyle\frac{\partial{\rho}_{11}^{(n)}}{\partial t}=$
$\displaystyle-\Gamma_{1}\rho_{11}^{(n)}+\frac{i}{2}\gamma\mu_{0}h_{x}\left(\rho_{01}^{(n-1)}-\rho_{10}^{(n-1)}\right),$
(15) $\displaystyle\frac{\partial{\rho}_{01}^{(n)}}{\partial t}=$
$\displaystyle(i\omega_{s}-\Gamma_{2})\rho_{01}^{(n)}-i\gamma\mu_{0}h_{z}\rho_{01}^{(n-1)}+$
$\displaystyle\frac{i}{2}\gamma\mu_{0}h_{x}\left(\rho_{11}^{(n-1)}-\rho_{00}^{(n-1)}\right).$

The first-order solution of $\rho$ is

$\displaystyle\rho_{11}^{(1)}$ $\displaystyle=0,$ (16)
$\displaystyle\rho_{01}^{(1)}$
$\displaystyle=-\frac{\Gamma_{\mathrm{p}}\gamma\mu_{0}}{2\Gamma_{1}}\sum_{\omega_{m}}\frac{\tilde{h}_{x}(\omega_{m})e^{i\omega_{m}t}}{\omega_{m}-\omega_{s}-i\Gamma_{2}}.$

Here, $\tilde{h}_{x}(\omega)$ and $\tilde{h}_{z}(\omega)$ are the Fourier
transforms of $h_{x}$ and $h_{z}$, which are only finite at
$\tilde{h}_{x}(\pm\omega_{j})=h_{j}\sin\theta e^{\pm i\varphi_{j}}/2$ and
$\tilde{h}_{z}(\pm\omega_{j})=h_{j}\cos\theta e^{\pm i\varphi_{j}}/2$ where
$j=1,2$. We define the first-order susceptibility
$\chi_{xx}^{(1)}(\omega_{j},\omega_{j})$ through
$\tilde{M}_{x}(\omega_{j})=\chi_{xx}^{(1)}(\omega_{j},\omega_{j})\tilde{h}_{x}(\omega_{j})$,
where
$M_{x}=\mathrm{Tr}(\rho\cdot\gamma\frac{\hbar}{2}\sigma_{x})=\frac{1}{2}\gamma\hbar(\rho_{01}+\rho_{10})$
is the $x$-axis component of magnetic moment. It can be expressed as

$\displaystyle\chi_{xx}^{(1)}(\omega_{j},\omega_{j})\approx-\frac{\Gamma_{\mathrm{p}}\gamma^{2}\mu_{0}\hbar}{4\Gamma_{1}(\omega_{j}-\omega_{s}-i\Gamma_{2})},$
(17)

when $|\omega_{j}-\omega_{s}|\ll|\omega_{j}+\omega_{s}|$ is satisfied
($j=1,2$).

The second-order solution of $\rho$ is

$\displaystyle\rho_{11}^{(2)}$
$\displaystyle=-\frac{\Gamma_{\mathrm{p}}\gamma^{2}\mu_{0}^{2}}{4\Gamma_{1}}\sum_{\omega_{m},\omega_{n}}\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{x}(\omega_{n})e^{i(\omega_{m}+\omega_{n})t}}{(\omega_{m}+\omega_{n}-i\Gamma_{1})(\omega_{m}-\omega_{s}-i\Gamma_{2})}+\mathrm{H.c.},$
(18) $\displaystyle\rho_{01}^{(2)}$
$\displaystyle=\frac{\Gamma_{\mathrm{p}}\gamma^{2}\mu_{0}^{2}}{2\Gamma_{1}}\sum_{\omega_{m},\omega_{n}}\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{z}(\omega_{n})e^{i(\omega_{m}+\omega_{n})t}}{(\omega_{m}+\omega_{n}-\omega_{s}-i\Gamma_{2})(\omega_{m}-\omega_{s}-i\Gamma_{2})}.$

We define $\chi_{xxz}^{(2)}(\omega_{2}\pm\omega_{1},\omega_{1},\omega_{2})$
through
$\tilde{M}_{x}(\omega_{2}\pm\omega_{1})=\chi_{xxz}^{(2)}(\omega_{2}\pm\omega_{1},\omega_{1},\omega_{2})\tilde{h}_{x}(\omega_{1})\tilde{h}_{z}(\omega_{2})$
and obtain

$\displaystyle\chi_{xxz}^{(2)}(\omega_{2}\pm\omega_{1},\omega_{1},\omega_{2})\approx$
$\displaystyle\mp\frac{\Gamma_{\mathrm{p}}\gamma^{3}\mu_{0}^{2}\hbar\omega_{s}}{4\Gamma_{1}(\omega_{2}\pm\omega_{1}-\omega_{s}-i\Gamma_{2})\omega_{1}\omega_{2}},$
(19)

when
$|\omega_{2}\pm\omega_{1}-\omega_{s}|\ll|\omega_{2}\pm\omega_{1}+\omega_{s}|$
and $\omega_{1,2}\gg\Gamma_{2}$ are satisfied. The third-order solution of
$\rho$ is

$\displaystyle\rho_{11}^{(3)}=$
$\displaystyle\frac{\Gamma_{\mathrm{p}}\gamma^{3}\mu_{0}^{3}}{4\Gamma_{1}}\sum_{\omega_{m},\omega_{n},\omega_{k}}\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{z}(\omega_{n})\tilde{h}_{x}(\omega_{k})e^{i(\omega_{m}+\omega_{n}+\omega_{k})t}}{(\omega_{m}+\omega_{n}-\omega_{s}-i\Gamma_{2})(\omega_{m}-\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}-i\Gamma_{1})}+\mathrm{H.c.},$
(20) $\displaystyle\rho_{01}^{(3)}=-$
$\displaystyle\frac{\Gamma_{\mathrm{p}}\gamma^{3}\mu_{0}^{3}}{4\Gamma_{1}}\sum_{\omega_{m},\omega_{n},\omega_{k}}\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{x}(\omega_{n})\tilde{h}_{x}(\omega_{k})e^{i(\omega_{m}+\omega_{n}+\omega_{k})t}}{(\omega_{m}+\omega_{n}-i\Gamma_{1})(\omega_{m}+\omega_{n}+\omega_{k}-\omega_{s}-i\Gamma_{2})}\left(\frac{1}{\omega_{m}-\omega_{s}-i\Gamma_{2}}+\frac{1}{\omega_{m}+\omega_{s}-i\Gamma_{2}}\right)$
$\displaystyle-$
$\displaystyle\frac{\Gamma_{\mathrm{p}}\gamma^{3}\mu_{0}^{3}}{2\Gamma_{1}}\sum_{\omega_{m},\omega_{n},\omega_{k}}\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{z}(\omega_{n})\tilde{h}_{z}(\omega_{k})e^{i(\omega_{m}+\omega_{n}+\omega_{k})t}}{(\omega_{m}+\omega_{n}-\omega_{s}-i\Gamma_{2})(\omega_{m}-\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}-\omega_{s}-i\Gamma_{2})}.$

We define $\chi_{xxxx}^{(3)}(\omega_{2},\omega_{2},-\omega_{1},\omega_{1})$
through
$\tilde{M}_{x}(\omega_{2})=\chi_{xxxx}^{(3)}(\omega_{2},\omega_{2},-\omega_{1},\omega_{1})\tilde{h}_{x}(\omega_{2})\tilde{h}_{x}(-\omega_{1})\tilde{h}_{x}(\omega_{1})$
and obtain

$\displaystyle\chi_{xxxx}^{(3)}(\omega_{2},\omega_{2},-\omega_{1},\omega_{1})\approx\frac{i\gamma^{4}\mu_{0}^{3}\hbar\Gamma_{\mathrm{p}}\Gamma_{2}}{4\Gamma_{1}(\omega_{2}-\omega_{1}-i\Gamma_{1})(\omega_{2}-\omega_{s}-i\Gamma_{2})^{2}(-\omega_{1}+\omega_{s}-i\Gamma_{2})},$
(21)

when $|\omega_{1,2}-\omega_{s}|\ll|\omega_{1,2}+\omega_{s}|$ and
$|\omega_{1}-\omega_{2}|\ll\Gamma_{2}$ are satisfied. The fourth-order
solution of $\rho_{11}^{(4)}$ is

$\displaystyle\rho_{11}^{(4)}=-$
$\displaystyle\frac{\Gamma_{\mathrm{p}}\gamma^{4}\mu_{0}^{4}}{8\Gamma_{1}}\sum_{\omega_{m},\omega_{n},\omega_{k},\omega_{l}}\left[\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{x}(\omega_{n})\tilde{h}_{x}(\omega_{k})\tilde{h}_{x}(\omega_{l})e^{i(\omega_{m}+\omega_{n}+\omega_{k}+\omega_{l})t}}{(\omega_{m}+\omega_{n}-i\Gamma_{1})(\omega_{m}-\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}-\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}+\omega_{l}-i\Gamma_{1})}+\mathrm{H.c.}\right]$
(22) $\displaystyle-$
$\displaystyle\frac{\Gamma_{\mathrm{p}}\gamma^{4}\mu_{0}^{4}}{8\Gamma_{1}}\sum_{\omega_{m},\omega_{n},\omega_{k},\omega_{l}}\left[\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{x}(\omega_{n})\tilde{h}_{x}(\omega_{k})\tilde{h}_{x}(\omega_{l})e^{i(\omega_{m}+\omega_{n}+\omega_{k}+\omega_{l})t}}{(\omega_{m}+\omega_{n}-i\Gamma_{1})(\omega_{m}+\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}-\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}+\omega_{l}-i\Gamma_{1})}+\mathrm{H.c.}\right]$
$\displaystyle-$
$\displaystyle\frac{\Gamma_{\mathrm{p}}\gamma^{4}\mu_{0}^{4}}{4\Gamma_{1}}\sum_{\omega_{m},\omega_{n},\omega_{k},\omega_{l}}\left[\frac{\tilde{h}_{x}(\omega_{m})\tilde{h}_{z}(\omega_{n})\tilde{h}_{z}(\omega_{k})\tilde{h}_{x}(\omega_{l})e^{i(\omega_{m}+\omega_{n}+\omega_{k}+\omega_{l})t}}{(\omega_{m}+\omega_{n}-\omega_{s}-i\Gamma_{2})(\omega_{m}-\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}-\omega_{s}-i\Gamma_{2})(\omega_{m}+\omega_{n}+\omega_{k}+\omega_{l}-i\Gamma_{1})}+\mathrm{H.c.}\right].$

##  Appendix C Origin of resonance at the sum or difference frequency of two
driving microwaves

The PL intensity is given by $I=I_{0}(1-\alpha\braket{\rho_{11}})$, where
$\braket{\rho_{11}}$ is the steady-state solution of $\rho_{11}$, $I_{0}$ is
the PL intensity in the case of full spin initialization (i.e.,
$\braket{\rho_{11}}=0$), and $\alpha$ is a phenomenological parameter to
account for the difference in the contribution of the $\ket{m_{s}=0}$
population and that of $\ket{m_{s}=\pm 1}$ population to PL intensity [44]. We
note that the fast-oscillating terms of $\rho_{11}$ will not contribute to the
PL intensity $I$.

For standard, linear magnetic resonance, the nonoscillating part in
$\rho_{11}^{(2)}$ in Eq. (18) is given by

$\displaystyle\braket{\rho_{11}^{(2)}}$
$\displaystyle\approx-\frac{\mu_{0}}{i\hbar\Gamma_{1}}\sum_{j=1,2}|\tilde{h}_{x}(\omega_{j})|^{2}\chi_{xx}^{(1)}(\omega_{j},\omega_{j})+\mathrm{H.c.}$
(23)
$\displaystyle=\frac{\gamma^{2}\mu_{0}^{2}\Gamma_{\mathrm{p}}\Gamma_{2}\sin^{2}\theta}{8\Gamma_{1}^{2}}\sum_{j=1,2}\frac{h_{j}^{2}}{(\omega_{j}-\omega_{s})^{2}+\Gamma_{2}^{2}},$

when $|\omega_{1,2}-\omega_{s}|\ll|\omega_{1,2}+\omega_{s}|$ are satisfied.
Equation (23) demonstrates the commonly observed resonance under a single
microwave frequency $\omega_{1}$ or $\omega_{2}$ close to $\omega_{s}$.

Next we consider nonlinear, two-photon magnetic resonance. The third term on
the right-hand side of Eq. (22) gives rise to the observed resonance when the
sum or difference frequency of $\omega_{1}$ and $\omega_{2}$ matches
$\omega_{\mathrm{A}}$. When $\omega_{1,2}$ are individually far away from
$\omega_{\mathrm{A}}$ but their sum or difference is around the resonance
condition of $\omega_{2}\pm\omega_{1}=\omega_{\mathrm{A}}$, we have

$\braket{\rho_{11}^{(4)}}\approx\frac{|\chi_{xxz}^{(2)}(\omega_{2}\pm\omega_{1},\omega_{1},\omega_{2})|^{2}h_{1}^{2}h_{2}^{2}\sin^{2}(2\theta)\Gamma_{\mathrm{2,A}}}{8\gamma^{2}\hbar^{2}\Gamma_{\mathrm{p}}},$
(24)

where $\chi_{xxz}^{(2)}$ is in the expression of Eq. (19), with $\omega_{s}$,
$\Gamma_{1}$, and $\Gamma_{2}$ replaced by $\omega_{\mathrm{A}}$,
$\Gamma_{\mathrm{1,A}}$, and $\Gamma_{\mathrm{2,A}}$, and we have utilized the
approximations of $\omega_{2}\pm\omega_{1}\approx\omega_{\mathrm{A}}$ and
$\omega_{1,2}\gg\Gamma_{\mathrm{2,A}}$. After averaging $\theta$ over
$[0,\pi]$ (note that the detected NV centers in the ensemble have randomly
oriented principal axes), the change of the PL intensity compared to that
without microwaves, defined as $\Delta
I=I_{0}-I_{0}(1-\alpha\braket{\rho_{11}^{(0)}})$, can be expressed as

$\displaystyle\Delta I$ $\displaystyle\approx-|\Delta
I|_{\mathrm{max}}\braket{\rho_{11}^{(4)}}$ (25)
$\displaystyle=-\eta|\chi_{xxz}^{(2)}(\omega_{2}\pm\omega_{1},\omega_{1},\omega_{2})|^{2}h_{1}^{2}h_{2}^{2},$

where $|\Delta I|_{\mathrm{max}}=I_{0}\alpha$ is the difference of the PL
intensity between the case of full spin initialization with
$\braket{\rho_{11}}=0$ and the case of full spin inversion with
$\braket{\rho_{11}}=1$ ($|\Delta I|_{\mathrm{max}}=I_{0}\alpha\epsilon$ with
modulation depth $\epsilon$ if we consider amplitude-modulated measurements)
and

$\eta=\frac{|\Delta
I|_{\mathrm{max}}\Gamma_{\mathrm{2,A}}}{16\gamma^{2}\hbar^{2}\Gamma_{\mathrm{p}}}$
(26)

is a constant factor. We see that $\Delta I$ in Eq. (25) is a Lorentzian with
a line width of $\Gamma_{\mathrm{2,A}}$ and a contrast of

$C=\frac{|\Delta
I|_{\mathrm{max}}\Gamma_{\mathrm{p}}\omega_{\mathrm{A}}^{2}\gamma^{4}\mu_{0}^{4}h_{1}^{2}h_{2}^{2}}{256\Gamma_{\mathrm{1,A}}^{2}\Gamma_{\mathrm{2,A}}\omega_{1}^{2}\omega_{2}^{2}}.$
(27)

Since $\Gamma_{\mathrm{2,E}}$ is much larger than $\Gamma_{\mathrm{2,A}}$, as
verified by the broad resonance dip at $\omega_{2}=\omega_{\mathrm{E}}$ [see
Figs. 2(a) and 7], the resonance at the sum or difference frequency of
$\omega_{1,2}$ is significant in ${}^{3}A_{2}$ but not observed in ${}^{3}E$.

Up to now, the line width of the sum- or difference-frequency resonance, i.e.,
the effective transverse spin relaxation rate, has been determined to be
$\Gamma_{\mathrm{2,A}}=\Gamma_{\mathrm{2,A}}^{0}+\Gamma_{\mathrm{p}}/2$ for a
single NV center or an ensemble of totally identical NV centers. In reality,
however, NV centers in the ensemble are different, since they are situated in
varying local environments, which leads to inhomogeneous broadening. We can
model this by assuming that the transition frequency $\omega_{\mathrm{A}}$ is
distributed based on the probability distribution function
$p(\omega_{\mathrm{A}})$. Then, the detected change of the PL intensity
averaged on the ensemble is $\langle\Delta I\rangle=\int_{0}^{\infty}\Delta
I\cdot p(\omega_{\mathrm{A}})d\omega_{\mathrm{A}}$, where $\Delta I$ is given
by Eq. (25). Typically, it is assumed that $p(\omega_{\mathrm{A}})$ is
Gaussian [52] and $\langle\Delta I\rangle$ is in the Voigt line shape, the
line width of which can be calculated by some complex error functions. Here,
we simply assume $p(\omega_{\mathrm{A}})$ is a Lorentzian, with a line width
of $\Gamma_{\mathrm{inh}}$. In this case, $\langle\Delta I\rangle$ is also a
Lorentzian, with a total line width of
$\Gamma_{\mathrm{2,A}}=\Gamma_{\mathrm{2,A}}^{0}+\Gamma_{\mathrm{p}}/2+\Gamma_{\mathrm{inh}}$.
Since $\Delta I$ around the resonance condition in our experiment is almost a
perfect Lorentzian [see Fig. 3(a)], this assumption is reasonable.

##  Appendix D Analysis on higher-order resonance dips

In Fig. 2(a), we observe additional bright lines (labeled from 4 to 7), which
can be traced to magnetic resonance excited by even higher-order processes.
Here, to examine whether our theory generally applies to higher-order
resonance dips, we select the three-photon resonance at
$2\omega_{2}-\omega_{1}=\omega_{\mathrm{A}}$ for a power- and frequency-
dependent test. From Eq. (20), the third-order susceptibility
$\chi_{xxxx}^{(3)}(2\omega_{2}-\omega_{1},\omega_{2},\omega_{2},-\omega_{1})$,
defined through

$\tilde{M}_{x}(2\omega_{2}-\omega_{1})=\chi_{xxxx}^{(3)}(2\omega_{2}-\omega_{1},\omega_{2},\omega_{2},-\omega_{1})\tilde{h}_{x}(\omega_{2})\tilde{h}_{x}(\omega_{2})\tilde{h}_{x}(-\omega_{1})$
(28)

and
$\chi_{xxzz}^{(3)}(2\omega_{2}-\omega_{1},\omega_{2},\omega_{2},-\omega_{1})$,
defined through

$\tilde{M}_{x}(2\omega_{2}-\omega_{1})=\chi_{xxzz}^{(3)}(2\omega_{2}-\omega_{1},\omega_{2},\omega_{2},-\omega_{1})\tilde{h}_{x}(\omega_{2})\tilde{h}_{z}(\omega_{2})\tilde{h}_{z}(-\omega_{1})$
(29)

are given by

$\displaystyle\chi_{xxxx}^{(3)}(2\omega_{2}-\omega_{1},\omega_{2},\omega_{2},-\omega_{1})$
$\displaystyle\approx\frac{i\gamma^{4}\mu_{0}^{3}\hbar\Gamma_{\mathrm{p}}\Gamma_{2}\omega_{s}}{8\Gamma_{1}\omega_{2}(\omega_{s}-\omega_{2})(\omega_{s}-\omega_{1})(2\omega_{2}-\omega_{1}-\omega_{s}-i\Gamma_{2})},$
(30)
$\displaystyle\chi_{xxzz}^{(3)}(2\omega_{2}-\omega_{1},\omega_{2},\omega_{2},-\omega_{1})$
$\displaystyle\approx\frac{i\gamma^{4}\mu_{0}^{3}\hbar\Gamma_{\mathrm{p}}\Gamma_{2}}{4\Gamma_{1}\omega_{1}\omega_{2}(2\omega_{2}-\omega_{1}-\omega_{s}-i\Gamma_{2})},$

respectively, when the near-resonance condition of
$2\omega_{2}-\omega_{1}\approx\omega_{s}$ is satisfied. It is noted that
$\chi_{xxxx}^{(3)}$ represents the spin transition from $\ket{m_{s}=0}$ to
$\ket{m_{s}=\pm 1}$ via $\sigma_{\pm}+\sigma_{\pm}-\sigma_{\pm}$. The
transition rate is proportional to
$(h_{2x}^{2}h_{1x}\omega_{s}/[\omega_{2}(\omega_{s}-\omega_{2})(\omega_{s}-\omega_{1})]$,
where $h_{jx}=h_{j}\cos\theta_{j}$ and $h_{jz}=h_{j}\sin\theta_{j}$ ($j=1,2$).
On the other hand, $\chi_{xxzz}^{(3)}$ represents the spin transition via
$\sigma^{\pm}+\pi+\pi$. The transition rate is proportional to
$h_{2x}h_{2z}h_{1x}/(\omega_{1}\omega_{2})$. For a general case with
$\theta\neq 0,\pi/2,\pi$, both of the two kinds of three-photon spin
transitions will contribute to the detected resonance signals. Note that in
the visited frequency ranges in Fig. 3(h), both of $\omega_{1,2}$ are close to
$\omega_{\mathrm{A}}$. Consequently,
$\omega_{\mathrm{A}}/[\omega_{2}(\omega_{\mathrm{A}}-\omega_{2})(\omega_{\mathrm{A}}-\omega_{1})]$,
the frequency dependence from $\chi_{xxxx}^{(3)}$, is more significant than
$1/(\omega_{1}\omega_{2})$, the frequency dependence from $\chi_{xxzz}^{(3)}$,
which is demonstrated by dashed lines in Fig. 3(h). The semiquantitative
consistency between the experimental data on the third-order resonance and our
theoretical model suggests the generality of the model to higher-order
resonances.

Figure 9: $\Delta I$ under $\omega_{1}$ and $\omega_{2}$ inputs when
$\omega_{1,2}$ are close to $\omega_{\mathrm{E}}$: (a) microwave powers
$P_{1}=5\ \mathrm{mW}$, $P_{2}=13\ \mathrm{mW}$; (b) microwave powers
$P_{1}=13\ \mathrm{mW}$, $P_{2}=5\ \mathrm{mW}$. Only (b) shows the dark line
satisfying $\omega_{1}=\omega_{2}$, which is a feature of EIT.

##  Appendix E Origin of electromagnetically induced transparency

In Figs. 2(a) and 7, on top of the series of the bright resonance lines, we
observe a dark line satisfying $\omega_{1}=\omega_{2}$ within the broad
resonance dip at $\omega_{2}=\omega_{\mathrm{E}}$. This feature of magnetic
resonance suppression under zero detuning of two waves is very similar to the
EIT phenomenon studied in nonlinear optics [45], where in the presence of a
strong pump wave, the interaction between the probe wave and the matter is
minimized due to interference effects. Mathematically, the first two terms on
the right-hand side of Eq. (22) play an important role. Combining Eqs. (22)
and (23), when $\omega_{1}\approx\omega_{2}$ and both of them are close to
$\omega_{\mathrm{E}}$, i.e.,
$|\omega_{1,2}-\omega_{\mathrm{E}}|\ll|\omega_{1,2}+\omega_{\mathrm{E}}|$ and
$|\omega_{1}-\omega_{2}|\ll\Gamma_{\mathrm{2,E}}$ are satisfied, we have

$\displaystyle\braket{\rho_{11}^{(2)}}+\braket{\rho_{11}^{(4)}}\approx-\frac{\mu_{0}|\tilde{h}_{x}(\omega_{2})|^{2}}{i\Gamma_{\mathrm{1,E}}\hbar}\left[\chi_{xx}^{(1)}(\omega_{2},\omega_{2})+|\tilde{h}_{x}(\omega_{1})|^{2}\chi_{xxxx}^{(3)}(\omega_{2},\omega_{2},-\omega_{1},\omega_{1})\right]+\mathrm{H.c.}+(\omega_{1}\Longleftrightarrow\omega_{2}).$
(31)

where $\chi_{xx}^{(1)}$ and $\chi_{xxxx}^{(3)}$ are in the expressions of Eqs.
(17) and (21), with $\omega_{s}$, $\Gamma_{1}$, and $\Gamma_{2}$ replaced by
$\omega_{\mathrm{E}}$, $\Gamma_{\mathrm{1,E}}$, and $\Gamma_{\mathrm{2,E}}$.
After averaging $\theta$ over $[0,\pi]$ and only keeping terms related to
$h_{2}$, since only the amplitude of the $\omega_{2}$ input is modulated, the
destructive interference between the spin transitions associated with
$\chi_{xx}^{(1)}$ and $\chi_{xxxx}^{(3)}$ results in

$\displaystyle\Delta I$ $\displaystyle\approx-|\Delta
I|_{\mathrm{max}}\left[\braket{\rho_{11}^{(2)}}+\braket{\rho_{11}^{(4)}}\right]$
(32)
$\displaystyle=\frac{iA}{\omega_{2}-\omega_{\mathrm{E}}-i\Gamma_{\mathrm{2,E}}-\frac{B}{\omega_{2}-\omega_{1}-i\Gamma_{\mathrm{1,E}}}}+\mathrm{H.c.}$

where $A=|\Delta
I|_{\mathrm{max}}\Gamma_{\mathrm{p}}\gamma^{2}\mu_{0}^{2}h_{2}^{2}/(32\Gamma_{\mathrm{1,E}}^{2})$
and
$B=i\Gamma_{\mathrm{2,E}}\gamma^{2}\mu_{0}^{2}h_{1}^{2}/[4(\omega_{1}-\omega_{\mathrm{E}}+i\Gamma_{\mathrm{2,E}})]\approx\gamma^{2}\mu_{0}^{2}h_{1}^{2}/4$.
We see that the peak value of $|\Delta I|$ is suppressed at zero detuning
$\omega_{2}-\omega_{1}=0$, and that the EIT effect is most significant when
$\Gamma_{\mathrm{2,E}}>\gamma\mu_{0}h_{1}/2>\Gamma_{\mathrm{1,E}}$, which is
satisfied in the spin transitions in ${}^{3}E$. Comparatively, the relatively
smaller $\Gamma_{\mathrm{2,A}}$ makes the EIT feature in ${}^{3}A_{2}$ less
noticeable. We note that here, we treat the $\omega_{1}$ signal as the probe
wave and the modulated $\omega_{2}$ signal as the pump wave. When using a weak
pump power and a high probe power [$P_{1}=5\ \mathrm{mW}$, $P_{2}=13\
\mathrm{mW}$; Fig. 9(a)], the transparency window is not observed. In
contrast, when using a high pump power and a weak probe power [$P_{1}=13\
\mathrm{mW}$, $P_{2}=5\ \mathrm{mW}$; Fig. 9(b)], the transparency window
clearly exists. This asymmetric result is consistent with Eq. (32), further
verifying the EIT origin of the resonance suppression feature.

The derivation from Appendixes B–E applies to the two-level model for the
transition between $\ket{m_{s}=0}$ and $\ket{m_{s}=+1}$. When we focus on the
transition between $\ket{m_{s}=0}$ and $\ket{m_{s}=-1}$, the treatment is
similar. The only difference in the derivation is to replace $h_{z}$ by
$-h_{z}$ due to the negative magnetic moment, which will add an extra negative
sign in the results of density-matrix elements that are odd functions of
$h_{z}$ (e.g., $\rho_{01}^{(2)}$ and $\rho_{11}^{(3)}$). The elements
$\rho_{11}^{(2)}$, $\rho_{11}^{(4)}$, and hence $\Delta I$ will not be
affected.

Figure 10: Numerically simulated results demonstrating (a) the sum- and
difference-frequency resonances at
$\omega_{2}\pm\omega_{1}=\omega_{\mathrm{A}}$ and (b) EIT effect when
$\omega_{1,2}$ are close to $\omega_{\mathrm{E}}$.

##  Appendix F Numerical simulations on the density-matrix master equation

Besides the theoretical derivation based on perturbation theory, we can also
solve Eq. (14) numerically. In Figs. 10(a) and 10(b), we show the numerically
simulated results for the sum- and difference-frequency resonances at
$\omega_{2}\pm\omega_{1}=\omega_{\mathrm{A}}$ and the EIT effect when
$\omega_{1,2}$ are close to $\omega_{\mathrm{E}}$, respectively. The
parameters used in the simulations are listed below.

Figure 10(a): $\omega_{\mathrm{A}}/(2\pi)=2.87\ \mathrm{GHz}$,
$\Gamma_{\mathrm{1,A}}^{0}/(2\pi)=1\ \mathrm{kHz}$,
$\Gamma_{\mathrm{2,A}}^{0}/(2\pi)=5\ \mathrm{MHz}$,
$\Gamma_{\mathrm{p}}/(2\pi)=3\ \mathrm{MHz}$, $h_{1,2}=20\ \mathrm{Oe}$.

Figure 10(b): $\omega_{\mathrm{E}}/(2\pi)=1.42\ \mathrm{GHz}$,
$\Gamma_{\mathrm{1,E}}^{0}/(2\pi)=1\ \mathrm{kHz}$,
$\Gamma_{\mathrm{2,E}}^{0}/(2\pi)=100\ \mathrm{MHz}$,
$\Gamma_{\mathrm{p}}/(2\pi)=3\ \mathrm{MHz}$, $h_{1,2}=20\ \mathrm{Oe}$.

Here, $\Gamma_{\mathrm{2,A}}^{0}$ and $\Gamma_{\mathrm{2,E}}^{0}$ include the
contribution of ensemble inhomogeneous broadening. $\theta$ is set as
$45^{\circ}$ to reflect the effect of the ensemble average. $|\Delta
I|_{\mathrm{max}}$ is set as 1. We take a time step of $\Delta t=2\times
10^{-13}\ \mathrm{s}$ and average $\rho_{11}(t)$ in $4\times 10^{-7}\
\mathrm{s}\leq t\leq 8\times 10^{-7}\ \mathrm{s}$ to obtain steady-state
solutions. The results with the same parameters except $h_{2}=0$ are
subtracted to imitate the modulation on $h_{2}$. In Figs. 10(a) and 10(b), we
see that the numerical simulations are consistent with the experimental data
in clearly demonstrating features induced by nonlinear multiphoton processes.

## References

  * Shen [1984] Y. R. Shen,  _The Principles of Nonlinear Optics_ (Wiley, New York, 1984).
  * Boyd [2008] R. W. Boyd,  _Nonlinear Optics_ , 3rd ed. (Academic Press, Boston, 2008).
  * Brabec and Krausz [2000] T. Brabec and F. Krausz, Intense few-cycle laser fields: Frontiers of nonlinear optics, Rev. Mod. Phys. 72, 545 (2000).
  * Reshef  _et al._ [2019] O. Reshef, I. De Leon, M. Z. Alam, and R. W. Boyd, Nonlinear optical effects in epsilon-near-zero media, Nat. Rev. Mater. 4, 535 (2019).
  * Fejer [1994] M. M. Fejer, Nonlinear Optical Frequency Conversion, Phys. Today 47, 25 (1994).
  * Caspani  _et al._ [2011] L. Caspani, D. Duchesne, K. Dolgaleva, S. J. Wagner, M. Ferrera, L. Razzari, A. Pasquazi, M. Peccianti, D. J. Moss, J. S. Aitchison, and R. Morandotti, Optical frequency conversion in integrated devices, J. Opt. Soc. Am. B 28, A67 (2011).
  * Guo  _et al._ [2016] X. Guo, C.-L. Zou, H. Jung, and H. X. Tang, On-Chip Strong Coupling and Efficient Frequency Conversion between Telecom and Visible Optical Modes, Phys. Rev. Lett. 117, 123902 (2016).
  * Li  _et al._ [2016] Q. Li, M. Davanço, and K. Srinivasan, Efficient and low-noise single-photon-level frequency conversion interfaces using silicon nanophotonics, Nat. Photon. 10, 406 (2016).
  * Wang  _et al._ [2021] J.-Q. Wang, Y.-H. Yang, M. Li, X.-X. Hu, J. B. Surya, X.-B. Xu, C.-H. Dong, G.-C. Guo, H. X. Tang, and C.-L. Zou, Efficient Frequency Conversion in a Degenerate ${\chi}^{(2)}$ Microresonator, Phys. Rev. Lett. 126, 133601 (2021).
  * Radic  _et al._ [2003] S. Radic, C. McKinstrie, R. Jopson, J. Centanni, and A. Chraplyvy, All-optical regeneration in one- and two-pump parametric amplifiers using highly nonlinear optical fiber, IEEE Photon. Technol. Lett. 15, 957 (2003).
  * Li  _et al._ [2015] G. Li, S. Chen, N. Pholchai, B. Reineke, P. W. H. Wong, E. B. Pun, K. W. Cheah, T. Zentgraf, and S. Zhang, Continuous control of the nonlinearity phase for harmonic generations, Nat. Mater. 14, 607 (2015).
  * Reimer  _et al._ [2015] C. Reimer, M. Kues, L. Caspani, B. Wetzel, P. Roztocki, M. Clerici, Y. Jestin, M. Ferrera, M. Peccianti, A. Pasquazi, B. E. Little, S. T. Chu, D. J. Moss, and R. Morandotti, Cross-polarized photon-pair generation and bi-chromatically pumped optical parametric oscillation on a chip, Nat. Commun. 6, 8236 (2015).
  * Mosca  _et al._ [2018] S. Mosca, M. Parisi, I. Ricciardi, F. Leo, T. Hansson, M. Erkintalo, P. Maddaloni, P. De Natale, S. Wabnitz, and M. De Rosa, Modulation Instability Induced Frequency Comb Generation in a Continuously Pumped Optical Parametric Oscillator, Phys. Rev. Lett. 121, 093903 (2018).
  * Loudon and Knight [1987] R. Loudon and P. Knight, Squeezed Light, J. Mod. Opt. 34, 709 (1987).
  * Suhara [2009] T. Suhara, Generation of quantum-entangled twin photons by waveguide nonlinear-optic devices, Laser Photonics Rev. 3, 370 (2009).
  * Quesada and Sipe [2015] N. Quesada and J. E. Sipe, Time-Ordering Effects in the Generation of Entangled Photons Using Nonlinear Optical Processes, Phys. Rev. Lett. 114, 093903 (2015).
  * Caspani  _et al._ [2017] L. Caspani, C. Xiong, B. J. Eggleton, D. Bajoni, M. Liscidini, M. Galli, R. Morandotti, and D. J. Moss, Integrated sources of photon quantum states based on nonlinear optics, Light Sci. Appl. 6, e17100 (2017).
  * Buckholtz and Yavuz [2020] Z. N. Buckholtz and D. D. Yavuz, Linear and nonlinear crystal optics using the magnetic field of light, Phys. Rev. A 101, 023831 (2020).
  * Orton  _et al._ [1960] J. W. Orton, P. Auzins, and J. E. Wertz, Double-Quantum Electron Spin Resonance Transitions of Nickel in Magnesium Oxide, Phys. Rev. Lett. 4, 128 (1960).
  * Clerjaud and Gelineau [1982] B. Clerjaud and A. Gelineau, Observation of Electron Paramagnetic Resonances at Multiples of the ”Classical” Resonance Magnetic Field, Phys. Rev. Lett. 48, 40 (1982).
  * Zur  _et al._ [1983] Y. Zur, M. H. Levitt, and S. Vega, Multiphoton NMR spectroscopy on a spin system with I=1/2, J. Chem. Phys. 78, 5293 (1983).
  * Kälin  _et al._ [2006] M. Kälin, M. Fedin, I. Gromov, and A. Schweiger, Multiple-Photon Transitions in EPR Spectroscopy, in  _Novel NMR and EPR techniques_, edited by J. Dolinšek, M. Vilfan, and S. Žumer (Springer, Berlin, Heidelberg, 2006) pp. 143–183.
  * Taylor  _et al._ [2008] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth, and M. D. Lukin, High-sensitivity diamond magnetometer with nanoscale resolution, Nat. Phys. 4, 810 (2008).
  * Balasubramanian  _et al._ [2008] G. Balasubramanian, I. Y. Chan, R. Kolesov, M. Al-Hmoud, J. Tisler, C. Shin, C. Kim, A. Wojcik, P. R. Hemmer, A. Krueger, T. Hanke, A. Leitenstorfer, R. Bratschitsch, F. Jelezko, and J. Wrachtrup, Nanoscale imaging magnetometry with diamond spins under ambient conditions, Nature (London) 455, 648 (2008).
  * Acosta  _et al._ [2009] V. M. Acosta, E. Bauch, M. P. Ledbetter, C. Santori, K.-M. C. Fu, P. E. Barclay, R. G. Beausoleil, H. Linget, J. F. Roch, F. Treussart, S. Chemerisov, W. Gawlik, and D. Budker, Diamonds with a high density of nitrogen-vacancy centers for magnetometry applications, Phys. Rev. B 80, 115202 (2009).
  * Maertz  _et al._ [2010] B. J. Maertz, A. P. Wijnheijmer, G. D. Fuchs, M. E. Nowakowski, and D. D. Awschalom, Vector magnetic field microscopy using nitrogen vacancy centers in diamond, Appl. Phys. Lett. 96, 092504 (2010).
  * Wolf  _et al._ [2015] T. Wolf, P. Neumann, K. Nakamura, H. Sumiya, T. Ohshima, J. Isoya, and J. Wrachtrup, Subpicotesla Diamond Magnetometry, Phys. Rev. X 5, 041001 (2015).
  * Barry  _et al._ [2020] J. F. Barry, J. M. Schloss, E. Bauch, M. J. Turner, C. A. Hart, L. M. Pham, and R. L. Walsworth, Sensitivity optimization for NV-diamond magnetometry, Rev. Mod. Phys. 92, 015004 (2020).
  * Dutt  _et al._ [2007] M. V. G. Dutt, L. Childress, L. Jiang, E. Togan, J. Maze, F. Jelezko, A. S. Zibrov, P. R. Hemmer, and M. D. Lukin, Quantum Register Based on Individual Electronic and Nuclear Spin Qubits in Diamond, Science 316, 1312 (2007).
  * Neumann  _et al._ [2010] P. Neumann, R. Kolesov, B. Naydenov, J. Beck, F. Rempp, M. Steiner, V. Jacques, G. Balasubramanian, M. L. Markham, D. J. Twitchen, S. Pezzagna, J. Meijer, J. Twamley, F. Jelezko, and J. Wrachtrup, Quantum register based on coupled electron spins in a room-temperature solid, Nat. Phys. 6, 249 (2010).
  * Togan  _et al._ [2010] E. Togan, Y. Chu, A. S. Trifonov, L. Jiang, J. Maze, L. Childress, M. V. G. Dutt, A. S. Sørensen, P. R. Hemmer, A. S. Zibrov, and M. D. Lukin, Quantum entanglement between an optical photon and a solid-state spin qubit, Nature (London) 466, 730 (2010).
  * Fuchs  _et al._ [2011] G. D. Fuchs, G. Burkard, P. V. Klimov, and D. D. Awschalom, A quantum memory intrinsic to single nitrogen–vacancy centres in diamond, Nat. Phys. 7, 789 (2011).
  * Childress and Hanson [2013] L. Childress and R. Hanson, Diamond NV centers for quantum computing and quantum networks, MRS Bull. 38, 134 (2013).
  * Bernien  _et al._ [2013] H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen, L. Childress, and R. Hanson, Heralded entanglement between solid-state qubits separated by three metres, Nature (London) 497, 86 (2013).
  * Humphreys  _et al._ [2018] P. C. Humphreys, N. Kalb, J. P. J. Morits, R. N. Schouten, R. F. L. Vermeulen, D. J. Twitchen, M. Markham, and R. Hanson, Deterministic delivery of remote entanglement on a quantum network, Nature (London) 558, 268 (2018).
  * Wang  _et al._ [2022] G. Wang, Y.-X. Liu, J. M. Schloss, S. T. Alsid, D. A. Braje, and P. Cappellaro, Sensing of Arbitrary-Frequency Fields Using a Quantum Mixer, Phys. Rev. X 12, 021061 (2022).
  * Doherty  _et al._ [2012] M. W. Doherty, F. Dolde, H. Fedder, F. Jelezko, J. Wrachtrup, N. B. Manson, and L. C. L. Hollenberg, Theory of the ground-state spin of the NV- center in diamond, Phys. Rev. B 85, 205203 (2012).
  * Rogers  _et al._ [2009] L. J. Rogers, R. L. McMurtrie, M. J. Sellars, and N. B. Manson, Time-averaging within the excited state of the nitrogen-vacancy centre in diamond, New J. Phys. 11, 063007 (2009).
  * Jensen  _et al._ [2013] K. Jensen, V. M. Acosta, A. Jarmola, and D. Budker, Light narrowing of magnetic resonances in ensembles of nitrogen-vacancy centers in diamond, Phys. Rev. B 87, 014115 (2013).
  * Rogers  _et al._ [2008] L. J. Rogers, S. Armstrong, M. J. Sellars, and N. B. Manson, Infrared emission of the NV centre in diamond: Zeeman and uniaxial stress studies, New J. Phys. 10, 103024 (2008).
  * Acosta  _et al._ [2010] V. M. Acosta, A. Jarmola, E. Bauch, and D. Budker, Optical properties of the nitrogen-vacancy singlet levels in diamond, Phys. Rev. B 82, 201202 (2010).
  * Choi  _et al._ [2012] S. Choi, M. Jain, and S. G. Louie, Mechanism for optical initialization of spin in NV- center in diamond, Phys. Rev. B 86, 041202 (2012).
  * Kamp  _et al._ [2018] E. J. Kamp, B. Carvajal, and N. Samarth, Continuous wave protocol for simultaneous polarization and optical detection of P1-center electron spin resonance, Phys. Rev. B 97, 045204 (2018).
  * Dréau  _et al._ [2011] A. Dréau, M. Lesik, L. Rondin, P. Spinicelli, O. Arcizet, J.-F. Roch, and V. Jacques, Avoiding power broadening in optically detected magnetic resonance of single nv defects for enhanced dc magnetic field sensitivity, Phys. Rev. B 84, 195204 (2011).
  * Fleischhauer  _et al._ [2005] M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Electromagnetically induced transparency: Optics in coherent media, Rev. Mod. Phys. 77, 633 (2005).
  * Jamonneau  _et al._ [2016] P. Jamonneau, G. Hétet, A. Dréau, J.-F. Roch, and V. Jacques, Coherent population trapping of a single nuclear spin under ambient conditions, Phys. Rev. Lett. 116, 043603 (2016).
  * Mrozek  _et al._ [2016] M. Mrozek, A. M. Wojciechowski, D. S. Rudnicki, J. Zachorowski, P. Kehayias, D. Budker, and W. Gawlik, Coherent population oscillations with nitrogen-vacancy color centers in diamond, Phys. Rev. B 94, 035204 (2016).
  * Sewani  _et al._ [2020] V. K. Sewani, H. H. Vallabhapurapu, Y. Yang, H. R. Firgau, C. Adambukulam, B. C. Johnson, J. J. Pla, and A. Laucht, Coherent control of NV- centers in diamond in a quantum teaching lab, Am. J. Phys. 88, 1156 (2020).
  * Du  _et al._ [2017] C. Du, T. van der Sar, T. X. Zhou, P. Upadhyaya, F. Casola, H. Zhang, M. C. Onbasli, C. A. Ross, R. L. Walsworth, Y. Tserkovnyak, and A. Yacoby, Control and local measurement of the spin chemical potential in a magnetic insulator, Science 357, 195 (2017).
  * Fukami  _et al._ [2021] M. Fukami, D. R. Candido, D. D. Awschalom, and M. E. Flatté, Opportunities for Long-Range Magnon-Mediated Entanglement of Spin Qubits via On- and Off-Resonant Coupling, PRX Quantum 2, 040314 (2021).
  * Koerner  _et al._ [2022] C. Koerner, R. Dreyer, M. Wagener, N. Liebing, H. G. Bauer, and G. Woltersdorf, Frequency multiplication by collective nanoscale spin-wave dynamics, Science 375, 1165 (2022).
  * Dobrovitski  _et al._ [2008] V. V. Dobrovitski, A. E. Feiguin, D. D. Awschalom, and R. Hanson, Decoherence dynamics of a single spin versus spin ensemble, Phys. Rev. B 77, 245212 (2008).