Correlated Spectroscopy of Electric Noise with Color Center Clusters

Experimental noise often contains information about the interactions of a system with its environment, but establishing a relation between the measured time fluctuations and the underlying physical observables is rarely apparent. Here, we leverage a multidimensional and multisensor analysis of spectral diffusion to investigate the dynamics of trapped carriers near subdiffraction clusters of nitrogen-vacancy (NV) centers in diamond. We establish statistical correlations in the spectral fluctuations we measure as we recursively probe the cluster optical resonances, which we then exploit to reveal proximal traps. Further, we deterministically induce Stark shifts in the cluster spectrum, ultimately allowing us to pinpoint the relative three-dimensional positions of interacting NVs as well as the location and charge sign of surrounding traps. Our results can be generalized to other color centers and provide opportunities for the characterization of photocarrier dynamics in semiconductors and the manipulation of nanoscale spin-qubit clusters connected via electric fields.

Often seen as detrimental, random fluctuations in the response of a probed system can nonetheless shed light on otherwise hidden physical processes.A paradigm example is Johnson noise, the intrinsic voltage fluctuations in a resistor, whose root mean square amplitude can be tied to the system temperature 1 .Another illustration can be found in the statistical fluctuations of a paramagnetic ensemble 2 , where random spin coherences forming and decaying with characteristic energies and time scales allow the experimenter to reconstruct the system's magnetic resonance spectrum in the absence of resonant drives 3,4 .Further, the ability to establish correlations in the observed fluctuations brings in the resolving power of multidimensional spectroscopy, broadly employed in magnetic resonance 5 .
While fundamental fluctuations are readily detectable in macroscopic systems 3,4 , experiments at the nanoscale are often better suited for noise spectroscopy because the fractional change on the observed random signal increases as the system size reduces [6][7][8] .In particular, simultaneous correlations in space and time -recently proposed as a strategy to enhance the information content in quantum sensing experiments 9 -are intrinsically easier to capture in nm-sized systems.
Here we monitor the optical transitions of individual NVs within clusters sharing the same diffraction-limited volume as we randomly alter the occupation of proximal charge traps.Comparing the spectra from multiple measurements -affected by reconfiguring electric environments -we extract single-and multi-NV spectral correlations that we then use to map out the NV relative positions in three-dimensions and co-locate proximal charge traps within the crystal host.Capitalizing on the dual role of NVs -alternatively serving as an electric field probe or a carrier trap -we further illustrate how controlled ionization of an individual center from a pair in the cluster allows us to deterministically change the optical resonances of the other.Related work has been reported recently, both for NVs in diamond 10 and organic color centers in carbon nanotubes 11 .

NV electronic structure and experimental design
Throughout our experiments, we use narrow-band, tunable laser excitation (637 nm) to monitor small sets of negatively charged nitrogen-vacancy (NV) centers in diamond via confocal microscopy.Formed by a substitutional nitrogen immediately adjacent to a vacancy 12 , these spin-active color centers are presently attracting broad attention as a platform for quantum information processing 13 and nanoscale sensing 14 .At low temperatures, the zero-phonon line features a fine structure, best captured through the energy diagram in Fig. 1a: Importantly, the 3 E excited state manifold is highly sensitive to the crystal environment 12 , resulting in a set of optical resonances whose exact frequencies depend on the transverse and longitudinal components of the local electric and strain fields (respectively,  ! and  ∥ ).Negative to most applications 13 , this NV susceptibility is instead valuable herein, as it heralds physical reconfigurations in the occupation of proximal charge traps (schematic in Fig. 1a).
We introduce our working strategy in Fig. 1b where we plot a series of photo-luminescence excitation (PLE) spectra from a first NV cluster.Combined use of selective NV -ionization and microwave excitation allows us to identify three NVs within the same diffraction-limited volume (hereafter denoted A, B, and C, see also Supplementary Material, Sections 1 through 4).Spectral 'diffusion' of the optical resonances is apparent as we intercalate pulses of green light (532 nm) between successive sweeps of the tunable laser.The power and duration of these pulses -adjusted to cycle the NV charge state between negative and neutral -also lead to concomitant changes in the occupation of proximal traps, hence resulting in a varying electric environment.
While fluctuations in the occupation of the many traps far from an NV induce a near-continuous diffusion of the optical spectrum, proximal charge state changes must arguably lead to discernible spectral jumps.Further, the energy diagram in Fig. 1a ties the optical spectra to  ! and  ∥ (respectively associated to shifts between optical resonances, and of the entire multiplet 15 ), which can be leveraged to extract information not apparent when considering a time-averaged spectrum (upper insert in Fig. 1b).We first demonstrate this notion in Fig. 1c where we build on the spectra measured for NVB to derive a twodimensional (2D) field histogram correlating the values derived for  ! and  ∥ .We find a highly non-uniform distribution, which hints at discrete jumps along a quasi-Gaussian electric noise background (see below).Crucially, data bunching is far less obvious if one only considers each field projection separately (upper and right inserts in Fig. 1b) or, more directly, the one-dimensional (1D) probability distributions of individual optical resonances, which underscores the need for a correlated, multivariate analysis.

Mapping out nearby charge traps
To better illustrate the advantages and limitations intrinsic to this class of noise spectroscopy we resort to NVC where the 2D field histogram of  ! and  ∥ reveals a more complex, though richer structure (Fig. 2a).To interpret our observations, we model the electric field environment as the stochastic sum of contributions from carriers at three proximal sites (Figs.2b and 2c).Assuming each trap intermittently hosts one fundamental charge, we converge to a set of four "primary" configurations corresponding to states where all traps are empty, or where a carrier occupies one of the three possible sites.
One can now combine these alternative scenarios to predict the transverse and longitudinal fields expected in cases where more than one trap is occupied.We demonstrate this notion in Fig. 2b where we add the fields produced by carriers in traps 1 and 2 (respectively, ( # , 0,0) and (0,  $ , 0) in Fig. 2b) to obtain the data subset in the upper right corner of the plot, ( # ,  $ , 0); it is easy to see the same applies to the two alternative cases (sites 1 and 3 or 2 and 3, faint vectors in Fig. 2b).Since transverse components -only determined in magnitude -follow correctly from a vector sum, we conclude that all traps must approximately lie on the same plane; further, because each individual trap hosts on average one carrier (13±2)% of the time (Fig. 2a), we expect to find two simultaneously full traps with a probability of only (2±1)%, in good agreement with our observations ((1.6±1)%, grey ellipses in Figs.2a and 2b).
Unfortunately, changes in the transverse and longitudinal fields as measured from a single NV are insufficient to completely determine the position and charge of a trap.In addition, our measurement uncertainty creates a spread in space for the probability distribution of the trap's position (see Supplementary Information, Sections 5 and 6).We plot in Fig 2c the isosurfaces of such probability distributions calculated for the trap set probed by NVC (Figs. 2a and 2b), assuming negative charges.The presence of a background bias field -combining the local strain and residual electric fields from traps not ionized by green illumination, see Supplemental Information, Section 7 -makes the shapes of these curves non-trivial.Interestingly, we can exclude sections of the solution set incompatible with observations where more than one trap is populated (red segments in each loop in Fig. 2c).We show below this form of "co-measurement" can be adapted to cases where more than one NV picks up the field from a single trap, an approach that dramatically increases the information on both the location and sign of the captured carrier.

Correlated spectroscopy across multiple color centers
Optical excitation is known to cycle the NV charge state between neutral and negatively charged 16 respectively, NV -and NV 0 -implying that NVs alternatively act as local probes or as point sources of electric field.A key implication is that a change in the charge state of a given NV must correspondingly lead to observable spectral shifts in their neighbors.We illustrate this idea in Figs.3a and 3b where we study a new subdiffraction-size cluster comprising four NVs (here labeled with letters D through G).The green laser pulse we apply prior to the optical frequency sweep randomly initializes the charge state of each NV into negative or neutral, hence making them observable or not upon resonant excitation (the zero-phonon lines of NV 0 lie far away from the scanning range of the tunable laser 17 ).Specifically, the yellow arrows on the right-hand side of Fig. 3b highlight instances where the NV % & resonances are missing, indicative of green-induced ionization (typically occurring with ~20% probability 16,18 ).We observe in each case a change in the relative peak amplitudes of NV ' & as well as a strong blueshift of its resonance frequencies (orange arrows); indeed, we conclude from the field histogram in Fig. 3c that the correspondence of these shifts with a direct optical read-out of charge state of NVG is nearly perfect.
The use of narrow-band excitation gives us the opportunity to gain control on the spectral response upon selective NV ionization.We validate this idea in Fig. 3d, where we first use weak red light to post-select an instance in which both NVF and NVG are negatively charged (blue trace).We subsequently tune the laser frequency and amplitude to ionize only NVG 18 ; following a second spectral sweep, the concomitant spectral changes we induce in NVF become apparent (red trace).Note that the converse experiment -where we probe NVG upon selectively ionizing NVF, Fig. 3e -yields analogous phenomenology, although the magnitude of the frequency shift is substantially different.This asymmetry reflects on the different orientations of the transverse background bias fields acting on each NV.
The longer range of Coulombic couplings hints at a fully interacting cluster, where altering the charge state of one of the traps reverberates on all others.Exposing the effect of more weakly coupled NVs, however, becomes challenging as the spectral diffusion caused by many coexisting traps becomes dominant.We circumvent this complication in Fig. 4a where we cross-correlate the dfields acting on NVD and NVF.We find that the charge state of NVG -here serving as a third classifier, blue and red dots -has an impact not only on  ∥ () !but, more importantly, on the mean values of  ∥ () " and  ! ()" , a response that reveals NVD's coupling to NVG (we find a similar response in the case where we exchange the roles of NVF and NVG, hence allowing us to conclude NVD also couples to NVF).Unlike green illumination, weak red excitation has a reduced effect on the charge states of NVs and most traps, hence suggesting the use of time correlations between consecutive spectral sweeps as a second alternative to mitigating spectral diffusion.Figure 4b shows the field histogram of NVD as well as the difference between the fields measured via successive spectral sweeps in the absence of NV charge initialization by green light.This time-correlation measurement effectively suppresses slow background spectral diffusion, and dramatically sharpens clusters created by light-induced changes in the charge states of NVF and NVG (Supplementary Material, Sections 7 and 8).The flipside is a longer experimental time and a partial loss of information as the population statistics of a trap cannot be determined.Worth noting, a similar analysis on NVE shows much weaker trends suggesting this color center is far from the rest (our diffraction limited volume stretches 1-2 µm along beam axis).
We can now combine the information gathered thus far to map the NV cluster in three dimensions (Fig. 4c).In short, we use a search algorithm that tests all positions of the NVs with respect to each other, using measurements from the field histograms of NVF and NVG to calculate the relevant probability distributions (Supplementary Information, Section 9). Figure 4c illustrates such a test performed for the most likely positions (spheres), showing perfect intersection between independent measurements of NVD (cyan and purple loops) and NVF (green and magenta).We estimate the relative distances of NVF and NVD relative to NVG respectively as (48 ± 1.6) nm and (150 ± 23) nm, with the greater uncertainty arising from the weaker couplings.
Recurrent optical spectroscopy over an extended data set shows simultaneous shifts of NVF and NVG approximately 2% of the time, which we interpret as a manifestation of intermittent carrier capture by a trap proximal to both NVs.Here too, we leverage the set of solutions extracted from either NV to co-locate the trap position, approximately (29 ± 2) nm from NVF ((60 ±3) nm from NVG, Fig. 4d).Note the probability distributions do not fully intersect, likely a consequence of underestimated error sources during data analysis.Interestingly, changing the nature of the trapped carrier from an electron to a hole leads to clearly disjoint solution sets and can be ruled out (fainter loops in Fig 4d).Colocalization of a trap, therefore, lifts ambiguities in the captured charge sign (and amplitude) inherent to single NV probe sensing, in the process providing clues on the physical nature of the trap: Potential candidates include a substitutional nitrogen impurity (here transitioning from N + to N 0 ) or a lattice vacancy (sporadically changing from V 0 to V -).
Figure 4e integrates the findings above into a combined 3D plot that includes all three interacting NVs in the cluster and dark trap DTA, as well as a second trap, DTB, observed to capture an electron with a probability -19% on average -that seems to depend on the charge state of NVF (Supplementary Material, Section 10); this occupancy is much higher than that of DTA (2%), and hints at a different physical nature.To compare the cluster geometry against that derived from optical imaging, we leverage the known NV orientations to determine the transformation matrix connecting the xyz frame in the figure to the x'y'z' reference frame in the lab.As shown in the upper inserts, projecting the calculated positions for NVD, NVF, and NVG onto the x'y'-plane -perpendicular to the incoming laser beam -yields an image consistent with that attained via super-resolution microscopy.

Conclusions and outlook
Statistical correlations in the spectral noise affecting the optical resonances of sub-diffraction NV clusters contain valuable information on the underlying electrostatic fluctuations of the crystal environment.Specifically, we exploit otherwise detrimental spectral diffusion in the optical response of an NV to single out proximal charge traps and set bounds on their physical locations.Further, we build on the ability to resonantly ionize NV -to expose Coulombic couplings between the NVs of a cluster -even in the presence of substantial environmental noise -which then allow us to determine their relative positions in three dimensions.Central to these findings is the notion of co-measurement by two (or more) atomic size probes, which we apply herein to map out the position of a remonte color center in a cluster of several, and to pinpoint the location of a non-fluorescent charge trap as well as the magnitude and sign of the carrier it captures.
While NV centers are key to the present findings, most ideas can be extended to other color centers in diamond or alternative material hosts, provided a known relation exists between the observed optical resonance spectrum and the electric environment 19,20 .Relevant examples include the silicon vacancy 21 (VSi) and carbon-silicon di-vacancy 22 (VSiVC) in SiC, as well as group-IV vacancy color centers in diamond [23][24][25][26] (which can be sensitive to proximal traps despite their symmetry-based first-order protection against electric fluctuations 27 ).From a methodological point of view, one can envision extensions in the form of protocols adapted to investigating the trap response under optical excitation not affecting the NV charge state or, alternatively, to probing the diffusion of electrons injected from proximal NVs serving as a source 28,29 .
Our approach paves the route to a deeper understanding of the microscopic mechanisms underlying spectral diffusion, and hence promises opportunities to developing novel schemes for electric field sensing 30 , or for tasks in quantum information processing relying on indistinguishable photons 31 .Along related lines, the comparatively long range of Coulombic couplings could be exploited to mediate interactions between spin qubits otherwise too far from each other to couple magnetically 32,33 .Besides applications in quantum science and technology, future material science studies will be required to shed light on the formation of NV clusters in bulk diamond -observed in our crystals with unanticipated frequency -as an intriguing alternative to implanting nitrogen-rich moieties 34 .grow as they are preferentially associated to  , = ±1 spin states.
We determine the spatial positions of all NVs in a subdiffraction cluster by implementing confocal imaging under resonant excitation 18 .To this end, we first tune the laser frequency to one of the optical transitions pertaining to the NV of interest (assumed known after proper assignment), and then implement the repump PLE protocol as we galvo-scan the laser across the field of view.We finally determine the NV position with sub-diffraction resolution from a two-dimensional point-spread Gaussian fit of the resulting image.

Automated analysis of multiple consecutive spectra
We process all optical spectra in a series via a customdesigned algorithm, the purpose of which is to determine throughout a given frequency scan the spectral positions of at least two optical resonances (including  0 and  1 ) for each NV in a cluster.A prerequisite in implementing this protocol is the manual assignment of all resonances (see Section on PLE assignment above), which we then leverage to properly configure our code for automated identification across a spectral series.In particular, we rely on the early assignments and thorough visual inspection of the spectral series to preset "search bands" for each resonance, i.e., spectral windows where we anticipate finding a PLE line (Fig. S3a).Note that for NVs known to experience large discrete jumps (e.g., NVF and NVG in the cluster of Fig. 3 in the main text), more than one band might be assigned to a given transition (see solid and dashed pink rectangles in Supplementary Fig. 4a).
Whenever two resonances (often belonging to two different NVs) overlap, we make use of other, more isolated optical transitions of the same NV to narrow down the likely frequencies (e.g., we use the strong correlation between the  $ and  0 resonances to identify the approximate frequency of transition  0 , even when proximal to a resonance from another NV).We extract the central frequency, width, and amplitude of each PLE line from a Gaussian fit, which we activate only when the observed PLE amplitude exceeds a predefined threshold (Supplementary Fig. 4b).
For a given set of NV resonances, we determine  !from the set of eigenenergies derived from diagonalizing  *+ .In principle, any one pair of two lines suffices, but some pairs are more sensitive (e.g.,  0 ,  1 ) than others (e.g.,  # ,  $ ).When relying on more than two resonances, we ensure optimal accuracy by first calculating the electric field from every combination of PLE lines, and subsequently determining an average weighted by the derivative of the extracted electric field relative to the frequency difference between the corresponding resonances.For each NV in a cluster, we use the same number of resonances -2 to 4 depending on the NV and experiment -to determine the electric fields for all spectra in the series.
While  !relates to the frequency differences between resonances, we extract  ∥ by comparing the average frequency in each NV spectrum to a fixed reference; note this reference can be arbitrary as our measurements are only susceptible to electric field changes (i.e., we cannot determine the absolute number of elementary charges in a  (') in the first spectral band (solid pink rectangle).For an amplitude threshold of 0.075 kcps (set by default in our analysis code), we find a PLE line in 61% of the spectra.
trap).Further, charges of opposite signs lead to different spatial probability distributions, hence creating an ambiguity that can only be lifted by a co-analysis of the field fluctuations experienced by a second probe NV.We return to this point below.
When analyzing the impact of a given NV (later referred to as the "source" NV) on another one in the cluster (the "probe" NV), we rely on the presence or absence of a reference optical line in the source NV spectrum to flag its charge state during a PLE sweep.Since spectral scans proceed from lower to higher frequencies, we choose this reference as the highest-energy optical transition in the source NV spectrum so as to mitigate errors arising from mid-scan ionization.
Throughout the manuscript, we express d-fields (i.e., strain/electric fields) in GHz.To convert between units, we use the proportionality factor  2 = 6.3 %34 5) 7 #$ , corresponding to a dipole moment difference Δ = 1.3 Debye between the NV -ground and excited states 43 .Interestingly, this value is not known with accuracy 12,[44][45][46] , and could be re-derived by combining the present technique with a super-resolution measurement of two NVs in the same plane.

Charge trap mapping via a single NV
Let us consider an individual charge trap near a probe NV and assume we can measure its charge state.By studying the NV field histogram, we can calculate the average magnitude of the electric field transverse and longitudinal components at the NV site depending on the charge state of the trap, which allows us to set boundaries on the trap's spatial location.Since we only determine the transverse field magnitude and not its orientation, the angle  describing the transverse electric field rotation upon an electric field change remains undetermined, implying the solution takes generically the form of a closed loop in 3D space (Supplementary Figure 5).
More formally, let ℰ ∥ , ℰ !(respectively,  ∥ ,  ! ) be the longitudinal and transverse fields with (without) an added charge .Assuming the NV sits at the origin, we find that for a given rotation  of the transverse field, the charge trap must be located at a position (, ℰ ∥ −  ∥ , ℰ !,  ! ) given by Where  8 is the dielectric constant of diamond and   denotes the field due to , here expressed as The above formula for the trap position  is only valid for exact measurements, and we must in practice deal with the uncertainty created by other fluctuations of the field or by experimental noise.We now focus on calculating the conditional probability distribution (|) of finding a charge at position  given a measurement  = (Δℰ ∥ , ℰ !,  ! ) with variance , where Δℰ ∥ = ℰ ∥ −  ∥ , and  @ $ denotes the variance for variable .We follow two alternative routes to calculate (|), each featuring complementary computational speed and accuracy.Our first strategy determines the most likely Assuming no prior on the charge position, () is a constant thus allowing us to cast (|) in the form where we have used the relation () = ∭ (|′)(′) A: and the integral extends over the diffraction volume.(|) can be calculated using the field histogram created by a charge at position  as well as the measurement variance.As an illustration, Supplementary Fig. 6a shows the probability distribution for the positions of NVG and NVD as seen by NVF.The Bayesian approach above is accurate but too resource intensive for cases where many iterations are necessary (the case in Section 7, below).We can gain computational speed via a simpler strategy that relies on summing many Gaussian distributions obtained by the impact of measurement errors at a given point.We can gain computational speed via a simpler strategy that relies on a Gaussian distribution for all errors.Assuming that the presence of a charge  rotates the transverse field by a known angle  B , we calculate the probability (|,  B ) from the error propagation formula as the 3D Gaussian distribution centered at   = ( B , ℰ ∥ −  ∥ , ℰ !,  ! ) with variances defined as and can slowly drift over timescale of several hours to days, as shown in Supplementary Fig. 7 for the NV cluster in Fig. 3 of the main text.Alternative mechanisms could be responsible for these drifts.As observed in Monge et al. 18 and consistent with works on bound membranes 48 , the diamond attachment method impacts the overall strain in the crystal, implying that minute, slow changes can create drifts.Changes could also be caused by a progressive rearrangement of a metastable electric environment (in the bulk, or on the surface).Indeed, the 100 MHz shift observed here could be induced by a change of one elementary charge at a 150-nm distance, 40 charges at 1 µm, or 1000 charges at 5 µm; such rearrangements have been observed 49 .All field histograms in the manuscript were obtained from spectral series showing no visible drifts.For Fig. 4 of the main text, two datasets obtained months apart and featuring spectral shift differing by hundreds of MHz were used separately to map out all three NVs (and the adjacent dark traps).Importantly, we found the resulting NV localization to be robust from one dataset to another provided we account for a rotation of the transverse field.While these drifts complicate the analysis of field histograms, one could likely apply a correction, or exploit the changing orientation of the transverse field to enhance sensitivity.The same multi-sensor spatial analysis can also be extended to the slower drift to determine their origin.

Time-correlation of the electric field
At low enough excitation power, we can perform several repump-free PLE sweeps in a row without green re-initialization of the charge environment.This opens the prospect of studying the evolution of the electric environment over time.Here, we recorded a series of 270 spectra under 3 nW of red excitation with one repump every three spectra.Note that under these parameters, the probabilistic ionization of the NVs we read out remains important.This could be overcome by using a combination of higher NA, lower power, lower detector dark count and higher acquisition time.For each three-PLE-sweep set, we track the electric field shift between consecutive spectra as well as the evolution of the charge state of all NVs.Fig. 4b of the main text shows the field shift histograms we measure.We observe this time-evolution of the electric environment under resonant excitation is slow such that the field shift histograms are dominated by NV ionization events.By measuring the field shifts instead of the absolute value of the electric field, we negate background contributions that do not vary, which improves the minimum shift we can measure and increases the sensitivity and range of our sensor.

Localization using a tri-partite interaction
In general, one can use the multipartite interaction in a cluster of more than two NVs to infer their relative positions and local bias fields.In the present case, we developed an algorithm adapted to 3 NVs, two of whom have the same crystalline axis (NVF and NVG in the cluster of Figs. 3 and 4 of the main text).As a brief overview, we first assume NVG is at the origin, and set its crystalline (bias d-field) axis along z (x).We then proceed to calculate the probability () for the relative positions of NVD and NVF, respectively  I and  ' , as well as the relative angle  J between the transverse bias fields of NVF and NVG (i.e., the angle between the average of  ! ()! and  ! ()* ).
Using the shortened notation  = ( I ,  ' ,  J ), we subsequently assess the probability () by sequentially combining measurements obtained from the field histograms of NVF and NVG.Throughout our analysis, we use each successive measurement to eliminate unlikely values of  and accelerate the search.The resulting probability density is then projected on the  I and  ' subspaces, yielding probability densities for those two NVs and an average value for  J = 80 ± 20°.
In the following, we refer to the NVs generating and sensing the electric field as the "source" and "probe" NVs, respectively.The steps below describe the search algorithm in further depth.
1 -For all combinations of source and probe NVs, respectively NV K and NV L , ,  = D, E, F, G, we calculate the probability distribution for NV K relative to NV L as the the transverse field is a vector, we have: We compare the calculated values for  ,! with that obtained experimentally from the field histogram  V,! *0W , recording the mismatch normalized by the experimental standard deviation, namely,  X,! u, , ,  B and subsequently weigh the probability densities of the remaining ones.
3 -We test the most likely relative position between two NVs of the same orientation, NV K and NV L .To do so, we iterate through all Gaussian distributions   K,L and   + L,K .
Note that crystallographically equivalent NVs do not necessarily experience local (transverse or longitudinal) d-fields of the same orientation.We therefore also test the two possible longitudinal strain directions and sweep the angle  Y K,L between the transverse strain of the two NVs. 4 -For every parameter set  B ,  B + ,  Y , we then calculate the position distribution of a third NV in two independent ways: directly relative to NV K , and using the relative positions of NV K and NV L ,  B,B + ,Z K,L , as well as the known distribution of position for NV O relative to NV L .We proceed similarly to step 3, i.e., we discard cases where the mismatch of the two calculated positions is large compared to the combined variances, and multiply the remaining Gaussian distributions to find a new one for the position of NV O relative to NV L .We determine a new renormalization factor for each distribution.
5 -We sum all Gaussian distributions for the positions of NV K and NV O relative to NV L , weighing each with their two renormalization factors.We finally renormalize each of the two distributions for NV K and NV O .

Co-localization of dark charge traps
As seen in Fig. 2, the NV cluster in Fig. 3 of the main text also contains non-fluorescent (i.e., "dark") charge traps.In an extended dataset partially shown in Supplementary Fig. 8, we observe discrete shifts of the resonances of both NVF and NVG, which we attribute to carrier capture by a dark trap DTA.These shifts occur 5 times in a series of 225 spectra (2.2±1%).Note that we observed a distinct, stronger shift with similar occurrences that could not be analyzed due to the limited range of our laser sweeps.In addition, the field histogram of NVG displays discrete jumps that do not correlate with the charge state of any other NV and that we attribute to carrier capture by a second dark trap, DTB.On average, these shifts occur in 19% of all spectra (13% and 25% when NVF is in its negative and neutral state, respectively).Once a trap is identified, we calculate the probability density for its location relative to NVF and NVG assuming a given charge (one electron or one hole).We then use our knowledge of the position and strain of the two NVs to combine that information by simply multiplying and renormalizing the two probability densities (Figs.4d and 4e in the main text).The relative quality of the intersection of the two loops informs us on the charge sign: By

Figure 1 |
Figure 1 | Optical spectroscopy of diffraction-limited NV clusters.(a) (Left) We use resonant confocal microscopy to individually address NVs from a small set sharing the same diffraction-limited volume.(Right) Under cryogenic conditions, the NV -excited state manifold 3 E splits into two triplets yielding a collection of optical resonances around 637 nm whose values depend on the electric and strain fields at each NV site.(b) Recursive photo-luminescence excitation (PLE) spectroscopy of a cluster featuring three different NVs.The bottom plot shows successive spectra under simultaneous MW (2.88 GHz) and tunable 637 nm light preceded by charge initialization with 532 nm light.The upper 1D plot is the integrated sum of all individual traces; spectral-diffusion-induced broadening of all resonances is apparent.(c) Correlated longitudinal and transverse fields on NVB as extracted from the spectra in (b).Data clustering (indicated by circles) is apparent in the 2D plot but not so much in the 1D projected histograms (upper and right inserts).In (b), the green (red) laser power during charge initialization (frequency sweep) is 1 µW (3 nW) and the reference frequency is 470.000THz.Unless otherwise noted, all experiments are carried out at 7 K in the absence of any externally applied magnetic field.

Figure 2 |
Figure 2 | Mapping out proximal charge traps.(a) Longitudinal and transverse fields for NVC as extracted from the spectroscopy set in Fig. 1b;  indicates the total number of observations, and  is the number of data points enclosed in each ellipse.(b) We model the electric field on NVC as the combined effect from point charges in three proximal traps and a fluctuating field environment of more distant carriers.The probability of a given charge configuration can be extracted from the fractional weight of each data cluster in the plot where not more than one trap is populated.This information is sufficient to predict the fields and occurrence probability of those configurations where more than one trap is occupied (see crosses within gray ellipses).(c) Probability isosurfaces for the point charge positions yielding the average fields shown in (b) for the case of electron-populated traps; percent values denote confidence intervals.Red lines indicate sections of otherwise possible solutions, here discarded based on observations with two simultaneously populated traps.

Figure 3 |
Figure 3 | Coulomb-field control of NV optical resonances.(a) Confocal image (532 nm excitation) of a second cluster comprising four different NVs.Superimposed circles indicate the in-plane positions as determined from confocal imaging using laser light resonant with the  !transition of each NV in the set; the radius in each circle indicates the uncertainty.The lower right insert is a zoomed-out image of the same cluster.(b) Recursive optical spectroscopy of the NV cluster in (a).Individual inspection of the spectra in the series reveals a correlation between NVG ionization (heralded by a missing optical resonance, yellow arrows) and the shift to longer wavelengths in the optical resonances of NVF (orange arrows).All conditions as in Fig. 1b.(c) Correlation spectroscopy between the longitudinal and transverse fields at NVG as derived from the results in (b).Red and blue dots respectively denote instances where NVD is negative or neutral; ellipses are guides to the eye.(d) Starting from a configuration where both NVF and NVG are negatively charged (blue trace), we probe the cluster response after resonant ionization of NVG (red trace).(e) Same as in (d) but for NVF ionization.The reference frequency in (b) and (d) is 470.000THz.

Figure 4 |
Figure 4 | 3D co-localization via multi-NV cross-correlated noise spectroscopy.(a) Correlations across NVD and NVF as derived from Fig. 3b; a charge state change of NVG (red or blue dots) leads to a shift in the mean electric field projections at NVD; ellipses are guides to the eye.(b) (Right) Histogram of the electric field as seen by NVD; blue (magenta) dots indicate instances where NVG (NVF) is neutral, while red indicates no ionization.Integers ,  highlight examples corresponding to time consecutive spectra.(Left) Same as before but calculated as d-field differences  =  "#$ −  " between two successive measurements; the suppression of slower temporal fluctuations leads to higher spectral resolution.(c) Co-localization of NVs in the cluster of Fig. 3.The magenta (green) loop represents the probability distributions for NVF (NVG) as seen by NVG (NVF).Similarly, the purple (cyan) areas represent probability distributions for NVD using NVG (NVF) as the probe; confidence intervals are 70% for NVD and 95% for the rest.(d) Colocalization of dark trap DTA, proximal to NVF and NVG; upper and lower loops show the 90% probability distributions for the position of the trap as seen by NVG and NVF, respectively.(e) (Main) Threedimensional spatial locations of NVD, NVF, NVG, and proximal traps DTA, DTB; confidence intervals from dark to light are 50% and 90%.(Upper inserts) Spatial locations of NVD, NVF, NVG projected on the xy optical plane as extracted from the 3D plot.The dashed circles show the NVs location measured from super-resolution (95% CI); all locations are relative to NVG.

Supplementary Figure 4 :
Automated analysis of a spectral series.(a) Spectral series for the cluster in Fig 3 of the main text (low frequency end only) illustrating some of the "bands" we use in our code during automatic searches of a given resonance.Colored boxes of indicate different NVs in the cluster.Dotted rectangles indicate secondary spectral bands in NVF chosen so as to take into account the shift introduced by sporadic ionization of NVG.Black horizontal lines every 3 spectra denote resets of the NV charge state with green light.The resonant laser power is 3 nW and the PL acquisition time per frequency step is 2 s; the number of points in each spectrum is 700 amounting to an acquisition time of 25 minutes per scan.Frequency shifts are referenced to 470.470 THz.(b) Cumulative distribution function (CDF) upon a Gaussian fit of  %

Supplementary Figure 5 :
Illustrating the impact of an underlying bias field.(a) Effect of the bias field (red arrow) on the measurement of the transverse field produced by an added charge.Since a measurement determines only the magnitude  ) (blue circle), different fields (e.g., purple and orange arrows) can both lead to the same end observation despite their different magnitude and orientation (blue and light-blue arrows).(b) Effect on the localization of a discrete charge trap.The purple and orange circles are the positions corresponding to the purple and orange arrows in (a), matched by color.The purple (black) line shows all possible positions for a negative charge to generate an end transverse field with the magnitude shown in (a) without (with) an identical change of longitudinal field (300 MHz).In this schematic, the green dot represents the NV at the origin, the red arrow points in the direction of the bias transverse field, and the blue arrow represents the direction of the NV crystal axis.position of the trap following a Bayesian approach.Let us assume a single charge is within a certain volume  of the NV center.A measurement  (with associated variance    ) has been carried out, which informs us on the charge position.Using Bayes theorem, we then write (|)() = ()(|).

2 -
with standard deviation  B K,L , and the sum extends over a discrete set of  positions along the solution loop in 3D space.We also calculate the electric field   u, ,  B K,L v that NV K at position   , would exert on NV L .For each probe NV L , and for every pair of source NVs, NV K and NV O , , ,  = D, E, F, G, we iterate through all Gaussian distributions of   K,L and   + O,L describing the positions of NV K and NV O to calculate the transverse electric fields created by only one NV,  P,! u, ,  B K,L v and  P,! u, ,  B + O,L v, and by the sum of both,  ,! .Since

For every degeneracy ofK
parameters  =  B K,L ,  B + O,L ,  Y K,L we calculate  B,B + ,Y K,L = | B K,L −  Z , -,/ ( B + L,K )|.u B the change in transverse strain axis between the two NVs.All remaining pairs of Gaussian distributions are multiplied, i.e., (). B + Z L,K (−), which results in a new normal Gaussian distribution of lower amplitude (with  B,B + ,Z the renormalization factor); note this distribution has, in general, a different center and standard deviation, which we denote  B,B + ,Z K,L and  B,B + ,Z K,L .