Confined nano-NMR spectroscopy using NV centers

Nano-NMR spectroscopy with nitrogen-vacancy centers holds the potential to provide high resolution spectra of minute samples. This is likely to have important implications for chemistry, medicine and pharmaceutical engineering. One of the main hurdles facing the technology is that diffusion of unpolarized liquid samples broadens the spectral lines thus limiting resolution. Experiments in the field are therefore impeded by the efforts involved in achieving high polarization of the sample which is a challenging endeavor. Here we examine a scenario where the liquid is confined to a small volume. We show that the confinement 'counteracts' the effect of diffusion, thus overcoming a major obstacle to the resolving abilities of the NV-NMR spectrometer.


INTRODUCTION
Nano Nuclear Magnetic Resonance (nano-NMR) spectroscopy is a new field that aimes to implement methods similar to classic NMR techniques on minute samples. The new technology may outperform microscopy methods and enable high resolution imaging at the single molecule level. Therefore, advances in the field will have major implications for fundamental science and pave the way for drug and medicine development.
Scaling down the sample size requires a different measurement device, because small samples generate magnetic fields that are undetectable by classical sensors. The Nitrogen-Vacancy (NV) center based spectrometer is a promising platform for nano-NMR, since the NV is a superb magnetometer of small fields [1,2]. Although many successful experiments using the NV spectrometer have been carried out in recent years [3][4][5][6][7], nearly all of them used a large sample size to generate a reasonable signal. The gain from using large volumes has to do with the dipolar interaction between the sensor and the nuclear spin ensemble. This creates an effective cut-off distance which is proportional to the NV's depth d, such that only nuclei within a radius of ∼ d from the surface contribute to the magnetic field. This limitation suggests that a resolution problem occurs when using shallow NVs since the diffusion time across the detection volume which sets the resolution capability is shorter (see Fig. 1a).
One way of avoiding the resolution problem is by using an ultra-polarized sample. While some experiments rely solely on statistical polarization, other types of polarizations; e.g. thermal polarization or hyperpolarized non-equilibrium states can be used to amplify the signal and circumvent the effects of diffusion [8,9]. It was also shown that given high polarization, the entanglement formed between the shallow sensors and the nuclear spin ensemble can be utilized to further increase precision [10]. However, despite considerable efforts to enhance the sample's polarization, current state-of-the art achieves about 1%. Going beyond this requires efficient hyperpolarization schemes which are technically hard to implement. Recent works suggest practical alternatives by showing that contrary to common assumptions, diffusion noise can be managed by implementing an appropriate measurement protocol. These methods rely on the long lasting correlation of magnetic field noise in liquid state nano-NMR with NV centers [11,12].
Another elegant and practical approach that is currently being explored by the community is to confine the nuclei in a small volume, thus limiting the effects of diffusion (see Fig. 1b). Here we provide a detailed theoretical analysis of confined nano-NMR spectroscopy. We calculate the expected average magnetic field and correlations in three different geometries: a cylinder, a hemisphere and a full sphere. We show that for any geometry of volume V the correlations decay to a constant value, resulting in unlimited resolution and an SNR which is a fraction of ∼ d 3 V from the SNR of the power spectrum at low frequency.
The advantage of confinement can easily be understood by the following argument. Denoting the liquid's diffusion coefficient by D, keep in mind that the system has two characteristic time scales: the time it takes a nucleus to diffuse out of the effective interaction region, τ D = d 2 D , and the time it takes a nucleus to reach the volume's boundary, τ V = V 2/3 D . Assuming that d < V 1/3 , we expect that for short times, t < τ D = d 2 D , when nuclei can rarely leave the interaction region, the magnetic field at the NV's position will be approximately constant. Therefore, the correlation is also constant. For intermidiate times τ D < t, as diffusion becomes dominant, the fluctuations in the magnetic field will cause the correlations to decay, while for τ D t τ V the correlation will decay as a power law [11]. However, for longer times t > τ V , nuclei are reflected from the boundary and therefore have a finite probability to return to the interaction region, which causes the correlation to decay to a constant value. The probability can be estimated by the ratio of the effective interaction region and the total volume ∼ d 3 V . The expected form of the correlation function is sketched in Fig. 2. In the following, we show in detail that these simple arguments indeed describe the correlation of confined arXiv:2006.12115v1 [quant-ph] 22 Jun 2020 (a) (b) Figure 1. (a) Nano-NMR using NV centers -the nuclear spin sample (red arrows) is placed on top of a diamond surface. The ensemble interacts with the NV center (blue arrow), which is situated at a distance d below the surface. The dioplar interaction creates an effective cutoff, such that only spins that are found in a hemisphere of radius ∼ d (blue dome) contribute to the magnetic field at the NV's position. While the generated field of the statistically polarized spins provides information about their Larmor frequency, after interrogation times t > τD ≡ d 2 D , the diffusion of spins outside the effective interaction region causes the field to fluctuate, which eventually leads to the decohernce of the sensor. This results in a resolution problem in nano-NMR using shallow sensors. (b) Confined nano-NMR -the nuclear spin ensemble is confined to a small space (gray cylinder). For interrogation times t < τV ≡ V 2/3 D the dynamics is similar to the non-confined scenario, since the interaction with the walls in negligible. For times t > τV , however, the nuclei have a finite probability of returning to the interaction volume, which leads to non-vanishing correlations of the magnetic field. These figures were created by Bar Soffer, Mirage studio.

Measurement protocol
We explore two measurement schemes: correlation spectroscopy [13][14][15][16][17] and a phase sensitive measurement (Qdyne) [8,18,19]. In the following we provide an outline of the two. We consider the NV center as a two level system (spin 1 2 ) with an energy gap of γ e B ext , where B ext is the externally applied magnetic field and γ e/n is the electron/nuclei gyromagntic ratio [2]. We denote S i (I i ) as the spin operator of the NV (nuclear spin) at the i direction, S ± (I ± ) as the NV's (nuclear spin's) raising/lowering operators and θ i as the angle between the NV's magnetization axis and the vector connecting the NV to the i-th nucleus.
In correlation spectroscopy, the NV interacts with nuclear spins for a period of time τ via the dipolar interaction with the ensemble at the end of which the NV state is written onto a memory qubit or alternatively the relevant information found in the relative phase is mapped to the population difference of the NV. After a waiting time t the experiment is carried out again for the same period of time τ . At the interrogation times the NV is driven with an external pulse sequence at frequency ω p , which is close to the Larmor frequency of the nuclei ω N , such that (ω N − ω p ) t ≡ δt 1. The Hamiltonian is effectively where φ is a random phase deriving from the unpolarized dynamics. Initializing the NV in the state |ψ 0 = 1 √ 2 (| ↑ + | ↓ ), after both interrogation periods, the probability of measuring the NV at the |ψ 0 state is where we took the weak back-action limit, namely Figure 2. A sketch of the expected behavior of the correlation function. Solid lines represent the approximate behavior in each regime, and dashed lines are extensions of the approximations to regimes where they are invalid. The value of the correlation at t = 0 is defined as B 2 rms . For times t τD the correlation decays exponentially, which is approximately linear for small arguments (blue line). As diffusion becomes dominant at times t > τD the exponential scaling is no longer valid (dashed blue line) and it changes to a power law ∝ t −3/2 (orange line). When the nuclei interact with the walls at times t > τV the correlation decays exponentially to γ e B(t)τ 1. Averaging (2) over realizations yields where C (1) (t) = B( t )B( 0 ) is the correlation function of the magnetic field at time t thus C (0) ≡ B 2 rms . The correlation is given by [11,20] (Ω 0 ) P (r,r 0 , t) , are the spherical harmonics, P (r,r 0 , t) is the diffusion propagator fromr tor 0 with time difference t andζ m are defined in [20]. For m = 1, for example, . For a detailed derivation of Eq.
(3) see [20]. We note that different values of m can be sampled by changing the dynamical decoupling pulse frequency. For example, m = 0 can be sensed by not applying any pulses (ω p = 0). In a phase sensitive measurement we employ a similar protocol, where the main difference stems from the fact that we measure multiple times in one realization; every period τ the NV is measured and then initialized at the state |ψ 0 . The evolution of the system between times t n ≡ nτ and t n+1 = (n + 1) τ for some integer n is where we assumed δτ 1. The probability of measuring with Φ (t n ) = 2γ e τ B (t n ) cos (δt n + φ). Assuming Φ (t n ) 1 we can approximate (6) by The joint probability is then When combined with the probability of both results being | ↓ y , an extra factor of 2 is gained. Averaging (8) over the different measurements gives This can be used to define a new Bernoulli variable by P = P |↑y (t n ) P |↑y (t k ) + P |↓y (t n ) P |↓y (t k ) (10) where the value of m in Eqs. (9) and (10) is determined by the external pulse frequency. In the following we present the calculation magnetic field correlations in three different geometries. We provide analytic calculations for the asymptotic behavior and numerical assessments for arbitrary times. Since the asymptotic behavior of the correlation functions is universal, we focus on the case m = 0. The calculation of the other correlation functions can be found in [20].

Cylinder
Given a statistically polarized ensemble which is confined to a cylinder with height L and radius R and an NV that is located at a depth d below the cylinder's base. The NV's position coincides with the cylinder's symmetry axis denoted asẑ. Assuming the NV's magnetization axis is along theẑ axis, the magnetic field correlation is given by (4), where the different m values can be probed by changing the pulse frequency mentioned in the previous section, see details in [11]. The instantaneous correlation, B 2 rms = C (0) (t = 0) can be calculated analytically as the propagator approaches a delta function, yielding: where we denote the physical coupling constant by J = hµ0γ N 4π . Eq. (11) reproduces the well-known scaling, , of the semi-infinite volume in the limit d R, L, as expected. In the other limit, d R, L, Eq. (11) can be approximated by B 2 rms ∝ J 2 V d 6 . This is because as the statistical polarization scales as √ V , the dipole-dipole interaction decays as 1 d 3 and B RM S scales with both.
The correlation (4) decays asymptotically to a constant as the propagator approaches the uniform distribution for times t τ V , For a sufficiently shallow sensor, d R, L Eq. (12) approaches The factor of d 3 V in this limit can be interpreted as the probability that a nucleus will be found in the effective interaction region, because at long times, diffusion spreads out the nuclei uniformly. Since the signal B 2 rms ∝ d −3 , remarkably, the correlation does not depend on the NV's depth. This regime (t τ V ), where the signal remains constant, is reminiscent of the situation of the unbounded region with a polarized sample, wherein the constant correlation is replaced by the average magnetic field.
At arbitrary times, (4) is given bỹ , J m is the mth Bessel function of the first kind, ν k,m is the kth zero of J m , δ i,j is Kronecker's delta and the characteristic decay time (4) might be very hard to estimate even numerically, Eq. (13) only contains a sum over two dimensional integrals (since the azimuthal integration is trivial) which are easy to calculate.
We evaluated (13) for k, s ∈ [1,25], the diffusion coefficient of immersion oil D = 0.5 nm 2 µs , d = 1 nm and different volumes [20]. In Fig 3 we present a representative case where R = L = 200 nm. The correlation at short times decays linearly, which agrees with an exponential decay for small arguments. In the intermediate regime the decays changes to a power law ∼ t −1.5 as predicted in previous works [11], while at long times the decay rate fits the slowest decaying mode of (13).
The numerical integration was done with a diffusion coefficient fitting immersion oil, since state-of-the-art experiments are conducted with high viscous fluids primarily for their long diffusion time, which results in a long correlation time and thus a high resolution capability. For water, for example, resolution is extremely limited, since the correlation vanishes rapidly. In a confined setting, however, the correlation decays rapidly to a constant value; therefore, diffusion poses no real limitation on the interrogation time thus eliminating the need to use high viscous fluids.

Hemisphere
The B rms of a hemisphere with radius R, given by the limit t → 0 of Eq. (4) is The limits of (14) are similar to those of the cylindrical case. The asymptotic behavior at long times is     V . In general the correlation can be expressed by the sum where j m is the mth spherical Bessel function of the first kind andν k,m is the kth zero that does not equal zero of j m (namely,ν 1,m = 0). The coordinatesr/r are the spherical coordinates where the origin is found at the center of the hemisphere's base at the NV's position. The results are similar to those of the cylindrical geometry and are found in [20].

Full sphere
Given a spherical geometry of radius R, the instantaneous correlation is The asymptotic behavior of (4) is then given by For arbitrary times, the correlation is We evaluated (19) for l ∈ [0, 30], k ∈ [1, 30], the diffusion coefficient of immersion oil D = 0.5 nm 2 µs , d = 1 nm and different volumes [20]. In Fig 4 we present a representative case where R = 200 nm. The correlation at short times and long times agrees with the results of the cylindrical geometry. In the intermediate regime the decay changes to a power law ∼ t −0.5 .

Sensitivity analysis
In our sensitivity analysis we used the tools of classic parameter estimation. We calculated the Fisher Information (FI), denoted by I; hence, the sensitivity can then be estimated with the Cramer-Rao bound [21] ∆δ We start by calculating the asymptotic behavior of the FI in either correlation spectroscopy or a phase sensitive    measurement. In the limits t ≤ τ D and t τ V the correlation is approximately a constant. Note that the regime t ≤ τ D is similar to the unconfined case. The FI of correlation spectroscopy is given by where T is the total time, and we used Eq. (3) for this derivation. Since C (t τ D ) ≈ B 2 rms we can deduce that for t ≈ τ D , Eq. (22) is the unconfined limit (UCL), since without confinement the measurement time is bounded by τ D .
where T m is the memory time. In the case where δT m > π 2 Eq. (23) can be simplified to Comparing Eq. (24) to Eq. (22) shows that we gained a factor of Since T m can be orders of magnitude larger than τ D and δτ D might be very small, this can lead to considerable enhancement in sensitivity. Numerical evaluations of Eq. (25) are found at the last section. In a Qdyne type measurement the probability (10) leads to For times t s τ V , the Fisher information (26) is Since the Fisher information is additive and the correlation is always larger than the asymptotic value we can bound the total Fisher information by where T is the total measurement time. The result (28) is a major improvement over (22), since it scales with T 3 which is not limited by the diffusion time. The bound (28) is quite crude, since for viscous fluids a great deal of information can be gained from the slow temporal decay. A comprehensive analysis can be found in [12].

Molecular dynamics
We corroborated the analytic prediction of the correlation function of the magnetic field as a result of the fluid motion in the closed geometry using molecular dynamics simulations, focusing on the case of the cylindrical geometry. The set-up is illustrated in figure 5. The system contains N particles, which are confined to a cylinder of radius R and height L. The particles interact via the Lennard-Jones (LJ) potential 4 [(σ/r) 12 − (σ/r) 6 ] with an interaction cutoff distance of r c = 2.5σ. Specular reflections are applied on the top and bottom walls of the cylinder and the Lennard-Jones 9/3 potential is applied at the curved walls to confine the particles to the interior of the cylinder. The system is initialized into a thermal state at temperature T by running a Langevin dynamics simulation until the system reaches thermal equilibrium. Each particle is assigned a random value of spin I z ∈ {−1, 1}. After initializing the system, we ran a deterministic molecular dynamics simulation, integrating Newton's laws using the Velocity-Verlet method with step size ∆t = 0.005 /(mσ 2 ). During the simulation we computed and stored the z component of the magnetic field induced by the particles at the position of the NV where r i is the distance between particle i and the NV center and θ i is the angle between r i and the NV quantization axis. The NV is placed along the axis of rotational symmetry of the cylinder, at a distance d below the bottom of the cylinder. Figure 6 shows the correlation function for the chosen cylindrical geometry at several values of d. At t = 0 the correlation decays as d −2.7 (shown at the inset), which is close to the d −3 scaling predicted by (11) at small values of d/R. The change in the decay is possibly due to the interaction between the particles, which is not taken into account in our model. At long times the the correlation decays to a constant value that increases with d/R as expected. A comparison between the asymptotic values of the MD simulation and the analytic prediction (12) is found at [20]. We compared the simulation results to the analytic correlation function (Eqs. (12) and (13)) by rescaling the asymptotic value of the former to match the one in Eq. (12) and choosing an appropriate diffusion coefficient. A representative case is shown in Fig. 7. We verified the results of the spherical geometry with a second set of simulations. The spherical geometry simulations follow a very similar protocol to the cylindrical geometry simulations. The number of particles was set to N = 28371 and the radius of the sphere was set to R = 20.47σ resulting in particle density of ρ = 0.79σ −3 . The temperature, integration timestep and the total integration time was identical to the cylindrical geometry simulation. The LJ 9/3 potential is applied across the whole surface of the sphere to confine the A comparison between the simulation results and the analytic predic-

Limits of the model
In the calculation of the correlation function we performed an average over realizations of the magnetic field noise. While this averaging holds for correlation spectroscopy, in a phase sensitive measurement which is based on a single realization it does not. The main difference is that in a single realization, the statistical polarization will cause the amplitude of the signal to fluctuate as different nuclei diffuse in and out of the effective interaction region. Assuming the phase of the signal remains constant, the FI is bounded by that of the average amplitude, which is given by the correlation function. As long as the spins retain their initial polarization, the phase will remain stable. The spin orientation will eventually fluctuate, mostly due to the interaction between the nuclei and walls. The total measurement time T in (28) is therefore limited by the characteristic time of the interaction, T walls . CONCLUSION We provided a theoretical basis for confined nano-NMR spectroscopy based on NV centers. We showed that the confinement indeed mitigates the deleterious effects of diffusion and therefore can be used to enhance the resolution abilities of both correlation spectroscopy and phase sensitive measurements. As advances are made with hyperpolarization [9,[22][23][24][25][26][27], it might prove more beneficial to combine polarization and confinement. The relevant calculations can be found in [20].  We start with the dipole-dipole interaction Hamiltonian: where we denote J = µ0γ N 4π for short. We write the inner product explicitly: The product of the two z components will yield (up to a minus sign): Where Y m l (Ω) is a spherical harmonic. This term corresponds to m = 0. The product of the other same axis terms will yield: We denote: Which is the second m = 0 term. The mixed XY terms will have similar contributions as the rest of the previous product: Summing both contributions we can define: These elements correspond to m = ±2. The remaining elements are: −3 |r| −3 [S z (I x cos(φ) + I y sin(φ)) + I z (S x cos(φ) + S y sin(φ))] cos(θ) sin(θ) Which leads us to our final definitions of the m = ±1 contributions: We define the prefactors

MEASUREMENT PROTOCOL
We explore two measurement schemes -correlation spectroscopy [1] and a Qdyne type measurement [2]. In the following we provide an outline of the two. We consider the NV center as a two level system (spin 1 2 ) with an energy gap of γ e B ext , where B ext is the externally applied magnetic field [3]. We denote S i (I i ) as the spin operator of the NV (nuclear spin) at the i direction, S ± (I ± ) as the NV's (nuclear spin's) raising/lowering operators and θ i as the angle between the NV's magnetization axis to the vector connecting the NV to the i-th nucleus. The free Hamiltonian of system in the presence of an external magnetic field is The NV interacts with the ensemble via the dipolar interaction, which is of the general form We assume that ω 0 is the largest frequency in the system. Therefore, in the interaction picture with respect to the (19), Eq. (20) can be written as where took the rotating wave approximation. Driving the NV center with N p π−pulses applied at a frequency of ω p , transforms (21) to the effective HamiltonianH where and τ p = π ωp . Assuming N p is odd, the Fourier series of (23) is Hence, if ω p is close to ω N , such that (ω N − ω p ) t ≡ δt 1, applying the rotating wave approximation to (22) leads toH where we used the explicit spatial formsg j ± = − 3 2 Jr −3 j cos θ j sin θ j e (−1) k iϕj and denoted g j ± = −J 3 2 r −3 j cos θ j sin θ j [4].
In correlation spectroscopy the NV interacts with nuclear spins for a period of time τ via (25) at the end of which the NV state is written onto a memory qubit. After a waiting time t the experiment is carried out again for the same period time τ . In the weak back-action limit, namely when N gτ 1 for a characteristic interaction strength g, the operators in (25) can be approximated to their classical values I i x = cos (φ i ) , I i y = sin (φ i ), where φ is a uniform random variable on [0, 2π]. Eq. (25) can be written accordinglỹ We now turn to the explicit derivation of the protocol. The NV is initialized to the state The state (27) is propagated with (26) for a period of time τ . Assuming δτ 1 this results in Note, that g i ± depends on time through the diffusion process. After a waiting time t, the state (28) evolves again under (26) yielding The probability of measuring |ψ 0 after a time t + 2τ is then In the weak back-action limit N gτ 1 Eq. (30) can be simplified further The averaging of φ i in (31) will result in Averaging (32) over realizations yields where C (t) is the correlation function at time t and C (0) ≡ γ 2 e B 2 rms .

THE MEAN FIELD CALCULATIONS
Assuming the nuclei have finite polarization p. The mean field of an NV whose magnetization axis coincides with theẑ axis is [4] Given the that the nuclei are confined to a cylinder of height L and radius R the integral (34) can be calculated by The same calculation carried out for a hemispherical volume of radius R yields Finally, we calculate the same integral for a spherical volume of radius R We notice that (36) recovers the half-space calculation of [4] in the limit R 0 → ∞ as expect. Eq. (37) recovers the half-space limit upto a factor of 2. Yet, Eq. (35) in the limit R → ∞ goes to zero. This is because the limit R → ∞ implies R L, d and specifically L − d R which is the 2D limit. Since 1 0 d cos θY 0 2 (Ω) = 0, the weighting of the r −3 of the dipolar interaction is crucial for a non-zero signal. But if we take the 2D limit Since ρ = d tan θ we have which is surprisingly the result of constant distance.

The diffusion propagator in a cylinder
In the following we calculate the diffusion propagator in a cylindrical geometry required for the calculation of the correlation function.
The propagator is the solution to the equation with the boundary conditions and the initial conditions We guess a solution to Eq. (40) of the form: Substituting (43) into (40) we arrive aṫ Since both sides of Eq. (44) depend on different independent variables, we deduce thaṫ where C 1 is a constant. The solution to the temporal equation of (45) is straight forward We are left with the right-hand side of Eq. (45) equation We guess a separated solution of the form Substituting (48) into (47) yields Using the same argument as before, we deduce that where C 2 is constant. The solution to the part of equation (50) on Z together with the boundary conditions (41) leads to and We continue solving Eq. (50) for g by guessing where n ∈ Z is required in order to satisfy the periodic boundary condition (41). Substituting Eq. (53) into (50) we arrive at We defineρ = αρ, such that Rewriting (54) in terms of (55) leads to where we chose α = C 2 + C1 D . This is the Bessel equation with solutions The propagator should be regular at ρ = 0, so ∀n B n = 0, reducing the solution to Imposing the boundary condition (41) on ρ leads to Denoting by ν m,n the m'th zero of the derivative of the n'th Bessel function of the first kind, Eq. (60) requires The radial dependence of the propagator is therefore Eq. (61) together with Eq. (51) determines the value of C 1 , We note that for n = 0 = C 1 = 0 = C 2 = 0 the solution Eq.  (41) is The initial condition (42) will set the coefficients A m,k,n : The sum over k is independent, hence, we can solve the z dependency of (65) separately: We definez Substituting Eq. (67) into (66) and using the delta function properties, By multiplying both sides of (68) by cos (sz) and integratingz on [0, π] we arrive at Therefore Eq. (64) simplifies to We continue with the rest of Eq. (65), We multiply both sides of Eq. (71) by ρe −ilφ J l ν c,l and integrate both sides. The integration of the right hand-side of (72) yields The integration on the left-hand side of (72) is equal to Substituting Eqs. (73) and (74) into (72) yields: To determine the value ofC we integrate both sides of (71) over the entire volumẽ However, d dx J 0 = −J 1 (x) so the entire series in Eq. (77) is zero, and The Green's function is therefore: where V = πR 2 L is the cylinder's volume.
The diffusion propagator in a hemisphere/sphere We now repeat the previous section for a hemispherical geometry. We are looking for a solution to the equation with the boundary conditions: and initial condition Using the method of images, we can impose the boundary condition ∂P ∂z z=0 = 0 by symmetrically extending the initial conditions and the boundary conditions.
We guess a solution to Eq. (80) of the form Substituting (83) into (80) leads tȯ (84) We deduceṪ The solution for the temporal part of (85) is We are left with the right-hand side of equation (85), (87) We guess and substitute Eq. (88) into Eq. (87), We can again assume 1 g ∂ ∂r r 2 ∂g ∂r The solutions for Y are the spherical harmonics: and The radial part of Eq. (90) is therefore if C 1 = 0 Eq. (93) simplifies to The solution to (93) that is regular at r = 0 is The boundary conditions (81) require from the solution (95) from which we deduce that the only valid solution is g = const.
Returning to (93), if C 1 = 0, we can define:r Rewriting Eq. (93) in terms of the variable (97), The two linearly independent solutions to this equation are the spherical Bessel functions j l and y l : The propagator should be regular at r=0, so for all l ∈ N 0 it follows that b l = 0, reducing the solution to: The radial boundary condition (81), now reads: We denote byν k,l the k'th zero of the derivative of the l'th spherical Bessel function of the first kind. Eq. (101), therefore, leads toν Hence, Note, that we assumed C 1 = 0, whereasν 1,l = 0 for all l = 1. We therefore change our definition ofν 1,l to the first root that does not equal zero. Summing up Eqs. (86), (91), (100) and (103), the general solution to (80) with the boundary conditions (81) is We now need to impose the symmetrically extended initial condition First, we find the angular contribution by equating (104) at t = 0 with (105), multiplying both sides by Y The radial contribution to the coefficients can be found similarly by Therefore, Substituting Eqs. (106) and (109) into (104) leads to (110) where V = 2π 3 R 3 . We find the constant A by equating (110) at t = 0 to (105) and integrating over the entire volume, which results in A = 1 V . So finally: The propagator for a full sphere is found in a similar fashion and is equal to:

Correlation calculation
In the following we provide general expressions for the magnetic field correlation function of a statistically polarized ensemble and then provide numeric results for the integrals.
The correlation is given by where P is the appropriate diffusion propagator. While solving (113) is generally hard, the asymptotic behavior of the correlation has a simplified expression. In short times t τ D the propagator approaches a delta function and therefore In long times t τ V , however, the propagator approaches the uniform distribution V −1 and therefor The correlation in long times is therefore proportional to the squared mean field. In the following we calculate (114) and (115) for three different geometries. In general, we expect the correlation to decay exponentially for times t τ D , which yields the well-known Lorenztian spectrum at high frequencies. In intermediate times, τ D t τ V the spectrum decays as a power law [5] as a result of diffusion, and at long times t τ V the spectrum decays exponentially to a constant value. In the following we offer some numerical results forC . which support these general arguments. cylinder We first solve (114) for each value of m individually: The limit d R, L recovers the half-plane B 2 rms of [5] as expected. We continue by calculating the long times limit (115).
Finally, we provide numerical results forC (m) . Substituting the diffusion propagator (79) into (113) we arrive at We evaluated the series (122) by calculating the first terms (k, s ∈ [1,25]). The integral was calculated numerically for d = nm and L = R = α for the values α 1 = 200 nm, α 2 = 100 nm, α 3 = 50 nm. The results where fitted in the regimes t τ D and t τ V to exponential functions and for τ D t τ V to a power law. The results for m = 0 are presented in Fig. 1, for m = 1 in Fig. 2 and for m = 2 in Fig. 3 Hemisphere We first solve (114) for each value of m individually:       The limit d R recovers the half-plane B 2 rms of [5]. We continue by calculating the long times limit (115).
We where not able to calculate the integral for m = 1 analytically. Finally, we provide numerical results forC (m) . Substituting the diffusion propagator (111) into (113) we arrive at We evaluated the series (129) for m = 0 by calculating the first terms (l ∈ [0, 30], k ∈ [1, 30]). The integral was calculated numerically for d = nm and R = 200 nm, 100 nm, & 50 nm. The The results where fitted in the regimes t τ D and t τ V to exponential functions and for τ D t τ V to a power law. The results for m = 0 are presented in Fig. 4. Note, that the series converges slowly for large volumes and in practice an increasing number of terms might be needed to get a reasonable estimation at times approaching t = 0. The zeroes of the derivative of the spherical Bessel function where calculated numerically using standard root finding methods with the analytical estimate [6] as the initial guess.

MD SIMULATIONS
We corroborated our analytical calculations by performing molecular dynamics (MD) simulations for two of the three geometries: cylindrical and spherical.

Cylinder
The simulation set-up contains N particles, which are confined to a cylinder of radius R and height L. The particles interact via the Lennard-Jones (LJ) potential 4 [(σ/r) 12 − (σ/r) 6 ] with an interaction cutoff distance of r c = 2.5σ. We performed two distinct simulation. In the first, specular reflections are applied on the top and bottom walls of the cylinder and the Lennard-Jones 9/3 potential is applied at the curved walls to confine the particles to the interior of the cylinder, while in the other the Lennard-Jones 9/3 potential is applied over the entire cylindrical boundary. The system is initialized into a thermal state at temperature T by running a Langevin dynamics simulation until the system reaches thermal equilibrium. Each particle is assigned a random value of spin I z ∈ {−1, 1}. After initializing the system, we ran a deterministic molecular dynamics simulation, integrating Newton's laws using the Velocity-Verlet method with step size ∆t = 0.005 /(mσ 2 ). During the simulation we computed and stored the z component of the magnetic field induced by the particles at the position of the NV where r i is the distance between particle i and the NV center and θ i is the angle between r i and the NV quantization axis. The NV is placed along the axis of rotational symmetry of the cylinder, at a distance d below the bottom of the cylinder. In both set-ups we took R = L = 16.44 in LJ units. Since our analytic model does not include interaction we first checked the dependency of B rms and C 0 (∞)/C 0 (0) with respect to d. The results are presented in Figs. 6 and 7 respectively. Both simulations present a scaling of d −2.7 at certain regions of d/R. Since the expected scaling is d −3 we expect our analytic results to to fit these regions. An interesting feature of the soft walls simulation is that the B rms does not increase as expected for shallow NVs (d/R 1). This is due to the interaction with the walls, which repels the nuclei outside of the effective interaction region causing the B rms to decrease. The normalized asymptotic constants of both simulations agree with the analytic results at large values of d/R. This is expected because in this regime the correlation decays rapidly to the constant value, allowing for a better sampling within the limited integration time. The disagreement for shallower NV's in the reflective walls simulation can be attributed to the underestimate of B rms discussed above. In the soft walls simulation it is more probably a result of the interaction with the walls that prevents nuclei from re-entering the effective interaction region, since the difference occurs at d < σ. Since the results of the reflective boundary better match the expected behavior of the correlation at t = 0 and t → ∞ we henceforth continue only with this simulation. We fitted the analytic results to the simulation by rescaling the simulation's results to match the analytic asymptotic constant and by choosing an appropriate diffusion coefficient. The result are presented in Fig. 8.

Sphere
The spherical geometry simulations follow a very similar protocol to the cylindrical geometry simulations. The number of particles was set to N = 28371 and the radius of the sphere was set to R = 20.47 resulting in particle density of ρ = 0.79σ −3 . The temperature, integration timestep and the total integration time was identical to the cylindrical geometry simulation. The LJ 9/3 potential is applied across the whole surface of the sphere to confine the particles to its interior. A comparison between the simulation results and the analytic prediction is presented at