††thanks: These authors contributed equally to this work††thanks: These
authors contributed equally to this work††thanks: These authors contributed
equally to this work††thanks: These authors contributed equally to this
work††thanks: These authors contributed equally to this work

# Nanoscale engineering and dynamical stabilization of mesoscopic spin
textures

Kieren Harkins  Department of Chemistry, University of California, Berkeley,
Berkeley, CA 94720, USA.  Christoph Fleckenstein  Department of Physics, KTH
Royal Institute of Technology, SE-106 91 Stockholm, Sweden.  Noella D’Souza
Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.  Paul M. Schindler  Max Planck Institute for the Physics of
Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany.  David Marchiori
Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.  Claudia Artiaco  Department of Physics, KTH Royal Institute of
Technology, SE-106 91 Stockholm, Sweden.  Quentin Reynard-Feytis  Department
of Chemistry, University of California, Berkeley, Berkeley, CA 94720, USA.
Ushoshi Basumallick  Department of Chemistry, University of California,
Berkeley, Berkeley, CA 94720, USA.  William Beatrez  Department of Chemistry,
University of California, Berkeley, Berkeley, CA 94720, USA.  Arjun Pillai
Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.  Matthias Hagn  Department of Chemistry, University of California,
Berkeley, Berkeley, CA 94720, USA.  Aniruddha Nayak  Department of Chemistry,
University of California, Berkeley, Berkeley, CA 94720, USA.  Samantha Breuer
Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.  Xudong Lv  Department of Chemistry, University of California,
Berkeley, Berkeley, CA 94720, USA.  Maxwell McAllister  Department of
Chemistry, University of California, Berkeley, Berkeley, CA 94720, USA.  Paul
Reshetikhin  Department of Chemistry, University of California, Berkeley,
Berkeley, CA 94720, USA.  Emanuel Druga  Department of Chemistry, University
of California, Berkeley, Berkeley, CA 94720, USA.  Marin Bukov  Max Planck
Institute for the Physics of Complex Systems, Nöthnitzer Str. 38, 01187
Dresden, Germany.  Ashok Ajoy  Department of Chemistry, University of
California, Berkeley, Berkeley, CA 94720, USA.  Chemical Sciences Division,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA.  CIFAR Azrieli
Global Scholars Program, 661 University Ave, Toronto, ON M5G 1M1, Canada.

###### Abstract

Thermalization phenomena, while ubiquitous in quantum systems, have
traditionally been viewed as obstacles to be mitigated. In this study, we
demonstrate the ability, instead, to harness thermalization to dynamically
engineer and stabilize structured quantum states in a mesoscopically large
ensemble of spins. Specifically, we showcase the capacity to generate,
control, stabilize, and read out “shell-like” spin texture with interacting
${}^{13}\textrm{C}$ nuclear spins in diamond, wherein spins are polarized
oppositely on either side of a critical radius. The texture spans several
nanometers and encompasses many hundred spins. We capitalize on the
thermalization process to impose a quasi-equilibrium upon the generated
texture; as a result, it is highly stable, immune to spin diffusion, and
endures over multiple-minute long periods — over a million times longer than
the intrinsic interaction scale of the spins. Additionally, the texture is
created and interrogated without locally controlling or probing the nuclear
spins. These features are accomplished using an electron spin as a nanoscale
injector of spin polarization, and employing it as a source of spatially
varying dissipation, allowing for serial readout of the emergent spin texture.
Long-time stabilization is achieved via prethermalization to a Floquet-induced
Hamiltonian under the electronic gradient field. Our work presents a new
approach to robust nanoscale spin state engineering and paves the way for new
applications in quantum simulation, quantum information science, and nanoscale
imaging.

##  I Introduction

Thermalization is a pervasive phenomenon in all of physics. The quest to
elucidate how isolated microscopic quantum systems approach equilibrium, has
spurred a thriving field at the intersection of contemporary theoretical and
experimental research. This has been marked by ongoing developments in
nonequilibrium dynamics, such as the Eigenstate Thermalization Hypothesis for
closed quantum systems [1, 2, 3], anomalous transport and emergent
hydrodynamics [4, 5, 6], or the creation of prethermal ordered states of
matter [7, 8].

Thermalization in quantum systems proceeds analogously to its classical
counterpart: entropy grows with time, erasing traces of the system’s prior
history. Physically, this translates into a gradual reduction of accessible
information that can be obtained from local measurements until, in thermal
equilibrium, a maximum entropy state is reached. In most cases, such loss of
information is irreversible, impeding the capacity to harness quantum systems.
Consequently, a broad endeavor has been underway to retard, even preclude,
these thermalization processes. Techniques range from physically isolating
quantum systems (e.g. in vacuum) or through quantum control [9], cooling them
to near absolute zero temperatures [10], or, in a complementary manner,
exploiting theoretical paradigms of many-body localization to inhibit the
onset of thermalization [11].

Though the thermalization process is often depicted as leading to mundane,
featureless states, in this work, we demonstrate its utility in preparing
structured mesoscopic quantum states. Specifically, in a system of interacting
nuclear spins at high temperature (${>}$100K), we exploit out-of-equilibrium
thermalizing dynamics to controllably engineer and stabilize mesoscopic shell-
like nuclear spin polarization textures that span several nanometers and
envelope hundreds of nuclear spins (Fig. 1A). Simultaneously, we continuously
observe the formation and stabilization of these textures with high temporal
resolution over prolonged, multiple-minute-long periods. This spatiotemporal
control, facilitated by thermalization, obviates the need for local spin
manipulation or readout, intrinsically protects against control errors, and
bypasses technical challenges of differentiating spins with near-identical
resonance frequencies. Given these methodological advantages, our approach has
direct implications for quantum memories [12, 13], spintronics [14, 15], and
nanoscale magnetic resonance imaging [16, 17].

Our experiments are in a model system of Nitrogen Vacancy (NV) center
electrons in diamond surrounded by ${}^{13}\textrm{C}$ nuclear spins (Fig. 1A)
[18]. The sparsely distributed NVs encompass ${\sim}10^{4}$
${}^{13}\textrm{C}$ nuclei, spanning a radius of ${\approx}12\,$nm [19]. Our
strategy is built upon three elements. Firstly, the optically polarizable NV
electron serves as a “polarization injector” and a nanoscale “antenna”. In its
role as an antenna, it generates a nanoscale gradient through the hyperfine
interaction, resulting in localized displacement of the nuclear spin resonance
frequencies [20]. Under this gradient field, a time-periodic (Floquet) drive
stabilizes the nuclei into a metastable state [21], characterized by a shell-
like polarization texture (Fig. 1B). Finally, spatially dependent dissipation
from the NV is exploited to observe these textures: ${}^{13}\textrm{C}$ nuclei
in the electronic proximity relax faster, providing a means to serially probe
spins farther away from the electron, without locally measuring them (Fig.
1B).

To demonstrate this new toolbox, we first demonstrate a simple means to
produce spin textures. Here the NV is exploited to successively inject
“hyperpolarization” – polarization that is orders of magnitude greater than
Boltzmann levels – of alternate sign into the ${}^{13}\textrm{C}$ nuclei.
However, the resulting State-Engineered spin texture lacks stability: its
domain boundary melts due to nuclear spin diffusion [22].

Exploiting thermalization, however, offers a solution — creating polarization
textures with domain boundaries protected against spin diffusion. We create an
effective inhomogeneous parent Hamiltonian with spatial characteristics
derived from the NV antenna field (Fig. 1C). Irrespective of the initial state
111Note that the initial state cannot be chosen fully arbitrarily but must
have a finite energy density with respect to $\mathcal{H}^{(0)}_{F}$ (3); for
more details see methods., the spin polarization naturally prethermalizes [24,
25, 26] to its drive-induced quasi-equilibrium state, forming spin shells with
constant-in-time domain boundaries. We observe the formation of this generated
Hamiltonian-Engineered [27, 28, 29] texture continuously with unprecedented
resolution over a lifetime spanning several minutes. Numerical simulations
corroborate these observations and illustrate a stable spin texture spanning
several nanometers.

Fig. 1: System and readout. (A) Spin-textures. Controllable “shell-like” spin
texture of positive and negative polarization (shaded red and blue,
respectively) is generated in a nanoscale ensemble of ${}^{13}\textrm{C}$
nuclear spins surrounding a central NV electron (manipulated by an external
laser, shown in green). The domains are separated by the critical radius
$r_{c}$ at the domain boundary (dashed white line). The texture encompasses
${\sim}\mathcal{O}(100)$ spins within the critical radius and remains stable
for minutes. The optically-pumped electron serves as a spin injector, and
produces a nanoscale magnetic field gradient $B$ (shaded light pink) that
stabilizes the spin texture in a prethermal state. (B) Experimental schematic
showing prethermalization caused by Floquet driving at high-magnetic field
(HF) $B_{0}{\geq}7$ T, with a train of spin-locking $\vartheta$-pulses of
length $t_{p}$. $\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ spin polarization
is interrogated in windows between the pulses for total readout times
$t{>}60$s (${>}0.5$M pulses). (C) ${}^{13}\textrm{C}$ spin texture is
generated and stabilized by a spatially-varying potential $\varphi(r)$ created
upon Floquet driving with $\vartheta{\approx}\pi$ in the presence of the NV
gradient (cf. Sec. III).  Fig. 2: Spin textures via hyperpolarization
injection (State Engineering). (A) Net ${}^{13}\textrm{C}$ spin polarization
under positive hyperpolarization buildup ($\tau_{+}$, red points), followed by
negative hyperpolarization injection after $\tau_{+}{=}60$s ($\tau_{-}$, blue
points) at room temperature. Solid lines represent a biexponential fit. Boxed
region around $t_{\textrm{pol}}{\approx}97$s corresponds to close to zero net
polarization. (B) Normalized spin-lock decays under Floquet drive with
$\vartheta{\approx}\pi/2$, showing (i) signal $S$ and (ii) rotating-frame
phase $\phi_{R}$. Dashed black line: spin-lock decay corresponding to
$\tau_{-}{=}0$ (a in A), displaying long lifetime $T_{2}’{\approx}$93.4s.
Colored lines: Data for different $\tau_{-}$ values (colorbar). For $\tau_{-}$
corresponding to the boxed region in (A), decay profiles exhibit a sharp zero-
crossing at $t{=}t_{\textrm{zc}}$ (II) and associated phase inversion (see
ii). Each trace has ${>}0.5{\times}10^{6}$ data points, with ${\sim}10^{4}$
points in the vicinity (1% variation) of $t_{\textrm{zc}}$. Dark orange line:
$\tau_{-}{=}44.8\,$s slice is emphasized for clarity. A full movie of this
dataset can be viewed at [30]. (iii) Movement of zero-crossing
$t_{\textrm{zc}}$ with $\tau_{-}$, demonstrates the ability to control spin
texture. (C) Schematic representation of spin textures with increasing
experiment progression time $t$ (arrow). Representative points (I-IV) are
marked on the $\tau_{-}{=}44.8$s slice in (B). NV is located at $r{=}0$.
Positive (negative) polarization is shaded red (blue). Electron-driven
dissipation during $t$ is indicated by a decrease in polarization (white
region) around $r{=}0$. Spin diffusion is indicated by growing texture size
and increasing white region at the boundary between polarization layers of
different sign, with dashed line as a guide to the eye for the zero
polarization regime. (D) ${}^{13}\textrm{C}$ NMR spectrum obtained via Fourier
transforms of full data in (B). Data is shown for varying $\tau_{-}$ (vertical
axis) and plotted on a logarithmic scale, comprising $91$ $\tau_{-}$ slices,
separated by $0.2\,$s. It displays a ${>}10^{9}$ variation in intensity
dynamic range (color bar). Signature of spin shell texture formation is the
narrowing of the FT spectrum with $\tau_{-}{\approx}37$s, corresponding to the
boxed region in (A).  Fig. 3:  Melting of spin-texture due to spin diffusion
(State Engineering). (A) Experiment schematic. Waiting period
$t_{\textrm{wait}}$ at high-field is introduced after successive
$\\{\tau_{+},\tau_{-}\\}$ spin injection, and prior to application of the
Floquet drive in Fig. 2. (B) Representative signal traces showing changes in
decay profiles with variable $t_{\textrm{wait}}$ for texture generated with
$\tau_{+}{=}60\,$s and $\tau_{-}{=}42\,$s (bolded trace in Fig. 2B). Decrease
in signal amplitude between $t_{\textrm{wait}}{=}0$ and
$t_{\textrm{wait}}{=}30\,$s evidences the instability of state-engineered spin
texture. A full movie of this data set, with phase information, can be found
at [31] (C) Schematic representation showing the melting of spin texture
during $t_{\textrm{wait}}$ due to spin diffusion at the boundary between
polarization layers, dashed line serves as guide to the eye for the domain
wall boundary. (D) Points show movement of the zero-crossing $t_{\textrm{zc}}$
position with $t_{\textrm{wait}}$ for data in (B). Solid line is a linear fit.
(E) Signal intensity at $t{=}50$s in (B) plotted for different values of
$t_{\textrm{wait}}$ normalized to its value corresponding to the
$t_{\textrm{wait}}{=}0$ case. Melting of the spin texture due to diffusion
manifests as a decrease in signal. Solid line is a biexponential fit. Light
blue line corresponds to case without spin texture ($\tau_{+}{=}60s$,
$\tau_{-}{=}0$). (F-G) Simulations corresponding to (D-E) showing zero-
crossing times and $|\mathrm{min}(\mathcal{I}_{x})|$ extracted from numerical
time-evolution in a one-dimensional short-range model using a similar shell-
like initial state (see Fig. 5 and Methods). Simulations show qualitative
agreement with experimental results.  Fig. 4: Robust spin textures by
Hamiltonian engineering. (A) Experiment schematic. Spins are hyperpolarized
for $\tau_{+}{=}90$s and subject to a spin-locking train with
$\vartheta{\approx}\pi$. (B) Measured signal $S$ undergoes sharp zero-crossing
and associated sign inversion (phase signal is analogous to Fig. 2B(ii) and
not shown). Colored lines show variation with $\vartheta$ (see colorbar). Each
trace has ${\approx}0.5$M points. Dark purple trace highlights representative
data at $\vartheta/\pi{=}0.94$ for clarity, with sharp zero-crossing at
$t_{\textrm{zc}}$. (C) Analogous signal to bolded purple trace in (B), with
small frequency offset (see SI S5). Signal is plotted on a logarithmic scale
and extends for $t{>}150$s. Dramatic signal zero-crossing is evident. (D)
Schematic representation of formed spin texture for key points of the bolded
trace in (B) (marked I-IV). Signal zero-crossing arises due to thermalization
to an effective Hamiltonian $\mathcal{H}_{\textrm{eff}}$ bearing spatial
texture arising from the NV-imposed gradient (Fig. 1C). Dashed black line
indicates the domain boundary at $r_{c}$. Spin texture remains robust against
spin diffusion with $t$ (see Fig. 5A). (E) Movement of zero-crossing with
$\vartheta$. 2D color plot showing logarithmic scale visualization of data in
(B) plotted with respect to $\vartheta$ (horizontal slices). Zero crossing
appears as an abrupt decrease in signal (colored blue). $\vartheta{=}\pi$
slice is marked and corresponds to a rapid signal decay. $t_{\textrm{zc}}$
occurs at later times for smaller $\vartheta$. Point III corresponding to
zero-crossing in (B) is marked. (F) Numerical simulations performed with LITE
for the simplified short-range model show good qualitative agreement with the
experimental data (see Fig. 5 and SI S8).  Fig. 5:  Simulations of spin
texture formation and evolution using the effective model Hamiltonians for (A)
State Engineering and (B-D) Hamiltonian Engineering. Panels (A,B) are obtained
from 1D short-range simulations for infinitely extended systems including
effects of dissipation using the LITE algorithm (see SI S8.1). System extends
in both directions, but only positive values are shown. Horizontal axis
displays interrogation time $t$ in $J^{-1}$ units, vertical axis displays
distance $r$ from NV in arbitrary units (a.u.). (A) Spin texture via State
Engineering. We imprint the spin texture in the initial state using 11
positively polarized spins close to the NV center at $r=0$ followed by 18
negatively polarized spins. (i) Colors display polarization
$\mathcal{I}_{r}^{x}(t)$ (colorbar) at different sites at time $t$ (vertical
slices). Panel displays spreading dynamics with $t$. Formed shells melt under
diffusion (see SI S8.3). (ii) Integrated polarization over the spin ensemble,
corresponding to signal $S$ measured in Fig. 2. Simulations reveal formation
of a sharp zero-crossing in either case ($t_{\textrm{zc}}$ marked), occurring
at the instance of zero total net polarization. (B) Spin texture via
Hamiltonian engineering. The initial state at $t{=}0$ contains $11$ positively
polarized spins within a background of spins in a fully mixed state. (i)
Polarization $\mathcal{I}_{r}^{x}(t)$ (colorbar) showing that spin texture
forms via thermalization. Late time behavior (dashed line) follows energy
diffusion ${\propto}\sqrt{t}$. Shell critical radius $r_{c}$ (horizontal
dashed line) is stabilized. (ii) Integrated ensemble polarization showing
formation of zero-crossing analogous to Fig. 4. (C, D) Formed shells in three-
dimensions via Hamiltonian Engineering, obtained from late-time dynamics of
classical three-dimensional long-range simulations (see SI Fig S28 for full
time traces). For the classical simulation dissipation is not considered and
spins within the frozen core (white region in D) are not simulated.

##  II Spin texturing via hyperpolarization injection

State-Engineered shells are created by alternately injecting positive and
negative hyperpolarization [32, 33] for periods $\tau_{+}$ and $\tau_{-}$,
respectively (Fig. 2A inset). Direct ${}^{13}\textrm{C}$ hyperpolarization
occurs rapidly over a short range, with distant nuclei polarized more slowly
via spin diffusion. Consequently, ${}^{13}\textrm{C}$ nuclei close to the NV
center are negatively polarized, while more distant nuclei are left positively
polarized, yielding the shell-like texture. Varying $\tau_{-}$ at a fixed
$\tau_{+}$ provides a means to tune the shell size.

Spin diffusion is driven by internuclear dipolar interactions with the
Hamiltonian

$\mathcal{H}_{\textrm{dd}}=\sum_{k<\ell}b_{k\ell}\left(3I_{k}^{z}I_{\ell}^{z}-\bm{I}_{k}\cdot\bm{I}_{\ell}\right),$
(1)

where
$b_{k\ell}{=}J_{\mathrm{exp}}(3\cos^{2}(\theta_{k\ell})-1)/{r}_{k\ell}^{3}$,
$J_{\mathrm{exp}}{=}\mu_{0}\hbar\gamma_{n}^{2}/4\pi$,
$\cos(\theta_{k\ell}){=}\bm{B}_{0}\cdot\mathbf{r}_{k\ell}/(|\bm{B}_{0}|r_{k\ell})$,
$\mathbf{r}_{k\ell}$ is the inter-spin vector,
$\bm{B}_{0}=B_{0}\hat{\textbf{z}}$ is the external field. The median
${}^{13}\textrm{C}$ coupling strength is $\left\langle
b_{kl}\right\rangle{=}J{\approx}0.6\,$kHz.

Fig. 2A depicts the net ${}^{13}\textrm{C}$ polarization during spin
injection, for $\tau_{+}{=}60\,$s, after which the sign of the injected
polarization is reversed (dashed vertical line). In the boxed region around
$t_{\textrm{pol}}{\approx}100\,$s ($\tau_{-}{=}40\,$s,) the total net
polarization approaches zero. Although there is limited net polarization, the
system can still exhibit locally large spin expectation at different sites. To
observe this, the ${}^{13}\textrm{C}$ nuclei are subject to a Floquet protocol
at $B_{0}{=}7\,$T (Fig. 2B(i) inset) involving a series of spin-locking
$\vartheta$-pulses of duration $t_{p}$, at Rabi frequency $\Omega$, and
separated by time interval $\tau$. Spin-locking with $\vartheta{\napprox}\pi$
yields the leading order effective Floquet Hamiltonian [29] (see SI S7.2),

$\mathcal{H}_{\mathrm{SL}}=-\frac{1}{2}\sum_{k,\ell}b_{k\ell}\left(3I^{x}_{k}I^{x}_{\ell}-\bm{I}_{k}\cdot\bm{I}_{\ell}\right)\,.$
(2)

Here ${>}560,000$ pulses are applied over a period of $t{\gtrsim}60\,$s. The
drive induced dynamics (quasi-)conserves $\hat{\textbf{x}}$-polarization as
evident from the leading order effective Hamiltonian Eq. (2) (see SI S7.2),
resulting in long transverse lifetimes $T_{2}^{\prime}$ [34].

This process can be quasi-continuously tracked in real-time. Between the
pulses, rapid and non-destructive (inductive) interrogation of
${}^{13}\textrm{C}$ Larmor precession occurs at a rate of
$\tau^{-1}{\sim}10\,$kHz. We record the nuclear polarization amplitude $S$ and
phase $\phi_{R}$ in the $\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane in
the rotating frame, where $\phi_{R}{=}0(\pi)$ refers to a vector along
$\hat{\textbf{x}}(-\hat{\textbf{x}})$.

Each point in Fig. 2A can therefore be expanded into a secondary dimension
$(S,\phi_{R})$, as shown in Fig. 2B, considerably increasing the information
content compared to previous approaches [35]. Consider first the point (a) at
the top of the polarization buildup curve in Fig. 2A ($\tau_{-}{=}0$).
Considering normalized signal $S$ in Fig. 2B (black dashed line), a long-lived
decay bearing lifetime $T_{2}^{\prime}{\approx}93.4\,$s is evident, reflective
of prethermalization to $\mathcal{H}_{\mathrm{SL}}$ and the resulting
$\hat{\textbf{x}}$ (quasi-)conservation [34]. In comparison to this slow,
featureless, decay, normalized signal profiles, the zero total net
polarization region (boxed in Fig. 2A) exhibits distinct differences. For
clarity, we consider first the representative trace (bold orange) in Fig.
2B(i) at $\tau_{-}{=}44.8\,$s. Normalized $S$ here features a sharp zero-
crossing at $t{=}t_{\textrm{zc}}$(${\approx}11.34\,\mathrm{s}$), accompanied
by a simultaneous reversal in total net polarization sign from
$-\hat{\textbf{x}}$ to $+\hat{\textbf{x}}$ (Fig. 2B(ii)). The region in the
zero-crossing vicinity (1% variation) itself encompasses ${\sim}10^{4}$
points, showcasing the rapid and non-invasive dynamics sampling —
$\tau^{-1}{>}14J$ in these experiments. Additionally, the data captures
dynamics up to very long times, $Jt{>}1.2{\times}10^{6}$, surpassing previous
experiments by several orders of magnitude [36, 37, 38].

The other colored traces in Fig. 2B(i) show corresponding signals $S$ for
different $\tau_{-}$ (see colorbar). The zero-crossing point $t_{\textrm{zc}}$
shifts to the right with increasing $\tau_{-}$. The extracted
$t_{\textrm{zc}}$ values are elaborated in Fig. 2B(iii). A movie showcasing
data for 151 changing values of $\tau_{-}$ is accessible at [30].

Quasi-conservation of the net $\hat{\textbf{x}}$-polarization naively suggests
constant-in-time signal curves, in contradiction to our observations. To
rationalize the data in Fig. 2, recall that the NV electron is strongly
coupled to a phonon bath while the ${}^{13}\textrm{C}$ spins are only weakly
coupled to it. Consequently, the NV acts as a dominant local relaxation source
for the ${}^{13}\textrm{C}$ nuclei. Proximal nuclei dissipate polarization
(${\sim}1/r^{6}$) faster than more distant ones, as can be derived with the
Lindblad master equation in the strong coupling limit (see Methods and Sec.
S6.2). Measurements at longer times $t$, therefore, serially probe nuclei
further away from the NV.

Fig. 2C schematically represents the polarization distribution during the
interrogation period $t$, focusing on specific points (I-IV) along the bold
orange trace ($\tau_{-}{=}44.8$ s) in Fig. 2B. For simplicity, we depict the
generated texture as spherical shells, although in reality, it possesses an
angular dependence inherited from the hyperfine interaction (see Fig. 5 C-D).
Starting from the initial texture ($t{=}0$), proximal ${}^{13}\textrm{C}$
polarization undergoes dissipation from the NV center (shown white), gradually
revealing polarization at greater distances as $t$ increases. Zero-crossing at
$t{=}t_{\textrm{zc}}$ corresponds to an equal distribution of positive and
negative polarization (Fig. 2C II). Further dissipation leads to a
polarization sign inversion (Fig. 2C III-IV). Overall, Fig. 2 illustrates the
ability to discriminate spins without relying on the electronic “frozen core”,
constituting a departure from previous work [39, 35].

The different stages in Fig. 2C exhibit distinct signatures in the
${}^{13}\textrm{C}$ NMR spectrum. In Fig. 2D we report the ${}^{13}\textrm{C}$
spectrum for 91 $\tau_{-}$ values, separated by $0.2\,$s intervals, obtained
by applying a Fourier transform (FT) to the sign-corrected data in Fig. 2B.
The Fourier intensity is shown on a logarithmic scale spanning over nine
orders of magnitude. The wide dynamic range reflects the high signal-to-noise
ratio in our experiments. ${}^{13}\textrm{C}$ spins closer to the NV produce
broader spectral lines because they experience faster relaxation, manifesting
as stronger contributions to the spectral wings around 0Hz. Conversely, more
distant ${}^{13}\textrm{C}$ nuclei are centrally located in the spectrum due
to longer relaxation times. With increasing $\tau_{-}$, depolarization
initially affects the spectral wings, resulting in an apparent narrowing of
the spectrum at $\tau_{-}{\approx}37\,$s. Subsequently, the inversion of the
central feature, corresponding to bulk ${}^{13}\textrm{C}$ nuclei, follows
suit.

Spin texture generated via State Engineering (Fig. 2) is not intrinsically
stable. Over time, any imprinted domain boundaries begin to dissolve due to
spin diffusion. Such “melting” of spin texture can be observed in experiments
where a delay period $t_{\textrm{wait}}$ (${<}30\,$s) is introduced before the
(shell-like) $\hat{\textbf{z}}$-polarized state is rotated into the
$\hat{\textbf{x}}$-$\hat{\textbf{y}}$ plane and the Floquet driving is started
(see Fig. 3A). During $t_{\textrm{wait}}$, NV-driven dissipation of the
$\hat{\textbf{z}}$-polarized texture is negligible. This is due to the
increased energy gap ($\gamma_{n}B_{0}$, as opposed to $\Omega$, see SI Sec.
S2.1), and is reflected in the extremely long nuclear lifetimes
$T_{1n}{\approx}1\,\mathrm{h}{\gg}T_{2}^{\prime}$ under these conditions.
Consequently, spin dynamics in the $t_{\textrm{wait}}$ period is primarily
affected by spin diffusion. Simulations based on a Lindblad master equation
confirm this picture (see Methods).

Fig. 3B shows the signal $S$ dependence on the delay period
$t_{\textrm{wait}}$, starting from a spin texture produced with
$\tau_{+}{=}60\,$s and $\tau_{-}{=}42\,$s. Spin diffusion gradually
homogenizes (flattens) the polarization distribution as $t_{\textrm{wait}}$
increases (Fig. 3B). This leads to a rightward shift in the zero-crossing
point $t_{\textrm{zc}}$ and a decrease in signal amplitude, as evident in Fig.
3B, and as schematically depicted in Fig. 3C.

The spin diffusion-mediated flattening of the polarization distribution can be
directly observed by considering the signal decreasing with
$t_{\textrm{wait}}$ for a fixed value of $t$. Fig. 3E shows
$S(t_{\textrm{wait}})$ at $t{=}50\,$s, normalized against the value at
$t_{\textrm{wait}}{=}0$ to highlight changes with $t_{\textrm{wait}}$. For a
homogeneously polarized state ($\tau_{-}{=}0$), the thus normalized signal
increases slightly with increasing $t_{\textrm{wait}}$ due to polarization
diffusing away from the NV, and hence being subjective to lower effective
relaxation from it. These observations are supported by numerical simulations
(Fig. 3F-G) obtained by solving the Lindblad equation for a simplified toy
model (see Fig. 5 and SI S8.1). All qualitative features and intuition from
the microscopic dynamics are well reproduced.

##  III Robust spin shells by Hamiltonian engineering

To stabilize the generated spin texture, we introduce a second approach of
Hamiltonian Engineering. This stems from a surprising observation: when
deploying the Floquet pulse train with flip angle $\vartheta{\approx}\pi$
(Fig. 4A) and starting with ${}^{13}\textrm{C}$ spins positively polarized
($\tau_{+}{=}60\,$s, $\tau_{-}{=}0\,$s), we observe that over time, the signal
$S$ exhibits a sharp zero-crossing and subsequent sign inversion that can
persist longer than $180\,$s. This is highlighted in representative traces in
Fig. 4B-C, on linear (B) and logarithmic (C) scales, respectively, both with
${>}10^{5}$ points. The emergence of spin shells here is counterintuitive.
Since for $\vartheta{\approx}\pi$ no conservation law shields the initial
state from rapid heat-death, one would expect the interacting nuclear spin
system to quickly relax to a featureless infinite-temperature state [40, 22].
Despite extensive use of $\pi$-trains (CPMG experiments [41, 42]) in various
contexts including dynamical decoupling and quantum sensing, to the best of
our knowledge, this phenomenon has hitherto not been previously reported.

We attribute the emergence of the long-lived signal to Hamiltonian engineering
facilitated by the simultaneous action of the nanoscale electronic field
gradient dressed by the $\pi$-train Floquet drive. In particular, the NV
electron spin, thermally polarized (${\approx}12.6\,\%$ polarized at $9.4\,$T
and $100\,$K), induces a hyperfine gradient $\eta({\mathbf{r}})$ on the
nuclear spins. As a consequence, the direction and magnitude of the Rabi field
$\Omega$ experienced by the nuclei depends on their proximity to the electron
(see Fig. 1A). When $\vartheta{\neq}\pi$ (cf. State Engineering, Fig. 2-Fig.
3), this merely alters the prethermalization axis
$\hat{\textbf{x}}{\rightarrow}\hat{\textbf{x}}^{\prime}(r)$ (see SI S7.2.1).
On the other hand, when $\vartheta{=}\pi$ exactly and we ignore the NV-induced
gradient field, the absence of $\hat{\textbf{x}}$-polarization conservation
leads to a rapid, yet trivial, decay of the spins within $T_{2}^{*}$ [34].

However, in the presence of the NV gradient field $\eta({\mathbf{r}})$, for
$\vartheta{=}\pi+\epsilon$ (with $\epsilon{\ll}\pi)$, spin dynamics is
governed by the effective Floquet Hamiltonian (see SI S7.2.2)

$\mathcal{H}_{\textrm{eff}}=\mathcal{H}_{\textrm{dd}}+\sum_{\mathbf{r}}\varphi_{\mathbf{r}}I_{\mathbf{r}}^{x}\,,$
(3)

where $\varphi_{\mathbf{r}}=\eta({\mathbf{r}})+(\vartheta-\pi)/\tau$ denotes
the effective spatially-varying on-site potential (Fig. 1C). We operate in the
regime of a weak electronic field gradient ($\varphi_{\mathbf{r}}{<}1\,$kHz),
far from the frozen-core limit — a significant departure from previous
experiments [16]. Notably, the spatial inhomogeneity induced by
$\varphi_{\mathbf{r}}$ separates into two regimes by a critical radius
$r_{c}$, where $\varphi_{\mathbf{r}}$ flips sign:
$\varphi_{\mathbf{r}}\lessgtr 0$ for
$r{=}\left|\left|{{\mathbf{r}}}\right|\right|\gtrless r_{c}$. The exact
position of $r_{c}$ is determined by the ${}^{13}\textrm{C}$ slice for which
$\varphi_{\bm{r}}{=}0$ (Fig. 1C). Note that here we ignore the angular
$\theta$ dependence for simplicity and return to it in Fig. 5C-D.

Crucially, the sign change of $\varphi_{\mathbf{r}}$ on either side of $r_{c}$
means that $\mathcal{H}_{\textrm{eff}}$ internally encodes spatial structure,
and $r_{c}$ serves as the domain boundary. Then, irrespective of the initial
state, applying the Eigenstate Thermalization Hypothesis [1, 2, 3], one
predicts that the spins thermalize to the quasi-equilibrium state described by
the density matrix
$\rho_{\textrm{eq}}{\propto}\exp(-\beta\mathcal{H}_{\textrm{eff}}){=}\mathds{1}-\beta\mathcal{H}_{\textrm{eff}}+\mathcal{O}(\beta^{2})$,
with an inverse temperature $\beta$ set by the energy density of the initial
state. The local expectation value of $\hat{\textbf{x}}$-polarization in the
thermal state,
$\mathcal{I}_{\mathbf{r}}^{x}{=}\mathrm{Tr}(I_{\mathbf{r}}^{x}\rho){\approx}-\beta\varphi_{\mathbf{r}}/2+\mathcal{O}(\beta^{2})$,
then attains opposite signs on either side of $r_{c}$ — imprinting a spatial
structure of a quasi-equilibrium shell stabilized within the prethermal
plateau (Fig. 1B).

The highlighted traces in Fig. 4B-C track such shell formation. The initial
polarization profile $\rho_{\textrm{i}}$ is far from the prethermal
equilibrium $\rho_{\textrm{eq}}$: it possesses a finite amount of energy with
respect to $\mathcal{H}_{\textrm{eff}}$ in the thermodynamic limit, which is
initially localized within the polarized region. Over time this energy
diffuses through the system due to energy (quasi-)conservation (see SI
S8.2.3). Consequently, spins in the region $r{>}r_{c}$ begin to sequentially
flip, driving the system, towards a quasi-equilibrium,
$\rho_{\textrm{i}}{\rightarrow}\rho_{\textrm{eq}}$ (see schematic depiction in
Fig. 4D). As the negative polarization in the outer region ($r{>}r_{c}$)
surpasses the positive polarization in the inner region ($r{<}r_{c}$) in
magnitude, the total polarization along the $\hat{\textbf{x}}$ direction
undergoes a sign inversion, resulting in a zero-crossing. Importantly, since
the spin texture emerges through thermalization and total polarization is not
conserved, the domain boundary $r_{c}$ remains stable (on prethermal
timescales). This stands in contrast to State Engineering (see Fig. 5).
Furthermore, the signal drop observed here is orders of magnitude larger
compared to State Engineering, since we do not operate in a regime of nearly
zero total net polarization, indicating a larger spin texture gradient
(Methods).

Experiments in Fig. 4B-C, therefore, uniquely enable real-time observation of
spin flipping dynamics during the prethermalization process. This capability,
combined with the large accessible values of $Jt\,{>}10^{6}$, allows us to
investigate thermalizing spin dynamics with exceptional resolution over
extended time periods, surpassing the capabilities of previous experiments by
many orders of magnitude [36, 37, 38].

Colored traces in Fig. 4B show similar experiments for different values of
$\vartheta$ in the vicinity of $\vartheta{=}\pi$ (see colorbar). The zero-
crossing position is observed to change with $\vartheta$. A movie of this data
is accessible in [43]. Fig. 4E recasts this movie for 120 different values of
$\vartheta$ in a 2D plot on a logarithmic scale in intensity, where the darker
blue colors highlight the zero crossing. This visualization clearly elucidates
the movement of $t_{\textrm{zc}}$ with $\vartheta$. Additionally, a distinct
slice in the plot displays rapid signal decay (indicated by the dashed line).
This is attributed to bulk ${}^{13}\textrm{C}$ nuclei, which are far removed
from the NV influence, corresponding for which $\vartheta{=}\pi$. The
qualitative behaviour can be reproduced in numerical simulations (see Methods
and Refs. [44, 45]). While the exact shape of the zero-crossing arc depends on
the details of the model (such as $\varphi_{\mathbf{r}}$), the simulations
capture the diverging behaviour of $t_{\textrm{zc}}$ found for
$\vartheta\to\pi^{-}$.

Based on data in Fig. 4, we estimate (see Methods) that the domain spans
${r_{c}{\approx}2.7\,\mathrm{nm}\,\cdot\,\sqrt[3]{|3\cos^{2}\theta-1|}}$
encompassing ${\approx}150$ spins for the representative line in Fig. 4B
($\vartheta{=}0.94\,\pi$) 222On the experimental timescale ($t{=}180\,$s),
total polarization spreading is expected to remain within a radius
${\lesssim}12\,$nm; the individual NV-${}^{13}\textrm{C}$ systems can be
considered isolated.. We are able to exert control over the $r_{c}$ domain
boundary by changing the Floquet engineered Hamiltonian (3). By adjusting the
kick-angle $\vartheta$ we can modify the effective on-site potential
$\varphi_{\bm{r}}$, and hence $r_{c}$, allowing us to tune the shell size
between tens (${\approx}50$) up to mesoscopic numbers (${\gtrsim}300$) of
spins (see SI Sec. S7.4) [47]. We observe that the movement of
$t_{\textrm{zc}}$ matches the theoretical prediction (Fig. 4E,F). Varying Rabi
frequency $\Omega$ or the frequency offset can modify the Hamiltonian and
allow furthur control over $r_{c}$ (SI Sec. S5). Furthermore,
thermalization to shell-like texture is found to be robust to lattice
orientation (SI Sec. S5.2) or initial states employed (SI Sec. S5.1).

##  IV Numerical Results

Quantum simulations here use Local-Information Time Evolution (LITE) [44, 45],
a recently developed algorithm suited to investigate diffusive and dissipative
quantum dynamics. Numerical tractability constrains us to work with a
simplified one-dimensional toy-model Hamiltonian, but which captures all
relevant features of $\mathcal{H}_{\mathrm{eff}}$ (Methods).

Fig. 5A(i) shows simulations corresponding to State Engineering experiments.
Time is represented on the horizontal axis in units of $J^{-1}$, the vertical
axis represents the distance from the NV (both axes are on a logarithmic
scale), and the colorbar shows the local $\hat{\textbf{x}}$-polarization
$\mathcal{I}_{r}^{x}$. Considering an initial shell-like texture produced by
spin injection, we observe that the subsequent dynamics are characterized by
melting domain walls, causing the magnetization gradient to diminish over
time. Fig. 5A(ii) shows the total net polarization
$S{=}|{\sum_{r}\mathcal{I}_{r}^{x}}|$ — analogous to that measured in
experiment. We observe a zero-crossing similar to Fig. 2B. Notably, in State
Engineering, we rely on dissipation to enable the observation of the sign-
inversion of the total polarization: in the absence of dissipation, total
$\hat{\textbf{x}}$ polarization is (quasi)-conserved and thus will appear
constant on prethermal timescales.

In contrast, Fig. 5B considers Hamiltonian Engineering, and polarization
initially confined in the vicinity of the NV. Over time, by virtue of energy
diffusion, the spins prethermalize to a polarization gradient under the NV-
induced potential. Peripheral spins become endowed with negative polarization
(blue region), and ultimately, over time, the total polarization inverts as
the majority of spins feature negative polarization. Importantly, the domain
$r_{c}$ (dashed horizontal line) separating the local positive and negative
polarization regions is solely set by the on-site potential
$\varphi_{\mathbf{r}}$ present in $\mathcal{H}_{\mathrm{eff}}$ and is,
therefore, stationary. As the spins thermalize, the (non-conserved)
polarization adapts to the stable gradient profile imposed by
$\mathcal{H}_{\mathrm{eff}}$. Irrespective of the initial state, the
diffusion-propelled dynamics therefore induce a polarization gradient with a
time-invariant domain boundary.

Considering the net polarization (Fig. 5B(ii)), we observe a steep zero-
crossing in qualitative agreement with experiments. We note that the zero-
crossing $t_{\textrm{zc}}$ arises even in the absence of dissipation (see SI
S8.2), as observed in the long-time slices in Fig. 5B(i). Dissipation only
serves to extinguish polarization at distances close to the NV (i.e.,
accelerating the inevitable diffusion-induced polarization inversion).

To demonstrate that these results hold beyond one-dimensional short-range
models, we additionally perform classical simulations based on a three-
dimensional long-range model comprising ${\approx}10^{3}$ dipolar-interacting
spins on a diamond lattice (Methods). In contrast to the one-dimensional
quantum simulations, here the system is finite and free of dissipation so that
equilibrium is reached in finite time. Notably, here, the system thermalizes
to a three-dimensional spin texture which exhibits the analytically predicted
$\theta$-dependence of the NV gradient field with $r_{c}$ spanning several
nanometers (cf. Sec. III).

##  V Discussion & Outlook

Our experiments introduce several novel features. The creation, stabilization,
and observation of spin texture are achieved through global spin control and
readout. Despite this, a good degree of local nuclear discrimination is shown
to be attainable by utilizing an electron as a controllable spin injector,
gradient source, and dissipator. Notably, unlike other experiments [16], our
approach is not confined to the constraints presented by diffusion barrier
limits, enabling examination of a mesoscopically large number of nuclear spins
around each electron. The injected hyperpolarization exceeds the Boltzmann
polarization levels by orders of magnitude compared to previous studies [39].

Our method of continuous interrogation in the rotating frame, as opposed to
point-by-point probes in the lab-frame [39], offer a distinct methodological
advantage. It facilitates a rapidly sampled ($J\tau{<}0.07$) visualization of
polarization dynamics while reaching into the very-long-time regime
($Jt{>}1.2{\times}10^{6}$). The latter is orders of magnitude beyond previous
experiments [36, 37, 38], permitting direct observation of emergent
stabilization. The dynamics observed complement recent experiments with cold
atoms [48, 49, 50, 5, 51], but occur in a distinct and novel context: at the
nanoscale, in the solid-state, and in the limit of strongly interacting
dipolar-coupled spins (for which $JT_{2}^{\ast}{\approx}1$).

Our dynamical stabilization protocol operates distinctly out of equilibrium.
The dynamics is governed by an emergent quasi-conservation of energy and
explicitly breaks polarization conservation, allowing us to inhibit diffusion
in the spin polarization channel. In conjunction with the electronic hyperfine
field, this protocol stabilizes localized, shell-like spin textures at a
finite energy density, even in the absence of local control. Note that no
obvious static protocol could yield similar outcomes in our system. The
underlying physics of this dynamical stabilization deviates from traditional
state engineering approaches, is versatile, and can encompass a variety of
models (short- or long-range, clean or weakly-disordered interacting systems,
quantum or classical models), across different dimensions.

The elucidated Hamiltonian Engineering protocol, therefore, highlights the
untapped capabilities of nonequilibrium control techniques in manipulating
quantum matter. For example, our approach can be utilized within quantum
simulators to induce a magnetization domain wall in a close to infinite-
temperature state – the starting point to investigate sub/superdiffusive spin
dynamics [52, 5, 53]. Depending on the structure of the effective Hamiltonian
and the initial state, one can also stabilize and investigate spin textures at
negative-temperature [54]. Our work, therefore, presents a first instance of a
far-reaching idea that demonstrates the applicability of concepts from the
emerging field of nonequilibrium quantum dynamics, to engineer stable many-
body states with tailored attributes [55].

Lastly, the manipulation and control of the orientation of spin polarization
at sub-nanometer length scales itself, opens up several promising future
directions. The use of ubiquitously occurring nuclear spins, as elucidated in
this work, broadens the application of spin texturing to diverse systems,
moving beyond previously considered magnetic materials [14, 56, 15, 57, 12,
58]. This suggests applications in quantum memories [12, 13], spintronics [14,
15], and spatiotemporal quantum sensing harnessing hyperpolarized nuclei as
sensors [59]. Controllable spin textures may also be applied to non-invasive
nanoscale chemical imaging in materials science and biology. We envision
employing targeted electron spin labels to map the radial spin densities of
different nuclei (e.g., ${}^{1}\textrm{H}$ and ${}^{13}\textrm{C}$) within
molecules.

##  VI Acknowledgements

We thank J. H. Bardarson, T. Klein Kvorning, R. Moessner, C. Ramanathan, J.
Reimer, D. Suter for insightful discussions, and J. Mercade (Tabor
Electronics) for technical assistance. This work was supported by ONR
(N00014-20-1-2806), AFOSR YIP (FA9550-23-1-0106), AFOSR DURIP
(FA9550-22-1-0156), a Google Faculty Research Award and the CIFAR Azrieli
Foundation (GS23-013). CF and CA acknowledge funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (Grant Agreement No. 101001902). The computations with the
LITE algorithm were enabled by resources provided by the National Academic
Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National
Infrastructure for Computing (SNIC) at Tetralith partially funded by the
Swedish Research Council through grant agreements no. 2022-06725 and no.
2018-05973. Classical simulations were performed on the MPI PKS HPC cluster.

##  VII Author Contributions

KH, CF, ND, PMS and DM contributed equally to this work. Order in which these
authors are listed was decided by dice roll.

DM, KH, ND, MM and ED built the experimental apparatus. ND, DM, KH, SB and AA
collected data and performed data analysis. KH, ND and DM conceived of
experimental design for Hamiltonian Engineering, and QRF, WB, XL and AA
conceived of experimental design for State Engineering. UB, AP, ND, KH and AN
made data visualizations.

PMS and CF developed the Floquet analysis. PMS did the numerically exact
quantum simulations and the classical simulations. CF did the approximate
quantum simulations using LITE and the diffusion analysis. CA performed the
open quantum systems analysis and contributed to LITE simulations. MB
supervised the theory work. PMS, CF, CA, and MB helped interpreting and
analyzing the numerical and experimental data. AA, MB, PMS, CF, CA, ND and KH
wrote the paper.

##  VIII Methods

EXPERIMENT

###  VIII.1 Sample

Experiments here employ a single-crystal diamond sample (3x3x0.3mm) from
Element6. It contains ${}^{13}\textrm{C}$ spins at natural abundance
$(1.1\,\%)$, and ${\approx}1\,$ppm NV- concentration, corresponding to inter-
NV spacing of ${\sim}25\,$nm, and each NV center has $10^{4}{}^{13}\textrm{C}$
surrounding it. The sample also hosts ${\gtrsim}20\,$ppm of substitutional
nitrogen impurities (P1 centers). For experiments in Fig. 2-Fig. 4, the sample
is oriented such that its [100] face is approximately parallel to $B_{0}$.
However, as elucidated in SI, Sec. S5.2, these results are qualitatively
independent of the sample orientation.

###  VIII.2 Experimental Apparatus

Fig. 6: Novel instrumentation for cyrogenic DNP. (A) Device construction. CAD
model showing instrument and highlighting main components. It consists of a
$9.4\,$T superconducting magnet, surrounded by an aluminum frame to which a 4K
cryostat is mounted on a belt-driven actuator, bearing a high-torque motor.
Laser illuminates sample from the bottom. (B) Photograph showing key features
of the instrument. Cryostat is shown at low-field position
$\textbf{B}_{\textrm{pol}}$; actuator truck and mounts are visible. Fig. 7:
Probe for ${}^{13}\textrm{C}$ hyperpolarization and readout (A) CAD model
showing the schematic of the NMR/DNP probe. It fits snugly within the
cryostat, providing additional RF shielding. ${>}3$ft long tuning rods extend
through the top of the probe and allow for RF cavity tuning and impedance
matching, even under vacuum and cryogenic conditions. (B) Photograph of
constructed probe. A close-up of coil arrangement is shown, highlighting the
RF saddle coil and centrally placed microwave loop employed for
hyperpolarization. Aperture at probe base enables optical access to the
sample.

State Engineering experiments (Fig. 2-Fig. 3), carried out at room
temperature, employ an apparatus described before in [60]. For Hamiltonian
Engineering (Fig. 4), on the other hand, we introduce a novel instrument (Fig.
6) for ${}^{13}\textrm{C}$ hyperpolarization and interrogation at cryogenic
temperatures. Low temperatures allow access to higher Boltzmann electronic
polarization, consequently stronger electronic gradient fields, and slower
electronic relaxation rates. The instrument supplies for the first time (to
our knowledge) a method for “cryogenic field cycling” for dynamic nuclear
polarization (DNP), allowing simultaneous operation at variable fields
(1mT-9.4T), and controllable cryogenic temperatures down to 4K (although we
restrict ourselves to 77K in this work).

The device utilizes an Oxford SpectrostatNMR cryostat under continuous flow
cryogenic cooling, which is mechanically moved (“shuttled”) from lower (few
mT) fields into a $B_{0}{=}9.4\,$T NMR magnet (Oxford). The low-field position
situated 640mm above the magnet center ($\textbf{B}_{\textrm{pol}}{=}36\,$mT)
is employed for hyperpolarization; the cryostat is then shuttled to high-field
($9.4\,$T) where the ${}^{13}\textrm{C}$ nuclei are interrogated. Shuttling
occurs via a belt-driven actuator (Parker) powered by a motor (ACS) fitted
with a high-torque gearbox for enhanced load-bearing capacity for the heavy
(25lb) cryostat. The actuator carries a movable stage to which two custom-
designed clamps secure the cryostat. A 1.6 m flexible transfer line allows for
continuous cooling during shuttling and operation over ${>}$1 week. Shuttling
occurs at 7 mm/s and takes ${\approx}90\,$s. Long ${}^{13}\textrm{C}$ $T_{1}$
lifetimes (${>}1$hr) at fields exceeding 0.1T, mean that there is minimal loss
of polarization during shuttling. We note that room temperature experiments
(Fig. 2-Fig. 3)), in contrast, involve shuttling in ${<}1$s.

The diamond sample is secured at the bottom portion of the cryostat in a
custom-built NMR/DNP probe (Fig. 7). The probe is top-mounted in the cryostat
and includes two coils: a loop through which microwaves (MWs) are applied for
DNP, and a saddle coil for NMR detection. A novel arrangement employing
O-rings at the probe top allows the NMR coil to be frequency-tuned and
impedance matched without breaking vacuum.

For optical DNP at $\textbf{B}_{\textrm{pol}}$, an optical window at the
cryostat bottom permits illumination by a laser ($\lambda{=}532$nm, Coherent)
directed into the bore using a 45∘ mirror. A pair of piezo-driven Zaber
mirrors ensure optimal alignment into the cryostat center, aided by a camera
mounted with a 637nm long-pass filter at the top of the cryostat. A TTL-
triggered mechanical shutter controls illumination timing to within 1$\mu$s.

Hyperpolarization employs MW chirps applied across the NV EPR spectrum. MWs
are generated by a Tabor Proteus arbitrary waveform transceiver (AWT) and
gated by a Mini-Circuits ZASWA-2-50DR+ switch. MWs are amplified in two stages
by ZHL 2W-63-S+ and ZHL-100W-63+ amplifiers. A Varian VNMRS console generates
the RF pulses at the ${}^{13}\textrm{C}$ Larmor frequency (${\approx}$100MHz),
while the NMR signal is detected in windows between the pulses, filtered, and
amplified by a Varian preamplifier, and digitized by the AWT. We apply
${>}0.5$M pulses, yielding a continuously interrogated NMR signal for up to
$t{=}180$s that is non-destructive since spins are only weakly coupled to the
RF coil. Timing of all events — lasers, MW application, mechanical shuttling,
triggering NMR detection, and signal digitization — is synchronized by the
Swabian pattern generator, controlled by MATLAB.

###  VIII.3 Hyperpolarization Methodology

Hyperpolarization follows [33]. Application of a continuous-wave 532nm green
laser (5W) induces polarization in the NV electrons to the $m_{s}{=}0$ state
at $\textbf{B}_{\textrm{pol}}$. This is transferred to ${}^{13}\textrm{C}$
nuclei through successive traversals of rotating frame Landau-Zener level
anti-crossings. Practically, this is accomplished by utilizing MW chirps
generated by a Tabor Proteus AWT, sweeping across the NV center EPR spectral
bandwidth (25 MHz) at a rate of 200 Hz.

###  VIII.4 Data Acquisition and Processing

Induced ${}^{13}\textrm{C}$ Larmor precession signal into the RF saddle coil
is captured by a Tabor Proteus AWT, and digitized every 1ns in
$t_{\textrm{acq}}$ windows between the pulses. It is decimated to preserve
memory and improve acquisition and processing speed (here 64-fold). A Fast
Fourier Transform (FFT) is applied to extract signal amplitude $S$ and phase
$\phi$ in each $t_{\textrm{acq}}$ window. Phase signals in each window are
linearly offset by the phase accrued during $t_{p}$ pulse periods, a phase
unwrapping algorithm is used to unfold this trivial phase development,
yielding phase $\phi_{R}$ of the spins in the rotating frame. $(S,\phi_{R})$
forms the basis for the experimental data in this work.

Fig. 8: Estimation of flip angle $\vartheta$. (A) Schematic of spin-locking
Rabi experiment, consisting of a train of $\pi/2$ pulses, where the angle of
the first pulse $\vartheta_{\textrm{init}}$, is varied. (B) Typical signals in
case of $\vartheta_{\textrm{init}}{=}\pi/2$ and
$\vartheta_{\textrm{init}}{=}\pi$, with total measurement time of $t{=}70$s,
corresponding to $>0.5$M pulses. (C) Variation of integrated signal
$S_{\textrm{int}}$ with length of the first pulse. Solid line represents
sinusoidal fit. Slight off-resonance in the pulses leads to a finite signal
even at $\vartheta_{\textrm{init}}{=}0$. Fig. 9: Comparison of State and
Hamiltonian Engineering signals. Comparison of signal $S$ for State
Engineering data in main text Fig. 2 (taken at room temperature) and
Hamiltonian Engineering data from main text Fig. 4 (taken at $100\,$K) for
comparable zero crossing times $t_{\textrm{zc}}{\approx}31\,$s. Data is shown
on a logarithmic scale for clarity. Measured signal is over an order of
magnitude greater in the Hamiltonian Engineering method.

###  VIII.5 Estimation of Kick Angle $\vartheta$

To estimate kick angle $\vartheta$, important in Fig. 4, we conduct
${}^{13}\textrm{C}$ Rabi experiments and fit the obtained nutation to a
sinusoidal curve. A spin-locking train of ${\sim}\pi/2$ pulses and continuous
readout (see Fig. 8B) for $t{=}70\,$s is employed to enhance measurement SNR.
The pulse length of the initial pulse $\vartheta_{\textrm{init}}$ is varied,
and the integrated signal is measured, as shown in Fig. 8C. High SNR yields a
high-precision estimate of $\vartheta$, with an error of within $1\,\%$.

###  VIII.6 Comparison of Signals in State and Hamiltonian Engineering

For clarity, we contrast in Fig. 9 the measured signal $S$ amplitude for
experiments of State engineering (Fig. 2-Fig. 3) and Hamiltonian engineering
(Fig. 4). Spin texture created by State Engineering starts from a regime of
low net polarization (Fig. 2A), and sign inversion (characterized by the
$t_{\textrm{zc}}$ zero-crossing), occurs due to electronic dissipation. On the
other hand, Hamiltonian engineering starts with large net polarization, and
sign inversion occurs on account of thermalization into the imposed potential.
As a result, the measured signal, and the magnitude of the zero-crossing drop
in the Hamiltonian engineering method is over one order of magnitude larger
than State engineering (see Fig. 9). The strong signal and robustness of the
domain boundaries highlight the advantages of Hamiltonian engineering for
stable spin texture generation and readout.

###  VIII.7 Additional Results on Hamiltonian Engineering

To further validate the physical picture in Fig. 4 and Fig. 5B for Hamiltonian
engineering, additional studies are conducted to probe the thermalization
process. A summary of these findings is provided here, with further details
available in the SI. First, to confirm that the zero-crossing observed in Fig.
4B does not arise from spins tipping towards the $\hat{\textbf{z}}$-axis where
they are unobservable, we interrupted the Floquet drive with a $\pi/2$ pulse
at $t{=}t_{\textrm{zc}}$. The data (SI Sec. S4) reveal no generation of
$\hat{\textbf{z}}$ during $t$, supporting the thermalization model above.
Next, we investigated the relative roles of diffusion and dissipation leading
up to $t_{\textrm{zc}}$ by studying profiles similar to Fig. 4B with varying
temperatures (Sec. S5.3). Temperature serves as a control parameter for
relaxation, as it strongly lengthens (${\gtrsim}$50-fold) the electronic
$T_{1e}$ while only linearly changing the strength of the electronic gradient.
Decreasing the temperature results in a rightward shift of $t_{\textrm{zc}}$
due to the slower rate of electron dissipation, consistent with theoretical
expectations (Fig. 5B).

THEORY

###  VIII.8 Effective Hamiltonians for State and Hamiltonian Engineering
protocols

As explained in the main text, the experimental system is subject to a Floquet
drive which consists of a periodic train of $\hat{\textbf{x}}$-pulses. When
the period of switching is small compared to the energy scales of the physical
system, this driving induces prethermalization wherein the dynamics is
governed by an effective Hamiltonian [21], before the system eventually heats
up to a featureless infinite-temperature state.

In the SI (Sec. S7), we provide a detailed derivation of the approximate
effective Hamiltonians $\mathcal{H}_{\mathrm{SL}}$, Eq. (2), and
$\mathcal{H}_{\textrm{eff}}$, Eq. (3), that govern the dynamics of the nuclear
spin system for the State and Hamiltonian Engineering approaches,
respectively. In particular, using exact simulations on system sizes up to
$N{=}16$ quantum spins, we obtain an excellent agreement between the dynamics
generated by the effective Hamiltonian and the exact Floquet system in all
temporal regimes of interest (see SI, Sec. S7). Therefore, for the theoretical
analysis in the main text, we work with static effective Hamiltonians. We
emphasize that their applicability is limited to the duration of the
prethermal plateau.

###  VIII.9 Dissipation induced by the NV center

In the strong coupling limit between the NV center and ${}^{13}\textrm{C}$
spins, the dissipation induced by the NV-center can be rigorously derived from
the Hamiltonian of the full system in three dimensions (see Sec. S6.2). Unlike
previous works [61], we integrate out the phonon bath and the NV electron
obtaining a Markovian master equation for the ${}^{13}\textrm{C}$ spins only.
Such a Markovian master equation involves on-site Lindblad jump operators with
coupling constants decaying as $\sim 1/r^{6}$ as a function of the distance
$r$ from the NV center (which is effectively short-ranged). Thus, as a
simplification, in our approximate one-dimensional quantum simulations, we
approximately solve the Lindblad equation with local jump operators $L$ acting
only on the site with index $r=0$ (representing the location of the NV) with
isotropic coupling constants,

$L_{+}=\frac{1}{2}\left(I_{0}^{x}+iI_{0}^{y}\right)\,,\leavevmode\nobreak\
\leavevmode\nobreak\
L_{-}=\frac{1}{2}\left(I_{0}^{x}-iI_{0}^{y}\right)\,,\leavevmode\nobreak\
\leavevmode\nobreak\ L_{z}=I_{0}^{z}.$ (4)

Such operators generate both dephasing and dissipation in the system. Both
effects are produced by the jump operators derived in Sec. S6.2 for the
Hamiltonian of the full system in three dimensions.

###  VIII.10 Numerical simulations

####  VIII.10.1 Quantum simulations

We use the novel algorithm LITE (local-information time evolution, see Refs.
[44, 45]) to simulate the dynamics of the ${}^{13}\textrm{C}$-spins with
respect to the effective Hamiltonians subject to dissipation. LITE is designed
to investigate the out-of-equilibrium dynamics of quantum systems, including
open systems governed by the Lindblad equation. Its adaptive system size
allows us to effectively simulate infinite systems. Motivated by the
experimental setup, the system is initialized in a spatially inhomogeneous
partially polarized state, where a relatively small number of spins in
proximity to the NV center ($r{=}0$) carry a finite partial
$\hat{\textbf{x}}$-polarization, while spins located far from $r{=}0$ are
fully mixed $\rho_{r\gg a}{=}\mathds{1}_{2}/2$.

To use the LITE toolbox effectively, we investigate numerically tractable one-
dimensional short-range toy models akin to the effective (three-dimensional)
Hamiltonians of the coupled ${}^{13}\textrm{C}$ system

$\displaystyle\mathcal{H}_{\pi/2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{r}J_{r}\left(3I^{x}_{r}I^{x}_{r+a}-\bm{I}_{r}\cdot\bm{I}_{r+a}\right),$
$\displaystyle\mathcal{H}_{\pi}$ $\displaystyle=$
$\displaystyle\sum_{r}J_{r}(3I^{z}_{r}I_{r+a}^{z}-\bm{I}_{r}\cdot\bm{I}_{r+a})+\varphi_{r}I_{r}^{x}.$
(5)

The subscripts indicate the corresponding partner Hamiltonian in the actual
system: $\pi/2$ refers to the system kicked with $\vartheta=\pi/2$ resulting
in the effective Hamiltonian $\mathcal{H}_{\mathrm{SL}}$; likewise, subscript
$\pi$ refers to a kick angle $\vartheta\approx\pi$ with the associated
Hamiltonian $\mathcal{H}_{\mathrm{eff}}$. In the simulations,
$J_{r}{=}J_{0}+W$ is taken weakly disordered with $W$ drawn uniformly at
random from the interval $[-0.3|J_{0}|,0.3|J_{0}|]$ and with lattice constant
$a$. The on-site potential follows
$\varphi_{r}{\sim}\frac{1}{r^{3}}-\delta\vartheta$ with $r=0$ corresponding to
the location of the NV center, and $\delta\vartheta$ is the deviation in the
kick angle from $\pi$ (see SI S7.2.2). For the simulations in Fig. 5 we use
$a=0.2/\sqrt[3]{\pi}$, $\delta\vartheta=0.05\pi$, $J_{0}=-0.025$, and
$\gamma{=}0.1|J_{0}|$

While these toy models may appear simplistic compared to the actual
experimental system, as they lack the complexity of higher dimensions and
long-range spin-spin couplings, they capture the essential physics of
diffusion (and dissipation), since they obey the same conservation laws
compared to their three-dimensional counterparts. Therefore, the conclusions
drawn from our numerical results can be qualitatively extended to the
experimental system.

####  VIII.10.2 Classical simulations

In addition, we also performed classical simulations analyzing the dynamics of
the system with hundreds of dipolar-interacting spins of the full three-
dimensional long-range interacting model. Spins are placed randomly on the
vertices of a 3D diamond lattice of a finite extent, with lattice constant
$a{=}0.356\,\mathrm{nm}$, and the classical evolution, generated under the
effective Hamiltonian $\mathcal{H}_{\textrm{eff}}$ (3), is studied. In the
classical limit, the evolution of the spins $\mathcal{I}_{k}$ is described by
(see SI, Sec. S10.1)

$\frac{\mathrm{d}}{\mathrm{d}t}\bm{\mathcal{I}}_{k}(t)=\bm{\mathcal{I}}_{k}(t)\times\bm{\nabla}_{\bm{\mathcal{I}}_{k}}\mathcal{H}_{\textrm{eff}}(\bm{\mathcal{I}}_{k}(t))\,.$
(6)

In particular, for the classical simulation we are only considering closed
system dynamics, i.e. we do not include dissipation as is done in the LITE
simulations. However, as we also comment on later dissipation is not an
essential ingredient for Hamiltonian Engineering.

As in the experiment, we average the simulation data over many (here $100$)
different lattice configurations. For each lattice configuration, we consider
an ensemble of $150$ partially polarized initial states drawn from a spatially
inhomogeneously polarized distribution, where all spins within a radius
$r{<}7\,\mathrm{nm}$ are polarized and all others unpolarized, without
correlations between the spins. We then evolve this ensemble using the
classical Hamilton equations (Eq. (6)), and compute the ensemble-averaged
polarization, see Sec. S10.

###  VIII.11 Details on the Theoretical Results for State Engineering

Fig. 10:  Zero crossing radius for State Engineering. Zoom-in to $Jt<200$
region in Fig. 5A(i) with adjusted color map scale to reveal polarization
dynamics close to NV ($r=0$). The dynamics of positive polarization close to
NV are dominated by dissipation such that it decays to zero before the
crossing radius can move outwards due to diffusion.

####  VIII.11.1 Zero crossing radius for State Engineering

In Fig. 5A(i) the polarization dynamics close to the NV center seem to show an
abrupt sign inversion from positive to negative polarization. However, this is
only an artefact originating from the logarithmic time scale. In Fig. 10 we
show the region close to $r{=}0$ and around the potential sign inversion on a
linear scale. Here, it becomes clear that the strong dissipation close to the
NV leads to a quick decay of positive polarization. Therefore, the positive
polarization cannot diffuse outwards. In fact, the negative polarization
diffuses symmetrically in both directions filling the polarization hole left
by the decayed positive polarization. Thus, the positive polarization does not
abruptly invert; rather it is absorbed by the NV and the negative polarization
can freely diffuse into what initially was a positively polarized regime. In
the absence of dissipation, both positive and negative polarization would
diffuse outwards leading to an increase in crossing radius (see SI Fig. S27).

####  VIII.11.2 Waiting time analysis

Typically, the dissipation strength is indirectly controlled by tuning the
temperature of the sample. Here, the peculiar form of the dissipation obtained
in the singular coupling limit (see SI, Sec. S6.2), where jump operators are
only the spin operators along the $\hat{\textbf{z}}$-axis, opens up another
possibility. While the strength of the dissipation itself remains immutable,
the effect of the Lindblad jump operators on the system strongly depends on
the direction of polarization: spins polarized along $\hat{\textbf{z}}$
experience only dephasing, whereas spins polarized in any other direction
undergo both dephasing and dissipation. Experimentally, this is reflected in
the extremely long $T_{1}$ lifetimes of $\hat{\textbf{z}}$-polarized states;
during this period, the system still evolves under diffusion and any initial
domain wall starts to melt over time (cf. SI, Fig. S25). Since dissipation
does not affect the spins uniformly, different signals are expected to emerge
depending on the waiting time $t_{\mathrm{wait}}$ after which the system is
finally rotated out of the $\hat{\textbf{z}}$-axis.

To simulate this effect we initialize the system in domain wall states akin to
Fig. S25 and evolve it with $\mathcal{H}_{\pi/2}$ up to times
$t_{\mathrm{wait}}$ after which we switch on dissipation. The results are
displayed in Fig. 11. As the waiting time increases, the zero-crossing time
grows accordingly. This is expected since the system diffuses polarization
away from $r=0$ during $t_{\mathrm{wait}}$ and becomes less sensitive to
dissipation (here acting only at $r=0$). A less obvious observation is the
decrease of the absolute value of the minimally reached total net polarization
$|\mathrm{min}(\mathcal{I}_{x})|$ as a function of the waiting time (Fig. 11
(c)). A possible explanation is that with increasing waiting time, diffusion
causes domain walls to melt, flattening the polarization profiles; hence, the
net amount of polarization left at large distances (i.e., probed at long
times) is reduced with increasing waiting time. This is also supported by
simulations performed with a uniformly polarized (within a small region around
the NV) initial state where no such decay appears as a function of waiting
time. In fact, in this scenario, even a slight increase can be observed.
Intuitively, this is expected as for increased $t_{\mathrm{wait}}$ more
polarization is able to diffuse away from the dissipation-active region around
the NV. All these results are qualitatively consistent with experimental
observations.

The numerical simulations shown in Fig. 3F-G have been performed by means of
the LITE algorithm. During the waiting time $Jt_{\textrm{wait}}$, the system
evolves under the one-dimensional short-range version of
$\mathcal{H}_{\textrm{dd}}$ and is subject to spin diffusion. When the readout
protocol is activated, dissipation comes into play (see above). Therefore, we
initialize the system with a different number of positively polarized spins
and activate dissipation after $Jt_{\textrm{wait}}$. To plot Fig. 3F we
measure the zero-crossing point $t_{\textrm{zc}}$, while for Fig. 3G we
compute the minimum of the total net polarization.

Fig. 11:  Waiting time simulations (a) Time-evolution curves of the total
$\hat{\textbf{x}}$-polarization for different waiting times. After
$t_{\mathrm{wait}}$ the dissipation acting on the site with
$r=r_{\mathrm{NV}}$ is switched on. For increased waiting time, the
corresponding zero-crossing time increases since polarization spreads
diffusively during $t_{\mathrm{wait}}$ diminishing the effect of dissipation
which acts most strongly at the position of the NV ($r_{\mathrm{NV}}=0$). We
use an initial domain wall state with a total of $N_{+}+N_{-}=31$ initially
polarized spins where the $N_{+}=13$ (panel $(i)$, blue) and $N_{+}=11$ (panel
$(ii)$, orange) spins closest to $r_{\mathrm{NV}}$ have $p_{r}=0.6$ and the
remaining $N_{-}=31-N_{+}$ spins have $p_{r}=-0.2$ (here, $\rho_{n}\propto
1+p_{r}I^{x}_{r}$). The dissipation parameters are
$\gamma_{z}=\gamma_{+}=\gamma_{-}=0.5J$, and the other simulation parameters
as in SI Fig. S15. (b) The zero-crossing time $t_{\mathrm{zc}}$ of the curves
in (a) measured from $t_{\mathrm{wait}}$ as a function of $t_{\mathrm{wait}}$.
(c) The absolute value of the minimum attained value of the total
$\hat{\textbf{x}}$-polarization of the curves in (a) as a function of the
waiting time, normalized with respect to $t_{\mathrm{wait}}=0$. In addition to
the domain wall initial states, we also perform a simulation with a uniformly
polarized initial state with $N_{-}=0$ ($N_{+}=31$) (red curve, (iii)).

###  VIII.12 Details on the Theoretical Results for Hamiltonian Engineering

####  VIII.12.1 Energy vs. Polarization diffusion

In the Hamiltonian Engineering simulations, we find that the variance of the
energy distribution grows as ${\sim}\sqrt{t}$ during the late-time dynamics
($Jt{>}10$), as expected for diffusive processes (see SI S8.2.3). Since the
energy operator has a large overlap with the $\hat{\textbf{x}}$-polarization
at long times, the diffusive scaling becomes visible in the spreading of the
polarizing front (dashed tilted line, Fig. 5B); however, we stress that in
Hamiltonian Engineering, proper diffusion of polarization is absent due to a
lack of $\hat{\textbf{x}}$-polarization conservation.

####  VIII.12.2 Effect of Dissipation

Let us emphasize that in contrast to State Engineering, the $t_{\textrm{zc}}$
zero-crossing for Hamiltonian Engineering arises even in the absence of
dissipation (see SI S8.2). As observed in the long-time slices in Fig. 5B,
dissipation only serves to extinguish polarization close to the NV (i.e.,
accelerating the inevitable diffusion-induced polarization inversion).

####  VIII.12.3 Critical radius $r_{c}$

In the main text, we argued that the polarization zero-crossing radius is
given by the zeroes, $\varphi(r_{c})=0$, of the effective on-site potential
$\varphi$. To leading order, this on-site potential is determined solely from
the microscopic gradient field induced by the NV together with the details of
the kick sequence (see SI S7.2). Thus, assuming that we have full knowledge
about the gradient field, we can estimate the crossing radius for the
experimental data, see Fig. 4C-D.

In particular, we assume that the gradient field induced by the NV electron
follows the dipole-dipole form $B=2PK_{\mathrm{exp}}(3\cos^{2}\theta-1)/r^{3}$
with $K_{\mathrm{exp}}=\mu_{0}\hbar\gamma_{n}\gamma_{e}/4\pi$ and electron
spin polarization of $P({\approx}10\,\%)$. Then, using the experimental
parameters from Fig. 4, i.e., $\Omega{\approx}50\,\mathrm{kHz}$,
$t_{\mathrm{acq}}=51\,\mu\mathrm{s}$ and, e.g., $\theta=0.94\,\pi$
($t_{\mathrm{p}}=95\,\mu\mathrm{s}$), we can estimate the crossing radius
$r_{c}{\approx}2.7\,\mathrm{nm}\times\sqrt[3]{|3\cos^{2}\theta-1|}$. For more
details see SI, Sec. S7.4.

Supplementary Information:
Nanoscale engineering and dynamical stabilization of mesoscopic spin textures

Kieren Harkins,1,∗ Christoph Fleckenstein,2,∗ Noella D’Souza,1,∗ Paul M.
Schindler,3,∗ David Marchiori,1,∗ Claudia Artiaco,2

Quentin Reynard-Feytis,1 Ushoshi Basumallick,1 William Beatrez,1 Arjun
Pillai,1 Matthias Hagn,1 Aniruddha Nayak,1

Samantha Breuer,1 Xudong Lv,1 Maxwell McAllister,1 Paul Reshetikhin,1 Emanuel
Druga,1 Marin Bukov,3 and Ashok Ajoy,1,4,5

_1 Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA._

_2 Department of Physics, KTH Royal Institute of Technology, SE-106 91
Stockholm, Sweden._

_3 Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str.
38, 01187 Dresden, Germany._

_4 Chemical Sciences Division, Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA._

_5 CIFAR Azrieli Global Scholars Program, 661 University Ave, Toronto, ON M5G
1M1, Canada._

###### Contents

  1. I Introduction
  2. II Spin texturing via hyperpolarization injection
  3. III Robust spin shells by Hamiltonian engineering
  4. IV Numerical Results
  5. V Discussion & Outlook
  6. VI Acknowledgements
  7. VII Author Contributions
  8. VIII Methods
    1. VIII.1 Sample
    2. VIII.2 Experimental Apparatus
    3. VIII.3 Hyperpolarization Methodology
    4. VIII.4 Data Acquisition and Processing
    5. VIII.5 Estimation of Kick Angle $\vartheta$
    6. VIII.6 Comparison of Signals in State and Hamiltonian Engineering
    7. VIII.7 Additional Results on Hamiltonian Engineering
    8. VIII.8 Effective Hamiltonians for State and Hamiltonian Engineering protocols
    9. VIII.9 Dissipation induced by the NV center
    10. VIII.10 Numerical simulations
      1. VIII.10.1 Quantum simulations
      2. VIII.10.2 Classical simulations
    11. VIII.11 Details on the Theoretical Results for State Engineering
      1. VIII.11.1 Zero crossing radius for State Engineering
      2. VIII.11.2 Waiting time analysis
    12. VIII.12 Details on the Theoretical Results for Hamiltonian Engineering
      1. VIII.12.1 Energy vs. Polarization diffusion
      2. VIII.12.2 Effect of Dissipation
      3. VIII.12.3 Critical radius $r_{c}$
  9. S1 Summary
  10. S2 Dissipation sources in the experiment
    1. S2.1 Nuclear relaxation in transverse and longitudinal directions
    2. S2.2 Indirect evidence for electron dissipation on the ${}^{13}\textrm{C}$ nuclear spins
  11. S3 Prethermal signal decays at a wide range of flip angles
  12. S4 Experiments probing $\hat{\textbf{z}}$ polarization
  13. S5 Hamiltonian engineering: variation of $t_{\textrm{zc}}$ with experimental parameters
    1. S5.1 Experiments with changing hyperpolarization time
    2. S5.2 Experiments with different orientations
    3. S5.3 Temperature dependence of spin texture
  14. S6 Theory model
    1. S6.1 Closed system
    2. S6.2 Open system
  15. S7  Effective Hamiltonian analysis
    1. S7.1  Exact diagonalization simulation
    2. S7.2  Derivation of Floquet Hamiltonians
      1. S7.2.1 Spin locking for state engineering ($\vartheta\napprox\pi$)
      2. S7.2.2 Spin locking for Hamiltonian engineering ($\vartheta\approx\pi$)
    3. S7.3  Analysis based on the Eigenstate Thermalization Hypothesis
    4. S7.4  Estimating the crossing radius $r_{c}$
    5. S7.5  Simplifications
  16. S8 One-dimensional approximate quantum dynamics
    1. S8.1 Introduction to the LITE algorithm
    2. S8.2 Energy diffusion around $\vartheta\approx\pi$: Hamiltonian engineering
      1. S8.2.1 Toy model Hamiltonian
      2. S8.2.2 Initial states
      3. S8.2.3 Energy diffusion
      4. S8.2.4  Constant on-site potential $\varphi$
      5. S8.2.5 Space-dependent on-site potential $\varphi_{n}$
      6. S8.2.6 Polarization gradients
      7. S8.2.7 Energy diffusion in inhomogeneous systems
      8. S8.2.8 Dissipation
    3. S8.3 Polarization diffusion around ($\vartheta\approx\pi/2$): state engineering
      1. S8.3.1 Toy model Hamiltonian
      2. S8.3.2 Initial states
      3. S8.3.3 Polarization diffusion
    4. S8.4 Comparison between energy diffusion and polarization diffusion
    5. S8.5 Summary
  17. S9 Role of dimensionality and interaction range
    1. S9.1 Floquet dynamics beyond the prethermal plateau
    2. S9.2 Equilibration and thermalization dynamics within the prethermal plateau
  18. S10 Three-dimensional classical simulations on a diamond lattice
    1. S10.1  Classical simulation algorithm
    2. S10.2  Spin inversion in classical simulations

##  S1 Summary

Given the length of this Supplementary Information, as an aid to the reader,
we present a summary of the main results here. Sec. S2, summarizes sources of
nuclear relaxation, with a focus on dissipation stemming from the central
electron. Sec. S2.1 discusses the dephasing time ($T_{2}^{\ast}$),
longitudinal relaxation time ($T_{1}$), and transverse relaxation under spin-
lock under Floquet driving ($T_{2}^{\prime}$) in our system. Sec. S2.2
elucidates indirect experimental signatures of NV-driven transverse relaxation
of the ${}^{13}\textrm{C}$ spins.

Section S3 describes spin-lock decays for a wide range of flip angles
$\vartheta\in[0,1.5\pi]$, showing slow featureless decays in all regions away
from $\vartheta{\approx}\pi$, and sharp zero-crossings in the latter due to
the formation of Hamiltonian engineering driven spin textures.

Section S4 elucidates complementary experimental studies focused on the
Hamiltonian engineering method. In Sec. S4, we describe experiments that
confirm that the loss of signal during the zero-crossing is not due to the
spin vector being trivially tilted onto the $\hat{\textbf{z}}$ axis. This was
demonstrated by inserting a $\pi/2$ pulse perpendicular to those in the train
of pulses. If the spin vector were in the $\hat{\textbf{z}}$ axis, this
sequence would immediately create observable polarization.

Section S5 examines changes in $t_{\textrm{zc}}$ to probe the sensitivity of
formed spin texture to experimental parameters. Increasing hyperpolarization
time in Fig. S7 does not significantly change the location of
$t_{\textrm{zc}}$, showcasing the robustness of Hamiltonian Engineering to
different initial states.

Decreasing the temperature in Fig. S9 diminishes the effect of electron
dissipation to the system, preserving the generated spin texture for longer
and moving $t_{\textrm{zc}}$ to later times in the experimental data. In Fig.
S5, changing the effective Floquet drive Rabi frequency, by changing the pulse
duty cycle while fixing $\vartheta$, leads to a movement in $t_{\textrm{zc}}$
as anticipated from theory.

Changing sample orientation creates composite behavior at kick angle $\pi$
qualitatively identical to other Hamiltonian Engineering experiments (compare
Fig. S8Aiii and Biii). The multi-dip behavior seen in Fig. S8Bii is consistent
with domain boundaries forming from NV centers inequivalently aligned with
respect to the applied magnetic field. Lastly, changing the pulse transmit
offset frequency (TOF) shifts the Rabi curve of the system by a finite phase,
yet zero crossing behavior can appear before or after the bulk $\pi$ value for
a given TOF (Fig. S6).

Section S6 details the theory model for dipolar-coupled ${}^{13}\textrm{C}$
nuclear spins and nitrogen-vacancy (NV) center electrons. The system is
described as a closed system with a hierarchical energy scale separation,
allowing for independent manipulation of the ${}^{13}\textrm{C}$ nuclei and NV
electrons. The Hamiltonian captures the couplings between spins and the
external magnetic field. Additionally, the open system nature of the
experimental setup is discussed in Sec. S6.2, where spins and NV centers
couple to a phonon bath, resulting in decoherence and dissipation. An
effective Lindblad master equation is derived, considering the singular
coupling limit, which describes the dissipative dynamics of the spins only.
The dissipative terms in the Lindblad equation account for dephasing and
dissipation effects operational in experiments.

In Sec. S7, we investigate in more detail the dynamics of the
${}^{13}\textrm{C}$ spin system coupled to the NV center. By considering the
fast dynamics of the NV spin, we replace the NV operators with their
expectation values. This leads to an effective on-site potential for the
${}^{13}\textrm{C}$ spins. We analyze the system’s evolution under a kick
sequence and derive the effective Hamiltonian using Floquet’s theorem. We
consider two cases: spin locking with a kick angle different from $\pi$ and
near $\pi$. In the former case, the effective Hamiltonian conserves the
$\hat{\textbf{x}}$-polarization, while it is not conserved in the latter case,
and the strong NV-induced potential leads to a spatially inhomogeneous
effective Hamiltonian. Exact numerical simulations confirm the validity of the
derived effective Hamiltonians. In Sec. S7.3, we discuss the dynamics of the
${}^{13}\textrm{C}$ system within the Eigenstate Thermalization Hypothesis
(ETH). The system is expected to (ultimately) thermalize to an infinite
temperature state due to the absence of energy conservation. However, at high
driving frequencies, a pre-thermal plateau emerges before full thermalization.
We find that in the $\vartheta{=}\pi$ case, the local polarization follows the
local on-site potentials induced by the NV center, which can lead to a sign-
inversion of the local polarization, and hence also of the integrated
polarization. In Sec. S7.4, we estimate the $r_{c}$ domain boundary where the
local polarization changes sign. Finally, in Sec. S7.5, we present two
simplified, effective Hamiltonians for $\vartheta{\neq}\pi$ and
$\vartheta{\approx}\pi$ that aid in a qualitative understanding of the
dynamics.

Section S8 first introduces the approximate time-evolution algorithm used to
perform the quantum dynamics of the simplified one-dimensional toy models used
in the rest of this Section. In particular, Sec. S8.1 summarizes the main
features of the local-information time-evolution (LITE) algorithm. LITE is
designed to investigate the out-of-equilibrium transport of short-range
systems. In contrast to similar algorithms [REF:], it preserves local
constants of motion. LITE decomposes the system into subsystems and solves the
von Neumann equation for each subsystem in parallel. Importantly, it can also
simulate open quantum systems described by the Lindblad master equation. In
Sec. S8.2, we apply LITE to simulate an effective one-dimensional short-range
toy model for the case $\vartheta{=}\pi$. Such a toy model is derived from the
experimental three-dimensional long-range model within the approximation of
sparse density of spins. Thus, the one-dimensional model retains only the
dominant nearest-neighbor mutual coupling and the space-dependent on-site
potential generated by the NV center. For $\vartheta{\approx}\pi$, the energy
diffusion and the behavior of the total $\hat{\textbf{x}}$-polarization are
analyzed, showing a sign inversion of polarization in the presence of a space-
dependent potential. The diffusion of energy in the inhomogeneous systems is
explored, revealing a slowing down of energy spread that can be explained in
terms of a spatially dependent diffusion constant. Dissipation effects are
then introduced. It is shown that dissipation favors energy diffusion,
inducing the emergence of Gaussian-like energy distributions. The presence of
dissipation does not alter the possibility of observing sign inversion in the
total $\hat{\textbf{x}}$-polarization. In Sec. S8.3, the dynamics of the one-
dimensional spin system obtained for a kick angle around $\pi/2$ is
investigated. Unlike the case of a kick angle close to $\pi$, the total (net)
$\hat{\textbf{x}}$-polarization is conserved in addition to energy. This leads
to the diffusion of polarization throughout the system, making it impossible
to detect polarization gradients using global measurements. Only by
introducing dissipation can a sign inversion of the integrated polarization be
induced. The effects are demonstrated using numerical simulations obtained
within the LITE algorithm. In Secs. S8.4 and S8.5, we highlight the main
differences between the two kick angles analyzed and summarize the results.

In Sec. S9, we discuss the role of dimensionality and interaction range in the
simplified one-dimensional versions of the long-range three-dimensional
quantum system used in the experiments. We argue that, although our simplified
models may not capture all the details of the out-of-equilibrium dynamics, we
expect them to provide qualitative insights into the three-dimensional system.
We discuss the effects of dimensionality on Floquet prethermalization and
equilibration dynamics within the prethermal plateau. We also consider the
possible effects of many-body localization and the impact of spin-spin
couplings and angular dependence.

In Sec. S10, we perform classical simulations on a three-dimensional diamond
lattice to complement the previous analysis. The classical simulation allows
us to study larger systems and reach longer time scales, although quantum
correlations are neglected. We analyze the spin locking regime near
$\vartheta\approx\pi$ and observe the formation of a spatially inhomogeneous
local polarization profile. The classical simulation confirms the predictions
from the analytical and numerical results obtained for one-dimensional toy
models. The simulations show that the local polarization at late times follows
the applied local potential up to an overall constant, which is related to the
inverse temperature. The classical simulation results provide further support
for the formation of a robust spatially inhomogeneous local polarization
profile in the experimental system.

##  S2 Dissipation sources in the experiment

###  S2.1 Nuclear relaxation in transverse and longitudinal directions

In this discussion, we elucidate the mechanisms for ${}^{13}\textrm{C}$
nuclear relaxation for different initial states. Let us begin by examining
Fig. S1A, which illustrates the ${}^{13}\textrm{C}$ dephasing when prepared in
a superposition state, ${\propto}I_{x}$. In the absence of the applied Floquet
drive, the measured signal rapidly decays in $T_{2}^{\ast}{\sim}1.5$ms. This
is primarily attributable to interspin interactions, specifically, the
formation of many-body spin states invisible under inductive readout. The
magnitude of the internuclear dipolar coupling can be estimated from this
decay as $\langle J\rangle{=}660$ Hz.

Moving on to Fig. S1B, we consider the spins prepared in the $\rho_{z}\propto
I_{z}$ state. In this case, the corresponding ($T_{1}$) relaxation is found to
be remarkably long, with $T_{1}{<}1$hr even at room temperature. Relaxation
arises from spin-flipping noise perpendicular to $\hat{\textbf{z}}$, stemming
from the spectral density component that matches the nuclear Larmor frequency
($\omega_{L}{=}\gamma_{n}B_{0}$). Phononic contributions likely play a
dominant role in these decay channels; their inherently weak nature results in
the long $T_{1}$. Conversely, to assess the contribution of electronic
dissipation to this relaxation, it is worth noting that the corresponding
spectral density is centered at $(T_{1e})^{-1}$, which is orders of magnitude
separated from $\omega_{L}$. Hence, this contribution can be considered to be
extremely weak. This is supported by experiments in Fig. 3 of the main paper.

Lastly, we direct our attention to $T_{2}^{\prime}$ relaxation, as shown in
Fig. S1C. Here spins are prepared in a superposition state
$\rho_{I}{\propto}I_{x}$ and subjected to Floquet driving, which preserves
them along $\hat{\textbf{x}}$. However, this effective Hamiltonian is only
accurate to leading order in a Magnus expansion, and higher-order terms can
induce heating to infinite temperature. Relaxation additionally arises from
$\hat{\textbf{z}}$-oriented noise spectral matched to the energy gap in the
rotating frame, $\Omega_{\textrm{eff}}{=}\Omega(t_{p}/\tau)$, where the latter
factor is the pulsing duty cycle. In our experiments,
$(T_{1e})^{-1}{\sim}100$kHz and exhibits a significant contribution at
$\Omega_{\textrm{eff}}$. Consequently, the electron can play a role as a
“dissipator” for the spins.

Note that while the lattice contains both NV and P1 center electrons, the NV
center is present at the same location $r{=}0$ for every ${}^{13}\textrm{C}$
ensemble considered while the P1s are randomly positioned. Thus, upon ensemble
averaging, the P1s contribute to only a background rate of $T_{2}^{\prime}$
relaxation, while the NV-driven relaxation is discernible in a spatially
dependent manner.

###  S2.2 Indirect evidence for electron dissipation on the
${}^{13}\textrm{C}$ nuclear spins

In previous work, we presented preliminary evidence supporting the role of the
electron spin as a nanoscale dissipator for neighboring nuclear spins [62].
Here, we provide a concise summary of these findings. Specifically, we
observed that the relaxation profiles of ${}^{13}\textrm{C}$ spins in the
$T_{2}^{\prime}$ regime appear to lengthen when a greater degree of
hyperpolarization (achieved by employing longer $\tau_{+}$ durations) is
injected into the spins. This is shown in Fig. S2, where the colorbar
represents $\tau_{+}$. We interpreted this behavior by considering that, with
longer $\tau_{+}$ durations, the polarization has more time to diffuse within
the lattice, thus reaching nuclei less affected by relaxation originating from
the NV center. Instead, the observation of sharp zero-crossings in Fig. 2 of
the paper provides more direct evidence for this, and allows it to be
exploited for the readout of spin texture.

Fig. S1: Nuclear dephasing and relaxation. (A) Panel illustrates
$T_{2}^{\ast}$ dephasing of ${}^{13}\textrm{C}$ spins during free induction
decay at room temperature. The fast decay observed,
$T_{2}^{\ast}{\approx}1.5$ms, is attributed to internuclear interactions.
Reproduced from [34] (B) Panel presents $T_{1}$ decay at room temperature,
lasting over twenty-five minutes at 7T. Reproduced from [62]. Fig. S2:
Indirect evidence for electron dissipation. Normalized spin lock curves
plotted against increasing polarization time $\tau_{+}$, reproduced from [62].
Data demonstrates that $T_{2}’$ values increase with $\tau_{+}$, indicating
that ${}^{13}\textrm{C}$ nuclei located further away are less affected by
relaxation from the NV, while proximal spins experience faster NV relaxation.
This indirect evidence is further supported by the formation of shells, as
depicted in Fig. 2. Fig. S3: Data for a wide range of spin-locking angles,
$\vartheta{\in}[0,1.5\pi]$. (A) Integrated signal $S_{\textrm{int}}$ over the
first $60\,$s of experiment time Shaded regions show representative regions
centered at $\vartheta{=}0.1\pi$, $\pi$, and $1.3\pi$, respectively. (B-D)
Corresponding decay profiles for the shaded regions in (A) are shown, with
$\vartheta$ values indicated by the colorbar. Prethermal decays are
featureless and follow stretched exponential decays on either side away from
$\vartheta{=}\pi$. Near $\pi$, however, there is a significant change in the
signal profile, associated with the formation of shell-like spin textures.
Full movie corresponding to this dataset is available at Ref. [63]

##  S3 Prethermal signal decays at a wide range of flip angles

We present here an analysis of spin-lock decay behavior for a wide range of
flip angles, $[0,1.5\pi]$, as a comparison to the results in the main text
which focused on the proximity of $\vartheta{=}\pi$ (Fig. 4) and
$\vartheta{=}\pi/2$ (Fig. 2). Experiments here (Fig. S3) were conducted at
115K, similar to conditions in Fig. 4. Fig. S3A shows the integrated signal
$S_{\textrm{int}}$ for the range of angles, and Fig. S3B-D shows
representative slices of signal $S$ at specific values of $\vartheta$ in three
ranges (shaded). A full movie of the dataset is available in the Youtube link
at Ref. [63].

Data in Fig. S3B,D illustrates that far from $\pi$, the spin-lock profiles
‘exhibit slow decays at $T_{2}’$ that follow stretched exponential behavior.
In the regime away from $\pi$, there is no inversion of the spin signal with
changing $\vartheta$. However, as observed in Fig. S3C (also Fig. 4), in the
proximity of $\pi$, there are sharp zero crossings corresponding to the
formation of spin-shell texture. Overall, these results demonstrate unexpected
sign inversions for the case of CPMG pulse trains, which to our knowledge, has
not been reported previously.

Fig. S4: “Interrupted” Floquet driving experiments investigating if zero-
crossing observed in Fig. 4 is caused by polarization tipping to the
$\hat{\textbf{z}}$ axis. (A) i) Schematic of the experiment shown in the
inset, consisting of a spin-locking train interrupted by an orthogonal $\pi/2$
pulse of phase x or y at $t=t_{int}$ that reorganizes population from the
$\hat{\textbf{z}}$ axis to the $\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$
plane. (A-B) i) Signal at $t_{int}{=}t_{zc}$ for case of
$\left.\pi/2\right|_{x}$ in (A) or $\left.\pi/2\right|_{y}$ in (B) interrupt
pulses. The measured $\hat{\textbf{z}}$ amplitude for $\left.\pi/2\right|_{x}$
pulses does not significantly change, indicating that zero-crossing in Fig. 4
does not result from spins tipping out of the
$\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane. Residual
$\hat{\textbf{z}}$-component can be attributed to slight off-resonance in the
pulses. (A-B) ii) Signal decays for the interrupted Floquet driving
experiments where $t_{int}\neq t_{zc}$. In (A), the $\left.\pi/2\right|_{x}$
pulse caused minimal disruption to the signal decays. However, in (B) the
$\left.\pi/2\right|_{y}$ pulse reorganized the population of spins from the
$\hat{\textbf{x}}$ axis to the unobservable $\hat{\textbf{z}}$ axis, causing
an abrupt loss of signal.

##  S4 Experiments probing $\hat{\textbf{z}}$ polarization

In the experiments describing shell formation through Hamiltonian Engineering
(Fig. 4 of the main paper), we observed an apparent change in the polarization
direction from $\hat{\textbf{x}}$ to $-\hat{\textbf{x}}$, which was attributed
to the thermalization of spins under the electron gradient induced potential.
In this section, we aim to provide additional evidence to support this.
Specifically, we conducted measurements probing the $\hat{\textbf{z}}$
component of the spin vector at various time points $t$ along experimental
traces in Fig. 4B. This allows us to examine whether the observed zero-
crossings in Fig. 4 are trivially a result of the spin vector tilting towards
the $\hat{\textbf{z}}$ axis, where the spins become unobservable.

To further investigate this, we performed additional experiments as depicted
in Fig. S4A. The spins are subject to the Floquet drive as before, but it is
interrupted at time $t{=}t_{c}$ with $\pi/2$ pulses applied in either the
$\hat{\textbf{x}}$ direction (Fig. S4A) or the $\hat{\textbf{y}}$ direction
(Fig. S4B). These pulses reorganize the spin populations from the
$\hat{\textbf{y}}{\text{-}}\hat{\textbf{z}}$ or
$\hat{\textbf{x}}-\hat{\textbf{z}}$ planes into the
$\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane at $t_{c}$. Setting
$t_{c}{=}t_{\textrm{zc}}$, we can precisely probe the $\hat{\textbf{z}}$
component at the zero-crossing point. Moreover, by varying $t_{c}$, we can
track the $\hat{\textbf{z}}$ component throughout the observed dynamics in
Fig. S4A-B.

Specifically, Fig. S4A(i) and Fig. S4B(i) display data obtained at
$t_{c}{=}t_{\textrm{zc}}$ slices for both cases. Notably, no significant
$\hat{\textbf{z}}$ signal was observed, with the residual signal primarily
attributed to the frequency offset of the pulses employed in these
experiments. To further support this observation, the inset Fig. S4B(ii)
illustrates the variation of the $\hat{\textbf{z}}$ amplitude as a function of
$t_{c}$. It reveals that the $\hat{\textbf{z}}$ component undergoes negligible
change with $t$. Consequently, Fig. S4 provides compelling evidence that the
zero-crossing observed in Fig. 4 arises from the thermalization dynamics of
the spins, rather than from the polarization vector deviating away from the
$\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane.

Fig. S5: Variation of spin-lock profiles with pulse duty cycle. Spin-lock
profiles for different effective Rabi frequencies (see colorbar) applied to
the nuclear spins and $\vartheta{\sim}\pi$, similar to Fig. 4 of the main
text. $\Omega_{\textrm{eff}}$ is controlled here by varying the pulse duty
cycle, while maintaining a fixed flip angle $\vartheta$. Inset shows pulse
sequence schematic, with definition of the duty cycle being the total time
between the start of one pulse and the next.  Fig. S6: Variation of spin-lock
profiles with frequency offset. (A) Spin-lock profiles at different frequency
offset values for the Floquet drive (see colorbar). (B) Movement of zero-
crossing with offset. The 2D color plot shows offset data from (A) plotted on
a logarithmic scale. Zero crossing appears as an abrupt decrease in signal
(colored blue). There is a wide range of offsets ($\sim$2 kHz) for which a
zero crossing is observed.

##  S5 Hamiltonian engineering: variation of $t_{\textrm{zc}}$ with
experimental parameters

Accurately estimating the domain boundary of the generated spin texture based
solely on the zero-crossing position $t_{\textrm{zc}}$ poses a significant
challenge. $t_{\textrm{zc}}$ relies on the complex interplay between electron
dissipation and spin diffusion, making it inadequate for providing an accurate
measurement of the critical radius $r_{c}$ associated with this boundary.
However, by considering the variations in experimental parameters, it is
possible to anticipate the trends in the movement of $r_{c}$ (and,
consequently, $t_{\textrm{zc}}$), which can then be validated through
experimental investigations.

Consider first the influence of the nuclear Rabi frequency
$\Omega_{\textrm{eff}}$. Since the effective flip angle $\vartheta$ is
determined as an interplay between the electronic gradient field and the
strength of the $\hat{\textbf{x}}$ Rabi field, changing
$\Omega_{\textrm{eff}}$ would directly affect $r_{c}$. We anticipate from
theory that decreasing $\Omega_{\textrm{eff}}$ should result in a leftward
shift in $t_{\textrm{zc}}$. In experiments (see Fig. S5), we manipulate
$\Omega_{\textrm{eff}}$ by modifying the duty cycle of the Floquet drive, as
schematically denoted in the inset of Fig. S5. Experimental data in Fig. S5
then confirms this theoretical prediction.

Similarly, changing pulse frequency offset at fixed $\Omega$ should result in
a shift of the ${}^{13}\textrm{C}$ spin layer that has a flip angle matching
$\pi$, thereby changing $r_{c}$ and $t_{\textrm{zc}}$. In the rotating frame
of the Floquet drive, changing offset refers to applying a global
$\hat{\textbf{z}}$-oriented field to the ${}^{13}\textrm{C}$ nuclei. Since
this field can either enhance or counteract the local field from the NV center
gradient, the precise magnitude of the shift to $r_{c}$ is challenging to
predict. However, qualitative trends in these shifts are easily observable.
For instance, Fig. S6A demonstrates the effect of changing the offset at a
fixed length $t_{p}$ of the applied pulses. A change in $t_{\textrm{zc}}$ is
observed, and for large offsets, the zero-crossing no longer appears in the
70s experiment time. This is consistent with effectively shifting the position
of the ${}^{13}\textrm{C}$ slice corresponding to $r_{c}$. We find that the
$\vartheta{=}\pi$ slice occurs after the parabola-like $t_{\textrm{zc}}$ curve
observed in Fig. 4C for off-resonant pulses. This again points to shifting
$r_{\mathrm{c}}$ via the resonance offset field. Nevertheless, the physics
underlying the formation and thermalization of the spinning shell remains
unchanged in this case.

Fig. S7: Effect of hyperpolarization injection $\tau_{+}$ on spin-shell
texture. (A) Schematic of the experiment sequence, where injection time
$\tau_{+}$ is varied. After shuttling time $t_{s}{=}90$s the sample is
subjected to a spin-locking train with $\vartheta{\sim}\pi$ for readout,
identical to main text Fig. 4A. (B) Normalized signal decays for different
values of $\tau_{+}$, as indicated by the colorbar. (C) Analysis of the zero-
crossing times $t_{\textrm{zc}}$ with $\tau_{+}$. Solid line is a spline guide
to the eye.  Fig. S8: Spin-texture for two different sample orientations.
(A-B)i EPR spectra. Indirectly obtained EPR spectrum via ${}^{13}\textrm{C}$
DNP, produced by sweeping MWs in narrow 25MHz windows. The integrated signal
of resulting NMR decay curves is plotted as a function of the center frequency
of the swept MWs. Solid line is a spline guide to the eye. (A-B)ii
Representative signal decays in the proximity of $\vartheta{=}\pi$ showing
zero crossings. (A-B)iii Color plots showing full data similar to Fig. 4E of
the main paper, highlighting the movement of the zero crossings with
$\vartheta$ for the two orientations. Panels (A-B)iii indicate physics
underlying Hamiltonian Engineering (Fig. 4) is qualitatively independent of
sample orientation. Fig. S9: Variation of spin-shell profiles with
temperature. (A) Representative spin-lock decays with the temperature (see
colorbar). Experiments are performed at a fixed flip angle
($\vartheta{\sim}\pi$), and identical Rabi frequency. (C) Variation of zero-
crossing time $t_{\textrm{zc}}$ with temperature.

###  S5.1 Experiments with changing hyperpolarization time

We examine the relationship between results displayed in Fig. 4 of the main
paper and the amount of hyperpolarization injected into the
${}^{13}\textrm{C}$ nuclear spins. Specifically Fig. S7 describes experiments
varying the amount of time $\tau_{+}$ for which hyperpolarization is injected.
Temperature is maintained consistent at 115K, flip angle here is 160∘, and
shuttling time $t_{s}{=}90\,$s.

The amplitude signal obtained in these experiments (Fig. S7B) shows the
expected zero crossing and subsequent phase inversion. The zero crossing time
$t_{\textrm{zc}}$ does not significantly change with increased
hyperpolarization time $\tau_{+}$. This stands in contrast to the movement of
the zero crossing with variations in the flip angle (Fig. 4) or temperature
(Fig. S9). These results support the fact that spin texture engineered via
Hamiltonian Engineering is qualitatively independent of the initial state of
nuclear spins.

###  S5.2 Experiments with different orientations

Here we examine the effect of changing the diamond sample orientation to
results described in Fig. 4 of the main paper. We compare two different sample
orientations, where the NV center families have different directions relative
to $B_{0}$. Fig. S8A-B shows the DNP-measured EPR spectra of the NV centers in
both cases. For experiments similar to Fig. 4B, we find qualitatively
identical development of spin-shell texture and zero-crossings in both cases,
as described in the colorplots in Fig. S8Aiii,Biii. The physics underlying
Fig. 4 can therefore be considered independent of the exact sample
orientation.

###  S5.3 Temperature dependence of spin texture

We now consider the effect of changing temperature on the observed spin
texture signal via Hamiltonian Engineering (Fig. 4 of the main paper). A key
advance enabled by the instrumentation (Fig. 6) introduced in this paper is
the ability to study the role of temperature (in 77K-RT range), which
previously has been shown to have a strong effect on rates of electronic
relaxation $T_{1e}$ (of both NV and P1 centers).

Fig. S9 presents experiments conducted at six representative temperatures. To
ensure consistency, we maintain the same nuclear Rabi frequency for all
experiments. Results indicate that the time at which the zero-crossing occurs
($t_{\textrm{zc}}$) increases as the temperature decreases. This can be seen
as a direct reflection of the increase in the NV center $T_{1e}$, which in
turn decreases the strength of the dissipation acting on the
${}^{13}\textrm{C}$ nuclear spins (see also Sec. S6.2). We note that Ref.
[64, 65] found that for samples with comparable NV center densities as the one
employed here, $T_{1e}$ increases sharply close to 100K and can exceed
$T_{1e}{>}1$s under these conditions.

##  S6 Theory model

###  S6.1 Closed system

As outlined in the main text, the system of interest is given by dipolar
coupled ${}^{13}\textrm{C}$ nuclear spins with the (rotating-frame)
Hamiltonian

$\mathcal{H}_{\textrm{dd}}=\sum_{k<\ell}b_{k\ell}\left(3I_{k}^{z}I_{\ell}^{z}-\bm{I}_{k}\cdot\bm{I}_{\ell}\right),$
(S 1)

where $\bm{I}_{k}=\left(I_{k}^{x},I_{k}^{y},I_{k}^{z}\right)^{T}$, the
spin-$1/2$ operators $I_{k}^{\eta}=\sigma_{k}^{\eta}/2$ describes the nuclear
spins, and $\sigma^{\eta}$ ($\eta\in\\{x,y,z\\}$) are the Pauli matrices. The
couplings between different ${}^{13}\textrm{C}$ spins strongly depends on
their distance and the angle with the $\hat{\textbf{z}}$-axis via
$b_{k\ell}=J_{\mathrm{exp}}(3\cos^{2}(\theta_{k\ell})-1)/{r}_{k\ell}^{3}$
where
$\cos(\theta_{k\ell}){=}\bm{B}\cdot\mathbf{r}_{k\ell}/(|\bm{B}||\mathbf{r}_{k\ell}|)$;
$\mathbf{r}_{k\ell}$ is the inter-spin vector and $\bm{B}=(0,0,B)^{T}$ is the
externally applied magnetic field and
$J_{\mathrm{exp}}{=}\mu_{0}\hbar\gamma_{n}^{2}/4\pi$, with the gyromagnetic
ratio of the neutron $\gamma_{n}$ and the vacuum permeability $\mu_{0}$. The
${}^{13}\textrm{C}$-spins relevant for the experiment are randomly placed on a
diamond lattice in the proximity of NV-centers. The NV-centers are
approximately described by an electronic two-level system, which is itself
coupled to all ${}^{13}\textrm{C}$ nuclear spins via dipole-dipole
interactions:

$\displaystyle\mathcal{H}_{\mathrm{NV}}$ $\displaystyle=$
$\displaystyle\epsilon_{\mathrm{NV}}S^{z}\,,$ (S 2)
$\displaystyle\mathcal{H}_{\text{NV-}^{\mathrm{13}}C}$ $\displaystyle=$
$\displaystyle\sum_{j}K_{j}\left(3S^{z}I_{j}^{z}-\bm{S}\cdot\bm{I}_{j}\right).$
(S 3)

Here $S^{\eta}=\sigma^{\eta}/2$ with $\alpha\in\\{x,y,z\\}$ are
(pseudo-)spin-$1/2$ operators describing the two-level system of the NV-
center. $K_{j}=K_{\mathrm{exp}}(3\cos^{2}(\theta_{j})-1)/{r}_{j}^{3}$ with the
vector $\mathbf{r}_{j}$ between the NV center and the $j$-th
${}^{13}\textrm{C}$ nuclear spin and the corresponding angle
$\cos(\theta_{j})=\bm{B}\cdot\mathbf{r}_{j}/(|\bm{B}||\mathbf{r}_{j}|)$
between $\mathbf{r}_{j}$ and the external magnetic field $\bm{B}$, and
$K_{\mathrm{exp}}{=}\mu_{0}\hbar\gamma_{n}\gamma_{e}/4\pi$, where $\gamma_{e}$
is the gyromagnetic ratio of the electron.

Note that the system obeys the hierarchy of energy scales
${\epsilon_{\mathrm{NV}}\gg K\gg J}$, where $J$ and $K$ are characteristic
energy scales associated with the Hamiltonians of Eq. (S 1) and (S 3). The
separation of energy scales is important as it allows to independently kick
(i.e., apply strong short pulses to) the ${}^{13}\textrm{C}$ spins while
leaving the NV two-level system unaffected. Moreover, it allows us to consider
the singular-coupling limit in order to derive an effective Lindblad master
equation for the ${}^{13}\textrm{C}$ spins coupled to the NV center (see
below). The kick Hamiltonian of interest, created from an external microwave
field at resonant Larmor-frequency of the external magnetic field, is given in
the rotating frame by

$\mathcal{H}_{x}(t)=\begin{cases}\Omega\sum_{j}I_{j}^{x}&n\tau<t<n\tau+t_{\mathrm{p}}\,,\hskip
5.69046ptn\in\mathbb{N}\\\ 0&\mathrm{otherwise}\end{cases}\,,$ (S 4)

with the driving period $\tau$ and pulse width duration $t_{\mathrm{p}}$.

For the combined system of NV and ${}^{13}\textrm{C}$ nuclear spins, it is
useful to perform an interaction-picture transformation with
$\mathcal{H}_{\mathrm{dd}}+\mathcal{H}_{\text{NV-}^{\mathrm{13}}\text{C}}$ as
the interaction Hamiltonian. Using the rotating-wave approximation, we readily
obtain

$\displaystyle\widetilde{\mathcal{H}}_{\mathrm{tot}}=\mathcal{H}_{\mathrm{dd}}+\sum_{j}2K_{j}S^{z}I_{j}^{z}\,,$
(S 5)

in the rotating frame with respect to $\exp(-it\mathcal{H}_{\mathrm{NV}})$.

###  S6.2 Open system

In the experimental platform, both the NV center and the ${}^{13}\textrm{C}$
spins are not completely isolated. They are coupled to a bath of phonons
generated by the diamond lattice; such a bath causes decoherence and
dissipation. Thus, a comprehensive theory of the system should take these
effects into consideration.

Following the rotating-wave approximation, Eq. (S 5), the microscopic
Hamiltonian of the total system composed of the NV center and the
${}^{13}\textrm{C}$ spins reads

$\displaystyle\mathcal{H}_{\mathrm{tot}}$
$\displaystyle=\mathcal{H}_{\mathrm{dd}}+\mathcal{H}_{\mathrm{NV}}+\sum_{j}2K_{j}S^{z}I_{j}^{z}$
(S 6)
$\displaystyle=\sum_{j<k}I_{jk}\left(3I_{j}^{z}J_{k}^{z}-\bm{I}_{j}\cdot\bm{I}_{k}\right)+\epsilon_{\mathrm{NV}}S^{z}+\sum_{j}2K_{j}S^{z}I_{j}^{z},$
(S 7)

in the lab frame. The first term represents the Hamiltonian of the
${}^{13}\textrm{C}$ spins, the second is the Hamiltonian of the NV center, and
the last is the effective interaction between the two. As mentioned, the total
system obeys the following hierarchy of energy scales:

$\epsilon_{\mathrm{NV}}\gg K\gg J,$ (S 8)

which implies that the timescale intrinsic to the NV center is much shorter
than any other timescale of the problem. This allows us to assume that the NV
center is strongly coupled to the thermal phonon bath while the
${}^{13}\textrm{C}$ spins are only weakly coupled to it. Therefore, we also
assume that the NV center is always in thermal equilibrium with the phonons.
For simplicity, we consider the NV at infinite temperature:
$\rho_{\mathrm{NV}}(t)=\frac{1}{2}\mathds{1}_{\mathrm{NV}},$ for all times
$t$; however, our results can be easily generalized to any temperature.

The thermal-equilibrium assumption for the NV allows us to derive an effective
quantum master equation for the reduced system of the ${}^{13}\textrm{C}$
spins by tracing out the NV center (and, indirectly, the phonons to which it
is coupled) in the form of a Lindblad master equation. The latter is the most
general Markovian master equation. It describes the equation of motion of a
reduced system in contact with a thermal bath which is memoryless, i.e., whose
timescale is much shorter than any other timescale in the problem, as the NV
center in our system. Moreover, the hierarchy of energy scales (S 8) indicates
the validity of the so-called singular coupling limit, that enormously
simplifies the microscopic derivation of the Lindblad equation [66, 67, 68].
Notice that previous works [61] have investigated the opposite scenario in
which one is interested in integrating out the nuclear spins to derive a
master equation for the NV electron only. From Eq. (S 8), it is clear that
this opposite scenario generates a highly non-Markovian master equation for
the NV electron.

The derivation of the Lindblad equation for the ${}^{13}\textrm{C}$ spins
proceeds as follows. Let us rewrite the total system Hamiltonian (S 6) in the
more general form as 333By considering a generic system Hamiltonian
$H_{\mathrm{{}^{13}C}}$ we emphasize that the Lindblad master equation derived
in the singular coupling limit does not depend on it, as will be clear in the
following.

$\mathcal{H}_{\mathrm{tot}}=\mathcal{H}_{\mathrm{{}^{13}C}}\otimes\mathds{1}_{\mathrm{NV}}+\alpha^{-2}\mathds{1}_{\mathrm{{}^{13}C}}\otimes\mathcal{H}^{\prime}_{\mathrm{NV}}+\alpha^{-1}V^{\prime},$
(S 9)

where we have renormalized the bath and interaction Hamiltonian as
$\mathcal{H}^{\prime}_{\mathrm{NV}}=\alpha^{2}\mathcal{H}_{\mathrm{NV}}$ and
$V^{\prime}=\sum_{j}2K^{\prime}_{j}S^{z}I_{j}^{z}$ with $K^{\prime}_{j}=\alpha
K_{j}$. The parameter $\alpha$ is the inverse of the coupling strength; in the
singular coupling limit we eventually take $\alpha\to 0$. Equation (S 9)
highlights the hierarchy of energy scales of our model (S 6).

We start by considering the Nakajima-Zwanzig equation [66, 68]

$\displaystyle\frac{d}{dt}\mathcal{P}\tilde{\rho}(t)=\alpha^{-1}\mathcal{P}\mathcal{V}(t)\mathcal{P}\tilde{\rho}(t)+\alpha^{-1}\mathcal{P}\mathcal{V}(t)\mathcal{G}(t,0)\mathcal{Q}\tilde{\rho}(0)$
$\displaystyle+\alpha^{-2}\int_{0}^{t}du\mathcal{P}\mathcal{V}(t)\mathcal{G}(t,u)\mathcal{Q}\mathcal{V}(u)\mathcal{P}\tilde{\rho}(u).$
(S 10)

where $\tilde{\cdot}$ denotes the time-evolved operator in the interaction
picture; $\mathcal{P}$ and $\mathcal{Q}$ are two orthogonal projector
operators (i.e., $\mathcal{P}^{2}=\mathcal{P}$, $\mathcal{Q}^{2}=\mathcal{Q}$,
and $\mathcal{PQ}=\mathcal{QP}=0$) given by
$\mathcal{P}\rho=\mathrm{Tr}_{\mathrm{NV}}(\rho)\otimes\rho_{\mathrm{NV}}$ and
$\mathcal{Q}\rho=(\mathds{1}-\mathcal{P})\rho$;
$\mathcal{V}(t)\,(\cdot)=-i[\tilde{V}^{\prime}(t),\cdot]$;
$\mathcal{G}(t,u)=\mathcal{T}e^{\int_{u}^{t}dt^{\prime}\mathcal{Q}\mathcal{V}(t^{\prime})}$
is the propagator with $\mathcal{T}$ the time-ordering operator.

The integro-differential equation (S 10) is exact and describes the equation
of motion of the relevant subspace $\mathcal{P}\tilde{\rho}(t)$ in the
interaction picture. Unfortunately, it is usually as difficult to solve as the
von Neumann equation describing the dynamics of the total system. Hence, we
need to make a few more assumptions to proceed further: First, we assume that

$\mathcal{P}\mathcal{V}(t)\mathcal{P}=0,$ (S 11)

which corresponds to
$\mathrm{Tr}_{\mathrm{NV}}\bigl{(}\tilde{V}^{\prime}(t)\rho_{\mathrm{NV}}\bigr{)}=0$.
The Hamiltonian (S 6) fulfils this condition since
$\mathrm{Tr}_{\mathrm{NV}}\bigl{(}S^{z}\rho_{\mathrm{NV}}\bigr{)}=\mathrm{Tr}_{\mathrm{NV}}\bigl{(}S^{z}/2\bigr{)}=0$.
Nonetheless, if the condition is not fulfilled one can shift the system
Hamiltonian $\mathcal{H}_{\mathrm{{}^{13}C}}$ such that it is satisfied [68].
Second, we assume that

$\rho(0)=\rho_{\mathrm{{}^{13}C}}(0)\otimes\rho_{\mathrm{NV}}(0),$ (S 12)

which implies that we have control over the system of the ${}^{13}\textrm{C}$
spins (for example, we can prepare them in a pure state at $t=0$). This
assumption is known as the Born approximation and it is necessary to obtain a
universal dynamical map [68].

Assumptions (S 11) and (S 12) make the first and second term on the right-hand
side of Eq. (S 10) vanish; the latter now becomes

$\frac{d}{dt}\mathcal{P}\tilde{\rho}(t)=\alpha^{-2}\int_{0}^{t}du\mathcal{K}(t,u)\mathcal{P}\tilde{\rho}(u),$
(S 13)

where
$\mathcal{K}(t,u)=\mathcal{P}\mathcal{V}(t)\mathcal{G}(t,u)\mathcal{Q}\mathcal{V}(u)$
is the memory kernel. Since the NV is always in thermal equilibrium, the
memory kernel becomes homogeneous: $\mathcal{K}(t,u)=\mathcal{K}(t-u)$.

Equation (S 13) is in general non-Markovian as the state of the system at time
$t$ depends on all the states at former times from 0 to $t$. By definition, a
Markovian master equation for the system is obtained if the kernel
$\mathcal{K}(t-u)$ behaves as a delta function with respect to
$\tilde{\rho}(u)$. This is verified in the singular coupling limit, where the
hierarchy of energy scales (S 8) implies that the typical timescale at which
the kernel vanish, $\tau_{\mathrm{NV}}$, is much smaller than the typical
variation timescale of the system, $\tau_{\mathrm{{}^{13}C}}$; thus,
$\tau_{\mathrm{NV}}/\tau_{\mathrm{{}^{13}C}}\to 0$ and
$\mathcal{K}(t-u)\propto\delta(t-u)$. This can be seen by integrating Eq. (S
13) and taking the zeroth-order expansion in $1/\alpha$ of $\mathcal{K}(t-u)$
444It can be proven that higher order terms vanish (see Refs. [107, 108,
109]).; we obtain

$\displaystyle\mathcal{P}\tilde{\rho}(t)=$
$\displaystyle\mathcal{P}\tilde{\rho}(0)$ (S 14)
$\displaystyle+\alpha^{-2}\int_{0}^{t}ds\int_{0}^{s}du\mathcal{P}\mathcal{V}(s)\mathcal{V}(u)\mathcal{P}\tilde{\rho}(u)+\mathcal{O}(\alpha^{-3}),$

which is equivalent to

$\tilde{\rho}_{\mathrm{{}^{13}C}}(t)=\tilde{\rho}_{\mathrm{{}^{13}C}}(0)\\\
+\alpha^{-2}\sum_{nm}\int_{0}^{t}ds\int_{0}^{s}du\frac{4K^{\prime}_{n}K^{\prime}_{m}}{r^{3}_{n}r^{3}_{m}}\biggl{(}C_{ij}(s-u)\bigl{[}\tilde{I}_{j}^{z}(u)\tilde{\rho}_{\mathrm{{}^{13}C}}(u),\tilde{I}^{z}_{i}(s)\bigr{]}\\\
+C^{*}_{nm}(s-u)\bigl{[}\tilde{I}^{z}_{n}(s),\tilde{\rho}_{\mathrm{{}^{13}C}}(u)\tilde{I}_{m}^{z}(u)\bigr{]}\biggr{)}+\mathcal{O}(\alpha^{-3}),$
(S 15)

where the correlation functions are defined as

$C_{nm}(s-u)=\mathrm{Tr}\bigl{(}\tilde{S}^{z}(s-u)S^{z}\rho_{\mathrm{NV}}\bigr{)}.$
(S 16)

By Fourier-decomposing the NV operator $\tilde{S}^{z}(s-u)$ and taking into
account the factor $\alpha^{-2}$ in the free NV evolution, we get

$C_{nm}(s-u)=\int_{-a}^{a}d\omega\,e^{-\frac{i\omega(s-u)}{\alpha^{2}}}\mathrm{Tr}\bigl{(}S^{z}(\omega)S^{z}\rho_{\mathrm{NV}}\bigr{)},$
(S 17)

where the integration limits $a$ denotes the spectral support. We see that,
when $\alpha\to 0$, the integral in Eq. (S 17) tends to a delta function:
$C_{nm}(s-u)\propto\delta(s-u)$. Thus, by making the change of variable
$w=\alpha^{-2}u$, and taking the limit $\alpha\to 0$, Eq. (S 15) becomes a
Markovian master equation for the ${}^{13}\textrm{C}$-spins; in the
Schrödinger picture, it reads

$\frac{d}{dt}\rho_{\mathrm{{}^{13}C}}(t)=-i\bigl{[}\mathcal{H}_{\mathrm{{}^{13}C}}+\mathcal{H}_{\mathrm{LS}},\rho_{\mathrm{{}^{13}C}}(t)\bigr{]}\\\
+\sum_{nm}\frac{4K^{\prime}_{n}K^{\prime}_{m}}{r^{3}_{n}r^{3}_{m}}\gamma_{nm}\biggl{(}I^{z}_{m}\rho_{\mathrm{{}^{13}C}}(t)I^{z}_{n}-\frac{1}{2}\\{I^{z}_{n}I^{z}_{m},\rho_{\mathrm{{}^{13}C}}(t)\\}\biggr{)},$
(S 18)

where the Lamb shift Hamiltonian is

$\mathcal{H}_{\mathrm{LS}}=\sum_{nm}L_{nm}I^{z}_{n}I^{z}_{m}.$ (S 19)

Here $\gamma_{nm}$ and $L_{nm}$ represents the Hermitian and anti-Hermitian
components of the integral over time of the bath correlation functions (S 17):

$\int_{0}^{\infty}ds\,C_{nm}(s)=\frac{\gamma_{nm}}{2}+iL_{nm}$ (S 20)

Hence, by inserting Eq. (S 17), we have $\gamma_{nm}\propto 1$ and $L_{nm}=0$.
By absorbing constants of order unity into $K^{\prime}_{i}K^{\prime}_{j}$, we
arrive at the final expression for the ${}^{13}\textrm{C}$ master equation:

$\frac{d}{dt}\rho_{\mathrm{{}^{13}C}}(t)=-i\bigl{[}\mathcal{H}_{\mathrm{{}^{13}C}},\rho_{\mathrm{{}^{13}C}}(t)\bigr{]}\\\
+\sum_{nm}\frac{K^{\prime}_{n}K^{\prime}_{m}}{r^{3}_{n}r^{3}_{m}}\biggl{(}I^{z}_{m}\rho_{\mathrm{{}^{13}C}}(t)I^{z}_{n}-\frac{1}{2}\\{I^{z}_{n}I^{z}_{m},\rho_{\mathrm{{}^{13}C}}(t)\\}\biggr{)}.$
(S 21)

A few comments are in order. The hierarchy of energy scales (S 8) and, thus,
the validity of the singular coupling limit approximation for this model
greatly simplify the derivation of the Lindblad master equation. In
particular, the Lindbladian term on the second line of Eq. (S 21) is
determined only by the bath Hamiltonian $\mathcal{H}_{\mathrm{NV}}$ and the
interaction Hamiltonian $\alpha V^{\prime}$ of Eq. (S 9), while it is
completely agnostic of the system Hamiltonian $\mathcal{H}_{{}^{13}C}$. This
is particularly useful if the system can be governed by different unitary
dynamics $\mathcal{H}_{{}^{13}C}$, e.g., in the prethermal regimes determined
by different kick angles (see Secs. S8.2 and S8.3). Moreover, we note that the
jump operators in Eq. (S 21) act as dephasing operators if the spins are
polarized in the $\hat{\textbf{z}}$ direction, while they generate both
dephasing and dissipation for a system that is polarized in the
$\hat{\textbf{x}}$ direction. This result is in agreement with the
experimental observations. In particular, it is at the origin of the non-
decaying behavior of the normalized signal $S$ as a function of
$t_{\textrm{wait}}$ for $\tau_{-}{=}0$ observed in Fig. 3E of the main text.
Indeed, when $\tau_{-}{=}0$ the whole system is positively polarized along the
$\hat{\textbf{z}}$ axis. Therefore, during $t_{\textrm{wait}}$ no diffusion
nor dissipation affect the system, giving a uniform behavior for any
$t_{\textrm{wait}}$. Finally, we note that the strength of the dissipative
terms $K^{\prime}_{n}K^{\prime}_{m}/(r^{3}_{n}r^{3}_{m})$ decays faster with
the distance from the NV center with respect to the dipole interaction. This
allows us to consider the simplified dissipative one-dimensional toy model
studied in Sec. S8.

##  S7  Effective Hamiltonian analysis

Fig. S10:  Effective Hamiltonian Summary. ${(\mathrm{A})}$ Energy scale
$E_{\mathrm{NV}}(r)$ of magnetic field induced by the NV spin as a function of
the distance from the NV, different colors correspond to different
polarization of the NV spin $P=0-1$. Gray regime corresponds to the
experimentally inaccessible frozen core. ${(\mathrm{B})}$ Spin locking at
$\vartheta{\napprox}\pi$: Effective Hamiltonian at small tilt of the
components of the conserved axis $\widetilde{I}$, Eq. (S 39), for spin locking
at $\vartheta{\neq}\pi$ as a function of distance from NV. ${(\mathrm{C})}$
spin locking at $\vartheta{\approx}\pi$: ${(\mathrm{i})}$ Effective on-site
potential $\varphi_{\mathrm{eff}}(r)$, Eq. (S 29), as a function of the
distance $r$ from NV and pulse width $t_{\mathrm{p}}$, with ${(\mathrm{ii})}$
Angle and radial distribution of effective on-site potential
$\varphi(r,\,\theta)$ for $t_{\mathrm{p}}=49.5\,\mu\mathrm{s}$. Zero crossing
radius $r_{c}$ is given by $r_{c}(\theta)\approx
3.4\,\mathrm{nm}\times\sqrt[3]{|3\cos(\theta)^{2}-1|}$. Summary.– In the case
of the spin locking sequence at $\vartheta{\neq}\pi$, the spatially
inhomogenous on-site potential induced by the NV spin only leads to a tilt in
the spin locking axis but no qualitative change in the physics compared to the
case of absent on-site potential. However, for spin locking at
$\vartheta{\approx}\pi$ the strong NV-induced on-site potential also leads to
spatially inhomogenous on-site potential in the effective Hamiltonian with a
possible sign change in space. Parameters.– We choose Rabi frequency
$\Omega=10\,\mathrm{kHz}$, interpulse time
$t_{\mathrm{acq}}=40\,\mu\mathrm{s}$ and for (B), (C) NV electron polarization
$P=0.1$.

Since the NV system experiences fast dynamics compared to the
${}^{13}\textrm{C}$-spins we neglect the NV-spin degree of freedom and replace
it by its mean value. Hence, we replace the operator $S^{z}$ in Eq. (S 5) by
its thermal expectation value $P=P(T)$, which depends on the temperature $T$,
thus, leading to an effective on-site potential $\eta_{j}=2K_{j}P$ for the
${}^{13}\textrm{C}$ nuclear spins. Thereby, we also neglect the effect of
dissipation.

The kick sequence (see main text Fig. 1B) at stroboscopic times
$n\tau{\approx}n(t_{\mathrm{acq}}+t_{\mathrm{p}})$, $n\in\mathds{N}$, can then
be written as

$U=e^{-i\mathcal{H}t_{\mathrm{acq}}}e^{-i\mathcal{H}_{\mathrm{SL}}t_{\mathrm{p}}}\,,$
(S 22)

with

$\displaystyle\mathcal{H}$
$\displaystyle=\mathcal{H}_{\mathrm{dd}}+\mathcal{H}_{\mathrm{pot}}\,,$ (S 23)
$\displaystyle\mathcal{H}_{\mathrm{SL}}$
$\displaystyle=\mathcal{H}_{x}+\mathcal{H}_{\mathrm{pot}}\,\,,$
$\displaystyle\mathcal{H}_{\mathrm{pot}}$
$\displaystyle=\sum_{j}\eta_{j}I_{j}^{z}\,.$

In Eq. (S 22) we additionally neglect the dipole-dipole coupling
$\mathcal{H}_{\mathrm{dd}}$ during the kicks, since its magnitude is much
smaller compared to $\mathcal{H}_{x}$, i.e.,
$\left|\left|{\mathcal{H}_{\mathrm{dd}}}\right|\right|{\ll}\left|\left|{\mathcal{H}_{x}}\right|\right|$.
However, we consider the effective on-site potential during the kicks, since
it can be of comparable strength at short distances from the NV, see also Fig.
S10 (A) for a schematic of the on-site potential and other energy scales.
Spins within a radius of $\leq 1.7\,\mathrm{nm}$, indicated as a shaded area
in Fig. S10 (A), are far detuned such that they are not contributing to the
bulk dynamics and also cannot be measured by NMR experiments.

The kicked evolution in Eq. (S 22) corresponds to a periodic (Floquet)
evolution with period $\tau$. Therefore, we can exploit Floquet’s theorem,
which states that stroboscopic dynamics ($t=n\tau$, $n\in\mathbb{N}$) are
described by an effective Hamiltonian

$U(n\tau)=U(\tau)^{n}=U_{F}^{n}$ (S 24)

with

$U_{F}=e^{-i\mathcal{H}t_{\mathrm{acq}}}e^{-i\mathcal{H}_{\mathrm{SL}}t_{\mathrm{p}}}\equiv
e^{-i\mathcal{H}_{F}\tau}\,.$ (S 25)

Although the exact effective Hamiltonian $\mathcal{H}_{F}$ is, in general, a
highly non-local object, analytical progress can be made by considering a
high-frequency (small $\tau$) expansion [71] up to some order $O(\tau^{n})$,

$\mathcal{H}_{F}=\sum_{n=0}^{\infty}\tau^{n}\mathcal{H}_{F}^{[n]}=\mathcal{H}_{F}^{(n)}+O(\tau^{n+1})\,,$
(S 26)

valid up to some finite number of drive cycles. For sufficiently small period
$\tau$ we expect that this Hamiltonian is a good description of the prethermal
plateau [71, 7, 8, 24, 34, 72, 25, 26] observed in the experiment. Any
deviations from the effective Hamiltonian are expected to result in Floquet
heating to a featureless infinite temperature state at long times [40, 73].

Before proceeding further, we outline the rest of this section and summarize
the main results. We start with a description in Sec. S7.1 of the numerical
algorithm used to support the analytical findings obtained in the rest of this
section. Then, in Sec. S7.2 we will derive the effective Hamiltonian
$\mathcal{H}^{(0)}_{F}$ to lowest order in the driving period $\tau$. The
derived effective Hamiltonians serve as the starting point of all other types
of numerical and theoretical analysis.

Using a toggling frame expansion, we derive two effective Hamiltonians
capturing the Floquet evolution at stroboscopic times for the two regimes of a
kick-angle $\vartheta\neq\pi$ and near $\vartheta\approx\pi$ in the presence
of the hyperfine splitting field $\eta_{j}$. In agreement with earlier
findings [74, 75], for $\vartheta\neq\pi$ we find that the effective
Hamiltonian conserves the $\hat{\textbf{x}}$-polarization, while the hyperfine
splitting $\eta_{j}$ only leads to a quantitative change in axis
$\hat{\textbf{x}}\to\tilde{\hat{\textbf{x}}}$. In the case,
$\vartheta\approx\pi$, this conservation law is broken and the hyperfine
splitting induces an effective spatially inhomogenous potential dominating the
relaxation dynamics of the polarization.

In Sec. S7.3 we analyze the long-time dynamics generated by the derived
effective Hamiltonians. We show that, for finite systems, the long-time
expectation values are well described by ETH-like arguments. In particular, we
find that, in the $\vartheta=\pi$ case, the local polarization follows the
local on-site potential induced by the electron which can lead to a sign-
inversion of the local polarization and hence – also of the integrated
polarization. We close this section by introducing simplified toy-model
effective Hamiltonians; they capture the main physics but, due to their
reduced complexity, make the qualitative analysis easier to understand.

###  S7.1  Exact diagonalization simulation

Throughout this section, we support our findings by numerically exact full
Floquet simulations. However, these simulations are limited to small system
sizes of $L{=}16$ spins. As we anticipate to find that the distance from the
NV plays a crucial role, we will consider in the following a one-dimensional
geometry to maximize the radial extension. Concretely, in this section we
consider spins, labeled by $j=1,\,\dots,\,L$, to have positions on a one-
dimensional chain described by $x_{j}=x_{0}+aj+\delta x_{j}$, with lattice
spacing $a=1$, offset distance $x_{0}=5$, some independent identically
distributed random displacements $\delta x_{j}$, drawn from a normal
distribution $\mathcal{N}(\mu{=}0,\,\sigma{=}0.01)$ with zero mean and
standard deviation $\sigma=0.01\,a$. However, if not explicitly stated
otherwise, we will consider the full Floquet evolution described by Eq. (S
23), i.e. the spins interact via long-range interactions which decay as
$1/r^{3}$ in the spin-spin distance $r$. In addition we choose
${\Omega/J=10}$, ${\delta\omega/J=0}$, ${Jt_{\mathrm{acq}}=0.1}$ and
$t_{\mathrm{p}}$ is chosen such that $\Omega t_{\mathrm{p}}=\vartheta$