For quantum sensing, it is vital to develop an efficient technique for determining the quantum state of the sensor. We optimize the weighting of the photoluminescence intensity for readout of the spin state of the nitrogen-vacancy (NV) center in diamond. We find that adopting a physical model that considers the optical transitions and relaxations of the NV center allows for an efficient readout. Our method improves the signal-to-noise ratio of the readout by 5.4% in a short time of 3 s, while the existing methods typically require 1 min of integration time. We also show that our technique enhances the readout of the nuclear spin memory. The demonstrated way is helpful for a wide range of measurements, from a few minutes to several days.

A nitrogen-vacancy (NV) center in diamond is a qubit that serves as a nanoscale sensor.^{1–7} The signal-to-noise ratio (SNR) of the spin-state readout is an essential factor directly determining the sensitivity.^{8,9} In most studies, the readout of NV spins relies on time-resolved photoluminescence intensity data (PL) obtained during about 500 ns until the electron spin is initialized by laser irradiation.^{10–12} Since the average number of photons obtained in a single readout is about 0.02, the same measurement must be repeated more than tens of thousands of times to obtain sufficient SNR.^{8,9} Data processing that maximizes the spin information from the limited photons is beneficial to reduce the measurement duration.

A method to improve the SNR of the spin-state readout by appropriately weighting the PL was reported previously.^{12} The weighting factors are obtained from the PL of the NV center electron spin *m*_{S} = 0 and *m*_{S} = −1 states as the reference data. This method requires accurate reference data, which is time-consuming to acquire. Another method has been proposed to determine weighting factors by machine learning that optimizes a unique cost function considering shot noise.^{13} However, a quantitative comparison of these methods is lacking, and it is still unclear which way is more effective.

In this article, we show that an adequate physical model is useful for obtaining weighting factors robust against shot noise. We focus on the five-level model (FLM),^{11,12} where the optical transitions and relaxations in the NV center are phenomenologically included. We prove that the weighting factors based on this model are optimal for improving the readout efficiency among the proposed weighting methods.^{12,13} Using the Rabi oscillation’s SNR as a benchmark, we quantitatively compare our model with the others. We also prove that our method applies to the repetitive readout of nuclear spin memory.^{14,15} Our method, which effectively reduces the measurement duration with low computational cost, is useful for various sensing measurements.

We start to explain a general weighting procedure for the readout of the NV electron spin state from the PL. We focus on the readout of *m*_{S} = 0 and *m*_{S} = −1 states only by considering a common situation in which the energy degeneracy of the *m*_{S} = ±1 state is resolved by applying a magnetic field parallel to the NV symmetry axis. Figure 1(a) shows the PL obtained by photon counting for 1250 ns with a time resolution of 10 ns. These work as reference data, i.e., PL {*a*_{i}} and {*b*_{i}} obtained for the electron spin states *m*_{S} = 0 and *m*_{S} = −1, respectively, where *i* = 1, 2, …, and 125 are the indices of the time bin. The PL continues to depend on the electron spin state for about 800 ns after the laser irradiation because of a unique optical transition of the NV center that initializes its electron spin.^{16,17} Since {*a*_{i}} and {*b*_{i}} depend on individual NV centers and laser intensity, it is desirable to measure them for each measurement. We define the probability *p* that the NV center is in the *m*_{S} = 0 state as follows:

where *x*_{i} is the measured intensity and *w*_{i} is the weighting factor. *w*_{i} should be set normalized appropriately to obtain *p* = 1 if *x*_{i} = *a*_{i} and *p* = 0 if *x*_{i} = *b*_{i}.

Next, we explain how to determine the set of *w*_{i}, {*w*_{i}}. Figure 1(b) shows {*w*_{i}} obtained by different weighting methods based on the PL shown in Fig. 1(a). We compare our method with the three existing methods; the most conventional method (CM), maximum likelihood estimation approximation (MEa),^{12} and machine learning (ML).^{13} In CM, {*w*_{i}} are a finite constant value for 0 ≤ *i* ≤ *ℓ* and otherwise 0, and *ℓ*(≤125) is determined to give the most accurate readout [the black line in Fig. 1(b)]. In MEa, {*w*_{i}} are obtained by normalizing the difference in the reference data {*a*_{i} − *b*_{i}} [the blue line in Fig. 1(b)]. When the SNR of {*a*_{i}} and {*b*_{i}} are sufficiently high, MEa is equivalent to maximum likelihood estimation.^{12} In ML, {*w*_{i}} are obtained by minimizing a certain cost function^{13} [the green line in Fig. 1(b)]. ML costs computationally because it should optimize {*w*_{i}}, which contains 125 parameters. We use MATLAB^{®} Optimization Toolbox™ to numerically optimize {*w*_{i}} using {*w*_{i}} of MEa as the initial value. When the SNR of {*a*_{i}} and {*b*_{i}} is low, {*w*_{i}} of MEa and ML are disturbed by photon shot noise as shown by the blue and green lines in Fig. 1(b), respectively.

As a new method unlike the above, we appropriately smooth the weighting factors based on the FLM for the optical transitions in the NV center; FLM can reproduce PL by adequately considering the ground states (*m*_{S} = 0, ±1), the excited states (*m*_{S} = 0, ±1), and the intermediate singlet state.^{11,12} The level diagram and the optical transition are illustrated in Fig. 1(c). There are only seven parameters: six transition coefficients (*γ*_{0}, *γ*_{1}, *S*_{0}, *S*_{1}, *D*_{0}, and *D*_{1}) and polarization *ζ* (see the Appendix). Here, *γ*_{0} and *γ*_{1} are the optical transition rate. *S*_{0} and *S*_{1} (*D*_{0} and *D*_{1}) are the transition rate between excited (ground) states and singlet states. Due to this simplicity, the computation time is negligible compared to the measurement duration. The three solid lines in Fig. 1(a) result from the fitting of the reference data using FLM. They agree with the experimental results within the uncertainty of the shot noise. By normalizing the difference between these fitted curves [the black solid line in Fig. 1(a)], we obtain *w*_{i}, as shown in the red line in Fig. 1(b). This smooth line is determined robustly against shot noise and effectively suppresses the effect of the noise on the analysis.

Under the common situation we are considering, it is possible to treat the optical transitions of the *m*_{S} = ±1 states together, and thus, the FLM works well.^{11,12} Although, strictly speaking, the optical transitions of the NV center are more complex, our method should work as long as it reproduces the spectral shape of the PL regardless of the details of the fitting parameters. Of course, we could also use different physical models depending on the experimental situation, for example, the seven-level model^{18} is used to account for the effects of mixed spin states due to the vertical component of the magnetic field.

We measure Rabi oscillations to quantitatively compare the efficiency of the four methods (CM, MEa, ML, and FLM). As an experimental platform, we use a single NV center in a ^{12}C-enriched diamond, which is grown by chemical vapor deposition^{19} with ^{15}N isotopically purified gas on a type-IIa diamond. A magnetic field of 233 mT is applied along the NV symmetry axis, and a sufficient microwave (MW) driving field is applied from a stripline nearby the present NV center.

Figure 2(a) shows the measurement protocol. As in previous studies,^{20,21} we acquire reference data {*a*_{i}} and {*b*_{i}} by repeating the sequence consisting of laser pulses and an adiabatic inversion pulse^{22,23} *N*_{ref} times and define the acquisition time as *T*_{ref}. We then measure the Rabi oscillation by sweeping the duration of the rectangular microwave pulse *t*_{j,burst} (*j* = 1, 2, …,71). We define “one Rabi measurement” as the repetition of the sequence for *N*_{acq} = 3.9 × 10^{5} times. We also define the measurement duration as *T*_{acq}, which corresponds to 86 s for the present Rabi measurement.

We now describe the SNR of the Rabi measurement as an evaluation index. Figure 2(b) shows an example of one Rabi measurement processed by MEa. Due to the high Rabi frequency compared to the hyperfine interaction with ^{15}N,^{24} *p* should exhibit a simple sinusoidal oscillation with a maximum value of 1 and a minimum of 0. When {*w*_{i}} is noisy, the range of obtained *p* fluctuates by a few percent. We normalize the obtained *p* with its amplitude deduced by fitting and get normalized data *p*_{norm} and the fitted curve *p*_{fit}. This normalization makes it possible to compare four methods with equal indices. Figure 2(c) shows the difference *δ*(= |*p*_{norm} − *p*_{fit}|) at each measurement point of the Rabi oscillation. We define the inverse of the root-mean-square of *δ* in the overall data as the SNR of the electron spin state estimation.

We examine the impact of the shot noise on the reference data by changing *N*_{ref}. Using the same data shown in Fig. 2(c) as an example, the red (blue) histogram in Fig. 2(d) shows the distribution of the SNR obtained by 100 Rabi measurements analyzed with *N*_{ref} = 3.9 × 10^{6}(*N*_{ref} = 3.9 × 10^{5}). The SNR with *N*_{ref} = 3.9 × 10^{6} is higher than that with *N*_{ref} = 3.9 × 10^{5}. This indicates that the analysis of the Rabi measurement is limited by the shot noise of the reference data {*a*_{i}} and {*b*_{i}}. For sufficiently large *N*_{ref}, the SNR is limited only by the shot noise of the Rabi oscillation data {*x*_{i}}, indicating that the analysis is optimal. To reduce the measurement duration, it is critical to obtain high SNR with less *N*_{ref}.

Figure 3(a) shows the *N*_{ref} dependence of SNR for different methods (CM, MEa, ML, and FLM). We make four remarks. First, the SNR of MEa (blue) is better than that of CM (black) for *N*_{ref} > 10^{6} and worse for *N*_{ref} < 10^{6}. *N*_{ref} < 10^{6} is the region where the estimation is limited by the shot noise on the reference data. Second, ML (green) gives almost the same result as MEa. This is reasonable because {*w*_{i}} are almost the same between ML and MEa, as shown in Fig. 1(b). Third, the SNRs of CM and our FLM (red) do not depend on *N*_{ref}. They are robust to the shot noise of the reference data as shown in Fig. 2(e), indicating that integration of *N*_{ref} < 3.9 × 10^{5} is sufficient. Fourth, our FLM always shows the highest SNR, even for small *N*_{ref}. This means that our method is the optimal analysis.

We now discuss how much our method can shorten the measurement. We compare each method with the SNR of CM at *T*_{ref} = 3 s and *T*_{acq} = 86 s, SNR_{CM}. The condition of SNR_{CM} is the best in CM for time efficiency and is appropriate for comparison. According to the central limit theorem, the SNR should improve in proportion to the square root of *T*_{acq}. In the analysis, the effective amount of PL acquired during *T*_{acq} varies with the weighting factors, resulting in a change in SNR as shown in Fig. 3(a). It is necessary to consider that the total measurement duration will be longer by *T*_{ref}. The total measurement duration required to achieve the equivalent SNR as SNR_{CM} is obtained as

Figure 3(b) shows the *T*_{ref} dependence of *T*_{tot} calculated with Eq. (2) for different methods. Clearly, for any *T*_{ref}, *T*_{tot} of MEa and ML is longer than the best of CM [black dashed line in Fig. 3(a)]. The reason why these methods are not beneficial is that they do not improve the SNR while *T*_{ref} is short compared to *T*_{acq}. On the other hand, our FLM only requires *T*_{tot} = 81.0 s at *T*_{ref} = 3 s. This indicates that FLM can reduce measurement duration by 9.0% compared to the best of CM. Therefore, our method is the only weighting method that is beneficial for short measurements.

The above-demonstrated efficient readout using physical models also applies to other cases. To exemplify this, in the rest of this article, we discuss the repetitive readout in NMR spectroscopy with NV centers,^{25–28} where the nuclear spin memory storing the electron spin state information is readout thousands of times. During this repetition, the nuclear spins relax in a time characterized by the relaxation time *T*_{1}. We introduce additional weighting factors for each readout, determined by exponential fitting.

Our repetitive readout consists of *N* readouts. We use the ^{15}N nuclear spin in the present NV center as the memory. Figure 4(a) schematically shows the pulse sequence and the resulting PL. We start to explain the first readout as an example (see “first” in Fig. 4(a)]. The MW pulse restores the nuclear spin state to the electron spin state, and the laser pulse is used to read it out. We get the PL as shown in the lower panel of Fig. 4(a). {*a*_{i,1}}({*b*_{i,1}}) is the PL in the case of *m*_{S} = 0(*m*_{S} = −1), which we obtain as reference data. Following Eq. (1), *p*_{1}, the probability that the electron spin is in *m*_{S} = 0, indicating that the nuclear spin state is *m*_{I} = +1/2, is obtained as $p1=\u2211i=1125wi(xi,1\u2212bi,1)$. Next, we repeat the readout for *N* times. Using the overall readouts {*p*_{n}} (*n* = 1, 2, …, and *N*), the probability *P* that *m*_{I} = +1/2 is stored in the memory is obtained as

where {*W*_{n}} are the weighting factors for each readout. Figure 4(b) shows $An=\u2211i=1125wiai,n$ and $Bn=\u2211i=1125wibi,n$, the sum of the reference data {*a*_{i,n}} and {*b*_{i,n}} in each readout. We obtain the weighting factors {*W*_{n}}: $Wn\u221de\u2212ntL/T1$ by applying exponential fitting of the difference {*A*_{n} − *B*_{n}}, as shown in Fig. 4(c). In the experiment, we set the duration of the single readout as *t*_{L} = 2.88 *µ*s and obtain the relaxation time of the memory as *T*_{1} = 1.3 ms.

We measure Rabi oscillations of the electron spin of the NV center with repetitive readouts. The obtained oscillation is analyzed to determine the SNR as in Fig. 2. We compare our weighting method with a conventional method. Our method is based on the combinations of FLM and memory relaxation (FLM + *T*_{1}) to determine {*w*_{i}} and {*W*_{n}}. We define the conventional method (CM) as the method of tuning only *N*, with {*w*_{i}} and {*W*_{n}} as constant. Figure 4(d) shows the *N* dependence of the SNR. For CM, as *N* increases, the SNR reaches a maximum value, SNR_{CM,max}, and then decreases. This reflects memory relaxation for large *N*.^{29} For our FLM + *T*_{1} method, the SNR remains constant within the error bars even for large *N*. Since our method properly includes the effect of memory relaxation, we can ignore the readout without spin information and, thus, maintain the SNR regardless of *N*.

We discuss the efficiency of our method. The maximum SNR for both methods is obtained at about *N* = 500. Our FLM + *T*_{1} method has an of 11.6% from CM for this condition. Typically, for sensing such as NMR spectroscopy, *T*_{acq} is ten times longer than *T*_{ref} or even more.^{26–28} Following Eq. (2), the time to obtain equivalent SNR with SNR_{CM,max} is obtained as $Ttot\u2248Tacq\xd7(SNR/SNRCM,max)\u22122$. Specifically, our method reduces the total measurement time by 19.8%. For example, when the sensing period is ten days with CM, our method can shorten it by two days. Additionally, our method will also be effective for repetitive readouts at low magnetic fields (below 100 mT) where the relaxation time of nuclear spin memory is short.^{15}

To conclude, we demonstrated that PL weighting based on physical models is optimal for improving the SNR of the optical spin readout of the NV center. Our technique requires only a low computational cost. Since the optimal analysis can be demonstrated with a reference data measurement of about 3 s, even a routine short measurement can be made more efficient. The repetitive readout of the nuclear spin can also improve SNR by 11.6%, which is helpful to shorten nuclear spin sensing time, typically taking several days.^{25–28}

This work was supported by Grants-in-Aid for Scientific Research (Grant Nos. JP18H01502, JP19H02547, JP19H00656 JP19H05826, JP20K22325, and JP22K03524), Q-LEAP (Grant No. JPMXS0118067395), and the Center for Spintronics Research Network, Keio University.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: DETERMINATION OF THE WEIGHTING FACTOR WITH FIVE LEVEL MODEL

The five-level model (FLM) is used as the model to give the photoluminescence (PL) from the NV center.^{11,12} We explain this model following the previous studies. FLM consists of the ground states (*m*_{S} = 0, ±1), the excited states (*m*_{S} = 0, ±1), and the intermediate singlet state. The level diagram and the optical transition are illustrated in Fig. 1(c). We define the occupancy at each level as

To reproduce the reference data, we define the occupancy after laser initialization (after applying an adiabatic inversion pulse) as *n*_{0}(*n*_{1}). Specifically, we set these as *n*_{0} = [1 − 2*ζ*,2*ζ*,0,0,0]^{T} and *n*_{1} = [*ζ*,1 − *ζ*,0,0,0]^{T}, where *ζ* is polarization. Here, we assume the occupancy at *m*_{S} = +1 of the ground state is equal to at *m*_{S} = −1, and define this as the polarization *ζ*. *n*_{0} and *n*_{1} are also the states corresponding to *p* = 1 and *p* = 0 in Rabi oscillation, respectively. FLM gives the following rate equation:

We obtain the spectral shape of PL from occupancy at the excited state,

We numerically calculate ** n**(

*t*) [Eq. (A2)] using the initial condition

**(**

*n**t*= 0) =

*n*_{0}(

**(**

*n**t*= 0) =

*n*_{1}) to obtain available reference data

*I*

_{0}(

*t*) (

*I*

_{1}(

*t*)).

We fit {*a*_{i}}({*b*_{i}}) with *I*_{0}(*t*) (*I*_{1}(*t*)). The fitting parameters are $F\u2261[\gamma 1,\gamma 0,S0,S1,D0,D1,\zeta ]T$. We use the nonlinear least square method in MATLAB\circledR Optimization Toolbox™ (lsqnonlin) to obtain the optimal fitting parameter *F*_{opt} as follows:

where *I*_{0,i} (*I*_{1,i}) is *I*_{0}(*t*) (*I*_{1}(*t*)) discretized by the index *i* of the time bin. The fitting result *I*_{0}(*t*) (*I*_{1}(*t*)) is shown as the blue (red) solid line in Fig. 1(a). We obtain the weighting factor with FLM as

where the denominator is the normalization constant.

In our experiment, a typical *F*_{opt} is *γ*_{1} = 13 MHz, *γ*_{0} = 40 MHz, *S*_{0} = 13 MHz, *S*_{1} = 60 MHz, *D*_{0} = 3.6 MHz, *D*_{1} = 3.4 MHz, and *ζ* = 0.20. The fitting parameters vary with the signal-to-noise ratio of the reference data and with each NV center. Since we use the FLM only to suppress noise in the weighting coefficients, whether the parameters of that model are realistic or not is outside the scope of this study. A more computationally tractable PL fitting function may exist.