#  Supporting Information
_to:_ 13C Hyperpolarization with Nitrogen-Vacancy Centers in Micro- and
Nanodiamonds for Sensitive Magnetic Resonance Applications

###### Contents

  1. S1 Material
  2. S2 Surface treatment characterization
    1. S2.1 Monitoring surface spin defects by EPR
    2. S2.2 13C $T_{1}$ relaxation
    3. S2.3 NV $T_{1}$ relaxation
    4. S2.4 XPS on smaller nanoparticles ($\sim$25\text{\,}\mathrm{nm}$$)
  3. S3 Sample illumination
    1. S3.1 Illumination protocol
    2. S3.2 Phenomenological description of the NV spectrum in X-band EPR under optical pumping
      1. S3.2.1 Redistribution of state populations
      2. S3.2.2 Spectral fits
    3. S3.3 Heating induced by illumination
    4. S3.4 Transient NV signal under pulsed illumination
  4. S4 Build-up and decay of 13C polarization
    1. S4.1 $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $100\text{\,}\mathrm{nm}$ samples, with sample rotation (25°/s)
    2. S4.2 Buildup for the $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample, without rotation
  5. S5 Evaluation of the hyperpolarization enhancement factor and the absolute polarization
  6. S6 PulsePol parameter tuning
    1. S6.1 Testing the robustness to detuning of the phase-offset PulsePol sequence
    2. S6.2 Comparison of resonance indices $n=3.5$ and $n=4.5$
    3. S6.3 The composite pulses
    4. S6.4 Other parameters adjustment
  7. S7 Potential adaptations for further increase of the bandwidth $\Delta_{\mathrm{pol}}$
  8. S8 NV spin Hamiltonian
    1. S8.1 Spin Hamiltonian diagonalization for NV centers under transverse $B$-field
    2. S8.2 Fraction of NV centers in the bandwidth $\Delta_{\mathrm{pol}}$ ($\approx$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$ to the magnetic field)
    3. S8.3 NV transition dipole and Rabi inhomogeneity
    4. S8.4 NV transition frequencies in diamond powder at different magnetic fields

##  S1 Material

The commercial samples are Pureon AG, type Ib, HPHT, MSY1.5-2.5 and Pureon AG,
type Ib, HPHT, MSY0-0.2 having about $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$
and $100\text{\,}\mathrm{nm}$ median size, respectively. Scanning electron
microscopy images are shown in Fig. S1. The size distribution parameters for
the batches employed in this work are given in Table 1.

Fig. S1: Scanning electron microscopy images of a) $100\text{\,}\mathrm{nm}$ and b) $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ diamond powders, respectively Pureon AG, MSY0-0.2 and MSY1.5-2.5. [HTML]FFFFFF sample | median size ($\text{\,}\mathrm{\SIUnitSymbolMicro m}$) | standard deviation ($\text{\,}\mathrm{\SIUnitSymbolMicro m}$)
---|---|---
Pureon AG, MSY1.5-2.5, type Ib | 1.95 | 0.31
Pureon AG, MSY0-0.2, type Ib | 0.115 | 0.028
Table 1: Size distribution properties for the batches of diamond powder
samples used in the present work.

##  S2 Surface treatment characterization

###  S2.1 Monitoring surface spin defects by EPR

EPR spectra were measured on the commercial $100\text{\,}\mathrm{nm}$ sample
(Pureon AG MSY0-0.2 Type Ib, labelled ‘base material’) and after the several
surface treatment steps summarized in Fig. S2.

In this study, the starting material is the commercial sample, with no NV
centers formed by irradiation and annealing. To characterize the surface
treatments, it is useful to perform two types of EPR measurements:

  * •

A narrow scan around $g=2$ (in a $\sim$100\text{\,}\mathrm{G}$$ range) at low
microwave power. This allows resolving bulk and subsurface impurities [1], as
well as dangling bonds with long relaxation times (e.g. in graphite). The
corresponding acquisitions are shown in Fig. S3.

  * •

An acquisition in a wide spectral range ($6000\text{\,}\mathrm{G}$), at higher
microwave power, in order to resolve impurities with a broad spectra and short
lifetimes. The spectra are shown in Fig. S4.

It is possible to observe the disappearance of a $g\sim 2$ contribution
following the air-oxidation (AO) treatment (Fig. S3a). We attribute this
effect to the removal of dangling bonds in graphite covering the nanodiamonds.
This is associated with a strong change in color of the nanoparticles,
initially black, which turn bright. in the color of the nanoparticles, which
are initially black and become bright after the treatment.

Fig. S2: Description of the surface treatment steps for the
$100\text{\,}\mathrm{nm}$ sample.

The wide-range, high-power spectrum in Fig. S4. reveals also the presence of
several broad contributions in the base material (covering a
$\sim$3000\text{\,}\mathrm{G}$$ range). Although one contribution (labelled
(1) in Fig. S4) is removed both by the AO and TAC treatments, another
important contribution (labelled (2) in Fig. S4) persists after air-oxidation,
but is removed by the triacid (TAC) treatment. Iron constitutes a likely
contaminant, as it is present in the beads used in the manufacturing process,
for milling the material down to the desired particle size. In smaller
nanodiamonds from the same manufacturer that underwent the AO treatment, the
presence of iron is confirmed by XPS (sect. S2.4). Consequently: 1) we
attribute the broad EPR contribution that is removed by the triacid cleaning,
to iron-related spin impurities 2) we hypothesize that the fast fluctuations
of these paramagnetic impurities explain the shortening effect on the $T_{1}$
of 13C and NV (detailed in sections S2.2 and S2.3, respectively).

Fig. S3: CW EPR spectra in a narrow spectral range ($100\text{\,}\mathrm{G}$)
around the $g=2$ resonance, at low microwave power
($P_{\mathrm{MW}}=$0.150\text{\,}\mathrm{\SIUnitSymbolMicro W}$$). The
acquisition frequency is $\nu=$9.84\text{\,}\mathrm{GHz}$$. The spectra are
normalized to the sample mass and $Q$-factor. a) Spectra showing the effect of
the first treatment steps (AO or TAC). The satellites at
$3480\text{\,}\mathrm{G}$ and $3540\text{\,}\mathrm{G}$ correspond to the
hyperfine lines of P1 centers. These satellites are unaffected by surface
treatments, as expected for bulk impurities. The signal in the center
comprises both the P1 central line, and $g\sim 2$ paramagnetic species.
Although a minor change is observed with the triacid (TAC), a strong
modification is observed following air-oxidation (AO), as the latter leads to
the disappearance of one central component that we attribute to the efficient
removal of graphite and associated dangling bonds. b) Effect of the successive
treatment combinations. A small change in the intensity of the central line
after the AO+TAC step, however, the spectrum remains essentially unaffected.
Fig. S4: CW EPR spectra in a broad spectral range
(0-$6000\text{\,}\mathrm{G}$), at high microwave power
($P_{\mathrm{MW}}=$23.7\text{\,}\mathrm{mW}$$). The acquisition frequency is
$\nu=$9.84\text{\,}\mathrm{GHz}$$ and the data is normalized to the sample
mass and $Q$-factor. The spectra are also trimmed vertically in order to
better visualize, specifically, the broad contributions. The first spectrum
shows an impurity contribution (1) removed both by the AO and TAC treatment.
Another contribution, labelled (2), is removed by the TAC step only, and is
associated to paramagnetic species of iron (see text). The contribution (3),
which appears in the half-field region (at $g\sim 4$), might be associated
with impurities of $S>1/2$ multiplicity within the nanodiamonds. After the
last treatment step (AO+TAC+AO), the contribution (3) is likely still present,
but not observed due to longer spectrum averaging and possible magnetic field
or frequency drifts.

###  S2.2 13C $T_{1}$ relaxation

In Fig. S5 are represented ${}^{{}^{13}\mathrm{C}}T_{1}$ acquisitions
performed with conventional NMR at $B=$7.05\text{\,}\mathrm{T}$$, using the
same batch of samples as in the EPR investigation in section S2.1. Acquisition
was made with the NMR saturation-recovery sequence, which consists in
observing the return to thermal polarization after a train of saturating
pulses (here corresponding to 20 $\pi/2$ pulses spaced by
$5\text{\,}\mathrm{ms}$). The experimental data points in Figure S5 consist of
the spectrum integral at each value of the time interval between saturation
and readout. Stretched exponential fits
$s(t)=A\left(1-\exp{-(\frac{t}{T_{\mathrm{1}}})^{\beta}}\right)$ are shown as
dashed lines in Fig. S5 and the fitted parameters are listed in Table 2.

Fig. S5: 13C relaxation for the $100\text{\,}\mathrm{nm}$ particles, at different stages of surface treatment (AO: air-oxidation, TAC: triacid-cleaning), acquired by saturation-recovery with conventional NMR at $B=$7.05\text{\,}\mathrm{T}$$. Dashed lines are stretched exponential recovery fits. [HTML]e6e6e6 $100\text{\,}\mathrm{nm}$ sample |  ${}^{{}^{13}\mathrm{C}}T_{1}$ (s) | stretch exponent $\beta$
---|---|---
Base material | 21(1) | 0.55(3)
TAC | 80(4) | 0.61(3)
AO | 25(2) | 0.54(3)
AO+TAC | 152(13) | 0.63(4)
AO+TAC+AO | 158(37) | 0.63(9)
Table 2:  Fitted parameters for the 13C $T_{1}$ relaxation of the
$100\text{\,}\mathrm{nm}$ sample at different stages of surface treatment
(samples are labelled as in Fig. S5).

###  S2.3 NV $T_{1}$ relaxation

The NV longitudinal relaxation times were measured on the
$100\text{\,}\mathrm{nm}$ samples, that were as well electron-irradiated and
annealed. Here, we compare two treatments (AO, AO+TAC). The air-oxidation (AO)
was done both before and after the electron irradiation and simultaneous
annealing. Triacid cleaning (for the AO+TAC sample) was performed only after
the second air-oxidation. Single exponential fits of the experimental data as
shown in main text, Fig. LABEL:fig:13C_NV_T1b provide the relaxation times
given in Table 3.

In contrast to our other measurements, the magnetic field was set to the
lowest at which intensity in the EPR spectrum could be detected considering
the field region in main text, Fig. LABEL:fig:optic_polar, that is, using
$\nu=$9.59\text{\,}\mathrm{GHz}$$, $B=$2403\text{\,}\mathrm{G}$$. This allows
being selective to NV centers that are aligned to the magnetic field, and
avoid state-mixing induced by the misalignement. The ${}^{\mathrm{NV}}T_{1}$
measured with this method reflects exclusively the rate of the single quantum
transitions ($m_{S}=0$ to $\pm 1$), and thus can be compared to the typical
${}^{\mathrm{NV}}T_{1}$ measurements reported in the literature. However, our
determination is performed at higher magnetic field than in most cases of
literature reports (using, e.g., zero-field). The increased magnetic field can
lead to further decoupling of the NVs from certain noise sources, and thus
favor a longer $T_{1}$.

[HTML]e6e6e6 $100\text{\,}\mathrm{nm}$ sample | ${}^{\mathrm{NV}}T_{1}$
---|---
AO | $2.20(3)\text{\,}\mathrm{ms}$
AO + TAC | $4.23(7)\text{\,}\mathrm{ms}$
Table 3: ${}^{\mathrm{NV}}T_{1}$ corresponding to the fits shown in main text,
Fig. LABEL:fig:13C_NV_T1b.

###  S2.4 XPS on smaller nanoparticles ($\sim$25\text{\,}\mathrm{nm}$$)

X-ray photoemission spectrum on a sample of $25\text{\,}\mathrm{nm}$ median
size is shown in Fig. S6. The sample is from the same manufacturer as the
$100\text{\,}\mathrm{nm}$ $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ samples
(Pureon AG), however contamination induced by the fabrication process
(milling) is more likely to present due to the smaller size. The Fe2p peak of
iron could be detected, as shown in inset of Fig. S6. In the
$100\text{\,}\mathrm{nm}$ sample (AO treatment), no iron could be detected by
XPS, which we attribute to the lower concentration at the surface of bigger
particles.

Fig. S6: X-ray photoemission spectroscopy spectrum of the sample Pureon AG,
MSY0-0.5, type Ib, of $25\text{\,}\mathrm{nm}$ median size, showing the carbon
(C1s) and oxygen (O1s) peaks. Inset: acquisition in the [700,730] eV range,
showing the signal characteristic of iron (Fe2p).

##  S3 Sample illumination

###  S3.1 Illumination protocol

Fig. S7: a) 3D printed micro-photonic structure consisting of more than a
hundred tips. b) Schematic illustration of the 3D structure inside the EPR
tube, which provides optical access to the contained diamond powder.

In contrast to single diamond crystals, the propagation of light inside of
diamond powder of micro- and nanosize is hindered by strong light scattering,
due to the high refractive index of diamond ($n=2.4$ at the wavelength of
$532\text{\,}\mathrm{nm}$). In the context of hyperpolarization, there is a
requirement to perform homogeneous optical pumping of large amounts of sample.
A microstructured waveguide array was developed to ensure an efficient
illumination of such quantities of diamond powder and circumvent the issue
caused by light scattering.

The waveguide arrays make it possible to directly access a larger volume of
micro- or nanodiamonds with light, leading to a better polarization of NV
centers upon illumination.

The micro-pillar array was manufactured using a commercially available direct
laser writing system, employing the 2-photon-polymerization technique
(Photonic Professional GT+, Nanoscribe GmbH, Eggenstein-Leopoldshafen,
Germany). We used the resin IP-S (Nanoscribe GmbH, Eggenstein-Leopoldshafen,
Germany) due to its excellent optical properties in the VIS and 10x lens with
NA=0.3 for fast prototyping progression. The structures were printed on a
silicon wafer using the dip-in approach, where the resin acts as its own
immersion medium. The slicing and hatching distances were chosen as 5 µm and 1
µm, respectively, which, although showing visible steps between layers, still
produced good quality results. Uncured resin was dissolved and removed by
submerging the structure in PGMEA and cleaning was done by a subsequent bath
and rinsing with isopropyl alcohol during which UV curing was applied to
ensure the complete polymerization of all components and to reduce shrinkage.
The structure itself comprises a base with a diameter of 2.9 mm on which a
hexagonal pattern of approximately 200 pillars, each around 1 mm in height, is
placed. The pillars are truncated cones which are topped by a hemispherical
cap (see Figure S7). The distance between the revolution axes of neighboring
pillars is $200\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

In the sample preparation, the diamond powder was mixed with ethylcinnamat
(C11H12O2), the latter, with a refractive index
($\lambda=$532\text{\,}\mathrm{nm}$$) $n=1.57$, has the following use for
index matching 1) it reduces losses by internal reflection inside the
structured array and 2) it mitigates further the effects of light scattering
by the diamond particles. The mixture of diamond powder and ethylcinnamat
(with a mass equal or lower to that of the diamond powder) provides a viscous
fluid, which was placed in a sapphire tube of ($4\text{\,}\mathrm{mm}$ OD,
Situs Technical GmbH, Wuppertal, Germany) and the waveguide structure placed
on top. As shown on Fig. S7b, an optical fiber ($1.5\text{\,}\mathrm{mm}$
diameter) delivers the illumination. In separate trials, glycerol was used as
the index-matching solution. As an alternative, water could be used, if
working below the melting point. While water is in general a convenient
solvent for room temperature operation, it has the issue of its strong
microwave absorption above the melting point. A specific resonator design or
sample geometry would be required for working with water-containing samples at
room temperature.

###  S3.2 Phenomenological description of the NV spectrum in X-band EPR under
optical pumping

####  S3.2.1 Redistribution of state populations

The EPR spectrum measured under illumination (as in main text, Fig.
LABEL:fig:optic_polar) shows the contribution from NV centers in randomly
oriented crystallites. Therefore, all possible NV center orientations must be
considered in the spectrum description.

Let us consider the case of a NV center aligned to the magnetic field. For
this NV, the eigenstates, defined by the quantum number $m_{s}$, can be
written as $\ket{-1},\ket{0},\ket{1}$. Illumination depletes equally the
$\ket{\pm 1}$ states, by population transfer to $\ket{0}$, so that the
occupation probability of each state ($\rho_{m_{S}}$) differs from that in
thermal equilibrium. It obeys

$\rho_{m_{S}}=(1-p_{\mathrm{NV}})\rho_{m_{s},\mathrm{thermal}}+p_{\mathrm{NV}}\lvert\bra{m_{s}}\ket{0}\rvert^{2},$
(1)

where $p_{\mathrm{NV}}$ is the degree of initialization of the NV center. The
$\lvert\bra{m_{s}}\ket{0}\rvert^{2}$ term equals to 1 for $m_{s}=0$, and 0
otherwise. For $p_{\mathrm{NV}}=1$, the probability of being into the
$m_{s}=0$ state is maximum (100%).

When NV centers are misaligned with respect to the magnetic field, $m_{s}$ no
longer remains a good quantum number, due to state mixing. In this case, one
can define new eigenstates $\ket{s_{i,\theta}}$, where $i=1,2,3$ and $\theta$
is the angle between the NV axis and the orientation of the magnetic field.
The misalignment-induced state mixing leads to different occupation
probabilities compared to the aligned case. In our present phenomenological
description, it is taken into account by performing the following replacement
in Eq. 1:

$\lvert\bra{m_{s}}\ket{0}\rvert^{2}\rightarrow\lvert\bra{s_{i,\theta}}\ket{0}\rvert^{2},$
(2)

so that instead of Eq. 1, the occupation probability ($\rho_{i,\theta}$) of NV
center to $\ket{s_{i,\theta}}$ state can be expressed as

$\rho_{i,\theta}=(1-p_{\mathrm{NV}})\rho_{i,\theta,\mathrm{thermal}}+p_{\mathrm{NV}}\lvert\bra{s_{i,\theta}}\ket{0}\rvert^{2}.$
(3)

In EPR, the signal is summed over each transition and scales with the
population difference between the two states involved in the transition. It
follows from Eq. 3 that the signal ($S_{\mathrm{EPR}}$), e.g. in the field-
swept spectrum, can be predicted as

$S_{\mathrm{EPR}}=(1-p_{\mathrm{NV}})S_{\mathrm{EPR,thermal}}+p_{\mathrm{NV}}S_{\mathrm{EPR,pol.}},$
(4)

where $S_{\mathrm{EPR,thermal}}$ and $S_{\mathrm{EPR,pol.}}$ correspond
respectively to thermally polarized NVs, and NVs with occupation probabilities
given by $\lvert\bra{s_{i,\theta}}\ket{0}\rvert^{2}$ (Eq. 2).

Remarkably, the occupation probabilities predicted by this equation are
strongly orientation-dependent [2]. The energy level plots in Fig. S8
illustrate the probabilities in the 0° and 90° case, calculated for
$p_{\mathrm{NV}}=0.3$.

Fig. S8: Energy levels of an NV center with a) 0° and b) 90° orientation to
the magnetic field. The lines represented in a) result from direct
diagonalization of the Hamiltonian. The ones in b) are obtained after the
derivation in section S8.1 (Eq. 26). The arrows indicate the position of the
transitions in a field-swept experiment (the microwave frequency,
$\nu\approx$9.6\text{\,}\mathrm{GHz}$$). The diameter of the circles on the
energy levels depicts the corresponding occupation probabilities under optical
pumping, according to the description in section S3.2, with the parameter
describing the optical initialization of the NV center in Eq. 3,
$p_{\mathrm{NV}}=0.3$.

Although this description is of convenient use, we insist on its
phenomenological character. In reality, the mechanism responsible for the
change in occupation probabilities under illumination involves the different
rates in the intersystem crossing (ISC) from the 3E excited state [3]. A
complete description of the EPR spectrum under illumination should thus likely
consider the change in polarization dynamics induced by state mixing within
the 3E excited state [4]. Nevertheless, the good fit of our model to the
experimental data, shows that the phenomenological description already
captures the main features of the spectral changes induced by illumination.

####  S3.2.2 Spectral fits

Simulations of the NV spectra (main text, Fig. 1 and Fig. S9) were implemented
using the pepper in the EasySpin package in Matlab [5]. The spin system is
built considering the Hamiltonian in Eq. 18. To obtain
$S_{\mathrm{EPR,thermal}}$ in Eq. 4, one can set the experiment setting
Temperature to the wanted temperature in Kelvins (for ambient conditions,
Temperature = 298). To obtain $S_{\mathrm{EPR,pol}}$, the Temperature setting
is set to [1 0 0] (Matlab list). When the signal is significantly enhanced
above thermal ($S_{\mathrm{EPR}}/S_{\mathrm{EPR,thermal}}\gg 1$), as an
approximation, one can keep only the second term in Eq. 4, which then yields
$S_{\mathrm{EPR}}\approx p_{\mathrm{NV}}S_{\mathrm{EPR,pol.}}$.

Besides the orientation-induced broadening, secondary contributions to the
inhomogeneous broadening of the NV center are the magnetic interaction to the
surrounding spin bath (consisting of P1, other NVs, and additional
paramagnetic defects), strain, and electric fields induced by local charges
[6]. The interaction to the paramagnetic spin bath is considered by defining a
lorentzian broadening, using the linewidth parameter lw. To take into account
strain and electric fields induced by local charges -both generating center-
to-center variation of $e_{x},e_{y}$ terms in Eq. 18\- one can in fact define
$e_{y}=0$, and add a distribution of $e_{x}$, characterized by its full width
at half maximum FWHM($e_{x}$). Indeed, since powder averaging is performed in
the simulation, the $e_{y}$ term need not to be considered. We approximate the
distribution of the term $e_{x}$ as a gaussian, which is done through the
parameter DStrain (DStrain=[0,FWHM($e_{x}$)]). We find:

  * •

For the $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ particles (main text, Fig.
1): $D=$2869\text{\,}\mathrm{MHz}$$ , $lw=0.4$ mT, FWHM($e_{x}$)=
$21\text{\,}\mathrm{MHz}$,

  * •

For the $100\text{\,}\mathrm{nm}$ particles (Fig. S9):
$D=$2869\text{\,}\mathrm{MHz}$$, $lw=$0.4\text{\,}\mathrm{mT}$$,
FWHM($e_{x}$)=$27\text{\,}\mathrm{MHz}$,

Fig. S9: NV EPR spectrum obtained after illumination with
$700\text{\,}\mathrm{\SIUnitSymbolMicro s}$ laser pulses for
$100\text{\,}\mathrm{nm}$-sized nanodiamonds with a sample mass of 12 mg
(black solid). The mass was subsequently decreased for hyperpolarization
experiments, to ensure better illuminations conditions. Contributions from the
two separate transitions fitting the experimental spectrum are represented in
red and blue dashed curves, respectively.  Fig. S10: Enhancement of the NV EPR
signal under continuous illumination at different laser powers for the samples
used for hyperpolarization (a. $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$, b.
$100\text{\,}\mathrm{nm}$). Dashed black lines are a guide to the eye. The
enhancement corresponds to the intensity of the low field EPR peak
corresponding to NV centers with 90$\degree$ to the magnetic field, normalized
to its intensity measured in the dark. Corresponding spectra are shown in c.
and d. respectively.  Fig. S11: Dynamics of the EPR-detected NV signal induced
by a $700\text{\,}\mathrm{\SIUnitSymbolMicro s}$ laser pulse for the samples
used in the experiments described in main text. The magnetic field is
$B=$0.286\text{\,}\mathrm{T}$$, and the microwave frequency
$\nu=$9.58\text{\,}\mathrm{GHz}$$, corresponding to the low-field peak of the
Pake doublet. The experimental trace gives the enhancement of NV signal, at
different times. The dotted lines indicate the laser pulse duration used in
the hyperpolarization experiments reported in main text
($400\text{\,}\mathrm{\SIUnitSymbolMicro s}$), at which the enhancement for
the $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $100\text{\,}\mathrm{nm}$
particles is 46 and 61 respectively, and the full length
($700\text{\,}\mathrm{\SIUnitSymbolMicro s}$) used in this experiment.

###  S3.3 Heating induced by illumination

The local temperature increase of the diamond powder sample under illumination
can be calculated by considering the temperature-induced change of NV zero-
field splitting $D$ using the experimentally determined relation from Acosta
et al. [7]:

$\frac{dD}{dT}=-$0.074\text{\,}\mathrm{MHz}\text{\,}{\mathrm{K}}^{-1}$.$ (5)

Combining Eq. 5 and 28, following a change in temperature $\Delta T$, the low
field NV EPR peak is expected to shift by magnetic field $\Delta B$ as
follows:

$\frac{\Delta B}{\Delta
T}=\frac{dB_{12}}{dD}\frac{dD}{dT}=$0.0158\text{\,}\mathrm{G}\mathrm{/}\mathrm{K}$.$
(6)

In this work, for precise determination of temperature changes induced by
illumination, we perform the fit as detailed in section S3.2 of the EPR
spectrum around the low-field peak for different values of continuous laser
power (Figure S10), which provides the value of $D$ at each laser power.
Comparing to that obtained at a chosen laser power (31 mW and 100 mW for
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $100\text{\,}\mathrm{nm}$,
respectively), we can extract the corresponding temperature change (Fig. S12).

From the heating rates determined by the linear fits in Fig. S12, we can
determine the temperature elevation of the samples during the application of
the NV initialization and polarization transfer cycle described in main text,
Fig. LABEL:fig:hp_protocolc. As an average laser power of
$420\text{\,}\mathrm{mW}$ was used, the temperature changes are
+$51\text{\,}\mathrm{K}$ and +$24\text{\,}\mathrm{K}$ for the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $100\text{\,}\mathrm{nm}$
respectively.

Here, the comparatively high temperature increase for the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample can be attributed to the two-
times higher mass employed and the higher NV concentration, both increasing
the heat deposition by photon absorption.

Fig. S12: Heating of the samples used in the experiments described in main
text (a. $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and b.
$100\text{\,}\mathrm{nm}$) induced by continuous illumination at different
laser power. Dashed lines are a linear fits, allowing to extract as heating
rates $0.122\text{\,}\mathrm{K}\text{\,}{\mathrm{mW}}^{-1}$ and
$0.058\text{\,}\mathrm{K}\text{\,}{\mathrm{mW}}^{-1}$ respectively.

###  S3.4 Transient NV signal under pulsed illumination

To investigate the NV initialization dynamics, the evolution of the NV signal
during a $700\text{\,}\mathrm{\SIUnitSymbolMicro s}$ laser pulse was acquired
with EPR (Fig. S11). The magnetic field and frequency correspond to the
transition between the two lower energy spin states, $s_{1,90\degree}$ and
$s_{2,90\degree}$ (depicted in Fig. S8).

In the hyperpolarization experiments, shorter,
$400\text{\,}\mathrm{\SIUnitSymbolMicro s}$ laser pulses were used. From Fig.
S8, corresponding NV signal enhancements for the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $100\text{\,}\mathrm{nm}$
particles, are $\epsilon_{\mathrm{NV}}=46$ and 61 respectively.

We now turn to the interpretation of the value of $\epsilon_{\mathrm{NV}}$, in
terms of state occupation probabilities. The EPR signal scales with the
population difference between the states involved in the transitions, and the
number of centers contributing to the signal. We note here that by
photoluminescence measurement, no conversion towards the neutrally charged NV0
was observed in the samples in the used laser power range, therefore we
consider the number of NV centers to be constant. In these conditions,

$\epsilon_{\mathrm{NV}}=\frac{(\Delta\rho)_{\mathrm{NV}}}{(\Delta\rho)_{\mathrm{NV,dark}}},$
(7)

where $(\Delta\rho)_{\mathrm{NV}}=\rho_{s_{1}}-\rho_{s_{2}}$, and
$(\Delta\rho)_{\mathrm{NV,dark}}=\rho_{s_{1},\mathrm{dark}}-\rho_{s_{2},\mathrm{dark}}$,
here $\rho_{s_{1,2}}$ are the occupation probabilities of states
$s_{1,90\degree}$ and $s_{2,90\degree}$ (‘$90\degree$’ is omitted in the
subscript for simplicity). The additional subscript ‘dark’ indicates the case
without illumination.

In the dark, the occupation probabilities are given by the Boltzmann
statistics:

$\rho_{s_{i},\mathrm{dark}}=e^{\frac{-u_{i}}{kT}}/\left(\sum\limits_{j=1}^{3}e^{\frac{-u_{j}}{kT}}\right),$
(8)

where $u_{i}$ ($i=1,2,3$) are the state energies given in Eq. 26. Using
$D=$2869\text{\,}\mathrm{MHz}$$, $B=$2860\text{\,}\mathrm{G}$$, at
$T=$293\text{\,}\mathrm{K}$$ one obtains:

$(\Delta\rho)_{\mathrm{NV,dark}}=\rho_{s_{1},\mathrm{dark}}-\rho_{s_{2},\mathrm{dark}}=5.234\times
10^{-4}.$ (9)

From Eq. 9, it is straightforward to estimate the population difference as
$(\Delta\rho)_{\mathrm{NV}}=\epsilon_{\mathrm{NV}}(\Delta\rho)_{\mathrm{NV,dark}}$.
We obtain $(\Delta\rho)_{\mathrm{NV}}=0.024$ and 0.032 for the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $100\text{\,}\mathrm{nm}$
particles, respectively. We note that CW illumination brings similar
populations in the steady-state, as comparable values of the enhancement
$\epsilon_{\mathrm{NV}}$ are observed at high laser power in Fig. S10. The
determined $(\Delta\rho)_{\mathrm{NV}}$ thus can be compared to the value
$\approx 0.07$ reported in NV-ensemble studies of single crystal diamonds
under continuous illumination [2].

##  S4 Build-up and decay of 13C polarization

###  S4.1 $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$100\text{\,}\mathrm{nm}$ samples, with sample rotation (25°/s)

The build-up time, obtained from fitting the data in Figure S13 with a
stretched exponential function
$s(t)=A\left(1-\exp{-(\frac{t}{T_{\mathrm{pol}}})^{\beta}}\right)$, is
$T_{\mathrm{pol}}=34\pm 2$ s (with $\beta=0.85\pm 0.06$) for the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample and $T_{\mathrm{pol}}=32\pm
5$ s (with $\beta=0.88\pm 0.12$) for the $100\text{\,}\mathrm{nm}$ sample. The
13 spectra after buildup ( $t=$300\text{\,}\mathrm{s}$$) are shown in main
text, Fig. LABEL:fig:Hyper_signal. The decay of polarization at 1 T, estimated
from the fitting the data in Figure S13 with a single exponential decay, is
$T_{1}=142\pm 16$ s and $T_{1}=116\pm 16$ for
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $100\text{\,}\mathrm{nm}$
diamond powders, respectively.

Fig. S13:  Build-up (left) and decay (right) of 13C polarization for a)
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and b) $100\text{\,}\mathrm{nm}$
diamond powders, respectively. The dashed lines are the fits of the data with
stretched exponential function for the build-up time and single exponential
for the decay of 13C polarization.

###  S4.2 Buildup for the $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample,
without rotation

The build-up time, obtained from fitting the data in Figure S14 with a
stretched exponential function is $T_{\mathrm{pol}}=29\pm 4$ s (with
$\beta=0.70\pm 0.08$).

Fig. S14:  Build-up of 13C polarization for the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample in the absence of sample
rotation and fit with a stretched exponential function giving
$T_{\mathrm{pol}}=29\pm 4$ s and $\beta=0.70\pm 0.08$.

##  S5 Evaluation of the hyperpolarization enhancement factor and the absolute
polarization

The signal after hyperpolarization $S_{\mathrm{{}^{13}C}}$, the enhancement
factor $\epsilon_{\mathrm{{}^{13}C}}$ and the absolute polarization
$p_{\mathrm{{}^{13}C}}$ (defined in main text) obey

$\epsilon_{\mathrm{{}^{13}C}}=\frac{S_{\mathrm{{}^{13}C}}}{S_{\mathrm{{}^{13}C,therm}}}=\frac{p_{\mathrm{{}^{13}C}}}{p_{\mathrm{{}^{13}C,therm}}},$
(10)

where $S_{\mathrm{{}^{13}C,therm}}$, $p_{\mathrm{{}^{13}C,therm}}$ are the
values of the thermal signal and the corresponding absolute nuclear
polarization at thermal equilibrium.

The value of the absolute polarization at thermal equilibrium can be
calculated as
$p_{\mathrm{{}^{13}\mathrm{C},therm}}=\frac{\hbar\gamma_{{}^{13}\mathrm{C}}B}{2k_{\mathrm{B}}T}$
, where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ is the temperature,
$\gamma_{{}^{13}\mathrm{C}}$ is the nuclear gyromagnetic ratio (from the
Boltzmann statistics as $k_{\mathrm{B}}T\gg\hbar\gamma_{{}^{13}\mathrm{C}}B$).
From Eq. 10 it can be seen that, from the hyperpolarized signal, both
$\epsilon_{\mathrm{{}^{13}C}}$ and $p_{\mathrm{{}^{13}C}}$ can be determined
if the thermal signal $S_{\mathrm{{}^{13}C,therm}}$ can be measured. However,
the sensitivity of the setup for in situ NMR detection is low and, as a
result, insufficient to detect the thermal 13C signal from the diamond powder
samples. To address this limitation, we employed a reference sample of
$40\text{\,}\mathrm{\SIUnitSymbolMicro L}$ water, in which the high thermal
signal from the 1H nuclei can easily be measured.

For the $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ particles sample, the
hyperpolarized spectrum seen in main text, Fig. LABEL:fig:Hyper_signala,
corresponds to a signal (observed value of the free induction decay at $t\sim
0$) of $1.407\text{\,}\mathrm{\SIUnitSymbolMicro V}$ for a mass
$m=$12\text{\,}\mathrm{mg}$$ of sample. For this sample first derive
$p_{\mathrm{{}^{13}C}}$ from the reference ${}^{1}{\mathrm{H}}$ measurement in
the water sample. Then we deduce $\epsilon_{\mathrm{{}^{13}C}}$ at the
polarization field ($0.287\text{\,}\mathrm{T}$).

Similar calculations apply for the $100\text{\,}\mathrm{nm}$ sample for which
the spectrum in main text, Figure LABEL:fig:Hyper_signalb corresponds to a
signal of $0.441\text{\,}\mathrm{\SIUnitSymbolMicro V}$ for a sample mass
$m=$6\text{\,}\mathrm{mg}$$. For each sample, both
$\epsilon_{\mathrm{{}^{13}C}}$ and $p_{\mathrm{{}^{13}C}}$ are reported in
main text.

##### Determination of $p_{\mathrm{{}^{13}C}}$

By measuring the intensity of the 1H thermal signal in the water sample
($S_{{}^{1}\mathrm{H,therm}}$), one can estimate the intensity of the 13C
thermal signal as follows:

$S_{\mathrm{{}^{13}C,therm}}=S_{\mathrm{{}^{1}\mathrm{H},therm}}\Gamma_{\mathrm{ref}}.$
(11)

We now focus on the derivation of $\Gamma_{\mathrm{ref}}$. The intensity of
thermal signal $S_{\mathrm{nucl}}$ in a magnetic field $B$ at temperature $T$
can be estimated as [8]

$S_{\mathrm{nucl}}=C\times(\gamma
B)\times(\frac{\gamma^{2}\hbar^{2}BN_{\mathrm{nucl}}}{4k_{\mathrm{B}}T}).$
(12)

where $C$ is a constant, determined by parameters of the detection channel,
$\hbar$ is the Planck constant, $N_{\mathrm{nucl}}$ is the number of detected
nuclear spins, and $\gamma$ is the nuclear gyromagnetic ratio. This expression
is valid in the limit $k_{\mathrm{B}}T>>\hbar\gamma B$.

The 13C signal was detected at $1.004\text{\,}\mathrm{T}$. Unlike
conventionally done in NMR, the detection of signals for 1H nuclear spins was
implemented at the same frequency, adjusting the magnetic field for the 1H
nucleus ($B=$0.252\text{\,}\mathrm{T}$$). Therefore, it is convenient to
replace the magnetic field $B$ in equation 12, as

$B=\frac{2\pi f}{\gamma},$ (13)

where $f$ is the detection frequency. Furthermore, given that the measurement
were both performed at ambient temperature, the ratio of thermal signals for
1H in water and 13C nuclear spins in the respective samples can be expressed
as

$\Gamma_{\mathrm{ref}}=\frac{S_{\mathrm{{}^{13}C,therm}}}{S_{\mathrm{{}^{1}H,therm}}}=\frac{\gamma_{\mathrm{{}^{13}C}}N_{\mathrm{{}^{13}C}}}{\gamma_{\mathrm{\mathrm{{}^{1}H}}}N_{\mathrm{{}^{1}H}}}.$
(14)

The number of hydrogen and carbon nuclear spins can be estimated as follows:

$\displaystyle
N_{\mathrm{{}^{13}C}}=\frac{m_{\mathrm{diam.}}N_{\mathrm{A}}C_{\mathrm{{}^{13}C}}}{M_{\mathrm{diam.}}},N_{{}^{1}\mathrm{H}}=\frac{n_{\mathrm{{}^{1}H/mol}}m_{\mathrm{water}}N_{\mathrm{A}}C_{\mathrm{{}^{1}H}}}{M_{\mathrm{water}}},$
(15)

where $m_{\mathrm{diam.}}$ and $m_{\mathrm{water}}$ are the masses of diamond
and water samples, respectively, $N_{\mathrm{A}}$ is the Avogadro constant,
$C_{\mathrm{{}^{13}C}}$ and $C_{\mathrm{{}^{1}H}}$ are the isotopic abundances
of 13C and 1H, respectively, $n_{\mathrm{{}^{1}H/mol}}$ is the number of
hydrogen atoms per molecule ($n_{\mathrm{{}^{1}H/mol}}=2$),
$M_{\mathrm{diam.}}$ and $M_{\mathrm{water}}$ are the molar masses of water
molecules and carbon in diamond at natural abundance, respectively. From the
parameters listed in Table 4 we determine $\Gamma_{\mathrm{ref}}=1/1652$ for
the $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample described in main text.

sample | mass |  $C_{{}^{13}\mathrm{C}}$ or $C_{{}^{1}\mathrm{H}}$ | molar mass |  $\gamma_{{}^{13}\mathrm{C}}$ or $\gamma_{{}^{1}\mathrm{H}}$
---|---|---|---|---
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ | $12\text{\,}\mathrm{mg}$ | 1.07% | $12.011\text{\,}\mathrm{g}\text{\,}{\mathrm{mol}}^{-1}$ |  $(2\pi)$$10.705\text{\,}\mathrm{MHz}\text{\,}{\mathrm{T}}^{-1}$
water reference | $40\text{\,}\mathrm{mg}$ | 99.98% | $18.015\text{\,}\mathrm{g}\text{\,}{\mathrm{mol}}^{-1}$ |  $(2\pi)$$42.576\text{\,}\mathrm{MHz}\text{\,}{\mathrm{T}}^{-1}$
Table 4: Sample properties used for determining the parameter
$\Gamma_{\mathrm{ref}}$, in Eq. 14, used for determination of the signal
enhancement (illustrated here for the $2\text{\,}\mathrm{\SIUnitSymbolMicro
m}$ sample).

Eq. 14 and Eq. 10 yield:

$p_{{}^{13}\mathrm{C}}=\frac{p_{{}^{13}\mathrm{C,therm}}\times
S_{{}^{13}\mathrm{C}}}{S_{{}^{13}\mathrm{C,therm}}}=\frac{p_{{}^{13}\mathrm{C,therm}}\times
S_{{}^{13}\mathrm{C}}}{\Gamma_{\mathrm{ref}}\times
S_{{}^{1}\mathrm{H,therm}}}$ (16)

We calculate, for the absolute thermal 13C polarization at the detection field
of $1.004\text{\,}\mathrm{T}$ ($T=$293\text{\,}\mathrm{K}$$),
$p_{{}^{13}\mathrm{C,thermal}}=$8.80\text{\times}{10}^{-7}\text{\,}$$. For the
reference water sample, we measure
$S_{{}^{1}\mathrm{H,therm}}=$5.420\text{\,}\mathrm{\SIUnitSymbolMicro V}$$.
Thus, for the $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample:

$p_{{}^{13}\mathrm{C}}=\frac{$8.80\text{\times}{10}^{-7}\text{\,}$\times$1.407\text{\,}\mathrm{\SIUnitSymbolMicro
V}$}{\frac{1}{1652}\times$5.420\text{\,}\mathrm{\SIUnitSymbolMicro
V}$}=$3.77\text{\times}{10}^{-4}\text{\,}$$ (17)

Similarly, for the $100\text{\,}\mathrm{nm}$ sample,
$p_{{}^{13}\mathrm{C}}=$2.36\text{\times}{10}^{-4}\text{\,}$$ is obtained.

##### Determination of $\epsilon_{\mathrm{{}^{13}C}}$ at the polarization
field ($0.287\text{\,}\mathrm{T}$)

For the absolute thermal 13C polarization at the _polarization_ field of
$0.287\text{\,}\mathrm{T}$ ($T=$293\text{\,}\mathrm{K}$$), we calculate
$p_{{}^{13}\mathrm{C,thermal}}=$2.52\text{\times}{10}^{-7}\text{\,}$$. For the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample, it follows from Eq. 17 that:
$\epsilon_{{}^{13}\mathrm{C}}=\frac{p_{{}^{13}\mathrm{C}}}{p_{{}^{13}\mathrm{C,thermal}}}\approx
1500$. Similarly, for the $100\text{\,}\mathrm{nm}$ sample, we obtain
$\epsilon_{{}^{13}\mathrm{C}}\approx 940$.

##  S6 PulsePol parameter tuning

###  S6.1 Testing the robustness to detuning of the phase-offset PulsePol
sequence

The original PulsePol sequence (main text, Fig LABEL:fig:hp_protocola) was
proposed by Schwartz _et al._ [9]. In our experiments, we however used the
‘phase-offset’ sequence described in (main text, Fig. LABEL:fig:hp_protocolb).
Tratzmiller _et al._ introduced a general set of sequences based on PulsePol,
in which the phase definition contains a free parameter $\phi$ [10]. In this
definition, the original PulsePol and the sequence presently labelled ‘phase-
offset’ are obtained by choosing $\phi=\pi/2$ and $\phi=3\pi/4$, respectively.

In order to characterize the robustness to spectral detuning, we employ a
single-crystal diamond, represented in Fig. S15a, that was utilized in an
earlier work [11]. Addressing NV in single-crystal diamond instead of powder
drastically reduces their inhomogeneous broadening, so that the sequence
robustness can be tested by varying the pulse carrier frequency with respect
to a group of NV transitions occupying a narrow spectral range.

Fig. S15: a) Single crystal diamond used for the sequence robustness tests. b)
Spectrum of NV around the resonator frequency $\nu=$9.7\text{\,}\mathrm{GHz}$$
at $B=$2436\text{\,}\mathrm{G}$$ aligned to the NV axis. The several lines
correspond to the hyperfine structure for the $m_{S}=0$ to +1 transition of
NV. Coupling to 14N nuclear spin ($I=1$) explains the three-lines pattern, and
nearby 13C provide weak additional satellites. The spectrum is obtained by
performing Fourier Transform of the Hahn echo signal, acquired in EPR using
pulses generated by the $1\text{\,}\mathrm{kW}$ TWT of Rabi frequency
$\Omega_{1}=(2\pi)$25\text{\,}\mathrm{MHz}$$, under continous illumination at
$P=$100\text{\,}\mathrm{mW}$$ laser power.  Fig. S16: Testing the robustness
to detuning of the original PulsePol sequence and its variant, represented in
main text Fig. LABEL:fig:hp_protocola and LABEL:fig:hp_protocolb,
respectively. 13C hyperpolarized signal is acquired at different detunings
$\Delta$ of the pulse carrier frequency with respect to the central line
($m_{I}=0$ for 14N) of the NV signal. The detuning is represented in units of
the ordinary Rabi frequency $\nu_{1}=\Omega_{1}/2\pi$. The top curve shows the
FFT spectrum from Fig. S15b on the same scale. A high-power amplifier (TWT,
$1\text{\,}\mathrm{kW}$) was used in order to ensure a high
$\nu_{1}=$25\text{\,}\mathrm{MHz}$$, driving the different hyperfine (14N)
transitions similarly. Resonance conditions $n=3$ and $n=3.5$ were used for
PulsePol and phase-offset PulsePol respectively. The number of repetitions was
$M=10$. Lines in the main plot are a guide to the eye.

The 13C signal was acquired after hyperpolarization with both the standard and
phase-offset PulsePol, at different magnetic fields positively and negatively
detuned from the resonance of the central hyperfine line ($m_{I}=0$ for 14N)
visible in Fig. S15b. A change in magnetic field $\delta B$ shifts the NV
spectrum frequency for the central line ($\nu_{\mathrm{NV}}$) creating for the
pulse carrier frequency ($\nu_{\mathrm{MW}})$ a detuning
$\Delta=2\pi(\nu_{\mathrm{MW}}-\nu_{\mathrm{NV}})=-\gamma\delta B$, where
$\Delta$ is given in angular frequency unit,
$\gamma=(2\pi)$2.8032\text{\,}\mathrm{M}\mathrm{Hz}\mathrm{/}\mathrm{G}$$ is
the gyromagnetic ratio of the NV center. Figure S16 shows the hyperpolarized
13C signal acquired with both the standard and phase-offset PulsePol. It is
well visible that the latter demonstrates increased robustness to detuning,
i.e. ensures the polarization transfer to be efficient in a wider spectral
range.

###  S6.2 Comparison of resonance indices $n=3.5$ and $n=4.5$

In Fig. S17 is shown the simulated polarization loss of NV when applying the
phase-offset PulsePol sequence for different $\tau$ intervals, considering in
one case the interaction to a 13C nuclear spin (Fig. S17a) and in the other
case the interaction to the strongly coupled 14N nuclear spin of NV (Fig.
S17b).

Fig. S17: Simulation of polarization dynamics for NV oriented at 90°to the
magnetic field, in the conditions of our experiment
($B=$0.287\text{\,}\mathrm{T}$$). The horizontal axis represents the $\tau$
delay in the sequence and the values on the vertical axis correspond to the
remaining polarization of the NV, after $M=20$ repetitions of the phase-offset
PulsePol sequence. a) Case of a 13C spin coupled to NV with hyperfine tensor
elements $A_{ZX}=(2\pi)$80\text{\,}\mathrm{kHz}$,A_{ZY}=0$, without
interaction to 14N b) Case of NV-14N interaction, where the 14N Hamiltonian is
taken from Felton _et al._ [12], in the absence of interaction to 13C.

The interaction to the 13C nuclear spin leads to the expected half-integer
resonances. The $n=3.5$ and $n=4.5$ resonances are the most intense in
comparison to the other half-integer ones in the range of represented $\tau$
values. Thus the $n=3.5$ and $n=4.5$ resonance dips corresponds to a strong
probability of transferring polarization to 13C. For the chosen number of
repetition of the 6-pulse block ($M=20$ in the figure), the $n=4.5$ resonance
dip is slightly more pronounced than the $n=3.5$ one, however this advantage
should disappear if we consider the same interaction times (as a given number
of repetitions is achieved in a time that is shorter by $1-3.5/4.5\approx
22\%$ with the $n=3.5$ resonance). From these considerations, i.e., assuming
similar interaction times, the $n=3.5$ resonance is expected to be
advantageous.

However, the presence of the 14N nuclear spin must be included and leads to
more complex dynamics. This dynamics is specific to the case of NV misaligned
to the magnetic field, which leads to microwave transitions affecting
simultaneously the state of the 14N and NV spin. These transitions, which can
be predicted from the diagonalization of the Hamiltonian, are illustrated in
Fig. S18. The simulation of PulsePol in Fig. S17b, which considers a NV
exclusively coupled to its 14N spin, shows that applying the sequence leads to
a polarization loss on the NV center due to this interaction. The polarization
loss on NV shows complex dependence on the parameter $\tau$.

To prevent the interaction to 14N from disturbing the dynamics of polarization
transfer to 13C, we choose a $\tau$ value for which the simulation predicts
the polarization loss on NV induced by the NV-14N interaction to be
negligible. It is visible from Fig. S17b that in the conditions of our
experiment the $n=4.5$ resonance fullfills this criterion better than the
$n=3.5$ one, therefore the former is chosen.

To confirm the benefit of using the $n=4.5$ resonance, we performed
polarization transfer at both the $n=3.5$ and $n=4.5$ resonances, which leads
to the hyperpolarized signal seen in Figure S19
($2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ particles). Using the $n=4.5$
resonance yields a 70% higher absolute signal than using the $n=3.5$ one. This
confirms the advantage of choosing the $\tau$ value among available resonances
so as to minimize the advert effects of the NV-14N interaction.

As the evolution dynamics of the coupled NV and 14N is affected by the
magnetic field, our conclusions do not hold at any value of $B$. At other
magnetic field values, it might be favorable to use yet other phases in
PulsePol (by adjusting the parameter $\phi$ in the general sequence definition
from ref. [10]). Indeed, the resulting change in the resonance positions might
allow finding conditions where the 14N-related dynamics is refocused, in the
same way as at our magnetic field, using the $\tau$ value corresponding to
$n=4.5$. However, the simulation of the NV-14N interaction dynamics at
different magnetic fields is outside the scope of this work.

Fig. S18: Energy diagram and transitions induced by microwave on the NV center
spin coupled to the 14N nuclear spin for a magnetic field perpendicular to the
$C_{3v}$ symmetry axis of NV. The diagram is restricted to the low energy spin
states $\ket{s_{1}}$ and $\ket{s_{2}}$ of the NV center, which were the one
driven in our experiment. In this configuration the $\ket{m_{I}=0}$ state of
14N (choosing the axis of the applied magnetic field for quantization) is the
only one that is unaffected by the hyperfine coupling to the NV. The other 14N
spin states are $\ket{\alpha_{1}},\ket{\alpha_{2}}$ in the $\ket{s_{1}}$
manifold of NV and $\ket{\alpha^{\prime}_{1}},\ket{\alpha^{\prime}_{2}}$ in
the $\ket{s_{2}}$ manifold. The arrows represent the microwave transitions:
solid for the allowed transitions, and dashed for the ‘forbidden’ transitions,
respectively.  Fig. S19: Hyperpolarized signal obtained using the $n=3.5$ and
$n=4.5$ resonances, after $M=20$ repetitions of the phase-offset PulsePol
($2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample). The change in sign of
polarization is an expected characteristic of the sequence. The signal in the
case $n=4.5$ is estimated to be 70% higher in absolute intensity than in the
case $n=3.5$.

###  S6.3 The composite pulses

Composite pulses consist of (quasi-)contiguous pulses of different phase
and/or amplitude. This provides degrees of freedom that can be used for
improving the robustness properties of sequences [13]. The composite pulses
used in this work require applying microwave only on a ‘on’ or ‘off’ basis
with no amplitude variation during the ‘on’ time. This makes their application
of experimental convenience as one can use any microwave amplifier in its
saturation regime, without the need to correct for response non-linearities.

The composite pulses used in the present work for increasing the bandwidth of
PulsePol are a variant of the symmetric phase-alternating composite pulses
introduced in ref. [13]: instead of analytically choosing the pulse durations
in order to cancel low-order contributions from an expansion of the pulse-
induced rotation in the detuning, we choose a numerical optimization which
instead aims to reach a minimum fidelity threshold of the intended $\pi$-pulse
in a maximal bandwidth. This can be advantageous as the additional error
correction intrinsic to PulsePol loosens the constraints on the individual
pulses. In order to increase the available degrees of freedom for the pulse,
we allowed for non-zero waiting times between sub-pulses. A pulse description
is defined via a set of durations $t_{i}=a_{i}\ \pi/\Omega$ written as
$[a_{n},\ldots,a_{1},a_{0}]$ which corresponds to a composite pulse with
subsections of duration $(t_{n},\ldots,t_{1},2t_{0},t_{1},t_{2},\ldots,t_{n})$
where the central $t_{0}$ durations correspond to positive-phase pulses, the
odd durations ${t_{2k+1}}_{k\in\mathbb{N}}$ correspond to waiting times, and
the even durations ${t_{2k}}_{k\in\mathbb{N}}$ have alternating positive and
negative phase. In this notation, a non-composite rectangular $\pi$ pulse is
defined as $[0.5]$

In Table 5, the composite pulses are named by the number of pulses to each
side of the central pulse (’sidebands’), the overall duration and power of one
composite pulse is measured in $\pi$ pulse multiples. As there are waiting
times involved, these values do not coincide. As robustness, we give the range
of $\Delta/\Omega$ for which PulsePol gives strong polarisation transfer
(approximate values). A visual representation of the pulses, as well as the
theoretical inversion profile and detuning- and robustness properties of the
pulses as part of PulsePol can be found in Fig. S21.

pulse | duration | power | robustness | definition
---|---|---|---|---
normal | 1 | 1 | $<0.4$ | [0.5]
2 | 3.7 | 2.7 | 0.7 | [0.187, 0.251, 0.408, 0.253, 0.721]
3 | 5.1 | 3.6 | 0.75 | [0.13 0.231 0.386 0.261 0.527 0.236 0.771]
5 | 8.1 | 5.6 | 1 | [0.167, 0.32, 0.427, 0.348,
|  |  |  | 0.658, 0.263, 0.72, 0.292, 0.822]
Table 5: Overview over definitions and properties of the numerically optimized
phase-alternating composite pulses of which the 2-sideband pulse was used in
this work to increase the bandwidth of PulsePol. The values for robustness
here refer to the range of detunings $\Delta$ in which polarisation transfer
remains successful, such that value of $1$ corresponds to the bandwidth
$\Delta\in[-1\Omega,+1\Omega]$ (cf. Fig. S21). Fig. S20: a) Illustration of
the PulsePol sequence, where composite pulses with two sidebands replace
rectangular pulses. b) Central pulse, consisting of two types of $\pi/2$
pulses, of symmetric shape. Fig. S21: Composite pulses (a–d), theoretical
(solid line) and experimental (dotted line where available) inversion
properties (e–h), both assuming the experimental Rabi frequency with
$\pi/\Omega=48\,$ns; and pulse-error robustness properties as part of PulsePol
using an increased Rabi frequency of $\Omega=(2\pi)\ 30\,$MHz to allow for
also using the longer composite pulses (i–l). Thus, (a,e,i) correspond to the
standard rectangular pulse, (b,f,j) to the 2-sideband pulse, (c,g,k) to the
3-sideband pulse, and (d,h,l) to the 4-sideband composite pulse. The
experimental profiles appearing as dotted lines in (e,f) were taken on NVs
aligned to the magnetic field, in the single crystal diamond also used for
experiments described in section S6.1.

In our experiments, the two-sideband pulse was used, with Rabi
$\Omega_{1}=(2\pi)$10.5\text{\,}\mathrm{MHz}$$. It can be seen from Table 5
that the width of the range of detunings for which PulsePol gives strong
polarization transfer is
$\Delta_{\mathrm{pol}}=1.4\Omega_{1}=(2\pi)$14.7\text{\,}\mathrm{MHz}$$.

As derived in [13], the symmetry of the composite $\pi$-pulses implies that
they are constant-phase pulses over their respective bandwidth, as well as
that the left and right halves of the pulses are broadband $z$-to-$x$ and
$x$-to-$z$ point-to-point rotations for the left and right half of the pulse,
respectively. These properties make these composite pulses well-suited for the
application in PulsePol, as the $(\pi/2)_{X}-(\pi)_{Y}-(\pi/2)_{X}$
combination can play out its full error correction for small deviations of the
optimal rotation. Furthermore, the oscillatory behavior of the sidebands of
phase-alternating pulses makes it possible to fill a relatively large fraction
of the waiting times of PulsePol without significantly degrading transfer
properties.

###  S6.4 Other parameters adjustment

Fig. S22: Hyperpolarized 13C NMR signals for the
$2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ sample a) as a function of the
number of repetitions $M$ in the PulsePol sequence, using sequence repetition
time $\tau_{\mathrm{SRT}}=$1.5\text{\,}\mathrm{ms}$$ b) as a function of the
sequence repetition time $\tau_{\mathrm{SRT}}$, using $M=68$. To mitigate
sample heating in b), the length of the laser initialization pulse was varied
to ensure constant laser duty cycle, as
$t_{\mathrm{laser}}=\frac{$400\text{\,}\mathrm{\SIUnitSymbolMicro
s}$}{$1500\text{\,}\mathrm{\SIUnitSymbolMicro s}$}\times\tau_{\mathrm{SRT}}$.

##  S7 Potential adaptations for further increase of the bandwidth
$\Delta_{\mathrm{pol}}$

##### Using a higher power amplifier.

In our experiments, we used a $40\text{\,}\mathrm{W}$ microwave amplifier for
applying the polarization transfer sequences, and obtained
$\Delta_{\mathrm{pol}}=(2\pi)$15\text{\,}\mathrm{MHz}$$. Using a higher power
amplifier ($500\text{\,}\mathrm{W}$) would ensure
$\Delta_{\mathrm{pol}}\approx(2\pi)$50\text{\,}\mathrm{MHz}$$ without change
in the protocol.

##### Sequence multiplexing.

For further improvement, one could also use pulses fulfilling two criteria,
that is not only performing efficient driving in a defined spectral range, but
also leaving the spin polarization essentially unaffected outside this region.
Although our approach based on composite pulses enabled us to meet the first
criterion, it did not fulfill the second. However, one could possibly resort
to envelope-modulated pulses [14, 15]. Despite their application being
experimentally more involved (as they requires good control over the non-
linearities in the amplifier and the excitation channel overall) compared to
our composite pulses, they would enable applying the polarization sequence at
nearby frequencies (this is in fact comparable to the idea of applying
parallel sweeps, i.e. frequency combs in the ISE [16]). The frequency
selectivity of the pulse would mitigate the advert effects of multiplexing,
i.e. cross-talks.

We expect that this approach could extend further the proportion of addressed
NVs, leading to $\Delta_{\mathrm{pol}}=(2\pi)$100\text{\,}\mathrm{MHz}$$. The
simulation shows that this corresponds to a 16% fraction of addressed NVs. To
make implementation more straightforward and for ensuring faster polarization
transfer, the frequency multiplexing and the use of the high power amplifier
could be combined.

##  S8 NV spin Hamiltonian

In this section, we consider the NV spin Hamiltonian and derive properties
related to experiments done in diamond powder, where all NV orientations are
present, under a magnetic field.

First, we recall that the spin Hamiltonian for NV in an external magnetic
field $\mathbf{B}$ is given by (neglecting the 14N and 13C hyperfine terms)

$\hat{H}/h=D\hat{S}_{z}^{2}+\bar{\gamma}\mathbf{B}\cdot\hat{\mathbf{S}}+e_{x}(\hat{S}_{y}^{2}-\hat{S}_{x}^{2})+e_{y}(\hat{S}_{x}\hat{S}_{y}+\hat{S}_{y}\hat{S}_{x}),$
(18)

where $D$ is the zero-field splitting of the NV center, $z$ is the defect
axis, $\bar{\gamma}=\gamma/2\pi$ ($\gamma$ is the gyromagnetic ratio of NV),
$e_{x}$ and $e_{y}$ are terms related to strain and local charges [6].

###  S8.1 Spin Hamiltonian diagonalization for NV centers under transverse
$B$-field

Hyperpolarization of nuclear spins was implemented by using NV centers
oriented with 90° to the magnetic field. In this section, we will describe the
diagonalization of the spin Hamiltonian of NV centers with 90° orientation to
the magnetic field, which allows deriving analytically the position of the
resonance peaks in the field-swept spectrum. The obtained results are used in
section S3.3, to predict the shift in the magnetic field of the 90° resonance
(low field) following a change in $D$, induced by heating.

For diagonalization, we follow steps that have already been described in the
literature [17, 18]. For simplicity, we now assume $e_{x},e_{y}\approx 0$ in
Eq. 18, as the variation of the EPR line position with the temperature,
originates dominantly from the change in $D$. The Hamiltonian in Eq. 18
becomes

$\hat{H}/h\approx
D\hat{S}_{z}^{2}+\bar{\gamma}\mathbf{B}\cdot\hat{\mathbf{S}}.$ (19)

Let’s now assume a magnetic field $B$ is applied along $x$ perpendicular to
the axis of the NV. In this scenario, the Hamiltonian in Eq. 19 transforms
into

$\hat{H}/h=D\hat{S}_{z}^{2}+\bar{\gamma}B\hat{S}_{x},$ (20)

where the operators $\hat{S}_{z}$ and $\hat{S}_{x}$ in Eq. 20 are the spin
matrices:

$\hat{S}_{z}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix},\>\hat{S}_{x}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1&0\\\
1&0&1\\\ 0&1&0\end{pmatrix}.$ (21)

In the basis of the zero-field eigenstates $\ket{-1}$, $\ket{0}$, $\ket{1}$,
the 90° NV Hamiltonian can be expressed as

$\hat{H}/h=\begin{pmatrix}D&\frac{\bar{\gamma}B}{\sqrt{2}}&0\\\
\frac{\bar{\gamma}B}{\sqrt{2}}&0&\frac{\bar{\gamma}B}{\sqrt{2}}\\\
0&\frac{\bar{\gamma}B}{\sqrt{2}}&D\end{pmatrix}.$ (22)

In order to find energies, we must diagonalize NV Hamiltonian:

$\det(\hat{H}/h-u)=0,$ (23)
$\begin{vmatrix}\begin{pmatrix}D-u&\frac{\bar{\gamma}B}{\sqrt{2}}&0\\\
\frac{\bar{\gamma}B}{\sqrt{2}}&-u&\frac{\bar{\gamma}B}{\sqrt{2}}\\\
0&\frac{\bar{\gamma}B}{\sqrt{2}}&D-u\end{pmatrix}\end{vmatrix}=0,$ (24)
$-u(D-u)^{2}-2(D-u)(\frac{\bar{\gamma}B}{\sqrt{2}})^{2}=0.$ (25)

The solutions of equation 25 are

$u_{1}=\frac{D-\sqrt{D^{2}+(2\bar{\gamma}B)^{2}}}{2},\>u_{2}=D,\\\
\>u_{3}=\frac{D+\sqrt{D^{2}+(2\bar{\gamma}B)^{2}}}{2}.$ (26)

The corresponding dependency of energies versus magnetic field is shown in
Fig. S8b. EPR spectra are typically obtained by sweeping a magnetic field at a
constant microwave frequency $\nu$. The resonance magnetic fields $B_{12}$ and
$B_{23}$, corresponding to the transitions $u_{2}\rightarrow u_{1}$ and
$u_{3}\rightarrow u_{2}$, are calculated as

$B_{12}=\frac{1}{\bar{\gamma}}\sqrt{\nu(\nu-D)},\>B_{23}=\frac{1}{\bar{\gamma}}\sqrt{\nu(\nu+D)}.$
(27)

The change of the low field $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ NV
resonance line position (for use in section S3.3) can be written as a function
of $D$ as

$\frac{dB_{12}}{dD}=\frac{-\nu}{2\bar{\gamma}\sqrt{\nu^{2}-D\nu}}=$-0.213\text{\,}\mathrm{G}\mathrm{/}\mathrm{M}\mathrm{Hz}$,$
(28)

taking $\nu=$9.6\text{\,}\mathrm{GHz}$$, $D=$2869\text{\,}\mathrm{MHz}$$,
$\bar{\gamma}=$2.8032\text{\,}\mathrm{M}\mathrm{Hz}\mathrm{/}\mathrm{G}$$.

The eigenstates corresponding to the energies in Eq. 26 are:

$\displaystyle\ket{s_{1,90\degree}}$
$\displaystyle=\frac{1}{\sqrt{1+c^{2}}}\left(\ket{0}+c\frac{\ket{-1}+\ket{+1}}{\sqrt{2}}\right)$
$\displaystyle\ket{s_{2,90\degree}}$
$\displaystyle=\frac{\ket{+1}-\ket{-1}}{\sqrt{2}}$ (29)
$\displaystyle\ket{s_{3,90\degree}}$ $\displaystyle=\frac{1}{\sqrt{1+c^{\prime
2}}}\left(\ket{0}+c^{\prime}\frac{\ket{-1}+\ket{+1}}{\sqrt{2}}\right),$

with
$c=\frac{\frac{D}{2}-\sqrt{(\frac{D}{2})^{2}+(\bar{\gamma}B)^{2}}}{\bar{\gamma}B}$
and
$c^{\prime}=\frac{\frac{D}{2}+\sqrt{(\frac{D}{2})^{2}+(\bar{\gamma}B)^{2}}}{\bar{\gamma}B}$.
At $B=$0.287\text{\,}\mathrm{T}$$, one has $c=-0.8375$ and
$c^{\prime}=1.1941$.

###  S8.2 Fraction of NV centers in the bandwidth $\Delta_{\mathrm{pol}}$
($\approx$90\text{\,}\mathrm{\SIUnitSymbolDegree}$$ to the magnetic field)

We derive the condition for NV with an axis that is slightly misaligned from
the perpendicular condition, i.e. at an angle $90\degree+\delta$ from the
magnetic field, to be present in the active bandwidth of the polarization
sequence, $\Delta_{\mathrm{pol}}$.

In the frame of the NV center $x,y,z$, the magnetic field can be written as
$\mathbf{B}=(B\cos(\delta),0,B\sin(\delta))$. In fact, the effect of the tilt
in the magnetic field can be treated considering the change in the magnetic
field vector $\mathbf{b}=(b_{\perp},0,b_{\parallel})\approx
B(-\delta^{2}/2,0,\delta)$. We solve the Hamiltonian in the basis of the
eigenstates for the exact, $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ tilt of
the magnetic field to the NV axis, which were determined previously in section
S8.1.

$\begin{vmatrix}\begin{pmatrix}u_{1}+\frac{2c}{1+c^{2}}\bar{\gamma}b_{\perp}-u^{\prime}&\frac{c}{\sqrt{1+c^{2}}}\bar{\gamma}b_{\parallel}&\frac{c+c^{\prime}}{\sqrt{(1+c^{2})(1+c^{\prime
2})}}\bar{\gamma}b_{\perp}\\\
\frac{c}{1+c^{2}}\bar{\gamma}b_{\parallel}&u_{2}-u^{\prime}&\frac{c^{\prime}}{1+c^{\prime
2}}\bar{\gamma}b_{\parallel}\\\
\frac{c+c^{\prime}}{\sqrt{(1+c^{2})(1+c^{\prime
2})}}\bar{\gamma}b_{\perp}&\frac{c^{\prime}}{1+c^{\prime
2}}\bar{\gamma}b_{\parallel}&u_{3}+\frac{2c^{\prime}}{1+c^{\prime
2}}\bar{\gamma}b_{\perp}-u^{\prime}\end{pmatrix}\end{vmatrix}=0,$

where $u_{i},i=1,2,3$ are the energies in Eq. 26 and $c,c^{\prime}$ are the
magnetic field-dependent terms defined in section S8.1, $u^{\prime}$ is the
variable to be determined.

Restricting to the terms up to the second order in $\delta$, the new energies
read

$\displaystyle u^{\prime}_{1}$
$\displaystyle=u_{1}-\bar{\gamma}B\left[\frac{c}{1+c^{2}}+\frac{c^{2}}{1+c^{2}}\frac{\bar{\gamma}B}{u_{2}-u_{1}}\right]\delta^{2}$
$\displaystyle=u_{1}+\frac{D}{2}\left(1-\frac{D}{\sqrt{D^{2}+(2\bar{\gamma}B)^{2}}}\right)\delta^{2},$
(30) $\displaystyle u^{\prime}_{2}$
$\displaystyle=u_{2}+\bar{\gamma}B\left[\frac{c^{2}}{1+c^{2}}\frac{\bar{\gamma}B}{u_{2}-u_{1}}-\frac{c^{\prime
2}}{1+c^{\prime 2}}\frac{\bar{\gamma}B}{u_{3}-u_{2}}\right]\delta^{2}$
$\displaystyle=u_{2}-D\delta^{2},$ (31) $\displaystyle u^{\prime}_{3}$
$\displaystyle\approx u_{3}+\bar{\gamma}B\left[+\frac{c^{\prime
2}}{1+c^{\prime
2}}\frac{\bar{\gamma}B}{u_{3}-u_{2}}-\frac{c^{\prime}}{1+c^{\prime
2}}\right]\delta^{2}$
$\displaystyle=u_{3}+\frac{D}{2}\left(1+\frac{D}{\sqrt{D^{2}+(2\bar{\gamma}B)^{2}}}\right)\delta^{2}.$
(32)

We now consider the transition $\ket{s_{1}}\rightarrow\ket{s_{2}}$ between the
states of energy $u^{\prime}_{1}$ and $u^{\prime}_{2}$ at
$B=$0.287\text{\,}\mathrm{T}$$. The transition frequency
$\nu_{12}^{\prime}=u^{\prime}_{2}-u^{\prime}_{1}$, varies with the tilt
$\delta$ so that

$\nu_{12}^{\prime}\approx\nu_{12}-($4052\text{\,}\mathrm{MHz}$)\times\delta^{2},$
(33)

where $\nu_{12}=u_{2}-u_{1}$. Now taking into account the width (here in
ordinary frequency) $\Delta_{\mathrm{pol}}/2\pi=$15\text{\,}\mathrm{MHz}$$, it
is possible to cover all the NV centers for which

$|\delta|<\delta_{\mathrm{M}}=\sqrt{\frac{15}{4052}}=$0.0608\text{\,}\mathrm{rad}$\hskip
20.0pt[3.48\degree].$ (34)

The fraction of NV centers that is addressed is the covered portion of the
solid angle, that is exactly $\delta_{\mathrm{M}}$. This suggests that 6.08 %
of NV centers are within the bandwidth $\Delta_{\mathrm{pol}}$.

However, this calculation neglects the inhomogeneous broadening from sources
other than the orientation-dependence of the transition frequency. Including
the additional sources of inhomogeneous broadening, determined from the
spectral fits (section S3.2), we find from a numerical simulation that the
fraction of NV centers is 4.15 % for the $\ket{s_{1}}\rightarrow\ket{s_{2}}$
transition. The other transition $\ket{s_{2}}\rightarrow\ket{s_{3}}$
contributes by a small amount, 0.21 %. Thus we keep $\mathbf{\approx 4\%}$ as
the number of NV addressed by the polarization sequence, in the
$15\text{\,}\mathrm{MHz}$ bandwidth, and correspondingly
$\delta_{\mathrm{M}}=$0.04\text{\,}\mathrm{rad}$\sim\mathbf{2.5\degree}$.

###  S8.3 NV transition dipole and Rabi inhomogeneity

The different orientations of the nanocrystals leads to a variation of the
magnetic transition dipole, even if the orientation of the NV center to the
external magnetic field is fixed ($\approx 90\degree$ in our experiments). It
is possible to derive the corresponding inhomogeneity of the Rabi frequency
for NV centers. For a linearly polarized microwave field and a transition
between two states $\ket{s_{i}}$, and $\ket{s_{j}}$, the transition dipole can
be expressed as

$d_{ij}=\bra{s_{i}}\hat{S}_{1}\ket{s_{j}},$ (35)

where $\hat{S}_{1}$ is the projection on the spin-operator in the direction of
the linearly polarized $\mathbf{B_{1}}$-field. In our configuration (EPR
resonator), the $\mathbf{B_{1}}$-field is perpendicular to the static magnetic
field $\mathbf{B}$. Using the coordinate system in Fig. S23a, it follows that
$\hat{S}_{1}=\cos\phi\hat{S}_{z}+\sin\phi\hat{S}_{y}$.

Fig. S23: a) Coordinate system of the NV frame, in the case of a perpendicular
magnetic field ($\parallel x$), and definition of the orientation of the
linearly polarized microwave field (in $yz$ plane). b) Rabi vs $\phi$ angle
for the employed $\ket{s_{1}}\rightarrow\ket{s_{2}}$ transition of NV at
$B=$0.287\text{\,}\mathrm{T}$$ after normalization to the Rabi frequency
$\Omega_{1,\infty}$ obtained in the high-field regime ($\bar{\gamma}B\gg D$).

Considering the 90°-oriented NV centers and the transitions between the states
$\ket{s_{1}}$ and $\ket{s_{2}}$, using Eq. 29 we obtain:
$d_{12}=(c\cos(\phi)+i\sin(\phi))/{\sqrt{1+c^{2}}}$, where $c$ is the
magnetic-field dependent term defined in section S8.1. As $c<1$, The magnitude
of $d_{12}$ is the weakest when the microwave field points in the direction of
the axis ($z$) of the NV center, and the strongest when it is directed along
$y$ (as $c^{\prime}>1$, the situation is however reverted for the other
transition). It is straightforward to derive the Rabi frequency, as it obeys
$\Omega_{1}\propto|d_{12}|$. One can write

$\frac{\Omega_{1}}{\Omega_{1,\infty}}=\frac{\sqrt{c^{2}\cos^{2}(\phi)+\sin^{2}(\phi)}}{\sqrt{1+c^{2}}}$
(36)

where $\Omega_{1,\infty}$ is the Rabi frequency assuming the Zeeman
eigenstates (obtained for $B\gg D$: $c\rightarrow 1$ and therefore
$d_{12}=1/\sqrt{2}$).

At $B=$0.287\text{\,}\mathrm{T}$$, the values of
$\Omega_{1}/\Omega_{1,\infty}$ vary in the range [0.91,1.08] (Fig. S23b).
Roughly speaking, the Rabi frequency has variations included within a $\pm
10\%$ range at this magnetic field.

###  S8.4 NV transition frequencies in diamond powder at different magnetic
fields

The distribution of transition frequencies in randomly oriented NV ensembles
is shown in Figure S24a for the main NV transitions. The simulation was
performed with Easyspin. In contrast with the simulation described in section
S3.2, the setting Temperature was not defined (to give equal positive weight
to both transitions) and the optional setting Intensity was set to Off (to
remove the effect of varying transition dipoles as a function of the
orientation). .

Fig. S24: a) Simulation of the distribution of the transition frequencies of
NV in diamond powder at magnetic fields
$B=$150\text{\,}\mathrm{G}$,$300\text{\,}\mathrm{G}$$ and
$2900\text{\,}\mathrm{G}$. The frequencies $\nu,\nu^{\prime}$ correspond to
the transitions as defined in b), where $B_{\mathrm{GSLAC}}$ is the magnetic
field of the ground state level anticrossing
($B_{\mathrm{GSLAC}}\sim$1024\text{\,}\mathrm{G}$$). The spectrum for the
transition at $\nu^{\prime\prime}$ (dotted arrows) is not represented in a)
(see text).

Besides the transitions represented in Figure S24a, we indicate a third
transition in Figure S24b (at frequency $\nu^{\prime\prime}$). In the case of
a magnetic field aligned to the NV axis, the transition is known as double
quantum ($-1\leftrightarrow+1$), and has vanishing dipole. In the case of a
magnetic field tilted from the NV axis, this third transition can be excited
as a single quantum transition, owing to spin state mixing. At the considered
magnetic fields, however, the corresponding dipole is typically lower than
that of the other transitions, thus it was not represented in the spectrum in
Figure S24a

## References

  * [1] B. V. Yavkin, G. V. Mamin, M. R. Gafurov, and S. B. Orlinskii.  Size-dependent concentration of N0 paramagnetic centres in HPHT nanodiamonds.  Magn. Reson. Solids, 17:15101, 2015.
  * [2] M. Drake, E. Scott, and J. A. Reimer.  Influence of magnetic field alignment and defect concentration on nitrogen-vacancy polarization in diamond.  New Journal of Physics, 18(1):013011, dec 2015.
  * [3] M. L. Goldman, M. W. Doherty, A. Sipahigil, N. Y. Yao, S. D. Bennett, N. B. Manson, A. Kubanek, and M. D. Lukin.  State-selective intersystem crossing in nitrogen-vacancy centers.  Physical Review B, 91(16), April 2015.
  * [4] J.-P. Tetienne, L. Rondin, P. Spinicelli, M. Chipaux, T. Debuisschert, J.-F. Roch, and V. Jacques.  Magnetic-field-dependent photodynamics of single NV defects in diamond: an application to qualitative all-optical magnetic imaging.  New Journal of Physics, 14(10):103033, October 2012.
  * [5] S. Stoll and A. Schweiger.  EasySpin, a Comprehensive Software Package for Spectral Simulation and Analysis in EPR.  J. Magn. Reson., 178(1):42–55, 2006.
  * [6] T. Mittiga, S. Hsieh, C. Zu, B. Kobrin, F. Machado, P. Bhattacharyya, N. Z. Rui, A. Jarmola, S. Choi, D. Budker, and N. Y. Yao.  Imaging the local charge environment of nitrogen-vacancy centers in diamond.  Physical Review Letters, 121(24), December 2018.
  * [7] V. M. Acosta, E. Bauch, M. P. Ledbetter, A. Waxman, L.-S. Bouchard, and D. Budker.  Temperature Dependence of the Nitrogen-Vacancy Magnetic Resonance in Diamond.  Phys. Rev. Lett., 104:070801, Feb 2010.
  * [8] Andrew Webb.  Increasing the Sensitivity of Magnetic Resonance Spectroscopy and Imaging.  Analytical Chemistry, 84(1):9–16, 2012.  PMID: 21978101.
  * [9] Ilai Schwartz, Jochen Scheuer, Benedikt Tratzmiller, Samuel Müller, Qiong Chen, Ish Dhand, Zhen-Yu Wang, Christoph Müller, Boris Naydenov, Fedor Jelezko, and Martin B. Plenio.  Robust optical polarization of nuclear spin baths using Hamiltonian engineering of nitrogen-vacancy center quantum dynamics.  Science Advances, 4(8), 2018.
  * [10] Benedikt Tratzmiller, Jan F. Haase, Zhenyu Wang, and Martin B. Plenio.  Parallel selective nuclear-spin addressing for fast high-fidelity quantum gates.  Physical Review A, 103(1), jan 2021.
  * [11] J. Scheuer, I. Schwartz, Q. Chen, D. Schulze-Sünninghausen, P. Carl, P. Höfer, A. Retzker, H. Sumiya, J. Isoya, B. Luy, M. B. Plenio, B. Naydenov, and F. Jelezko.  Optically induced dynamic nuclear spin polarisation in diamond.  New Journal of Physics, 18(1):013040, January 2016.
  * [12] S. Felton, A. M. Edmonds, M. E. Newton, P. M. Martineau, D. Fisher, D. J. Twitchen, and J. M. Baker.  Hyperfine interaction in the ground state of the negatively charged nitrogen vacancy center in diamond.  Phys. Rev. B, 79:075203, Feb 2009.
  * [13] A.J Shaka and A Pines.  Symmetric phase-alternating composite pulses.  Journal of Magnetic Resonance (1969), 71(3):495–503, 1987.
  * [14] D. J. Lurie.  A systematic design procedure for selective pulses in NMR imaging.  Magnetic Resonance Imaging, 3(3):235–243, 1985.
  * [15] Kyryl Kobzar, Sebastian Ehni, Thomas E. Skinner, Steffen J. Glaser, and Burkhard Luy.  Exploring the limits of broadband 90° and 180° universal rotation pulses.  Journal of Magnetic Resonance, 225:142–160, 2012.
  * [16] A. Ajoy, R. Nazaryan, K. Liu, X. Lv, B. Safvati, G. Wang, E. Druga, J. A. Reimer, D. Suter, C. Ramanathan, C. A. Meriles, and A. Pines.  Enhanced dynamic nuclear polarization via swept microwave frequency combs.  Proceedings of the National Academy of Sciences, 115(42):10576–10581, oct 2018.
  * [17] F. E. Mabbs and D. Collison.  Electron paramagnetic resonance of d transition metal compounds.  Elsevier, 2013.
  * [18] Boris Yavkin, Marat Gafurov, Mikhail Volodin, Georgy Mamin, and Sergei B. Orlinskii.  EPR and double resonances in study of diamonds and nanodiamonds.  In Experimental Methods in the Physical Sciences, pages 83–113. Elsevier, 2019.