Efficient continuous wave noise spectroscopy beyond weak coupling
Abstract
The optimization of quantum control for physical qubits relies on accurate noise characterization. Probing the spectral density of semiclassical phase noise using a spin interacting with a continuouswave (CW) resonant excitation field has recently gained attention. CW noise spectroscopy protocols have been based on the generalized Bloch equations (GBE) or the filter function formalism, assuming weak coupling to a Markovian bath. However, this standard protocol can substantially underestimate at low frequencies when the CW pulse amplitude becomes comparable to . Here, we derive the coherence decay function more generally by extending it to higher orders in the noise strength and discarding the Markov approximation. Numerical simulations show that this provides a more accurate description of the spin dynamics compared to a simple exponential decay, especially on short timescales. Exploiting these results, we devise a protocol that uses an experiment at a single CW pulse amplitude to extend the spectral range over which can be reliably determined to .
I Introduction
The problem of a qubit interacting with a noisy environment (bath) is of fundamental importance in the field of quantum information processing. Choosing the optimal strategy to fight decoherence depends on the noise characteristics of a particular qubit implementation. For many solidstate qubits, singleaxis phase noise is dominant, and treating the environment in a stochastic semiclassical approximation suffices to describe the dephasing process. For example, a systemenvironment Hamiltonian of the form , where () is the Pauli operator for the system (environment) and is the coupling strength, is approximated as the semiclassical by tracing over the environmental degrees of freedom suterPRA (). In the limit of many environment qubits forming a spin bath, with intrabath couplings strong compared to , can be treated as a stationary, Gaussiandistributed function with zero mean, i.e. . These properties will be assumed throughout the remainder of the paper.
In the context outlined above, knowledge of the spectral density function , the Fourier transform of the twopoint correlation function for , can be used to optimize quantum control, such as dynamical decoupling (DD) and dynamically corrected gates PhysRevLett.102.080501 (). The spectral information can also be used in decoherence suppression techniques such as holeburning revRef3.PhysRevLett.115.033601 (). One way to estimate is to monitor the response of the qubit as it undergoes DD pulse sequences with certain spectral properties PRL.107.170504 (); suterPRL (); AlmogJPhysB (); N.BarGill_nature (); PhysRevB.72.134519 (); bylander (); fei_yan_nature (); PhysRevLett.109.153601 (); 095389842933333001 (). This can be understood intuitively using the overlap integral approach PhysRevLett.109.020501 (); Biercuk_NJP_2013 (); PhysRevLett.93.130406 (); PhysRevLett.87.270405 (); Cywinski (); suterPRA (); JPB_Biercuk (); uhrigNJP (); PhysRevLett.98.100504 (); PhysRevLett.113.250501 (); 095389842933333001 (). For example, under a series of equally spaced, instantaneous pulses, the bathtraced Hamiltonian becomes where alternates between and at a period corresponding to the pulse spacing , and the frequency of the decoupling cycle is . In this case, an exponential decay of qubit coherence is predicted, , where the timedependent decay rate is determined by the overlap integral of the noise spectral density and the frequencydomain filter function (the Fourier transform of the timedomain filter ), :
(1) 
In other words, the external control sequence acts as a bandpass filter and can be tailored such that the qubit is most sensitive to certain spectral bands of the noise power. The same formalism can be applied in quantum sensing such as spinbased magnetometry for oscillating fields using nitrogenvacancy centers in diamond NaturePhysics.4.810.2008 (); PhysRevLett.106.080802 (); PhysRevA.86.062320 (); Staudacher561 (); NatureNano.10.129.2015 (); NatureNano.10.125.2015 (), further motivating the development of accurate DDbased spectroscopy over an extended bandwidth. If the function is spectrally broad or peaked at many frequencies, extracting becomes challenging, particularly if the functional form of is not known a priori. can be made spectrally narrow if there are sufficiently many decoupling cycles, i.e., where is the decoupling period. The filter function then approaches a deltafunction at the decoupling cycle frequency and its harmonics , where is a nonzero integer. When is much longer than the time scale of the bath correlation decay, can be regarded as constant within the probed spectral width, allowing equation 1 to be written in a discrete form. A protocol for estimating using the discrete form of equation 1 and by taking the harmonics into account was designed and implemented experimentally in Ref. suterPRL ().
The pulsed method becomes disadvantageous at high probe frequencies such that finite pulse width effects cannot be ignored and limit the minimum pulse spacing. Moreover, the lowest frequency (given by the maximum pulse delay) dictates the frequency resolution with which the spectral density function is probed suterPRL (). This makes the protocol inefficient in certain situations. For example, probing an unknown noise spectrum over a wide frequency window requires a very large number of experiments.
An alternative approach is to monitor coherence decay under a continuous wave (CW) “spinlocking” pulse, also known in the NMR literature as a measurement. In this case, the qubit dynamics can only be studied perturbatively due to the noncommutativity of the effective Hamiltonian. experiments have been used in NMR to probe slow atomic motions that give rise to fluctuations in the dipolar field PhysRev.135.A1099 (); PhysRev.137.A235 (); ailionAMR (); LookAndLowe (). The NMR literature, however, has not directly addressed the problem of extracting an unknown and arbitrary from a series of measurements. This was first addressed in the context of the generalized Bloch equations (GBE) formalism GBE (); PhysRevB.72.134519 (). The generalized Bloch equations were derived to describe the relaxation dynamics of a system simultaneously interacting with a heat bath and an arbitrarily strong excitation field. The derivation is based on the following assumptions: (1) the system and the bath are weakly coupled, and are initially in a product state; (2) the time scale of the relaxation of the system is much slower than that associated with the decay of the bath correlation functions and the period of the driving field; (3) the bathinduced coherent system dynamics are negligible compared to that induced by the system Hamiltonian; (4) the rotating wave approximation (RWA). The weak coupling assumption means that we keep terms only up to second order in the systembath coupling strength . The second assumption leads to the aforementioned deltafunction approximation. For the noise model described earlier, the GBE predicts an exponential decay of coherence in the experiment. The decay rate is directly proportional to where is the pulse amplitude (Rabi frequency). The steadystate coherence is negligible as long as the high temperature limit is satisfied. Note that the decay rate here is time independent and thus cannot capture nonMarkovian dynamics. The CW approach can often perform well to higher frequencies than the pulsed method, since finite pulse width effects tend to appear before the maximum excitation power is reached or before the RWA is violated. Moreover, the CW protocol can be more efficient than pulsed methods since a single coherence decay measurement yields the spectral density of noise at the target frequency. noise spectroscopy was demonstrated experimentally in Ref. AlmogJPhysB () for opticallytrapped ultracold atoms coupled to a collisional bath, and in Ref. fei_yan_nature () in the context of superconducting qubit decoherence. In the latter case, the analysis was based on the GBE but included more general noise (relaxation) than considered here.
Neither the CW or pulsed methods can probe to arbitrarily low frequencies using the standard analyses above. In these analyses, the number of drive field periods (decoupling cycles) should be large enough to justify the delta function approximation, i. e. . Since the signal decay timescale is , the minimum probe frequency is limited by the condition (where in the pulsed method). The main goal of this paper is to study spin dynamics under CW excitation beyond approximations (1) and (2) above, so that the signal decay at low frequencies can be better modeled. This information is then used to increase the spectral range over which noise spectroscopy produces valid results. We describe the state evolution in the Liouville representation liouville () and apply the cumulant expansion method cumulant1 (); cumulant2 () to calculate the ensemble average, finding the functional form of the coherence decay up to fourth order in (or second order in ). The resulting equations are derived without any assumptions about the CW pulse length or the bath correlation time, in order to capture nonMarkovian behaviour. These results are used to design a CW noise spectroscopy protocol that extends the range for which can be accurately determined down to .
The remainder of the paper is organized as follows. Section II derives the coherence decay function in the experiment up to fourth order in (i.e., second order in ). Section III compares our results to the standard exponential decay function and shows that our model predicts the signal decay significantly better in the short time regime. Section IV.1 presents the noise spectroscopy protocol exploiting the coherence decay function derived in Sec. II, and the improved accuracy in the determination is demonstrated in Sec. IV.2 via numerical simulations. Section V concludes.
Ii Coherence decay function
In this section we derive the coherence decay function of a system under CW driving as a function of the spectral density, , of Gaussian, zeromean semiclassical phase noise as introduced above. In the interaction frame of the CW pulse of amplitude along in the rotating frame, the semiclassical stochastic Hamiltonian transforms in time as
(2) 
The derivation of the coherence decay function under this Hamiltonian must involve perturbation series, since . Using stochastic Liouville theory and super operator formalism liouville (), the ensemble averaged qubit evolution can be described as , where is the density matrix decribing the qubit at time , denotes ensembles averaging over noise realizations, denotes vectorization that stacks the rows of a matrix into a vector, and
(3) 
where is the unit matrix, and is the Dyson time ordering operator cumulant2 (). The ensemble average of the noisy operator can be evaluated with the cumulant expansion
(4) 
where is called the cumulant function and is called the th cumulant cumulant1 (); cumulant2 (). By Taylorexpanding and comparing equations 3 and 4, the cumulants can be found.
For the Hamiltonian in equation 2, the powers of are linear combinations of the operators from following set:
(5) 
Moreover, is proportional to the th power of . As a result, is a linear combination of the operators in . Thus, can be expressed as
(6) 
where is the th element in . In the spinlocking experiment, the normalized signal is
(7)  
(8) 
where is the element of an operator at row and column . Combining equations 6 and 7,
(9) 
It is convenient to write above equation as
(10)  
(11) 
where , and the index indicates that the contribution is linked to the th cumulant. The above equation has several significant implications. First, the average signal in the measurement is an exponential function whose argument (decay rate) is given as a perturbation series. Second, the th order decay rate is proportional to , and therefore proportional to the ensemble average of products of Liouvillians, , and hence proportional to the average over products of the classical Gaussian distributed random variable, ). Third, the th order decay rate can simply be calculated from the coefficients of the , , , and terms corresponding to the th order cumulant, independently from other cumulants. Moreover, since is Gaussiandistributed with zero mean, for an integer , according to Isserlisâ Gaussian moment theorem 10.2307/2331932 (). Subsequently, only even order cumulants are nonzero, and the first few terms are
(12)  
(13)  
(14) 
In the following, we present analytical calculation of the first two nonzero decay terms. We also show that in the limit of , the nd order decay rate is identical to the GBE result. To go beyond the GBEbased analysis, we omit the large approximation, and describe the qubit dynamics including the th order decay rate.
ii.1 Secondorder decay rate
By expanding equation 12 (the nd order cumulant) using the operators in equation 5, we find . Then from equation 11 we acquire the decay rate attributed to the nd order cumulant
(15)  
(16) 
where
The imaginary part of is an odd function with respect to , and is an even function, thus equation 16 can be further simplified as
(17) 
Figure 1 shows the real parts of for fixed values of (a) and (b) the number of Rabi cycles . The even function behaves like as (i.e., as ). Thus, for a fixed value of and in the limit of , the nd order decay rate can be approximated as
(18) 
Therefore, the normalized spin signal attributed to the nd order cumulant in the limit of large is a simple exponential decay:
(19) 
This expression is equal to the result derived from the GBE in the high temperature limit .
ii.2 Fourthorder decay rate
Following the same steps as for the nd order term evaluation, the th order cumulant can be expressed in terms of the operators in equation 5, and the coefficients can be found. The decay rate attributed to the th order cumulant can be expressed as
(20) 
We can use Isserlis Gaussian moment theorem again, and write the th order correlation function as the products of the nd order correlation functions: The product of the nd order correlation function can be expressed in terms of the spectral density, as before. Then, the th order decay rate can be rewritten compactly as
(21) 
where
(22) 
The above time integration can be carried out analytically. Symmetries in the filter function result in the imaginary component of the integral going to zero as was the case for 2nd order decay rate. Finally, using the symmetry of the expression for the th order decay rate can be expressed as
(23) 
where
Iii Accuracy of coherence decay
To evaluate the accuracy of the cumulant expansion decays, we simulate a set of experiments with known input noise spectra, . The simulated signal decays are generated by timediscretized unitary evolution of an initial state density matrix, using randomly generated noise realizations. A cosine series representation, as described in reference Shinozuka1991 (), is used to generate the noise realizations, in equation 2. This simulation method accurately represents stationary, Gaussian noise matching the input noise spectrum, and converges to the correct spectrum at a rate of 1/N, where N is number of noise samples used. For all noise spectra in this paper, plateaus below , i.e. for . These simulated signal decays are used in the following sections to represent experimental data. Since the number of noise realizations is finite, we recalculate based on the simulated signal decays, and this is the that appears in the plots below.
Using the simple exponential expression from equation 19, can be accurately determined when the assumptions listed in section 1 are valid, i.e. when the relaxation timescale is long compared to the drive field period (). When this condition is not met, the signal decay can be nonexponential and the standard analysis can give inaccurate values for . Figure 2a shows the result of applying a leastsquares exponential fit and using equation 19 to determine for a input noise spectrum. In this example, and was obtained by fitting simulated CW noise spectroscopy experiments at . The insets show examples of the normalized signal decay , which becomes nonexponential at low , leading to inaccurate values. Figures 2b and c show this same procedure applied to (b) a noise spectrum (c) and a noise spectrum with 100 times higher noise power than that shown in (a). Simulated CW experiments at were used for figure 2b, and experiments at were used for figure 2c.
The nonexponential signal decay displayed in the inset of figure 2a can be understood based on the shape of the filter function with respect to when the function approximation no longer holds. Figure 2d shows a nonexponential signal decay, with displayed in insets at certain time points in the decay, along with the noise spectrum . The finite width of the main filter function peak, as well as its satellite peaks, overlaps with large values of at low frequencies, producing oscillations in the signal decay.
Figure 3 shows the error, as a function of decay time and drive field amplitude, between the simulated “experimental” signal decays and theoretical decays calculated using one of three different methods. Each calculation uses the same that generated the simulated decays. Figure 3a shows the result of applying the simple exponential from equation 19. Figures 3b and c show the decays calculated via the cumulant expansion from equation 11, with (b) using only , and (c) using both and .
The cumulant expansion method provides much better matching to signal decays at short times (t 0.5 s in the example in figure 3) compared to the standard exponential decay of equation 19. However, using only introduce errors in the intermediate (2035 rad/s) region. Using the cumulant expansion up to fourth order (), the mismatch is reduced at all times over a large portion of the drive amplitudes ( for the parameters chosen in this example).
Iv Noise spectroscopy based on the cumulant expansion
The cumulant expansion method can be used to improve CW noise spectroscopy when the experimental signal decays become nonexponential. The accuracy of a given noise spectrum estimate can be tested by comparing the cumulant expansion signal decay, calculated using , and the experimental decay. Furthermore, the nonexponential, oscillatory behaviour observed at short timescales is the result of a wide frequency filter that overlaps with across a range of frequencies, sometimes extending to . This shorttime behaviour thus contains broad spectral information and can be used to extend the range over which can be determined. In particular, one can choose a drive frequency for which the signal decay is wellmatched by the calculation (e.g. in figure 3) and extract information about for from detailed fitting of the shorttime behaviour.
To illustrate this, figure 4 shows the shorttime behaviour of a calculated signal decay for two different noise spectra. One spectrum is labelled ‘correct’, while the other represents an error in which at low frequencies has been changed. Here, the error is introduced for , while the drive amplitude . This shows that the decay is sensitive to variation in at frequencies far below the probing frequency .
iv.1 Noise spectroscopy protocol
To take advantage of the accuracy of the cumulant expansion for improving noise spectroscopy, we propose a gradient ascent protocol based on matching to the experimental signal decay using a single chosen pulse amplitude, . An initial estimate of the noise spectrum, , is obtained from the standard approach of fitting to exponential decays for a group of probe frequencies. Then, by accurately fitting a detailed signal decay at using the cumulant expansion method, particularly in the shorttime regime, the full noise spectrum can be determined.
The signal decay given by the cumulant expansion method is
We discretize the above expression, and define a fitness function as the rootmean square error between the experimentally measured decay, , and the calculated decay for a given :
(24) 
We can then calculate the gradient of the fitness function, , for any target frequency . The gradient is used to update the estimate of towards a closer matching of the experimental and calculated decays. The full protocol is:

To obtain an initial estimate, , use exponential fits of decays over a range of , and matching to .

Select a drive amplitude, , for detailed matching of the decay curve. should be low enough to display nonexponential features at short times, but not so low that the calculation is inaccurate. Otherwise, the following steps will not converge to a high fitness function .

Calculate the fitness function, , for the decay at using the initial estimate .

While for some target threshold ,

Calculate for the current estimate

Update the estimate: for all , for some fixed

Calculate the updated fitness and increment .

To improve the speed of calculation/convergence, we can use the knowledge that the simple exponential decay fitting is accurate when the signal decays are smooth exponentials, and only update for where that condition is not satisfied.
iv.2 Demonstration of protocol
The cumulant expansion noise spectroscopy protocol described above was applied to the simulated experiments presented in section 3 corresponding to three different noise spectra. Figure 5 shows the result obtained using the simulation with the input noise (with a plateau as ). The initial estimate, , uses the standard exponential fitting method at 11 pulse amplitudes in the range of 20125 rad/s. The final estimate was obtained by detailed fitting to a single decay curve as described above. For comparison, this final step was done with three different choices for the parameters , where is the total pulse duration. The final noise spectrum estimate is a much better fit in the low frequency regime to the correct (input) spectrum.
Some artifacts are introduced in the form of oscillations in the intermediate frequency range (). These artifacts have characteristic periods of order and are a consequence of remaining error between the decay and the true decay, such as contributions from and higher terms. These oscillations are not a consequence of the gradient optimization and we have not found a straightforward way to remove them. However, given that is typically a smooth function, and that these oscillations are confined to a certain band of frequencies, smoothing or sliding window averaging can be used to suppress the oscillations in the final estimate. Alternately, they can be fully removed if the noise spectrum can be fit to a certain functional form, such as . Figure 6 shows the results of applying the cumulant expansion noise spectroscopy protocol to the same three simulated experiments shown in figure 2. The upper panels show the input, initial, and final determined by fitting the signal decay at pulse amplitude and total time . The lower panels of figure 6 show the result of fitting from the upper panel with a general power law , where and are free parameters. Note that this power law fit is applied in a frequency range that excludes the plateau region in . To obtain the initial estimate, , figure 6a uses the same experiment set as figure 5, figure 6b uses experiments at 15 pulse amplitudes in the range 1120 rad/s, and figure 6c uses 9 pulse amplitudes in the range 101250 rad/s.
V Conclusion
In summary, we have treated the problem of spin evolution in the presence of singleaxis phase noise during an experiment with CW excitation, with a goal to improve the determination of an arbitrary noise spectral density. By retaining cumulant expansion terms up to fourth order in the systembath coupling, we can more accurately match coherence decay dynamics that exhibit nonexponential and oscillatory behaviour, and thereby extract more accurate spectral information, especially at low frequencies. We present a twostep protocol: (1) estimate using the standard exponential fitting approach by probing over a set of frequencies (lowresolution signal decay experiments); (2) refine based on fitting a single, highresolution signal decay using the fourth order cumulant expansion calculation. Since this second step consists of probing at a single frequency, it is efficient in terms of experimental resources. For the cases of and noise (with low frequency cutoff plateau), we have shown that this protocol allows for accurate determination of to zero frequency, i.e. the low frequency regime where standard CW and pulsed noise spectroscopy fail. While the examples given were noise spectra of the form , the theoretical analysis and protocol are applicable to arbitrary spectra, and in future work we plan to test this applicability in simulations and real experiments. In addition, inhomogeneous broadening is typical in physical spin systems, and should also be included. This can be expressed as an additional Hamiltonian component , where is a static random variable. Thus, inhomogeneous broadening yields a peak in at , which should enhance the oscillations in signal decay at short timescales for low probing frequencies. Our protocol should therefore reveal such broadening. Additional work could also extend the cumulant expansion noise spectroscopy protocol to include multiaxis noise and/or higher order cumulants () for more general applications.
Acknowledgements.
This work was supported by Natural Sciences and Engineering Research Council (NSERC). D.K.P. was supported by the National Research Foundation of Korea (Grants No. 2015R1A2A2A01006251 and No. 2016R1A5A1008184).References
 Ashok Ajoy, Gonzalo A. Álvarez, and Dieter Suter. Optimal pulse spacing for dynamical decoupling in the presence of a purely dephasing spin bath. Phys. Rev. A, 83:032303, Mar 2011.
 Kaveh Khodjasteh and Lorenza Viola. Dynamically errorcorrected gates for universal quantum computation. Phys. Rev. Lett., 102:080501, Feb 2009.
 Dmitry O. Krimer, Benedikt Hartl, and Stefan Rotter. Hybrid quantum systems with collectively coupled spin states: Suppression of decoherence through spectral hole burning. Phys. Rev. Lett., 115:033601, Jul 2015.
 Tatsuro Yuge, Susumu Sasaki, and Yoshiro Hirayama. Measurement of the noise spectrum using a multiplepulse sequence. Phys. Rev. Lett., 107:170504, 2011.
 Gonzalo A. Álvarez and Dieter Suter. Measuring the spectrum of colored noise by dynamical decoupling. Phys. Rev. Lett., 107:230501, 2011.
 Ido Almog, Yoav Sagi, Goren Gordon, Guy Bensky, Gershon Kurizki, and Nir Davidson. Direct measurement of the system environment coupling as a tool for understanding decoherence and dynamical decoupling. J. Phys. B: At. Mol. Opt. Phys., 44(15):154006, 2011.
 N. BarGill, L.M. Pham, C. Belthangady, D. Le Sage, P. Cappellaro, J.R. Maze, M.D. Lukin, A. Yacoby, and R. Walsworth. Suppression of spinbath dynamics for improved coherence of multispinqubit systems. Nat. Commun., 3:858, 2012.
 G. Ithier, E. Collin, P. Joyez, P. J. Meeson, D. Vion, D. Esteve, F. Chiarello, A. Shnirman, Y. Makhlin, J. Schriefl, and G. Schön. Decoherence in a superconducting quantum bit circuit. Phys. Rev. B, 72:134519, Oct 2005.
 Jonas Bylander, Simon Gustavsson, Fei Yan, Fumiki Yoshihara, Khalil Harrabi, George Fitch, David G. Cory, Yasunobu Nakamura, JawShen Tsai, and William D. Oliver. Noise spectroscopy through dynamical decoupling with a superconducting flux qubit. Nat Phys, 7(7):565–70, 2011.
 Fei Yan, Simon Gustavsson, Jonas Bylander, Xiaoyue Jin, Fumiki Yoshihara, David G. Cory, Yasunobu Nakamura, Terry P. Orlando, and William D. Oliver. Rotatingframe relaxation as a noise spectrum analyser of a superconducting qubit undergoing driven evolution. Nat. Commun., 4:2337, 2013.
 D. H. Slichter, R. Vijay, S. J. Weber, S. Boutin, M. Boissonneault, J. M. Gambetta, A. Blais, and I. Siddiqi. Measurementinduced qubit state mixing in circuit qed from upconverted dephasing noise. Phys. Rev. Lett., 109:153601, Oct 2012.
 P SzaÅkowski, G Ramon, J Krzywda, D Kwiatkowski, and Å CywiÅski. Environmental noise spectroscopy with qubits subjected to dynamical decoupling. Journal of Physics: Condensed Matter, 29(33):333001, 2017.
 Todd Green, Hermann Uys, and Michael J. Biercuk. Highorder noise filtering in nontrivial quantum logic gates. Phys. Rev. Lett., 109:020501, Jul 2012.
 Todd J Green, Jarrah Sastrawan, Hermann Uys, and Michael J Biercuk. Arbitrary quantum control of qubits in the presence of universal noise. New J. Phys., 15:095004, 2013.
 A. G. Kofman and G. Kurizki. Unified theory of dynamically suppressed qubit decoherence in thermal baths. Phys. Rev. Lett., 93:130406, Sep 2004.
 A. G. Kofman and G. Kurizki. Universal dynamical control of quantum mechanical decay: Modulation of the coupling to the continuum. Phys. Rev. Lett., 87:270405, Dec 2001.
 Łukasz Cywiński, Roman M. Lutchyn, Cody P. Nave, and S. Das Sarma. How to enhance dephasing time in superconducting qubits. Phys. Rev. B, 77:174509, May 2008.
 M J Biercuk, A C Doherty, and H Uys. Dynamical decoupling sequence construction as a filterdesign problem. J. Phys. B: At. Mol. Opt. Phys, 44(15):154002, 2011.
 Götz S. Uhrig. Exact results on dynamical decoupling by pulses in quantum information processes. New J. Phys., 10(8):083024, 2008.
 Götz S. Uhrig. Keeping a quantum bit alive by optimized pulse sequences. Phys. Rev. Lett., 98:100504, 2007.
 Gerardo A. PazSilva and Lorenza Viola. General transferfunction approach to noise filtering in openloop quantum control. Phys. Rev. Lett., 113:250501, Dec 2014.
 J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth, and M. D. Lukin. Highsensitivity diamond magnetometer with nanoscale resolution. Nature Physics, 4:810, 09 2008.
 G. de Lange, D. Ristè, V. V. Dobrovitski, and R. Hanson. Singlespin magnetometry with multipulse sensing sequences. Phys. Rev. Lett., 106:080802, Feb 2011.
 Masashi Hirose, Clarice D. Aiello, and Paola Cappellaro. Continuous dynamical decoupling magnetometry. Phys. Rev. A, 86:062320, Dec 2012.
 T. Staudacher, F. Shi, S. Pezzagna, J. Meijer, J. Du, C. A. Meriles, F. Reinhard, and J. Wrachtrup. Nuclear magnetic resonance spectroscopy on a (5nanometer) sample volume. Science, 339(6119):561–563, 2013.
 Stephen J. DeVience, Linh M. Pham, Igor Lovchinsky, Alexander O. Sushkov, Nir BarGill, Chinmay Belthangady, Francesco Casola, Madeleine Corbett, Huiliang Zhang, Mikhail Lukin, Hongkun Park, Amir Yacoby, and Ronald L. Walsworth. Nanoscale nmr spectroscopy and imaging of multiple nuclear species. Nature Nanotechnology, 10:129, 01 2015.
 T. Häberle, D. SchmidLorch, F. Reinhard, and J. Wrachtrup. Nanoscale nuclear magnetic imaging with chemical contrast. Nature Nanotechnology, 10:125, 01 2015.
 Charles P. Slichter and David Ailion. Lowfield relaxation and the study of ultraslow atomic motions by magnetic resonance. Phys. Rev., 135:A1099, Aug 1964.
 David C. Ailion and Charles P. Slichter. Observation of ultraslow translational diffusion in metallic lithium by magnetic resonance. Phys. Rev., 137:A235, Jan 1965.
 D. C. Ailion. Nmr and ultraslow motions. Adv. Mag. Res., 5:177–227, 1971.
 D. C. Look and I. J. Lowe. Nuclear magnetic dipoleâdipole relaxation along the static and rotating magnetic fields: Application to gypsum. The Journal of Chemical Physics, 44(8):2995, 1966.
 Eitan Geva, Ronnie Kosloff, and J. L. Skinner. On the relaxation of a twolevel system driven by a strong electromagnetic field. J. Chem. Phys, 102(21):8541, Jun 1995.
 K. Blum. Density Matrix Theory and Application. Physics of Atoms and Molecules. Springer US, 1996.
 Ryogo Kubo. Stochastic liouville equations. J. Math. Physics, 4(2):174, 1963.
 Ryogo Kubo. Generalized cumulant expansion method. J. Phys. Soc. Jap., 17(7):1100–1120, 1962.
 L. Isserlis. On a formula for the productmoment coefficient of any order of a normal frequency distribution in any number of variables. Biometrika, 12(1/2):134–139, 1918.
 Masanobu Shinozuka and George Deodatis. Simulation of stochastic processes by spectral representation. Applied Mechanics Reviews, 44(4):191–204, 1991.