Mapping a 50-spin-qubit network through correlated sensing

Spins associated to optically accessible solid-state defects have emerged as a versatile platform for exploring quantum simulation, quantum sensing and quantum communication. Pioneering experiments have shown the sensing, imaging, and control of multiple nuclear spins surrounding a single electron spin defect. However, the accessible size of these spin networks has been constrained by the spectral resolution of current methods. Here, we map a network of 50 coupled spins through high-resolution correlated sensing schemes, using a single nitrogen-vacancy center in diamond. We develop concatenated double-resonance sequences that identify spin-chains through the network. These chains reveal the characteristic spin frequencies and their interconnections with high spectral resolution, and can be fused together to map out the network. Our results provide new opportunities for quantum simulations by increasing the number of available spin qubits. Additionally, our methods might find applications in nano-scale imaging of complex spin systems external to the host crystal.

State-of-the-art experiments have demonstrated the imaging of spin networks containing up to 27 nuclear spins [1,19,[33][34][35].The ability to map larger spin networks can be a precursor for quantum simulations that are currently intractable, would provide a precise understanding of the noise environment of spin-qubit registers [32,36,37], and might contribute towards efforts to image complex spin systems outside of the host material [10,20,21,23,38,39].A key open challenge for mapping larger networks is spectral crowding, which causes overlapping signals and introduces ambiguity in the assignment of signals to individual spins and the interactions between them.
Here, we develop correlated sensing sequences that measure both the network connectivity as well as the characteristic spin frequencies with high spectral resolution.We ap-ply these sequences to map a 50-nuclear-spin network comprised of 1225 spin-spin interactions in the vicinity of a nitrogen vacancy (NV) center in diamond.The key concept of our method is to concatenate double-resonance sequences to measure chains of coupled spins through the network.The mapping of spin chains removes ambiguity about how the spins are connected and enables the sensing of spins that are farther away from the electron-spin sensor in spectrally crowded regions.These results significantly increase the size and complexity of the accessible spin network.Additionally, our methods are applicable to a wide variety of systems, and might inspire future methods to magnetically image complex samples such as individual molecules or proteins [10,33,39].

Spin-network mapping
We consider a network of N coupled nuclear spins in the vicinity of a single electron spin that acts as a quantum sensor [1,19].The effective dynamics of the nuclear-spin network, with an external magnetic field along the z-axis, are described by the Hamiltonian (see Supplementary Note 1): where Î(i) z denotes the nuclear Pauli spin-1 2 operator for spin i, A i are the precession frequencies associated with each spin, and C ij denotes the nuclear-nuclear coupling between spin i and j.The frequencies A i might differ due to differences in species (gyromagnetic ratio), the local magnetic field and spin environment, and due to coupling to the sensor electron spin.Our goal is to extract the characteristic spin frequencies FIG. 1. Mapping spin networks.Graph representing a spin network, where vertices denote spins and edges denote spin-spin interactions (Cij).Spins are distributed among spectral regions (coloured disks) by their precession frequency (Ai).a) If all spin frequencies are unique (one spin in each disk), the network can be mapped by measuring only pairwise interactions (C12) between frequencies (A1, A2).b) If spins spectrally overlap (e.g.spins 2 and 5 with A2 ≈ A5) due to the finite line width set by the dephasing time T ⋆ 2 , pairwise measurements alone are ambiguous when assigning interactions to specific spins.By measuring chains (e.g. through A1, A2, A3) we directly retrieve the connectivity of the network.c) We also exploit spin chains to measure interactions between spins that are otherwise challenging to access.As an example, couplings belonging to Spin 4 are not directly accessible from the spin at A1 -due to spectral crowding or negligible couplings -but can be obtained through a chain.d) Finally, we complement the spin chains with a correlated double-resonance method that enhances the spectral resolution for the spin-frequency shifts (∆i) from ∼ 1/T ⋆ 2 to ∼ 1/T2, so that spectrally overlapping spins can also be resolved directly.This figure shows a conceptual network with vertices organized in frequency space.In Supplementary Fig. 1, we discuss the specific relations between frequency and spatial position for the experimental system considered here: an NV center in diamond and surrounding 13 C-spin network, for which increased spectral crowing (panels b-d) naturally occurs for 13 C spins that are farther away from the NV center.
A i and spin-spin couplings C ij that capture the structure of the network.
Fig. 1 shows an example network, with coloured disks denoting frequency regions, and numbered dots inside signifying spins at these frequencies.Although in principle all spins are coupled to all spins, we draw edges only for strong, resolvable, spin-spin couplings, defined by: C ij ≳ 1/T 2 , where T 2 is the nuclear Hahn-echo coherence time (∼ 0.5 s) [5].The network connectivity constitutes the presence (or absence) of such resolvable couplings.In general, the number of frequency disks is smaller than the number of spins, as multiple spins might occupy the same frequency region (i.e.overlap in frequency).
State-of-the-art spin-network mapping relies on isolating individual nuclear-nuclear interactions through spin-echo double resonance (SEDOR) [1].Applying simultaneous echo pulses at frequencies A i and A j preserves the interaction C ij between spins at A i and A j , while decoupling them from (quasi-static) environmental noise and the rest of the network, so that the coupling C ij is encoded in the nuclear-spin polarisation with high spectral resolution (set by the nuclear T 2 -time rather than T ⋆ 2 -time).The signal is acquired by mapping the resulting nuclear spin polarisation, for example at frequency A i , on the NV electron spin and reading it out optically [5].Such a measurement yields a correlated list of three frequencies {A i , C ij , A j } (Fig. 1a).If all spins are spectrally isolated, so that the A i do not overlap, these pairwise measurements completely characterise the network.
However, due to their finite spectral line widths (set by 1/T ⋆ 2 ), multiple spin frequencies A i may overlap (indicated by multiple spins occupying a disk).This introduces ambiguity when assigning measured couplings to specific spins in the network, and causes complex overlapping signals, which are difficult to resolve and interpret [1,19].Figure 1b shows an example where pairwise measurements break down; spins 2 and 5 overlap in frequency (A 2 ≈ A 5 ).Applying pairwise SEDOR between frequencies A 1 , A 3 , A 4 and a frequency that overlaps with A 2 and A 5 returns three independent pairwise correlations: {A 1 , C 12 , A 2 }, {A 3 , C 23 , A 2 } and {A 4 , C 45 , A 5 }.Crucially, however, such measurements cannot distinguish this uncoupled 2-spin and 3-spin chain (Fig 1b) from a single 4-spin network (with a single central spin at A 2 ), nor from a network of 3 uncoupled 2-spin chains (three spectrally overlapping spins).Without introducing additional a-priori knowledge or assumptions about the system, pairwise measurements cannot be assigned to specific spins and are thus insufficient to reconstruct the network [1].
Our approach is to measure connected chains through the network, and combine these with high-resolution spin frequency measurements.First, spin-chain sensing (detailed in Section ) correlates multiple frequencies and spin-spin couplings, directly accessing the underlying network connectivity, and thus reducing ambiguity due to (potential) spectral overlap.Consider the previous example: by probing the correlation between the three frequencies A 1 , A 2 and A 3 in a single measurement, we directly reveal that Spin 1 and Spin 3 are connected to the same spin at A 2 (Spin 2).Such a spinchain measurement yields a correlated list of 5 frequencies: {A 1 , C 12 , A 2 , C 23 , A 3 }, characterising the 3-spin chain.Applying the same method but now with spin 4 (A 3 ← A 4 ) reveals that it is not connected to Spin 2, but couples to another spin (spin 5) that overlaps in frequency with Spin 2.
Second, spin-chain sensing enables measuring couplings that are otherwise challenging to access, enabling exploration further into the network.Consider the case where starting from some spin (e.g.Spin 1 in Fig. 1c) it is challenging to probe a part of the network, either because the couplings to Spin 1 are too weak to be observed or spectral crowding causes signals to overlap.The desired interactions (e.g.those belonging to Spin 4 in Fig. 1c) can be reached by constructing a spin chain, in which each link is formed by a strong and resolvable spin-spin interaction.The chain iteratively unlocks new spins that can be used as sensors of their own local spatial environment.
Finally, we combine the spin-chain measurements with a correlated double-echo spectroscopy scheme that increases the resolution with which different A i are distinguished from ∼ 1/T ⋆ 2 to ∼ 1/T 2 (Fig. 1d).This directly reduces spectral overlap of spin frequencies, further removing ambiguity.
In principle, the entire network can be mapped by expanding and looping a single chain.In practice, measuring limited-size chains is sufficient.A N-spin chain measurement yields a correlated list of N spin frequencies A, alongside N − 1 coupling frequencies C, which quickly becomes uniquely identifiable, even when some spin frequencies in the network are degenerate.This allows for the merging of chains that share a common section to reconstruct the network (Methods).

Experimental system
We demonstrate these methods on a network of 50 13 C spins surrounding a single NV center in diamond at 4 K.The NV electron spin is initialized and measured optically and is used as the sensor spin [1].We employ dynamical decoupling sequences to sense nuclear spins at selected frequency bands, using sequences with and without radio-frequency driving (DDRF) of the nuclear spins to ensure sensitivity in all directions from the NV (Methods) [1,5].The nuclear spins are polarized via the electron spin, using global dynamicalnuclear-polarisation techniques (PulsePol sequence [2,4]), or by selective projective measurements or SWAP gates [1,5].
The 13 C nuclear spin frequencies are given by A i = ω L + m s ∆ i , with ω L the global Larmor frequency and ∆ i a local shift due to the hyperfine interaction with the NV center (see for example Ref. [3] and Supplementary Note 1).Here, we neglected corrections due to the anisotropy of the hyperfine interaction, which are treated in Supplementary Note 4. The experiments are performed with the electronic spin in the m s = ±1 states.Because, for the spins considered, ∆ i is typically two to three orders of magnitude larger than the nuclear-nuclear couplings C ij , nuclear-spin flip-flop interactions are largely frozen, and Eq. 1 applies (Supplementary Note 1).
In the NV-nuclear system, spectral crowding forms a natural challenge for determining the spin-network structure.The spin frequencies are broadened by the inhomogeneous linewidth ∼ 1/T ⋆ 2 , which is mainly set by the coupling to all other nuclear spins.A limited number of nuclear spins close to the NV center are spectrally isolated (defined as: , making them directly accessible with electron-nuclear gates [1,5], and making pairwise measurements sufficient to map the interactions.However, the hyperfine interaction, and thus ∆ i , decreases with distance (∼ r −3 ), resulting in an increasing spectral density for lower ∆ i (larger distance).Interestingly, there exists a spectrally crowded region 2 ) for which nuclear spins still do not couple strongly to other spins in the same spectral region (C ij ≲ 1/T 2 ∀ j ), for example when they are on opposite sides of the NV center.Contrary to previous work [1], the methods outlined in Section allow us to measure interactions between spins in the spectrally crowded region (see Supplementary Note 2), unlocking a part of the network that was previously not accessible.

Spin-chain sensing
We experimentally demonstrate the correlated sensing of spin chains up to five nuclear spins (Fig. 2), by sweeping a multi-dimensional parameter space (set by 5 spin frequencies and 4 spin-spin couplings).We start by polarizing the spin network [2,4] and use the electron spin to sense a nuclear spin (Spin 1) at frequency A 1 , which marks the start of the chain.
First, we perform a double-resonance sensing sequence (Fig. 2b) consisting of a spin-echo sequence at frequency A 1 and an additional π-pulse at frequency RF 2 .The free evolution time t 12 is set to 50ms, to optimise sensitivity to nuclear-nuclear couplings (typically ∼ 10Hz).By sweeping RF 2 , strong connections (C 1j ≫ 1/T 2 ) are revealed through dips in the coherence signal of Spin 1 (Fig. 2d, left).We select a connection to a spin at RF 2 = A 2 (Spin 2) and determine (C 12 ) by sweeping t 12 (Fig. 2d, right).
Next, we extend the chain.To map the state of Spin 2 back to the electron sensor through Spin 1, we change the phase of the first π 2 -pulse (labelled 'yx') and set t 12 = 1/(2C 12 ) to maximise signal transfer (see Supplementary Note 3).We then insert a double-resonance block for frequencies RF 2 = A 2 and RF 3 in front of the sequence (Fig. 2c and e,left) to explore the couplings of Spin 2 to the network.This concatenating procedure can be continued to extend the chain, with up to 5 nuclear spins shown in Fig. 2. In general, the signal strength decreases with increasing chain length, as it is set by a combination of the degree of polarisation and decoherence (T 2 relative to C ij ) of all spins in the chain (See Supplementary Note 3).This limits the chain lengths that can be effectively used.
By mapping back the signal through the spin chain, the five spin frequencies and the 4 coupling frequencies are directly correlated: they are found to originate from the same branch of the network.As spins are now characterized by their connection to the chains, rather than by their individual, generally degenerate, frequencies (Fig. 1b), they can be uniquely identified.Additionally, the chains enable measuring individual spin-spin couplings in spectrally crowded regions (Fig. 1c).As an example, the expected density of spins at frequency A 4 is around 30 spins per kHz (Supplementary Fig. 6), making Spin 4 challenging to access directly from the electron spin.However, because Spin 3 probes only a small part of space, Spin 4 can be accessed through the chain, as demonstrated by the single-frequency oscillation in Fig. 2f.Another advantage over previous methods [1] is that our sequences are sensitive to both the magnitude and the sign of the couplings, at the cost of requiring observable polarisation of the spins in the chain.The sign of the couplings provides additional information for reconstructing the network (Fig. 2g).

High-resolution measurement of spin frequencies
While the sensing of spin chains unlocks new parts of the network and reduces ambiguity by directly mapping the network connections, the spectral resolution for the spin frequencies (A i ) remains limited by the nuclear inhomogeneous dephasing time T ⋆ 2 ∼ 5 ms [5].Next, we demonstrate highresolution (T 2 -limited) measurements of the characteristic spin frequency shifts ∆ i .These frequencies provide a way to label spins, and thus further reduce ambiguity regarding which spins participate in the measured chains, particularly when a spectral region in the chain contains multiple spins (see Fig. 1d).
We isolate the interaction of nuclear spins with the electron spin through an electron-nuclear double-resonance block acting at a selected nuclear-spin frequency region.The key idea is that the frequency shift imprinted by the electronspin sensor can be recoupled by controlling the electron spin state.We use microwave pulses that transfer the electron population from the |−1⟩ to the |+1⟩ state (Fig. 3b, Methods).The nuclear spin is decoupled from quasi-static noise and the rest of the spins, extending its coherence time, while the interaction of interest (∆ i ) is retained.Figure 3 shows an example for a nuclear spin at A 1 , for which we measure a hyperfine shift ∆ 1 = 14549.91(5) Hz and a spectral linewidth of 1.8 Hz (Fig. 3d and e).Besides a tool to distinguish individual spins in the network with high spectral resolution, this method has the potential for improved characterisation of the hyperfine interaction in electron-nuclear spin systems.
The observed coherence time T 2 = 0.36(2)s is slightly shorter than the bare nuclear spin-echo time T 2,SE = 0.62 (5) s.This reduction is caused by a perturbitive component of the hyperfine tensor in combination with the finite magnetic field strength (see Supplementary Note 4).Flipping the electron spin between m s = ±1 changes the quantization axes of the nuclear spins, which causes a change of the nuclearnuclear interactions [1], which is not decoupled by the spinecho sequence (see Supplementary Fig. 1).The effect is strongest for spins near the NV center.For larger fields or for spins with weak hyperfine couplings, we expect that further resolution enhancement is possible by applying multiple refocusing pulses (see Supplementary Note 4).
Finally, we combine spin-chain sensing and electronnuclear double resonance to correlate high-resolution spin frequencies (∆ i ) with specific spin-spin couplings (C ij ), even when a chain contains multiple spins with overlapping frequencies.We illustrate this scheme on a chain of spins, where two spins (2 and 3) have a similar frequency (A 2 ≈ A 3 ) and both couple to A 1 and A 4 (Fig. 4a).The goal is to extract ∆ 2 , ∆ 3 and the couplings to Spin 4 (C 24 , C 34 ).As a reference, standard double-resonance shows a quickly decaying time-domain signal, indicating couplings to multiple spins that are spectrally unresolved (Fig. 4b).
Figure 4c shows how the electron-nuclear double resonance sequence (mint green) is inserted in the spin-chain sequence to perform high-resolution spectroscopy of the A 2 , A 3 frequency region.Sweeping the interaction time t 1 shows multiple frequencies (Fig. 4e), hinting at the existence of multiple spins with approximate frequency A 2 .The result is consistent with two spins at frequencies ∆ 2 = 8019.5(2)Hz and ∆ 3 = 7695.2(1)Hz, split by an internal coupling of C 23 = 7.6(1) Hz (Fig. 4a and Supplementary Fig. 2e,f).
Next, we add a nuclear-nuclear block (pink block in Fig. 4d) and sweep both electron-nuclear (t 1 ) and nuclear-nuclear (t 2 ) double-resonance times to correlate ∆ 2 and ∆ 3 with nuclearnuclear couplings C 24 and C 34 .After the t 1 evolution, the hyperfine shifts ∆ i are imprinted in the z-expectation value of each spin, effectively modulating the nuclear-nuclear couplings observed in t 2 .The 2D power spectral density (PSD) shows signals in two distinct frequency regions along the f 1 -axis, corresponding to ∆ 2 and ∆ 3 (Fig. 4f).Analysing the nuclear-nuclear (f 2 ) signal at these frequencies (Fig. 4g), we find C 24 = −11.8(2)Hz and C 34 = −0.2(5)Hz.We attribute the splitting to the coupling C 23 between Spins 2 and 3 (Methods, Supplementary Fig. 2g,h).Varying RF 4 enables the measurement of the interactions of spins 2 and 3 to other parts of the network (for example to determine C 12 , C 13 ).Beyond the examples shown here, the electronnuclear block can be inserted at specific positions in the spinchain sequence (Fig. 2c) to extract ∆ i of all spins in the chain (Supplementary Fig. 9).

Reconstruction of a 50-spin network
Finally, we apply these methods to map a 50-spin network.The problem resembles a graph search (Methods) [42].By identifying a number of spin chains in the system, and fusing them together based on overlapping sections, we reconstruct the connectivity (Fig. 5).Limited sized chains are sufficient because the couplings are highly non-uniform, so that a few overlapping vertices and edges enable fusing chains with high confidence.We use a total of 249 measured interactions through pairwise and chained measurements.Fusing these together provides a hypothesis for the network connectivity (Fig. 5b).
To validate our solution for the network we use the additional information that the nuclear-nuclear couplings can be modeled as dipolar and attempt to reconstruct the spatial distribution of the spins.Compared to work based on pairwise measurements [1], our spin-chain measurements provide additional information on the connectivity and coupling signs, reducing the complexity of the numerical reconstruction.Additionally, we constrain the position using the measured hyperfine shift ∆ i (Methods).Because the problem is highly overdetermined [1], the fact that a spatial solution is found that closely matches the measured frequencies and assignments validates the obtained network connectivity.Additionally, the reconstruction yields a spatial image of the spin network and predicts the remaining unmeasured 976 spin-spin interactions, most of which are weak (< 1Hz).An overview of the complete 50 spin cluster, characterized by 50 spin frequencies and 1225 spin-spin couplings can be found in Supplementary Table 1 and in Fig. 5b.FIG. 5. Mapping a 50-spin network.a) Schematic illustrating the procedure for mapping large networks.1: Separate high-resolution chains through the network are measured (two example chains shown here).2: We merge chains that share a common section of the network.3: Optionally, an algorithm adapted from [1] estimates the most likely spin positions (Methods), which predicts all unmeasured nuclear-nuclear couplings (dotted lines) and provides a validation for the assignment and merging of step 2. b) Graph of the 50-spin network mapped in this work, with edges indicating spin-spin interactions above 2 Hz and vertex colors denoting spin frequencies Ai.The spins are labelled according to Supplementary Table 1.Black circles indicate the 23 newly mapped spins compared to previous work [1].A 3D spatial image of the network is presented in Supplementary Fig. 3.

DISCUSSION
In conclusion, we developed correlated double-resonance sensing that can map the structure of large networks of coupled spins, with high spectral resolution.We applied these methods to reconstruct a 50-spin network in the vicinity of an NV center in diamond.The methods can be applied to a variety of systems in different platforms, including electron-electron spin networks [7,11,[13][14][15][16][17][18][19]43].Mapping larger spin systems might be in reach using machine-learning-enhanced protocols and sparse or adaptive sampling techniques, which can further reduce acquisition times [44,45].Combined with control fields [2,5,32], the methods developed here provide a basis for universal quantum control and readout of the network, which has applications in quantum simulations of many-body physics [2].Furthermore, the precise characterisation of a 50-spin network provides new opportunities for optimizing quantum control gates in spin qubit registers [5,32,36], for testing theoretical predictions for defect spin systems [13], and for studying coherence of solid-state spins on the microscopic level, including quantitative tests of open quantum systems and approximation of the central spin model [47].Finally, these results might inspire high-resolution nano-MRI of quantum materials and biologically relevant samples outside the host crystal.

Sample and setup
All experiments are performed on a naturally occurring NV center at a temperature of 3.7K (Montana S50 Cryostation), using a home-built confocal microscopy setup.The diamond sample was homoepitaxially grown using chemical vapor deposition and cleaved along the ⟨111⟩ crystal direction (Element Six).The sample has a natural abundance of 13 C (1.1 %).The NV center has been selected on the absence of couplings to 13 C stronger than ≈ 500 kHz.No selection was made on other properties of the 13 C nuclei distribution.A solid immersion lens (SIL) that enhances photon collection efficiency is fabricated around the NV center.A gold stripline is deposited close to the edge of the SIL for applying microwave (MW) and radio-frequency (RF) pulses.An external magnetic field of B z = 403.553G is applied along the symmetry axis of the NV center, using a (temperature-stabilized) permanent neodymium magnet mounted on a piezo stage outside the cryostat [5].The field is aligned to within 0.1 degrees using a thermal echo sequence [1].

Electron and nuclear spins
The sample was previously characterized in Abobeih et al. [1] and the 27 nuclear spins imaged in that work are a subset of the 50 nuclear-spin network presented here.The NV electron spin has a dephasing time of T ⋆ 2 = 4.9(2) µs, a Hahn spin echo time of T 2 = 1.182(5) ms, and a relaxation time of T 1 > 1 hr [1].The spin state is initialized via spin-pumping and read out in a single shot through spin-selective resonant excitation, with fidelities F 0 = 89.3(2)(F 1 = 98.2(1)) for the the m s = 0 (m s = −1) state, resulting in an aver-age fidelity of F avg = 0.938 (2).The readout is corrected for these numbers to obtain a best estimate of the electronic spin state.The nuclear spins have typical dephasing times of T 2 = 5 − 10 ms and Hahn echo T 2,n up to 0.77(4) s [5].T 2 -times for spins with frequencies closer to the nuclear Larmor frequency (∆ i ≲ 5 kHz) typically decrease to below 100 ms (see e.g.Fig. 2g, right panel), as the spin echo simultaneously drives other nuclear spins at these frequencies which are re-coupled to the target (instantaneous diffusion).

Pulse sequences
We drive the electronic m s = 0 ↔ m s = −1 (m s = 0 ↔ m s = +1) spin transitions at 1.746666 (4.008650) GHz with Hermite-shaped pulses.For transferring the electron population from the m s = −1 to the m s = +1 state (Fig. 3 and  4), we apply two consecutive π-pulses at the two MW transitions, spaced by a waiting time of 3 µs.For all experiments, we apply RF pulses with an error-function envelope in the frequency range 400 − 500 kHz.Details on the electronics to generate these pulses can be found in Ref. [2].
For most experiments described in this work, the measurable signal is dependent on the degree of nuclear spin polarisation.We use a dynamical nuclear polarisation sequence, PulsePol,to transfer polarisation from the electron spin to the nuclear spin bath [2,4].The number of repetitions of the sequence is dependent on the specific polarisation dynamics of the spins being used in the given experiment, but ranges from 500-10000.The PulsePol sequence is indicated by the 'Init' block in the sequence schematics.All double resonance sequences follow the convention illustrated in the dotted boxes in Fig. 2b,c and Fig. 3b, where the horizontal grey lines denote different RF frequencies and the top line the electronic MW frequency.The two letters in the double resonance blocks ('xx' or 'yx') denote the rotation axes of the first and final π/2-pulses.The π-pulses (along the x-axis) are applied sequentially (following Ref. [1]).The lengths of all RF pulses are taken into account for calculating the total evolution time.Nuclear spins are read out via the electron by phase-sensitive ('yx') dynamical decoupling; DD or DDRF sequences [5], indicated by the 'RO'-marked block in the sequence schematics.Typically, the spin that is read out with the electron is reinitialised via a SWAP gate before the final SEDOR block in order to maximise its polarisation.However, all experiments presented here can be performed by using just the DNP initialisation, albeit with a slightly lower signal to noise ratio.

2D spectroscopy experiments
For the 2D measurement we concatenate an electronnuclear double resonance with a nuclear-nuclear SEDOR.For every t 1 -point, we acquire 20 t 2 points, ranging from 10 to 260 ms.The final π/2-pulse of the electron double resonance and the first of the SEDOR are not executed, as they can be compiled away.To correct for any slow magnetic field drifts that lead to miscalibration of the two-qubit gate used for read-out, causing a small offset in the measured signal, we set our signal baseline to the mean of the final five points (≈ 200 − 260 ms), where we expect the signal to be mostly decayed.Note that these field drifts do not affect any of the double resonance blocks in which the quantities to be measured are encoded (due to the spin-echo).
Both the 1D (Fig. 4d) and 2D (Fig. 4e) signals are undersampled to reduce the required bandwidth.To extract ∆A 2 , ∆A 3 , we fit a sum of cosines to the time domain signal of Fig. 4d.To extract the frequencies along the f 2 -axis, which encode the nuclear-nuclear couplings (C 24 , C 34 ), we take an (extended) line-cut at f 1 = ∆A 2 and f 1 = ∆A 3 .To increase the signal, we sum over the four bins indicated by the dotted lines.We fit two independent Gaussians to the f 2data to extract C 24 and C 34 .We find splittings of 7.8(2) Hz and 10.2(5) Hz, respectively, whose deviation with respect to measurements in Fig. 4d is unexplained.The skewed configuration of the two peaks (lower left, upper right) is a result of the correlation of the neighbouring spin state between the t 1 and t 2 evolution times.The different ratio of signal amplitudes belonging to Spin 2 and Spin 3, between the 1D and 2D electron-nuclear measurements are due to using different settings for the chained readout (evolution time, RF power).As we are only interested in extracting frequencies, we can tolerate such deviations.Supplementary Figure 2 shows numerical simulations of the experiments presented in Fig. 4.These are generated by evaluating the Hamiltonian in Supplementary Eq. 3, taking into account the two spins at A 2 , A 3 , the spin at A 1 and the electron spin.

Network reconstruction
Here, we outline a general procedure for mapping the network by performing specific spin-chain and high-resolution ∆ i measurements.The mapping-problem resembles a graph search, with the NV electron spin used as root [42].We base the protocol on a breadth-first like search, which yields a spanning tree as output, completely characterising the network.The following pseudocode describes the protocol: New vertices that are detected by chained measurements are iteratively added, once we verify that a vertex was not characterised before (i.e. has a duplicate in the spanning tree T ).The function MeasureCoupling(v, f ) performs a spinecho double resonance sequence between vertex v and a frequency f , (a spin chain of length i − 1 is used to access v) and checks whether a single, resolvable coupling is present (stored in the boolean variable 'singlecoupling').In the case that v is the electron spin (el) an electron-nuclear DD(RF) sequence is performed [3,5].The function CheckVertex(w, T ) instructs the experimenter to perform a number of spin-chain and electron-nuclear double resonance measurements, comparing the vertex w and its position in the network with that of the (possibly duplicate) vertex k (see Supplementary Note 5).If one of these measurements is not consistent with our knowledge of k, we conclude w is a unique vertex and add it to T .If all measurements coincide with our knowledge of k, we conclude it is the same vertex and merge w and k.If the CheckVertex(w, T ) is inconclusive (e.g.due to limited measurement resolution), we do not add w to T .Note that the measurement resolution, determined by the nuclear T 2 -time, is expected to decrease for spins further away from the NV center (See Supplementary Note 2).This eventually limits the number of unique spins that can be identified and added to the network map.
The platform-independent procedure outlined above can be complemented by logic based on the 3D spatial structure of the system [1].For example, when the CheckVertex(w, T ) function is inconclusive, one can sometimes still conclude that w must be unique (or vice versa equal to k), based on the restricted number of possible physical positions of these two spins in 3D space [1].In practice, we alternate the graph search procedure with calls to a positioning algorithm [1], which continuously checks whether the spanning tree T is physical and aides in the identification of possible duplicates.

3D spatial image
For the 3D reconstruction of the network, we use the positioning algorithm developed in Ref. [1].To limit the experimental time we re-use the data of Ref. [1] and add the new measurements to it in an iterative way.We set the tolerance for the difference between measured and calculated couplings to 1 Hz.Although we only measure the new spin-spin couplings and chains when the electron is in the m s = −1 state, we can assume this is within tolerance to the average value of the coupling if the perpendicular hyperfine component is small (< 10kHz) [1].The spin positions are restricted by the diamond lattice.Spins that belong to the same chain are always added in the same iteration and up to 10000 possible configurations are kept.Chains starting from different parts of the known cluster can be positioned in a parallel fashion if they share no spins, reducing computational time.For spins that are relatively far away from the NV, we also make use of the interaction with the electron spin, approximating the hyperfine shift ∆ i to be of dipolar form within a tolerance of 1 kHz (neglecting the Fermi contact term [13]).For those cases, we model the electron spin as a point dipole with origin at the center of mass, as computed by density functional theory [13].If multiple solutions are found, we report the standard deviation of the possible solutions as a measure of the spatial uncertainty (see Supplementary Table 1).

Error model and fitting
Confidence intervals assume the measurement of the electron state is limited by photon shot-noise.The shot-noiselimited model is propagated in an absolute sense, meaning the uncertainty on fit parameters is not rescaled to match the sample variance of the residuals after the fit.For all quoted numbers, the number between brackets indicates one standard deviation or error indicated by the fitting procedure.We calculate the error on the PSD according to Ref. [14], assuming normally distributed errors.
Supplementary Table II.Spin labels of spins featuring in the experiments in the main text.The bands indicate spin frequencies Ai, shifted by the hyperfine interaction ∆i with the electron spin (mint green).This interaction diminishes with distance from the NV center, leading to spectral crowding (multiple spins per band).For simplicity we do not visualize the angular dependence of the hyperfine interaction.Frequencies A1, A2 and A3 contain only a single spin, allowing for a direct readout with the electron spin.Additionally, the measured couplings between these frequencies (e.g.C12) can be assigned unambiguously from pairwise measurements (see Fig. 1a).b) When multiple spins occupy the same frequency band (e.g. at A4), spin-chain measurements resolve ambiguity by retrieving the connectivity of the network (see also Fig. 1b).c) Spins in spectrally crowded areas (e.g. at A6, see also Supplementary Fig. 1), can still be accessed via a spin-chain, starting from a directly accessible spin (e.g. at A2) (see Fig. 1c).d) High-resolution measurements of ∆i (indicated by narrow pink and blue bands) allow for directly distinguishing multiple spins in a single frequency band (e.g. at A4) (see Fig. 1d).Note that we use different frequency labels compared to Fig. 1. ρ0 = ρe,0 with Î the (two-dimensional) identity matrix, the polarisation degree p i ∈ [−1, 1] for spin i and ρe,0 the initial state of the NV electron spin.The p i can be obtained by independent state preparation and measurement characterisations [2].

Readout
For all experiments presented in this work, the signal is read out via dynamical-decoupling sensing sequences (DD or DDRF [3,5]) with the electron spin.Only a limited number of spins can be read out selectively (see Supplementary Table I), due to spectral overlap between nuclear spins.Hence, we choose the first spin in all sensing chains to be one that is directly accessible to the electron spin.Combined initialisation and readout fidelity varies between spins (0.44(2) − 0.95(2), corrected for electron readout fidelity).
To quantify the transition between the spectrally crowded and spatially crowded regions, we compute the expected drop in spin-echo coherence, taking into account instantaneous diffusion [11] (Supplementary Fig 6b).Assuming linear addition of the dephasing rates [12]: with 1/T 2,ID the instantaneous diffusion dephasing rate and T 2 ≈ 500 ms the isolated spin-echo coherence time [5].We compute T 2,ID by examining the mean T ⋆ 2 of the subsystem of spins that occupy a frequency bin of size 100 Hz for all ∆.We mark the start of the spatially crowded region at the point T 2,SE ≈ 0.5 T 2 , which implies |∆| ≲ 1 kHz (Supplementary Fig. 6b).For the 23 newly characterised spins in this work, most belong to the spectrally crowded or spatially crowded regions (8 satisfy 3.5 < |∆ i | < 7.5 kHz, 10 satisfy 1 < |∆ i | < 3.5 kHz, and 4 satisfy |∆ i | < 1 kHz, see Supplementary Table I).
The discussion so far has focused on the degree of spectral crowding using the numerical values for a natural abundance sample (1.1% 13 C).We now examine how these relations depend on the isotope concentration.In the basic picture of a point-dipole electron spin, surrounded by a bath of nuclear spins with spatial density ρ 0 , the system exhibits a form of scale invariance.That is, we can define the (dimensionless) spectral density: describing the number of spins per line width.For dipolar interactions, both the electron-nuclear hyperfine shift ∆ as well as the nuclear line width (∼ 1/T ⋆ 2 ) scale linearly with the density, so that the dimensionless spectral density (Supplementary Eq. 11) is independent of ρ 0 .Intuitively, this can by understood as both the nuclear-spin line widths and the spacing between nuclear-spin frequencies scaling with ρ 0 , keeping their ratio (i.e. the degree of spectral crowding) constant.In principle, in this elementary system, the physics and quantities like the number of spins that can be mapped are independent of concentration, with only absolute time and distances being rescaled.
In practice, however, this scale invariance breaks down for both for low and high isotope concentrations.At low isotope concentrations, other noise sources will start to limit both the nuclear line widths (via reduced T ⋆ 2 ), as well as the spin-echo coherence times (via reduced T 2 ).At high concentrations -including around the natural 1.1% abundance -the discreteness of the lattice and the contact hyperfine interaction due to the finite NV electron wave function need to be taken into account.Predictions for the optimal concentration for a given goal, such as controlling the largest network, likely need detailed numerical simulations taking these details into account, which we do not pursue here.

FIG. 2 .
FIG.2.Sensing spin chains.a) Schematic of a N = 5 nuclear-spin chain through different spectral regions {A1, A2, A3, A4, A5} (coloured disks), starting from the NV electron spin ('el').Even though there might be multiple spins at each of the nuclear frequencies, only a single one is connected to this chain.b) Pulse sequence (Methods) for the prototypical N = 2 sequence (SEDOR)[1].c) Pulse sequence for sensing a chain of N = 5 nuclear spins, correlating 5 spin frequencies and 4 spin-spin interactions.In this case, the RF frequency (RF5) and freeevolution time (t45) are varied to probe the connections of the spin at A4 to other spins.The resulting signal is mapped back via concatenated SEDOR sequences and finally read out ('RO') through the electron spin (Methods).d-g) Experimental data, sweeping the frequency RFN of the recoupling pulse (left) to detect the frequencies of spins coupled to Spin 1, and varying the free evolution time tN−1,N (right) to extract their coupling strengths (for N = {2, 3, 4, 5}).For the frequency sweeps, evolution times tij are selected a-priori (annotated).Colored highlights denote the signals due to the spins in the chain and solid lines are fits to the data (see Supplementary Note 3).The signal in the bottom panel is inverted, due to the coupling C34 being negative.

FIG. 3 .
FIG. 3. Electron-nuclear double resonance.a) The nuclear-spin frequencies Ai are shifted by the hyperfine interaction with the electron spin ∆i (blue dotted line).Performing double resonance between the nuclear spin and the electron spin (mint green) retains this interaction while decoupling from quasi-static noise.b) Pulse sequence for measuring ∆1 for a spin at ≈ A1.The nuclear spin undergoes a double resonance sequence, picking up a phase (downward arrow) from the interaction with the electron spin, whose population is synchronously transferred from the |−1⟩ to the |+1⟩ state.Finally, the signal is read out (denoted 'RO') via the electron spin (Methods) c) Time domain signal of ∆1 = 14549.91(5)Hz (undersampled), with a coherence time of T2 = 0.36(2) s, fitted by a sinusoid with Gaussian decay.d) Zoom-in of spectroscopy data as in Fig. 2d, showing a broad resonance (135 Hz FWHM), limited by the nuclear T ⋆ 2 -time.e) Power spectral density (PSD) of (d), showing a line width that is ∼ 75 times improved compared to (c).

FIG. 4 .
FIG.4.Two-dimensional spectroscopy of spectrally crowded spins.a) Schematic of the studied system, which contains two spins with overlapping frequencies A2 ≈ A3 (grey circle), with slightly different hyperfine shifts (∆2, ∆3).Both are coupled to the spin at A1, which is used for to transfer the signal to the electron spin for readout ('RO').b) SEDOR spectroscopy (as in Fig.2d) of the frequency region A2, A3, with the estimated spin frequencies indicated.Sweeping the t12 evolution time results in a quick decay.c) Pulse sequence for the electron-nuclear double resonance sequence used in (e), where the ∆i are extracted by sweeping t1.d) Pulse sequence combining electronnuclear and nuclear-nuclear double resonance, used in (f).Adding a nuclear-nuclear block (pink) and sweeping both t1 and t2 reveals the correlation between ∆i and spin-spin couplings.e) Sweeping t1 yields a high-resolution PSD of the A2, A3 frequency region, showing two (split) frequencies ∆2 and ∆3.The solid curve is a four-frequency fit to the data.f) Signal (PSD) for the two-dimensional sequence, revealing two distinct regions along the f1-axis at ∆2 and ∆3.g) Binned line cut of (f) along the f2-axis at frequencies ∆2, ∆3 (region indicated by dotted lines).The positions of the (split) peaks indicate the coupling to the spin at A4 (C24 = −11.8(2)Hz, C34 = −0.2(5)Hz).The solid line is a fit of two Gaussians to extract the couplings.

. 1 .
FigureSpin number Spin label and frequency f if singlecoupling then ▷ Checks if MeasureCoupling returned a single, resolvable coupling create vertex w A w = f C vw = C unique, duplicate = CheckVertex(w, T ) ▷ Checks if w was already mapped in T if unique then ▷ w was not yet mapped add w to V i+1 in T ▷ w is added to T as a new vertex end if if not unique and duplicate == k then ▷ w is the same vertex as k in T add C vk = C vw in T ▷ The measured coupling is assigned to k

Table I :
Information on the nuclear spins mapped in this work.A dash denotes that no ∆i signal could be obtained due to low polarisation, coherence, or readout contrast.Value in parentheses denotes standard deviation on the last digit.Spins C1-C27 were previously characterised in Ref.[1]and the used labels are consistent with their work.These data are also available online at: https://doi.org/10.4121/aba1cc84-0aea-4cdc-93ca-68b0db38bd81.v1