{'fraction_non_alphanumeric': 0.03526096320395913, 'fraction_numerical': 0.01878751775649544, 'mean_word_length': 4.967596390484004, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Template-Assisted Self Assembly of Fluorescent Nanodiamonds for Scalable Quantum Technologies

Milled nanodiamonds containing nitrogen-vacancy (NV) centers provide an excellent platform for sensing applications as they are optically robust, have nanoscale quantum sensitivity, and form colloidal dispersions which enable bottom-up assembly techniques for device integration. However, variations in their size, shape, and surface chemistry limit the ability to position individual nanodiamonds and statistically study properties that affect their optical and quantum characteristics. Here, we present a scalable strategy to form ordered arrays of nanodiamonds using capillary-driven, template-assisted self assembly. This method enables the precise spatial arrangement of isolated nanodiamonds with diameters below 50 nm across millimeter-scale areas. Measurements of over 200 assembled nanodiamonds yield a statistical understanding of their structural, optical, and quantum properties. The NV centers' spin and charge properties are uncorrelated with nanodiamond size, but rather are consistent with heterogeneity in their nanoscale environment. This flexible assembly method, together with improved understanding of the material, will enable the integration of nanodiamonds into future quantum photonic and electronic devices.


I. INTRODUCTION
The nitrogen-vacancy (NV) center in diamond functions as an optically-addressable spin qubit with roomtemperature spin coherence and sensitivity to environmental perturbations such as magnetic and electric fields, strain, temperature, and pressure [1][2][3][4]. Although the best spin coherence properties are achieved in bulk diamond crystals [5,6], nanodiamond particles with diameters < 100 nm still shield the spin qubit from decoherence [7], and unlike bulk diamond, they can be dispersed in solvents to form colloids and placed within nanoscale distances of biological or inorganic structures. As an optically robust, biocompatible, quantum sensitive colloid, nanodiamonds provide a platform for sub-diffraction limited imaging, in vivo and in vitro sensing, and integration with complex heterojunction devices [3,8,9].
Two methods are typically used to prepare nanodiamonds: detonation synthesis and milling bulk crystals. Detonation nanodiamonds are useful as dyes and biomarkers; they offer good size uniformity, but the outer layer of sp 2 carbon and high background impurity levels associated with their synthesis inhibit the formation and stability of NV centers [3,10]. On the other hand, milled nanodiamonds inherit the chemical purity of their source crystal, at the cost of uniformity [11]. Milled nanodiamonds vary widely in size, shape, surface chemistry, and number of NV centers [3]. These chemical and structural variations are reported to give rise to inhomogeneity in the nanodiamonds' optical and quantum properties [7,12,13]. Statistical studies of individual, isolated nanodiamonds are therefore necessary to correlate their structure and properties and ultimately to optimize nanodiamond quantum devices.
The variation in nanodiamond size typically exceeds 10-50% depending on the milling parameters [12][13][14]. This variation, in addition to heterogeneity in shape and surface chemistry hinders the placement or assembly of individual particles and their incorporation into devices. In past studies, researchers have used atomic force microscopes (AFMs) to serially position individual nanodiamonds [15][16][17][18], or chosen to forgo ordered assemblies and search for randomly dispersed particles [12,13,19,20]. Self-assembly techniques have been employed to assemble clusters of several nanodiamonds [21]. Nanodiamond placement has been achieved through electrostatic interactions between the particles and specific substrates [22]. These approaches however place limits on the type of templates that can be used and hinder large statistical studies of individual nanodiamonds.
Here, we employ capillary-driven, template-assisted, self assembly (TASA) to form millimeter-scale area arrays of individual, milled nanodiamonds. TASA combines the flexibility and precision of nanolithography with the scalability of colloidal assembly. It has been used to position individual nanoparticles, as well as to direct the orientation of anisotropic particles, of various materials [23][24][25][26][27][28][29][30]. By using these arrays and an automated optical measurement system, we study the statistical het-erogeneity of the nanodiamonds' structural, optical, and spin properties. Our studies uncover statistical variations in nanodiamond properties correlated with characteristics of the underlying diamond material, the nanodiamond creation process, and the NV centers' local environment. This understanding will inform the realization of improved nanodiamond materials, and the versatile TASA method can be utilized in fabricating devices for quantum sensing and quantum photonics.

A. Assembly Process
We utilize commercially-available, low-fluorescence, milled nanodiamonds from Adámas Nanotechnologies. The nanodiamonds are specified to contain an average of 1-4 NV centers (with approximately 13% containing a single NV center), have an average size of ∼ 40 nm, and have surface carboxylate anions that enable their stable colloidal dispersion [14,31]. Figure 1 depicts the TASA process. Electron-beam lithography is used to fabricate templates with cylindrical traps, ranging in design diameter from 32-200 nm and having heights of 62 nm, in poly-(methyl methacrylate) (PMMA) thin films deposited on 250 nm SiO 2 on Si substrates. Larger trap diameters (> 50 nm) result in multiple particle assembly (Supporting Information Figure  S1) and smaller templates (< 35 nm) show limited particle assembly. Here, we focus on 35-48 nm trap diameters with a 2.6 µm pitch to simplify single-nanodiamond spectroscopic and morphological characterization, and we fabricate 26×26 nanodiamond arrays.
To assemble the nanodiamonds, we use a home-built apparatus in which an aqueous nanodiamond dispersion is deposited between a glass slide and the template surface ( Figure 1). A motororized stage translates the nanodiamond dispersion at a velocity of 3.5 µm/s across the template surface. An accumulation of particles forms FIG. 1. The TASA method. Schematic of the capillarydriven, template-assisted self assembly (TASA) of milled, fluorescent nanodiamonds into lithographically defined, PMMA templates on a Si/SiO2 substrate.

FIG. 2.
Assembly yield and size distribution. (A) Stitched AFM images of an assembly based on a 26×26 array of 42-nm-diameter trap sites. (B) Higher resolution AFM image of the 5×5 subarray indicated by a white box in (A). The lateral dimensions of the individual nanodiamonds are enlarged due to convolution with the AFM probe tip. (C) Assembly yield vs. trap diameter. The box plot statistics are calculated by dividing the array into 5×5 nanodiamond subarrays to account for regional variations. (D) Histogram of AFM-measured, assembled nanodiamond heights over all templates (258 individual nanodiamonds) and a DLS measurement (red line) of the parent nanodiamond dispersion. (E) Nanodiamond size distribution as a function of template trap diameter. Outliers (red crosses) in (C) and (E) indicate measurements whose values are >1.5 times the interquartile range from the bottom or top of the box.
at the meniscus and capillary forces drive the nanodiamonds into the trap sites [23,[26][27][28]32]. The assembly apparatus is enclosed in an environmental chamber to maintain a dew point of 8.6-9.5 • C. A recirculating chiller and heater is used to pump water through a copper block that serves as a sample holder and maintains a substrate temperature of 24-25 • C. Drier chamber conditions yield assemblies with excess nanodiamond deposition, overfilling the trap sites. A wetter chamber results in limited assembly. After assembly, the PMMA template is removed by sequential immersion of the samples in N-Methyl-2-pyrrolidone (NMP) and acetone baths followed by an isopropanol wash.
We employ AFM measurements to characterize the yield of our large-area nanodiamond assemblies. The piezo-driven stage on the AFM has a maximum 40 µm × 40 µm scan area. We stitch together sequential scans to view the full arrays [33]. Figure 2A shows a representative AFM image of a TASA assembly of nanodiamonds in a 26×26 array. Although assembly yield is high (see below), we still observe instances of multi-particle and noparticle assembly. These imperfections are expected in part due to the polydispersity in size, shape, and surface chemistry of the nanodiamonds which can both impede assembly and lead to particle agglomeration. We also hy-pothesize that regional variations in assembly arise from fluctuations in chamber conditions, which alter the capillary forces, and incomplete liftoff, which may fail to remove larger nanodiamond deposits. Dividing the larger array ( Figure 2A) into subarrays, as shown in Figure 2B, accounts for these regional variations and allows for statistical characterization of the assembly yield.
We repeat the assembly and AFM characterization for five different template trap diameters (Supporting Information Figure S2). Collectively, we observe a median yield of single nanodiamonds of 76%, with individual regions in most arrays reaching 100% yield ( Figure 2B). Figure 2B represents a lower bound on assembly yield, since higher resolution AFM scans occasionally reveal previously unseen nanodiamonds. No significant statistical variation in yield is observed for trap diameters ranging from 38-48 nm, however we see a ∼15% drop in yield for the 35 nm trap diameter template compared to those for the larger traps.

B. Nanodiamond Size Distribution
In addition to detecting the presence of assembled nanodiamonds for yield estimates, AFM measurements quantify the nanodiamond heights. Tip convolution prevents reliable measurement of the lateral dimensions of individual nanodiamonds. The distribution of particle heights (histogram, Figure 2D), with average size and standard deviation of 39 ± 12 nm, is consistent with dynamic light scattering (DLS) measurements (red curve, Figure 2D; 36 ± 12 nm) and with the manufacturer's specifications [14,31]. Although the nanodiamonds are irregular in shape, we posit that the AFM height serves as a suitable measure for individual nanodiamond size.
Each template design yields arrays of nanodiamonds with a similar median size ( Figure 2E). The smallest 35 nm template displays a tighter particle size range, while the larger templates have extended tails, especially towards larger particle sizes. The observation of particles larger in height than the trap diameter is consistent with the likelihood of irregularly-shaped nanodiamonds assembling with their long axis perpendicular to the substrate surface.

C. Optical Characterization
We use a custom-built, automated, confocal microscope (see Methods) to probe the optical and quantum properties of the nanodiamond arrays. We stitch together multiple photoluminescence (PL) images to map the fluorescence of nanodiamonds in large arrays ( Figure  3A). Using an automated procedure, we characterize each fluorescent nanodiamond within a subarray ( Figure 3B and Supporting Information Figure S3 and S4), recording measurements of the particles' PL spectrum, photon autocorrelation function, PL saturation as a function of excitation power, and spin lifetime. Figures 3C-F depict the results of these measurements for four sites (I-IV) that exemplify the variation of observed optical properties. All optical measurements are performed using 532 nm excitation. The excitation power is 0.425 mW, measured before the objective lens for all measurements except for saturation where the power varies. AFM scans confirm that each of the four sites contains a single fluorescent nanodiamond (see Figure 2C and Supporting Information Figure S3). Figure 3C shows single nanodiamond PL spectra. Wavelengths below 550 nm are cut off by the long pass filter used to exclude the 532 nm excitation. The NV center's PL spectrum consists of a linear combination of spectra associated with the NV 0 and NV − charge states. Under optical excitation, the NV center cycles between these two charge states through processes of ionization and recombination [34], resulting in different spectral weightings as shown by the examples of nanodiamonds I-IV. We use a nonnegative matrix factorization method to decompose the spectra of individual nanodiamonds into their corresponding N V 0 (blue) and N V − (red) charge components [13]. From these decomposed spectra we calculate the charge ratio, is the integrated spectral intensity associated with the NV − (NV 0 ) charge state. We graphically represent F NV − using a blue (predominantly NV 0 ) to red (predominantly NV − ) color scale throughout Figure 3. Numeric results for nanodiamonds I-IV can be found in Table I. Observations of mixed charge states, 0 < F NV − < 1, can result from dynamic ionization and recombination of individual NV centers or from ensemble averaging over multiple NV centers within a nanodiamond.
We measure the photon autocorrelation function, g (2) (τ ), which probes the likelihood of two fluorescence photons being detected at varying temporal delays, τ [35,36]. Figure 3D shows examples of g (2) (τ ) for individual nanodiamonds. g (2) (τ ) is corrected for background fluorescence, detector dark counts, and detector timing jitter. The zero-delay value, g (2) (0), is used to quantify the number of NV centers in each nanodiamond. Specifically, we calculate the effective number (N Eff ) of NV centers using the equation: The progressively shallower dips seen in nanodiamonds I-IV for small delays, where g (2) (τ ) < 1 as τ → 0, indicate quantum fluorescence from a progressively larger number of discrete emitters. Nanodiamond I exhibits g (2) (0) < 0.5 which indicates that the emission is dominated by a single NV center. This relationship assumes identical PL brightness from each emitter, which is not strictly accurate for nanodiamonds given their morphological inhomogeneities, varying crystal orientations, and different charge state populations. Furthermore, measurements of g (2) (0) can be biased by systematic errors such as incomplete background correction, particularly for dim emitters (typical of small N Eff , e.g., for nanodiamond I) or when g (2) (0) → 1 (as N Eff → ∞, e.g., for nanodiamond IV). As described further in the Statistical Characterization Section, we find that N Eff is linearly related to the saturation PL brightness of the nanodiamonds ( Figure 3E), especially for the range 1 < N Eff < 10 and we therefore interpret N Eff as a reasonable quantitative estimate for the number of NV centers at each location.
While each nanodiamond varies in brightness ( Figure  3B), by varying the illumination power, P , and recording the PL intensity, I PL , we find that the nanodiamond I PL displays a similar characteristic power dependence (Figure 3E). As expected, the observed brightness increases with g (2) (0) and the number of emitters, N Eff . We perform this PL saturation measurement multiple times to check for stability, subtract a background measurement from a location on the sample nearby each nanodiamond, and fit the resulting data to the empirical model, from which we quantify the nanodiamond's saturation PL brightness (C Sat ) and saturation power (P Sat ). Finally, to probe the quantum features of the nanodiamonds, we measure the electron spin lifetime, T 1 ( Figure  3F), using a time-resolved fluorescence technique. We initialize the NV centers into their m s = 0 spin-triplet sublevel using a 1.6 µs, 532 nm laser excitation pulse and then count PL photons during a 300 ns readout laser pulse at delay times ranging from 50 ns to 900 µs. Due to the NV center's spin-dependent optical dynamics, the PL amplitude represents the probability that the NV centers remain in the m s = 0 spin sublevel of the N V − ground state at the time of the readout pulse. Longer T 1 lifetimes indicate higher levels of quantum isolation and enable improved quantum sensitivity [3,7]. Variations at short ( 200 ns) times can result from relaxation of the metastable spin singlet, whereas longer-time variations can reflect charge or spin relaxation [3,7,20]. We fit the data using a set of mathematical functions that account empirically for the rates of these disparate processes as well as for the potential presence of multiple NV centers within a nanodiamond (see Methods). In some cases, such as nanodiamond II, the analysis is inconclusive; in these cases we do not report a T 1 lifetime. Nanodiamond II's small PL contrast is consistent with its low F NV − , since spin relaxation is only observed for NV centers in the negative charge state.

D. Statistical Characterization
Ordered arrays of nanodiamonds significantly simplify experiments, allowing us to study the statistical properties of large numbers of nanodiamonds, correlating optical and structural measurement modalities. We collect AFM scans and the suite of optical measurements, as described for the example nanodiamonds in Figure 3, but now for 219 nanodiamonds assembled in multiple arrays, to build statistical data sets of nanodiamond size and emitter number, brightness, charge state, and spin lifetime ( Figure 4). The Supporting Information also includes an analysis of correlations between these variables (Supporting Information Figures S5-S9).

Emitter Number Distribution
We study the number of NV centers hosted within each nanodiamod ( Figure 4A). A significant fraction (68 out of 219, or 31%) of nanodiamonds are observed in AFM measurements but invisible in optical measurements. We set N Eff = 0 for these non-fluorescent nanodiamonds, as indicated by the purple bar in Figure 4A. The nonfluorescent nanodiamonds are excluded from the other datasets in Figure 4. For the fluorescent nanodiamonds, we calculate N Eff from autocorrelation measurements of g (2) (0), as in Figure 3D, using eq. (1). Throughout Figure 4, we indicate the nanodiamonds likely to contain one (0 < N Eff < 1.5) or two (1.5 < N Eff < 2.5) NV centers using separately colored bars. Of the 151 fluorescent nanodiamonds, four are too dim and sixteen are too bright to reliably determine N Eff from g (2) (0); we estimate N Eff for these nanodiamonds using the observed linear relationship between N Eff and C Sat , as described below (Figure 4B, inset) and add the results to Figure 4A as grey bars.
The observed fractions of non-fluorescent nanodiamonds (31%) and single emitters (12%) in this sample are in close agreement with the manufacturer's product specifications (25% and 13%, respectively) [31]. However, we show below that the observed N Eff distribution in Figure 4A is not consistent with a simple stochastic model assuming an isotropic distribution of NV centers within the nanodiamonds. We develop a quantitative model for the observed distribution, motivated by the mechanisms of diamond growth, irradiation, and milling.
A model for the N Eff distribution must account for the variance in particle size and the stochastic incorporation of NV centers within each particle. First, we consider a model in which the probability of observing N emitters in a spherical nanodiamond of diameter d obeys a Poissonian distribution, P (N |d) = Poiss (N ; N ), where the mean, N = ρV , is determined by the NV density, ρ, and the particle volume, V = πd 3 /6. The probability distribution for N is subsequently given by is the probability of finding a particle with diameter d, as characterized by AFM measurements (Figure 2D). This model fails to match the observed distribution, with a marked discrepancy especially at low N Eff . (The best fit is rejected on statistical grounds with reduced chi-squared (χ 2 r ) = 3.05, p < 10 −10 , as shown in Supporting Information Figure S10.) However, when we exclude the non-fluorescent nanodiamonds, this model can explain the distribution of remaining fluorescent particles (χ 2 r = 1.2; p = 0.29). We conclude that the relatively large fraction of non-fluorescent particles in our data cannot be explained using this simple model.
To reconcile this discrepancy, we hypothesize that the nanodiamonds can be divided into two groups: one where the NV centers are distributed isotropically with nonnegligible density, resulting in the observed N Eff distribution of fluorescent nanodiamonds, and another group where the density is so low that there is a negligible probability of finding an NV center. Nitrogen is incorporated inhomogeneously during high-pressure, hightemperature growth of bulk diamond crystals [37][38][39]. When the crystals are milled, the resulting nanodiamonds have an inhomogeneous nitrogen concentration. Subsequent electron irradiation creates a homogeneous distribution of carbon vacancies, which become mobile upon annealing to form NV centers. In regions with high nitrogen content, the NV-center formation will be vacancy limited and hence the NV-center density is proportional to the radiation dose. In nitrogen-depleted regions, the NV-center formation will be nitrogen limited, and the density will be lower. This hypothesis is supported by the observation that Adámas Nanotechnology produces "highly-fluorescent" nanodiamonds with the same size distribution, using an identical synthesis and milling pro-cedure but significantly higher irradiation levels. These samples also contain a 30% fraction of non-fluorescent particles. Furthermore, the height distribution of nonfluorescent particles ( Figure 4A, inset) is similar to the full distribution ( Figure 2D), whereas an exclusively isotropic NV-center creation model would predict a narrowed distribution of primarily non-fluorescent smaller particles. The invariance of the non-fluorescent fraction with respect to milling and irradiation conditions supports the conclusion that the heterogeneity in NV-center creation derives from a nonuniform nitrogen distribution.
To account for this effect in our model, we adopt a bimodal distribution of NV-center density, where ρ = ρ 1 with probability α, and ρ = ρ 2 with probability 1 − α; see Methods, Equation (11). A fit to the data using this model yields an acceptable result as shown by the curve in Figure 4A (χ 2 r = 1.3, p = 0.17), where ρ 1 = 14 ± 1.4 × 10 16 cm −3 , ρ 2 = 0.001 ± 0.45 × 10 16 cm −3 , and α = 0.77 ± 0.05. The implication of this analysis is that the observed 31% of dark nanodiamonds arise from a combination of the underlying material inhomogeneity (accounting for 23% of dark nanodiamonds in the sample) and the stochastic creation of NV centers (accounting for the remaining 8% of dark nanodiamonds). Accounting for the expected ranges of NV-center creation efficiency, these fits place lower (upper) limits on the nitrogen content in the nitrogen-rich (nitrogen-poor) particles of 7.1 ppm (2.5 ppm). See the Supporting Information section III for more details on the nitrogen content and NV-center creation estimates.

Nanodiamond Brightness
The empirical saturation brightness distribution is shown in Figure 4B. We observe a clear linear relationship between C Sat and N Eff (Figure 4B inset), with a slope of 22.8 ± 2.2 kCts/s per NV center. The linear relationship implies that NV centers in these nanodiamonds have nearly identical apparent brightness, despite wide variations in charge stability and spin stability, as discussed in the following sections. The strong linear relationship between C Sat and N Eff implies that more timeefficient measurements of C Sat in comparison to g (2) (0) can serve as a robust indicator for N Eff . We use the relationship to estimate N Eff values for the 20 nanodiamonds with an ambiguous g (2) (0) in our dataset (Figure 4, grey bars and Supporting Information Figure S11-S12). The distribution for P Sat can be found in Supporting Information Figure S13. Figure 4C shows the F NV − distribution for fluorescent nanodiamonds in our assemblies. The mean value ( F NV − = 0.6 ± 0.1; see Table I) agrees with previous measurements of randomly dispersed, ∼150-nm-diameter, milled nanodiamonds [13]. However, it is significantly lower than the measurement of F NV − = 0.75 for the parent nanodiamond dispersion used in assembly (Supporting Information Figure S14). The dispersion and individual nanoparticle measurements are performed at similar optical excitation power densities, however we note that dispersed particles can exchange freely in and out of the excitation path during the measurement, whereas assembled, individual particles are continuously probed. Although ionization and recombination rates depend on excitation power density, the ratio F NV − generally does not vary strongly with power under 532 nm excitation [34,40,41]. We did not observe a power dependence in the qualitative shape of the PL spectra of individual nanodiamonds for excitation powers ranging from 0.1 to 1.6 mW before the objective.

NV-Center Charge-State Distribution
Accordingly, the nanoparticles' environment may play the dominant role in determining F NV − . A bias towards the negative charge state is often explained by a combination of core nitrogen impurities, which are reported to donate electrons to NV centers [41,42], together with the presence of carboxylate anions on the nanodiamond surface [13,43,44]. Once assembled, we hypothesize that silanol groups or water of hydration on the SiO 2 surface, known to act as electron traps [45][46][47], can neutralize the negative surface charge of the nanodiamonds and thereby decrease F NV − .
A closer examination of the F NV − distribution of individual nanodiamonds containing small numbers of NV centers reveals a bimodal shape ( Figure 4C). Nanodiamonds with one or two NV centers exhibit F NV − values clustered either around 0.4 or around 0.9. A few particles exhibit F NV − > 0.95, essentially a pure negative charge-state configuration. The bimodal distribution for nanodiamonds with N Eff ∼1-2 suggests the NV centers fall into two distinct categories based on F NV − . The more mixed charge states (F NV − ∼ 0.6) seen in particles with larger N Eff therefore arise from averaging the signals from multiple NV centers. The bimodal distribution could reflect distinct local environments, for example, distinguishing NV centers near the surface from those within the nanodiamond core. We checked for a correlation between F NV − and particle size, which would support this interpretation, but any such correlation is not significant in our data (see Supporting Information Figures S9-10). However, if environmental effects are limited to NV centers within a few nanometers from a particle's surface, we would not expect to observe a strong correlation for this range of particles sizes. To quantify this relationship, we compared the fraction of nanodiamonds with F NV − < 0.5 with the predicted shell-volume fraction as a function of particle size (see Supporting Information Figure S15). We find that our data are consistent with this interpretation, where NV centers located within 2-3 nm of the surface have reduced F NV − , however even larger sample sizes would be required to establish a direct correlation between F NV − and particle size.
The observation of a bimodal charge-state distribution has important practical implications. The pure negative charge state is typically desired for quantum spin sensing applications, whereas mixed charge states may be more sensitive to variations in electrochemical potential [48]. Hence, deeper understanding of this phenomena that leads to manufacturers' ability to promote one category over the other, or to efficiently sort the particles based on F NV − , would improve the material's performance in these key applications.

Spin Lifetimes
Spin relaxation measurements are performed as described in the Optical Characterization section, and those nanodiamonds whose measurements are adequately fit by our models are included in Figure 4D; we are unable to resolve T 1 values for 48 nanodiamonds. For the remaining nanodiamonds, we observe a wide distribution of T 1 values, ranging from 1.0 ± 0.5 µs to 608.0 ± 211.5 µs; see Figure 4D. We also observe a positive correlation between T 1 and F NV − (see Supporting Information Figures S5 and  S9). This correlation is consistent with the interpretation proposed in the previous section that NV centers close to the nanodiamond's surface display higher environmental sensitivity, as seen by lower F NV − and T 1 , according to a core-shell relationship. The lack of a correlation with height, which has been observed for studies considering much larger size ranges [7,12], further stresses the importance of NV center placement rather than particle size for these nanodiamonds. However, other interpretations are possible. For example, proximal fluctuating charges that lower F NV − for an NV center in the nanodiamond core could also limit T 1 . Charge fluctuations have been observed for NV centers in nanodiamonds in the dark following a laser pulse [20]; such fluctuations could systematically bias the T 1 measurements.
As for the case of charge-state stability, understanding and controlling the spin lifetime has important implications for applications. Particularly for spin relaxometry sensing and imaging, an increased T 1 translates directly to improved sensitivity [2,3,7,20]. We note that several of the longest T 1 values observed in our sample occurred in nanodimaonds with only one or two NV centers, boding well for the prospects of single-spin applications.

III. CONCLUSION
We describe the use of the TASA method to assemble large arrays of single nanodiamonds with high yield. We exploit these arrays to perform automated measurements of the morphological, optical, and quantum properties of nanodiamonds from 219 isolated particles. We observe a wide distribution of emitters in the particles, ranging from zero to >20 NV centers per particle, with 31% dark nanodiamonds and 12% single NV centers. We attribute the distribution to spatial heterogeneity in the nitrogen incorporation during growth of the parent, bulk crystal as well as stochastic creation of NV centers within nanodiamonds. Greater control over the nitrogen incorporation rates during crystal growth, vacancy generation during electron irradiation, and the particle size uniformity could provide tighter control over the NV center distribution. The fluorescent nanodiamonds show a wide distribution in their charge states and spin lifetimes. We find that F NV − and T 1 are not correlated with nanodiamond size in this sample, but instead indicate a heterogeneity in the NV-center environment, which could arise from variations in the internal crystal structure and composition or the NV center distance from the nanodiamond surface. Controlling the surface conditions should therefore enable more uniform charge states and longer spin lifetimes.
Building on the TASA method and leveraging our understating of the NV center's properties, we see a path for integrating nanodiamonds and other nanoscale quantum materials into photonic and electronic quantum devices. We envision extending the TASA method for the placement of individual nanodiamonds to that of other irregularly shaped nanoparticles. Templates can also be fabricated on functional substrates to better optimize the nanodiamonds' surface environment. We can adapt the TASA method to create multiparticle heterostructure assemblies [18,19,23], where plasmonic, magnetic, or dielectric nanoparticles could act as optical or magnetic antennas and enhance the NV-center's optical signal and quantum sensitivity.

Template Fabrication and Nanodiamond Assembly
The starting material for these studies is an aqueous suspension of 40 nm, milled, fluorescent nanodiamonds (Adámas Nanotechnology, NDNV40nmLw10ml).
We fabricate templates on single-side-polished Si wafers with a 250 nm SiO 2 thermal oxide layer (Virginia Semiconductor). Poly(methyl methacrylate) (PMMA) positive resist, 950 K molecular weight PMMA in anisole (MicroChem), is spun onto the substrates at 3000 rpm for 60 s followed by a soft bake at 180 o C for 90 s. Electronbeam lithography (Elionix ELS-7500EX) is used to pattern the template at a 50 kV accelerating voltage and a 50 pA beam current. Circular traps are designed with diameters ranging from 35-200 nm. The PMMA is developed in a 1:3 water:isopropanol solution for 90 s. Prior to assembly we perform a 6 s, 75 W, 20 sccm, O 2 plasma (Gatan Solarus) descum resulting in an approximately 62-nm-thick template.
The assembly is performed using a custom-built capillary assembly apparatus [23]. A 35 µL volume of a 0.001 mg/mL nanodiamond dispersion in a 1% Sodium dodecyl sulfate (Sigma-Aldrich) aqueous dispersion is dropped between the template and a glass slide. The substrate is translated by a motorized linear stage (New Focus Picomotor) at a rate of 3.5 µm/s. The assembly apparatus is enclosed in a home-built, humidity-controlled environmental chamber, and the substrate temperature is regulated by a water-cooled stage and monitored by a thermocouple. During assembly, the substrate temperature is 24 • C, and the ambient dew point is maintained between 8.5 and 9.5 • C. These temperature and humidity conditions maintain a contact angle >24 • throughout the assembly. Lift off is then performed by soaking the assembled samples in N-Methyl-2-pyrrolidone (NMP) for 1 min, followed by 5 min acetone and IPA washes.

AFM Measurements
Atomic force microscopy (AFM) measurements are taken using a MFP-3D-BIO AFM from Asylum Research. All measurements are collected using a standard Si probe (AC240TS-R3, Olympus). For quantifying yields, large-area scans (40 µm × 40 µm) are taken at a scan rate of 0.3 Hz. Mid-area scans (20 µm × 20 µm) for mapping subarrays are taken at a scan rate of 0.27 Hz and individual nanodiamond measurements are collected using a 1.5 µm × 1.5 µm scan area at a scan rate of 0.3 Hz. Igor Pro version 6.38B01 is used to analyze AFM data. Composite images are formed from the AFM scans using ImageJ and the program described by Preibisch et al. [33] Total yield counts and are analyzed by hand from the resulting images.

Optical Measurements
Optical measurements are taken using two custombuilt confocal microscopes. In both microscopes a fast steering mirror (FSM, Optics in motion, OIM101) projects a beam onto to the back of a 100x objective (0.9 NA 100x objective Olympus, MPI Plan Fluor/OL, or Nikon, Plan Flour x100/0.5-1.3), allowing us to raster a 532 nm (green) continuous-wave laser (Coherent, Compass 315M-150) over the sample and collect PL emission on single-photon avalanche diode detectors. We employ two microscopes to allow multiple measurements to run in parallel and to extend the spectral bandwidth. The first microscope, described by Exarhos et al. [49], is used to collect saturation and autocorrelation measurements. The second microscope, described by Huang et al. [50], is modified with a 532 nm low pass dichroic mirror (Semrock), to extend the spectral bandwidth and collect PL maps, spectra, and T 1 measurements. Throughout each measurement, a periodic tracking PL scan is used to maintain alignment for optimal collection from individual nanodiamonds. During tracking, PL line scans in the x and y directions (controlled by the FSM) are fit using a Gaussian peak plus a constant background. The FSM is then adjusted to align the microscope to the Gaussian's mean position. A similar tracking scan and Gaussian fit in the z (focus) direction is also performed utilizing a piezoelectric stage (Thorlabs MZS500-E) connected to either the sample mount or the objective, depending on the microscope, to optimize focus.
Spectral measurements are taken on a Princeton Instruments IsoPlane-160 spectrometer with a 700-nm blaze, 300 G/mm grating and a thermoelectrically cooled charge-coupled device (CCD) camera (Princeton Instruments PIXIS 100BX) yielding a spectral resolution of 0.7 nm. Measurements consist of one 90 s background measurement combined with multiple 90 s signal measurements. We then perform cosmic ray rejection and background subtraction on each scan before combining the signals and correcting for wavelength-dependent photon detection efficiencies.
The power-dependent PL measurements are achieved by placing a variable neutral density filter (Thorlabs NDC-50C-2-A) in the excitation path, with pre-objective powers calibrated with a Thorlabs PM100D power meter.
T 1 lifetime measurements are implemented by programming the required pulse sequences (described in the Optical Characterization subsection) onto an arbitrarywaveform generator (AWG; AWG520 Tektronix). The AWG control signals are passed to the microscope's acousto-optic modulator (AOM) for generating optical pulses, and to three high isolation switches (ZASWA-250DR Mini-Circuits) for time-gating photon detection events recorded by counters in a data acquisition card (National Instruments PCIe-6323).
DLS measurements are conducted on a Malvern Instruments Zetasizer Nano-s and analyzed using Malvern's software.

Autocorrelation analysis
Autocorrelation data is collected using a Hanbury Brown-Twiss setup with a time correlated single-photon counting module (PicoQuant PicoHarp 300) in timetagged, time-resolved collection mode. Background correction is performed by determining ρ = S/(S + B) from the signal (S) and background (B) levels for each spot (both determined from Gaussian fits during tracking scans) and then calculating the corrected autocorrelation function: The background corrected autocorrelation data are then fit using the empirical model: Here, A and D represent the bunching and antibunching amplitudes, respectively, t 0 accounts for signal delays, and τ 1 represents the antibunching timescale [36].
The signal-to-noise ratio of autocorrelation measurements depends on S, ρ, and on the acquisition time.
Therefore, we dynamically adjust the acquisition time for each nanodiamond in order to achieve a desired uncertainty in g (2) (τ ) (∆g), given a particular time resolution (∆τ ): For this study, we use ∆τ = 2 ns, and we set ∆g = 0.025 in order to achieve sufficiently high resolution, however we limit the total measurement time for individual spots to a maximum of 1 h. For most instances, this results in a resolvable g (2) (0) value. For some exceptionally bright or dim nanodiamonds, the antibunching dip is smaller than ∆g, and we cannot reliably determine N Eff values from g (2) (τ ).

Fitting Spin Relaxation Data
Recorded spin lifetime measurements for single NV centers are typically fit to a four-level rate model [7]: where T m accounts for the relaxation of the metastable singlet spin state. To account for the presence of multiple NV centers, we expand our empirical fitting model to account for a second spin relaxation rate: Due to signal-to-noise constraints and the potential overlap with T 1 , we may not always observe the T m term. Accordingly, we add two more potential models to describe the spin relaxation dynamics of the nanodiamonds.
We fit the data for each nanodiamond using all four models. We exclude fits where uncertainties for the T 1 parameter exceed the T 1 result, where parameters are inconsistent with known physical properties (T m > 10 µs or C 1 < 0), or where χ 2 r is greater than 1 +2 2 DoF . Here DoF is the number of degrees of freedom in the model. If no valid fits remain, we exclude that data point from later statistical analysis. In the event of multiple acceptable fits, χ 2 r is within 1 ± 2 DoF , we choose the model with the smallest number of parameters.

Modeling the Emitter Number Distribution
We fit the N Eff data by modeling the total probability of observing N emitters given a spherical nanodiamond of diameter d weighted by the likelihood of observing each d value: Here, the conditional probability of N given d is modeled as a Poissonian distribution: P (N |d) = Poiss (N ; N ), where the mean, N = ρV , is determined by the NV density, ρ, and the particle volume, V = πd 3 /6. The unconditional probability of finding a particle with diameter d, φ(d), is determined using AFM measurements as shown in Figure 2D. To account for a variation in ni-trogen content, we assume two different ρ values with their own probabilities α and 1 − α. The model of P(N ) then becomes: The results of these models can be found in Supporting Information Figure S10.
Supporting Information Figure 7. Plot of Height vsF NV − .

III. ESTIMATING NV CENTER DENSITY
Supporting Information Figure 10. N Eff histogram with three fits based on different models. The red curve represents a model with a constant NV density, ρ, and fits N Eff data for all nanodiamonds.
The magenta curve is a fit using the same model but only fitting data for fluorescent nanodiamonds.
The green curve assumes a bimodal distribution of ρ and utilizes all N Eff data.
We estimate a rough NV center density in the final nanodiamonds using the specifications provided by Adaḿas Nanotechnology. The nanodiamonds were irradiated with 2-3 MeV electrons to a dose of 10 18 cm −2 [1]. For this energy and dose, we assume a vacancy creation efficiency of order 1 cm −1 [2,3]. While NV center creation efficiency depends on annealing temperature, nitrogen content, and diamond size, recent studies employing similar synthesis conditions demonstrate conversions efficiencies on the order of 1-10% [4,5]. Based on these assumptions, we estimate a final NV-center density of ∼ 10 17 cm −3 , which is consistent with the value found in the vacancy limited region (Supporting Information Figure S10). The dispersion spectra are acquired by placing 300 µL of 1 mg/ml nanodiamond disperion loaded in a glass cuvette (SpectroCell R-403) in the confocal excitation path and illuminating with a 10x objective (0.3 NA 10X objective Nikon) at power of 4.7 mW before the objective. We perform a nonnegative matrix factorization analysis to compare this data with the spectral data of the assembled nanodiamonds.
To better understand the correlations between F NV − , T 1 and, size, we model the nanodiamonds as spherical particles with a core radius (R C ) and a shell thickness (S) such that R C + S is equal to the total radius of the sphere (R). The shell volume ratio (V S/R ) is then given by: This relationship is plotted in Supporting Information Figure S15 for a series of S and R values. We can compare these volume ratios to the distribution of observed F NV − and T 1 values, assuming that NV centers residing in the shell region are more affected by the environment then defects in the core region, and hence have a lower F NV − and smaller T 1 .
In this model, the likelihood of observing an NV center in the shell takes the form of a binomial distribution with a success probability V S/R . To perform this analysis, we sort the nanodiamonds by height and calculate the observed fraction of measurements that fall below a given threshold value for F NV − or T 1 . From this observed fraction and the number of measurements, we calculate the Wilson score interval as an estimate for the underlying binomial probability and corresponding uncertainty. In Supporting Information Figure S15, we plot the Wilson score intervals alongside the shell volume fraction for particular threshold values. Although the uncertainties are large, the results are consistent with the interpretation that both F NV − and T 1 are reduced for NV centers located within a few nanometers of the nanodiamond surface.