Sensing Electrochemical Signals Using a Nitrogen-Vacancy Center in Diamond

Chemical sensors with high sensitivity that can be used under extreme conditions and can be miniaturized are of high interest in science and industry. The nitrogen-vacancy (NV) center in diamond is an ideal candidate as a nanosensor due to the long coherence time of its electron spin and its optical accessibility. In this theoretical work, we propose the use of an NV center to detect electrochemical signals emerging from an electrolyte solution, thus obtaining a concentration sensor. For this purpose, we propose the use of the inhomogeneous dephasing rate of the electron spin of the NV center (1/T2★) as a signal. We show that for a range of mean ionic concentrations in the bulk of the electrolyte solution, the electric field fluctuations produced by the diffusional fluctuations in the local concentration of ions result in dephasing rates that can be inferred from free induction decay measurements. Moreover, we show that for a range of concentrations, the electric field generated at the position of the NV center can be used to estimate the concentration of ions.


Introduction
Accurate detection of ionic concentrations in an electrolyte solution is of fundamental interest for process control and optimization in the chemical, pharmaceutical, and food industries. It is also important in environmental applications as a means to detect pollutant levels, as well as in fundamental studies in chemistry and biology. Examples include measuring the concentration of copper ions in a copper refinery or measuring the pH of a solution, which is defined as the negative logarithm of the concentration of hydrogen ions. In copper electro-refining processes, electrolysis is used to increase the purity of copper in combination with an electrolyte Cu 2+ /SO 2− 4 solution, and monitoring of the ionic concentration is required to control the purity of copper with high precision.
Carbon-based materials have been widely used as sensors [1][2][3][4]. Among them, the negatively charged nitrogen-vacancy (NV) center in diamond has attracted much attention. Diamond is a bio-compatible crystal that is resilient to extreme conditions. Moreover, the electron spin of the NV center can be prepared and read out through optical excitation [5]. The zero field splitting between the spin triplet sub-levels in the ground state of the NV center is sensitive to various physical parameters, such as crystal strain [6], temperature [7], electric [8] and magnetic fields [9][10][11][12]. The charge state of the NV center can also be used as a sensor. In Refs. [13,14], the effect of the pH of an electrolyte solution on the charge state of the NV center in a bulk diamond and in a functionalized nanodiamond were investigated.
In addition, the coherence time of the electron spin of the NV center is sensitive to its surrounding environment. The decoherence time of the NV center is used to estimate fluctuating magnetic fields from paramagnetic impurities [15], a liquid's diffusion coefficient [16], the drift velocity inside a microfluid channel [17], and the concentration of ions in a cell membrane [18]. In Ref. [19], the decoherence time of the NV center was monitored during a change in the pH of an electrolyte solution.
Here, we consider an electrolyte-diamond interface, as shown in Figure 1. We propose the use of the dephasing rate of the electron spin of an NV center in diamond in order to estimate the concentration of ions in the electrolyte solution. The rationale behind this choice, as we will explain in detail in the following sections, is that the diffusional fluctuations in the local concentration of ions in the electrolyte solution generate an electrical noise signal, which results in an additional dephasing rate, 1/T 2 , that competes with the intrinsic dephasing rate of the NV center. As we will show, for a range of mean ionic concentrations in the bulk of the electrolyte, c b > 0.04 mol/m 3 , the induced 1/T 2 is larger than 10 kHz, and for c b > 100 mol/m 3 , it is larger than 300 kHz. Therefore, depending on the intrinsic 1/T 2 of the electron spin of the NV center, the range of c b that can be estimated with this method can be larger than 0.04 mol/m 3 . We also show that for a lower range of ion concentrations, c b < 0.1 mol/m 3 , the gradient of the electric field generated at the position of the NV center is large enough that it may be used to estimate the concentration of ions.

Diamond
Electrolyte NV Figure 1. The proposed experimental setup, which contains a Cu 2+ /SO 2− 4 electrolyte solution in contact with the surface of a bulk diamond crystal with a nitrogen-vacancy (NV) center (red arrow). The upper surface of the diamond is taken to be at z = 0 with potential φ 0 , while z = ∆ is far into the bulk of the electrolyte solution with potential φ b,e . The potential in the opposite surface of the diamond is shown with φ b,d , which is taken to be zero in our simulations.

Electrolyte Solution
In this section, we present a brief analysis of the basic physico-chemical mechanisms governing the diffusion of ionic chemical species in a liquid electrolyte solution, which is displayed as the sample to be analyzed in the proposed experimental setup (Figure 1). Since the ionic species and their concentration profiles determine a local charge density, by solving the Poisson equation, we obtain the electric potential and the electric field inside the electrolyte solution. Through a further statistical analysis, we then obtain mathematical expressions for the fluctuations of the electric field inside the electrolyte solution.

Diffusion of Ionic Species in Electrolyte Solutions
A typical model for the flux N s (mol m −2 s −1 ) of ionic species in a liquid electrolyte solution, taking into account the effect of the concentration gradient, Fick's law, and the effect of an electric field, is given by [20] In this equation, known as the Nernst-Planck equation, D is the diffusion constant, z s is the valence of each ion species in the electrolyte solution, and c s is the concentration of each species. The Faraday constant, F = 96, 485.3365 C/mol, is a measure of charge per mole of substance and is the product of Avogadro's number and the charge of an electron, R = 8.314 J/(mol K) is the universal constant of gases, T is the temperature in K, and φ is the electric potential in volts. By requiring mass balance, we obtain the following linear partial differential equation [20]: The above equation accounts for how the diffusion process involves electrical forces on each ionic species in an electrolyte solution. To account for the electric field due to the presence of the ionic species, we consider the Poisson equation: Here, e is the electric permittivity of the electrolyte solution (approximately the same as in water for concentrations below 10 3 (mol/m 3 ) [21]) and ρ(x) is the local charge density determined by the ionic concentrations: Due to the small thickness of the fluid layer of the sample compared to its transverse dimensions, we assume that the system presents a concentration gradient and potential gradient only along the direction z normal to the interface (see Figure 1). This approximation is equivalent to considering an infinitely large interface surface between the liquid solution and the diamond crystal. Therefore, we can write N s =ẑN s with the boundary condition N s | z=0 = 0, which expresses the fact that the ionic species in the liquid solution cannot penetrate into the solid crystal. As a result, in the steady state, i.e., ∂c s /∂t = 0, from Equation (2), we obtain N s = 0 for s = ±. With the previous considerations, the simplified system of coupled non-linear differential equations to be solved is given by Taking the integral of Equation (5) from z to ∆, with ∆ being the distance from the interface far into the bulk of the solution, we obtain where c b,s is the bulk concentration of each species. Taking into account that the electrolyte should be electrically neutral in the bulk, due to charge conservation, we obtain the electric potential as (see Appendix A): Here is the potential in the bulk of the electrolyte (see Figure 1). We have defined κ, the inverse screening length, as Note that Equation (8) is valid for κ∆ 1 (see Appendix A). From this potential, we can calculate the electric field, which, on the surface of the diamond inside the electrolyte, is given by Here, the subscript eq emphasizes that this corresponds to the local electric field in thermal equilibrium.

Electric Field Fluctuations at the Interface between the Solution and the Diamond
Due to thermal noise, the concentration of ions will have small fluctuations around its equilibrium value, i.e., δc(z, t) = c(z, t) − c eq (z). These fluctuations in the concentration will induce fluctuations in the electric field, i.e., δE(z, t) = E(z, t) − E eq (z). The electric field fluctuations are directly linked to the fluctuations in concentration through the Poisson equation: We can integrate this equation over z in the range [z, ∆]. For κ∆ 1, we have δE(z = ∆, t) ∼ 0. On the other hand, as is shown in Appendix B, the correlation of the fluctuations of concentration is given by Therefore, the correlation of the fluctuations of the electric field at the surface of the diamond from the liquid solution side, i.e., at z = 0 + , is obtained as Using the equilibrium concentration given in Equation (7), in the limiting case κ∆ 1, the electric field fluctuations at z = 0 + can be simplified to It is clear that the correlation of the fluctuations is directly proportional to the concentration in the bulk c b . Figure  In the following section, we obtain the fluctuations of the electric field at the position of the NV center in the diamond. For doing so, we first obtain the electric potential and the electric field inside the diamond.

Diamond and the NV Center
The electron spin of the NV center in the diamond is sensitive to the electric field and to its fluctuations. The fluctuations of the electric field result in dephasing of the NV electron spin. In this section, we consider an NV center in a bulk diamond in contact with the electrolyte, as depicted in Figure 1. The electric field inside the diamond is screened by the dielectric response of the diamond and ionization of donors and acceptors inside the diamond [22]. We first find the electric field and its fluctuations, resulting from the concentration of ions in the electrolyte and their fluctuations, at the position of the NV center. We then find the dephasing rate of the electron spin of the NV center as a result of the fluctuations of the electric field. Finally, we show that for a range of c b , the gradient of the electric field induced at the position of the NV center is large enough to be resolved by a Ramsey measurement.

Potential and Electric Field Inside the Diamond
As is shown in Figure 1, we consider a bulk diamond in contact with the electrolyte with its interface at z = 0. We assume that the bulk diamond is implanted with nitrogen ions, forming mainly substitutional nitrogen, N s , (≈96%) and NV defects (≈4%), and that the concentration of vacancy defects is negligible. The substitutional nitrogen and the negatively charged NV center, NV − , act as donors with ionization energies E N d = 1.7 eV and E NV − d = 2.7 − 2.8 eV, respectively [23]. The positively charged NV center, NV + , acts as an acceptor with ionization energy E NV + a = 0.9 − 1.1 eV (see Figure 3) [23]. We take the areal density of implanted nitrogens to be D s = 10 12 cm −2 . Thus, the volume density of substitutional nitrogen and NV are given by N N d = 0.96D s /d max and N NV a = 0.04D s /d max , where we take d max = 14 nm as the maximum implantation depth of nitrogen ions.
The charge density inside the diamond is given by the density of electrons n(z), density of holes p(z), and density of ionized donors and ionized acceptors where e > 0 is the magnitude of the electron charge. Here, we have used the subscript d for the charge density to avoid confusion with the charge density inside the electrolyte solution. Assuming that we can use the Boltzmann approximation, the density of electrons and holes in the presence of a potential φ(z) can be written as [24] n(z) = N c exp Here, k is the Boltzmann constant, T is the temperature, and N c and N v are the effective density of states in the conduction and valence bands, respectively, given as and in whichh is the reduced Planck constant, E v is the valence band maximum, E c is the conduction band minimum, and m * n = 0.57m 0 and m * p = 0.8m 0 are the effective masses of conduction and valence bands, respectively [23]. The corresponding densities of ionized donors and acceptors are given by To find the potential φ(z) inside the diamond, we use the Poisson equation where d = 5.8 0 is the dielectric constant of the diamond, with 0 being the vacuum permittivity. At the interface of the electrolyte and diamond, we impose the continuity of the electric potential and displacement, i.e., Here, d and e stand for diamond and electrolyte, and e is the dielectric constant of the electrolyte. Solving the Poisson equation, we obtain the electric field as Here, we have assumed that ρ d depends on φ and does not depend explicitly on z. Integrating the above equation gives the potential φ at the position z as The above equation has to be solved numerically for a specific value of z. Replacing the obtained potential in Equation (23), we obtain the electric field at the position z. The electric field and potential inside the bulk diamond are shown in Figure 4 for c b = 1 mol/m 3 . In this figure, the potential in the bulk of the electrolyte is taken to be φ b,e = 1.5 V at ∆ = 1 mm. The potential on the surface of the diamond is then found by numerically solving Equation (24) by setting the potential inside the bulk of the diamond at 100 nm depth to φ bd = 0 V. Note that we have taken the direction of the z axis towards the bulk of the electrolyte. The other parameters are taken as T = 298 K, A = 4 mm 2 , and D = 2.3 × 10 −9 m 2 /s.
3 × 10 −9 m 2 /s, and T = 298 K. Figure 5 shows the electric field at the position of NV at 10 nm as a function of c b , as calculated from the system of equations comprising Equations (15)- (24). The inset of this figure shows that for values of c b < 0.1 mol/m 3 , the gradient of the electric field is higher. In Section 2.2.3, we will show that this feature may be used to sense the concentration of ions in the electrolyte solution through the Stark effect on the electron spin of the NV center.

NV Center Inhomogenous Dephasing Rate
The electron spin of the NV center has a spin triplet ground state. Its ground state Hamiltonian in the presence of a magnetic and an electric field is given by (h = 1) [25] where D = 2.87 GHz is the ground state zero field splitting, γ e = 2.8 MHz/G is the electron gyromagnetic ratio, and H E i are the terms due to the electric field that cause transitions between the spin states with the difference in the spin projection ∆m s = i, Here, E NV i are the components of the electric field in the NV reference frame. The coupling parameters are experimentally found to be d || = 0.35 Hz cm/V and d ⊥ = 17 Hz cm/V [26]. On the other hand, the coefficient d ⊥ is expected to be of the same order of magnitude as d ⊥ [27]. We assume that the NV symmetry axis sits in the x, z plane of the laboratory frame. Therefore, for an electric field that is in the z direction of the lab frame, in the NV frame, we have Due to the large value of the zero field splitting D and the small values of d || and d ⊥ , we can neglect the H E 0 and H E 1 terms in the Hamiltonian. Moreover, for weak magnetic fields, B D/γ e , we can neglect the perpendicular component of the magnetic field. Therefore, the eigenstates of the Hamiltonian in terms of the eigenstates of S z are given by |0 and [27] Here, The energy splitting between the |0 and |± states is given by With our choice of the NV frame, we have E NV y = 0 and, therefore, It is shown in Ref. [28] that for nonzero values of B z , the electron spin is protected from electric field noise. Therefore, to be able to detect electric field fluctuations, we consider B z = 0. For the case of B z = 0, the fluctuations in ν can be written as Assuming that δξ ⊥ , δβ z << ξ ⊥ , we can expand the square root and obtain Therefore, where we have assumed < (δβ z ) 2 δξ ⊥ >= 0. The fluctuations in the energy splitting ν result in dephasing of the electron spin of the NV center, which can be measured through free induction decay measurement [29]. For this purpose, the spin is prepared in |0 , followed by a Ramsey sequence that consists of a π/2 microwave pulse, a free evolution for time τ, and another π/2 microwave pulse followed by a projective measurement on the electron spin (see Figure 6). The probability of detecting the electron spin in the state |0 after the Ramsey sequence is given by  Figure 6. The sequence for free induction decay measurement. The electron spin is initially prepared in |0 using a green laser. With the application of a π/2 microwave (MW) pulse, the electron is prepared in the superposition 1 √ 2 (|0 + |1 ). The electron spin then goes under a free evolution for time τ. The accumulated phase during the free evolution is projected to the populations using another π/2 MW pulse. At the end, the state of the electron is read out using a green laser.
The exponential decay factor determines the dephasing rate of the NV center. To calculate δψ 2 , we write We assume that the fluctuations are stationary, i.e, their correlation is a function of the time difference, δν(t 1 )δν(t 2 ) = f (t 1 − t 2 ). We define τ = t − t and T = (t + t )/2. Thus, From Equation (35), the correlation function of the fluctuations in ν due to the electric field can be written as The electric field at the position of the NV center depends on the electric field at the interface between the electrolyte liquid solution and the diamond, i.e., at z = 0 + in our coordinate system (see Figure 1). Therefore, the fluctuations of the electric field at the position of the NV center can be written as As a result, We numerically calculate the partial derivative of E NV with respect to E e (0) using Equation (23). Note that φ, which appears as the upper limit of the integral, is also a function of E e . For the correlation of the fluctuations of the electric field on the surface, we use Equation (14). For ∆ ≈ 1 mm and for times up to a few miliseconds, we have ∆/ √ 4Dt >> 1. Therefore, the error functions can be approximated by 1 and the second term in the summation can be approximated by zero, i.e., (44) As a result, the exponential decay factor in the free induction decay signal, Equation (38), scales as t 2 . The factor of t 2 in the exponential gives 1/(T 2 ) 2 , i.e., In Figure 7, we have plotted 1/T 2 as a function of c b for an NV center at a depth of 10 nm in the diamond. For the numerical simulations (solid line), we considered φ b,e = 1.5 V and φ b,d = 0 V. The solid line was obtained by numerically finding the potential on the surface of the diamond using Equation (24). The dashed red line is a fit to the simulations.  For an NV center in a bulk diamond with an abundance of 13 C isotopes ranging from very low to natural, 1/T 2 due to the surrounding nuclear spins is of the order of a few kHz to a few hundred kHz, respectively [30][31][32]. To be able to detect fluctuations of the concentration, the corresponding 1/T 2 should be larger than the intrinsic 1/T 2 of the NV center. For c b > 100 mol/m 3 , the electric field fluctuations produced by the diffusional fluctuations in the local concentration of ions result in 1/T 2 > 300 kHz, and for c b > 0.04 mol/m 3 , they result in 1/T 2 > 10 kHz. Therefore, the range of concentrations that can be estimated with this method depends on the properties of the diamond sample. We note that, for an NV that is positioned at a smaller distance from the surface of the diamond, the induced electric field at the position of the NV is larger and the induced T 2 is larger.
In the next section, we show that for a suitable range of concentrations, the Stark effect may also be used to estimate c b .

Electric Field Sensing
To be able to measure the electric field, the electron spin needs to be prepared in a state that is susceptible to the term H E 2 of the Hamiltonian (Equation (28)) [8,33]. As the H E 2 term causes a transition between the m s = ±1 states, such a state is a superposition of |m s = ±1 states. A perpendicular magnetic field much smaller than the zero field splitting, D, in the absence of an axial magnetic field, can prepare the states [33] where ϕ B is the azimuthal angle of the magnetic field. The energy shifts of these states due to H E 2 are where we have taken E NV y = 0. As can be seen from Figure 5, for c b < 0.1 (mol/m 3 ), a change of one order of magnitude in c b results in a change of about 230 kV/m in the electric field at the position of the NV center, E NV x . This change in the electric field results in an energy shift of about 39 kHz for ϕ B = 0. This change in the electric field can be resolved in a Ramsey sequence by preparing the NV center in the superposition state of |0 and either of the |± states and measuring the energy shift.

Conclusions
We considered an NV center in a bulk diamond in contact with an electrolyte. The fluctuations in the concentration of ions induce a fluctuating electric field at the position of the NV center. We first showed that the dephasing rate of the electron spin of the NV center, 1/T 2 , can be used to estimate the concentration of ions in the electrolyte for a range of concentrations c b > 0.04 mol/m 3 , depending on the intrinsic 1/T 2 of the electron spin of the NV. For this range of c b , the induced 1/T 2 resulting from the fluctuations of the concentration of charged species is larger than 10 kHz. We also showed that for c b < 0.1 mol/m 3 , the gradient of the electric field induced at the position of the NV center is large enough to be resolved through Ramsey spectroscopy for a change of one order of magnitude in c b . Therefore, through the estimation of the electric field, c b can be estimated.
Using a dynamical decoupling sequence, such as the spin echo or Carr-Purcell-Meiboom-Gill sequence, the 1/T 2 of the electron spin of the NV center due to surrounding nuclear spins can be decreased to a few kilohertz. In that case, the NV center would be sensitive to the fluctuations in the concentration over a wider range of c b .

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: NV nitrogen vacancy

Appendix A. Electric Potential and Electric Field Inside the Electrolyte
In this appendix, we give details of the derivation of the electric potential and the electric field inside the electrolyte. First, we define a dimensionless potential as where φ(∆) is the potential in the bulk of the electrolyte (see Figure 1). In the bulk, the solution should be electrically neutral due to charge conservation, i.e., ∑ s=± z s c b,s = 0, with z s being the valence of the ion species and c b,s being the bulk concentration of the species. Assuming z + = −z − = z s , we can write c b,+ = c b,− = c b . Therefore, in terms of the dimensionless potential, the Poisson equation given in Equation (6) can be written as We first consider the limit where c b → 0. Since Φ is proportional to c b , because the second derivative of Φ is proportional to c b , in the limit of c b → 0, we can linearize the above equation and obtain where κ is defined in Equation (9). The solution of this differential equation can be written as The constants A and B can be obtained by considering the boundary conditions Φ(0) = z s F V 0 /(RT) and Φ(∆) = 0. Thus, we obtain Therefore, the electric field can be written as In the limit of c b → 0, i.e., κ → 0, we obtain This is the expected electric field between two parallel plates with the distance ∆ and the potential difference V 0 .

Appendix B. Thermal Fluctuations in the Local Concentration of Chemical Species
In this appendix, we derive the fluctuations of the concentration of ions in the electrolyte. From the statistical point of view, the local concentration of chemical species is the average number of particles per unit volume n(x, t) , where n(x, t) is a random variable subject to stochastic dynamics due to Brownian motion. The corresponding average then satisfies the diffusion equation where D is the diffusion constant. This equation possesses a Green function G(x, t; x 0 , t 0 ), which represents the transition probability for a particle to diffuse from a point x 0 at time t 0 to a point x at time t. This transition probability satisfies the differential equation with the initial condition Equation (A17) has an explicit solution (in d dimensions) given by Following Ref. [34], we define the two-particle correlator It can be shown that this two-particle correlator satisfies the following equation: We now define the particle number fluctuations δn(x, t) = n(x, t) − n(x, t) , which, by definition, satisfy the following properties: δn(x, t) = 0, δn(x 1 , t 1 )δn(x 2 , t 2 ) = n(x 1 , t)n(x 2 , t) − n(x 1 , t) n(x 2 , t) = [n(x 1 , t)n(x 2 , t)] + δ( Combining Equations (A21) and (A22), after some algebraic manipulations, one obtains the differential equation for the correlator of fluctuations in the particle number The steady-state (equilibrium) solution to this equation is We can now consider the time evolution of the fluctuation correlator by using the conditional probability (Green's function) as follows: Inserting Green's function from Equation (A19) and calculating the integral, we obtain Here, Θ(t) is the step function. This equation clearly shows an exponential decay of the correlation of the fluctuations with the relative distance |x 2 − x 1 |. For the particular case of a three-dimensional system that, due to physical constraints, is subjected to concentration gradients and fluctuations along only one dimension (for instance, the z-direction), i.e., n(x, t) = n(z, t), the particle concentration is still defined per unit volume, i.e., in units of 1/m 3 . In this case, the equilibrium correlator, Equation (A24), must be modified as where the Dirac delta function is one-dimensional (with units of 1/m) and A is the area transversal to the direction of the gradients (the surface area of the bulk diamond). With this argument, Equation (A25) is now given by a one-dimensional integral with Green's function for d = 1 as δn(z 1 , t)δn(z 2 , 0) = dz 3 G(z 2 , t; z 3 , 0) δn(z 3 )δn(z 1 ) eq = n(z 1 ) eq A(4πD t) 1/2 exp − (z 2 − z 1 ) 2 4Dt Θ(t).
From the above equation, the concentration fluctuations on a molar basis for two species s, s = 1, 2 correspond to δc s (z 1 , t)δc s (z 2 , 0) = δ s,s c eq s (z 1 ) N A A(4πD s t) 1/2 exp − where N A = 6.02 × 10 23 (particles/mole) is Avogadro's number, and D s the diffusion constant for each species.

Appendix C. Sensitivity
In this appendix, we derive the sensitivity of the proposed scheme in estimating c b . The sensitivity, η, is defined as the square root of the product of the variance of c b and the measurement time As c b is estimated by estimating the dephasing time T 2 , from the error propagation formula, we can write (∆c b ) 2 = (∆T 2 ) 2 |∂T 2 /∂c b | 2 . (A31) From the fit to the numerical simulations of 1/T 2 versus c b (see Figure 7), we have , where A ≈ 39, 295 Hz (mol/m 3 ) −B and B ≈ 0.417 (unitless). On the other hand, T 2 is estimated by measuring the electron spin of the NV center through the detection of its photoluminescence. Thus, the observable that is used to estimate T 2 can be taken as M = a|0 0| + b|1 1| [35]. The parameters a and b are random variables that determine the number of photons collected during the readout interval (about 300 ns) if the spin state is in the |0 and |1 states, respectively, with their averages over many measurements given by α = a and β = b . The variance of T 2 can then be written as where (∆M) 2 = M 2 avg − M 2 avg , with the average being over many measurements. For a free induction decay measurement with evolution time t, we obtain [35] and where ψ is the accumulated phase during the evolution time t given in Equation (37). Note that, in Section 2.2.2, we found that the exponential decay factor scales with t 2 . Taking the derivative of Equation (A33) and substituting into Equation (A32), we obtain (∆T 2 ) 2 = 15(T 2 ) 6 e 2(t/T 2 ) 2 2α t 4 cos 2 (ψ) , where we have used α = 3β/2, and since α ≈ 0.03, therefore, 1/α >> 1, and we only considered the first term in Equation (A34) [35]. Thus, the sensitivity is obtained as where, to obtain a limit, we have taken cos(ψ) ≈ 1. As an example, we obtain η for c b = 10 mol/m 3 , for which we have found T 2 ≈ 10 µs. Taking t = 10 µs, we obtain η = 3.27 mol m −3 Hz −1/2 . This sensitivity can be improved by a protocol in which t is changed during the measurement repetitions. An extra enhancement may be achieved by optimizing the measurement sequence using a Bayesian approach [12]. In addition, the sensitivity can also be improved by increasing the photon collection efficiency, which can be achieved by, for instance, coupling the NV to a photonic waveguide [36] or improving the spin readout of the NV center through spin-to-charge conversion [37].