Sensing phases of water via nitrogen-vacancy centres in diamond

Ultra-thin layers of liquids on a surface behave differently from bulk liquids due to liquid-surface interactions. Some examples are significant changes in diffusion properties and the temperature at which the liquid-solid phase transition takes place. Indeed, molecular dynamics simulations suggest that thin layers of water on a diamond surface may remain solid even well above room temperature. However, because of the small volumes that are involved, it is exceedingly difficult to examine these phenomena experimentally with current technologies. In this context, shallow NV centres promise a highly sensitive tool for the investigation of magnetic signals emanating from liquids and solids that are deposited on the surface of a diamond. Moreover, NV centres are non-invasive sensors with extraordinary performance even at room-temperature. To that end, we present here a theoretical work, complemented with numerical evidence based on bosonization techniques, that predicts the measurable signal from a single NV centre when interacting with large spin baths in different configurations. In fact, by means of continuous dynamical decoupling, the polarization exchange between a single NV centre and the hydrogen nuclear spins from the water molecules is enhanced, leading to differences in the coherent dynamics of the NV centre that are interpreted as an unambiguous trace of the molecular structure. We therefore propose single NV centres as sensors capable to resolve structural water features at the nanoscale and even sensitive to phase transitions.


Master equation for Liquid Water
In this section we derive the master equation for the liquid water scenario that will lead to Eq. (2) in the main text. Starting from the general Hamiltonian in Eq.(1) of the main text, and under the assumption that for fast-moving particles the internuclear interaction averaged out, we can rewrite the Hamiltonian as In the last expression and for the sake of clarity, we have make use of the ladder operators to indicate the possible interactions. Also, we have renamed g i (t) ≡ 1 4 A i x (t) + iA i y (t) and h i (t) ≡ 1 2 A i z (t) and B i eff = B − 1 γ H 1 2 A i I i is an effective field created on the nuclear spins by the NV centre. First, we consider that we work in high-fields such that γ N |B| A i ∀i, and thus consider that all the nuclei rotate with identical Larmor frequency ω N .
Second, in order to treat such a Hamiltonian one may split it into a time independent part, H ω ≡ ΩS z + ∑ ω N I i z , and a time-dependent stochastic part The system will evolve according to the Von-Neumann master equation, This is a stochastic differential equation. Nonetheless, we want to look at the deterministic master equation, that is, averaged over all possible stochastic trajectories. There exist several methods for solving such an equation, in particular we follow 1 , which have been used previously in 2, 3 , obtaining accurate results. For solving the dynamics we go in a interaction picture with respect to H ω and split the Hamiltonian in a time dependent and time independent part obtaining where · represent the average over all possible stochastic trajectories. We have introduced the detuning ∆ = Ω − ω and the hyperfine coupling fluctuations around the average Also, it is assumed all the particles posses identical average properties, g i = g ∀i, h i = h ∀i, which is true for homogeneous motion. We note that in the assumption that ∆ ≈ 0, and Ω ∆, we can neglect fast rotating terms formH 0 obtaining which is a time independent Hamiltonian where all the fast-rotating terms have been neglected. Moreover, in our scenario g = 0.
Always in the interaction picture with respect to H ω , the master equation may be written aṡ In the last expression we have made use of the Liouville superoperator, L α (t)· ≡ −i H α , · . The presented master equation determines the behaviour ofρ(t) for a specific stochastic trajectory. Yet, we are interested in the average value of the density matrix over all possible stochastic trajectories since it is the solely measurable quantity. Thus, performing a cumulant expansion and following straightforwardly reference 1, 3 , when |H 0 |τ c 1 and G 2 τ c 1, which is true for moderate coupling and short τ c , the average solution of such a system iṡ where ρ (t) is the average of the total density matrix ρ(t) over all possible stochastic trajectories. This is an approximation up to second order in terms of G n i τ n−1 c . Up to this order, the effects of the stochastic motion can be split into three different contributions. A shift to the energy levels with the frequencies defined as A dissipator originated from flip-flop and flip-flip interactions And the last term which will cause pure dephasing on the NV centre We remark that for symmetry reasons, in our set-up it is verified g = 0. That means, there is not a net coherent coupling between NV centre and nuclei, allowing us to express the system state as ρ (t) = ρ NV (t) ⊗ ρ B . The evolution of the NV population n(t) ≡ 1 2 + Tr (S z ρ NV (t)), may be calculated obtaininġ With n B the population of a single nuclear spin. When thermal baths are considered, n B = 1/2 and we arrive to a final expressioṅ With the depolarization rate α(t) ≡ 1 4 N (γ x (∆,t) + γ x (2 Ω + ∆,t) + γ z (Ω,t)).

Equivalence between Master Eq. and dynamical equation for covariance matrix
Now we focus on the derivation of the evolution equation for the covariance matrix that will describe the evolution of the system when a mixture of phases is present. In fact, when the two phases coexist the description of the NV centre dynamics becomes more challenging since we need to combine the evolution given by Eq.(7) in the main text, with the depolarization created by a liquid bath Eq.(2) in the main text. While the former is a equation for bosons, since we have applied the HPA, the latter is a equation for spins. Since we are interested in obtain numerical results, we prefer to translate the equation obtained for moving spins to an equivalent expression for bosons.

2/4
We start tracing out the liquid bath in S.Eq.(9) in order to obtain a master equation solely for the NV For the sake of simplicity we have introduced the new auxiliary variables,Ω, γ + , γ − , we omit its explicit form for a moment since it is not relevant for our porpoise here. Under the influence of the former master equation, the population on the NV evolves aṡ which correspond to S.Eq.(13), with n B = 1 2 and γ + = γ − = 1 2 α(t). If now we look at a bosonic system, we have that a master equation of the forṁ leads to a equation for the population aṡ Therefore, the bosonic system will evolve equivalently to spin system as far as γ a † = γ + and γ a = 2 γ + + γ − . For n B = 1 2 we can again express these quantities in terms of the depolarization rate, γ a † = 1 2 α(t), γ a = 3 2 α(t). Thus, with appropriate rates, we can reproduce on a bosonic system the dynamics generated by a Lindblad dissipator on a spin system.
Provided now that for bosons we write the total Hamiltonian (NV and solid spins) as H, the master equation of the system iṡ This master equation leads to a evolution equation for the covariance matrix of the form of Eq.(8) of the main text. On the other hand, this reflects that a direct application of the HPA on the master equation does not always lead to faithful result. In our case, this forbids us to use the HPA for obtaining results in the case of fast moving nuclei.

NV as magnetometer
In this section we present the derivation of Eq.(9) of the main text. As we are principally interested in the solid part, we use the HPA and treat a system of bosons. Under the action of a high magnetic field gradient parallel to the NV centre quantization axis, the Larmor frequency of a nucleus at a position z i is just ω i = γ H (B 0 + λ z i ), where B 0 is a constant magnetic field and λ is the gradient strength λ ≡ dB/dz. The Hamiltonian for a bosonic system expressed in Eq.(6) of the main text, can be straightforwardly modified resulting in For a crystal such as Ice Ih, the 1 H nuclei are distributed in layers separated by certain distance d. Therefore, if a given layer seats at a distance z L = n d + z 0 from the NV centre, all the nuclei in this layer will precess with the same Larmor frequency ω i (nd + z 0 ). Moreover, they are detuned from nuclei belonging to adjacent layers,∆ = ω i (nd + z 0 ) − ω j ((n + 1)d + z 0 ) = γ H λ d. If our field gradient is high enough such that∆ ∑ z i =z L |g 2 i |, then setting Ω =ω(z L ), the effective Hamiltonian may be written as