409 lines
16 KiB
Rust
409 lines
16 KiB
Rust
//! NV-ensemble sensor model — Pass 4 of the implementation plan.
|
||
//!
|
||
//! Linear-readout proxy for ODMR ensemble magnetometry. Per plan §2.3, the
|
||
//! full Hamiltonian + Lindblad solver is *out of scope* (plan §6); we
|
||
//! implement the leading-order ensemble sensitivity formula that Barry et al.
|
||
//! *Rev. Mod. Phys.* 92, 015004 (2020) §III.A validates as adequate for
|
||
//! ensemble magnetometers operated in the linear regime.
|
||
//!
|
||
//! # What this module models
|
||
//!
|
||
//! - **ODMR transition**: `ν± = D ± γ_e |B_∥|` per Doherty 2013 §3.
|
||
//! - **Lorentzian lineshape** at FWHM Γ ≈ 1 MHz (Barry 2020 Fig. 4).
|
||
//! - **T₂ decay envelope**: `exp(−t/T₂)` (Jarmola PRL 108, 2012; Barry 2020).
|
||
//! - **Shot-noise floor**: `δB ∝ 1/(γ_e · C · √(N · t · T₂*))` —
|
||
//! leading-order projection-noise-limited sensitivity (Barry 2020 Eq. 35).
|
||
//! - **4-axis crystallographic projection**: `[1,1,1]/√3`, `[1,-1,-1]/√3`,
|
||
//! `[-1,1,-1]/√3`, `[-1,-1,1]/√3` (Doherty 2013 §3).
|
||
//! - **Least-squares 3-vector recovery** from the 4 projection scalars.
|
||
//!
|
||
//! # What this module does NOT model
|
||
//!
|
||
//! Strain broadening, hyperfine coupling, magnetic-resonance saturation,
|
||
//! pulsed dynamical decoupling, photon shot noise vs spin projection noise
|
||
//! distinction, microwave power broadening. These are flagged in plan §6 as
|
||
//! out-of-scope; if any matters for a future use case, the simulator
|
||
//! escalates to the QuTiP path.
|
||
//!
|
||
//! # Determinism
|
||
//!
|
||
//! Shot noise is sampled from a ChaCha20 PRNG seeded explicitly per `sample`
|
||
//! call. Same `(seed, B_in, dt)` produces byte-identical [`NvReading`] —
|
||
//! the foundation of the proof-bundle commitment in plan §5.
|
||
|
||
use crate::{D_GS, GAMMA_E};
|
||
use rand::SeedableRng;
|
||
use rand_chacha::ChaCha20Rng;
|
||
use serde::{Deserialize, Serialize};
|
||
|
||
/// Default ODMR linewidth (FWHM, Hz). 1 MHz typical for COTS bulk diamond
|
||
/// (Barry 2020 Fig. 4). Strain-free lab samples can be narrower; CW-ODMR
|
||
/// power broadening can widen this in production hardware.
|
||
pub const DEFAULT_GAMMA_FWHM_HZ: f64 = 1.0e6;
|
||
|
||
/// Default T₁ (s). 5 ms at room temperature (Jarmola PRL 108, 2012;
|
||
/// Barry 2020 Table III).
|
||
pub const DEFAULT_T1_S: f64 = 5.0e-3;
|
||
|
||
/// Default T₂ (s). 1 µs for COTS bulk (Barry 2020 Table III).
|
||
pub const DEFAULT_T2_S: f64 = 1.0e-6;
|
||
|
||
/// Default T₂* (s). 200 ns for COTS bulk (Barry 2020 Table III).
|
||
pub const DEFAULT_T2_STAR_S: f64 = 200.0e-9;
|
||
|
||
/// Default ODMR contrast `C`. 0.03 = 3% for COTS bulk diamond
|
||
/// (Barry 2020 Table III).
|
||
pub const DEFAULT_CONTRAST: f64 = 0.03;
|
||
|
||
/// Default sensing spin count `N`. ~10¹² spins per ~1 mm³ DNV-B-class
|
||
/// diamond (Barry 2020 §IV.A).
|
||
pub const DEFAULT_N_SPINS: f64 = 1.0e12;
|
||
|
||
/// NV crystallographic axes (4 of them, normalised). Doherty 2013 §3.
|
||
/// Tetrahedral 〈111〉 family in the diamond lattice.
|
||
pub fn nv_axes() -> [[f64; 3]; 4] {
|
||
let s = 1.0 / 3.0_f64.sqrt();
|
||
[[s, s, s], [s, -s, -s], [-s, s, -s], [-s, -s, s]]
|
||
}
|
||
|
||
/// Sensor configuration. All defaults match plan §2.3 / Barry 2020 Table III
|
||
/// for COTS-grade bulk diamond at room temperature.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
|
||
pub struct NvSensorConfig {
|
||
/// ODMR FWHM (Hz). Default 1 MHz.
|
||
pub gamma_fwhm_hz: f64,
|
||
/// T₁ (s). Default 5 ms.
|
||
pub t1_s: f64,
|
||
/// T₂ (s). Default 1 µs.
|
||
pub t2_s: f64,
|
||
/// T₂* (s). Default 200 ns.
|
||
pub t2_star_s: f64,
|
||
/// ODMR contrast `C`. Default 0.03.
|
||
pub contrast: f64,
|
||
/// Sensing spin count `N`. Default 1e12.
|
||
pub n_spins: f64,
|
||
/// Disable shot noise (analytic mode). Default `false`.
|
||
pub shot_noise_disabled: bool,
|
||
}
|
||
|
||
impl Default for NvSensorConfig {
|
||
fn default() -> Self {
|
||
Self {
|
||
gamma_fwhm_hz: DEFAULT_GAMMA_FWHM_HZ,
|
||
t1_s: DEFAULT_T1_S,
|
||
t2_s: DEFAULT_T2_S,
|
||
t2_star_s: DEFAULT_T2_STAR_S,
|
||
contrast: DEFAULT_CONTRAST,
|
||
n_spins: DEFAULT_N_SPINS,
|
||
shot_noise_disabled: false,
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Output of one sensor sample.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
|
||
pub struct NvReading {
|
||
/// Recovered 3-vector field (T) — LSQ inversion of 4 noisy axis
|
||
/// projections back to xyz.
|
||
pub b_recovered: [f64; 3],
|
||
/// Per-axis 1σ noise estimate (T).
|
||
pub sigma_per_axis: [f64; 3],
|
||
/// Shot-noise floor for this integration window (T/√Hz).
|
||
pub noise_floor_t_sqrt_hz: f64,
|
||
/// Effective ODMR transition frequencies (Hz) for the higher branch
|
||
/// `ν+ = D + γ_e · |B_∥|` of each NV axis. Useful for downstream lockin
|
||
/// demod cross-checks; not required by the basic pipeline.
|
||
pub odmr_nu_plus_hz: [f64; 4],
|
||
}
|
||
|
||
/// NV-ensemble sensor.
|
||
#[derive(Debug, Clone, Copy)]
|
||
pub struct NvSensor {
|
||
/// Active configuration.
|
||
pub config: NvSensorConfig,
|
||
}
|
||
|
||
impl NvSensor {
|
||
/// Construct a sensor with the supplied config.
|
||
pub fn new(config: NvSensorConfig) -> Self {
|
||
Self { config }
|
||
}
|
||
|
||
/// Construct a sensor with COTS-grade defaults (Barry 2020 Table III).
|
||
pub fn cots_defaults() -> Self {
|
||
Self::new(NvSensorConfig::default())
|
||
}
|
||
|
||
/// Lorentzian normalised at peak: `L(δν) = (Γ/2)² / [(δν)² + (Γ/2)²]`,
|
||
/// returning 1.0 on resonance and falling to 0.5 at the half-width.
|
||
/// `delta_nu_hz` is the offset from line centre.
|
||
pub fn lorentzian(&self, delta_nu_hz: f64) -> f64 {
|
||
let half = self.config.gamma_fwhm_hz * 0.5;
|
||
let half_sq = half * half;
|
||
half_sq / (delta_nu_hz * delta_nu_hz + half_sq)
|
||
}
|
||
|
||
/// T₂ decay envelope: `exp(-t/T₂)`. Used to model coherence loss at
|
||
/// long integration times.
|
||
pub fn t2_envelope(&self, t_s: f64) -> f64 {
|
||
if t_s <= 0.0 {
|
||
return 1.0;
|
||
}
|
||
(-t_s / self.config.t2_s).exp()
|
||
}
|
||
|
||
/// Photon-shot-noise-limited sensitivity floor for the chosen
|
||
/// integration time. Plan §2.3: `δB ∝ 1/(γ_e · C · √(N · t · T₂*))`.
|
||
/// Returns T/√Hz at the BW=1 Hz reference; multiply by √BW to get the
|
||
/// per-sample noise σ in T.
|
||
pub fn shot_noise_floor_t_sqrt_hz(&self, integration_s: f64) -> f64 {
|
||
let t = integration_s.max(self.config.t2_star_s);
|
||
let denom = GAMMA_E
|
||
* self.config.contrast
|
||
* (self.config.n_spins * t * self.config.t2_star_s).sqrt();
|
||
if denom <= 0.0 {
|
||
f64::INFINITY
|
||
} else {
|
||
1.0 / denom
|
||
}
|
||
}
|
||
|
||
/// Sample the sensor — projects `b_in` onto each of the 4 NV axes,
|
||
/// applies shot noise, and recovers an LSQ 3-vector estimate. `dt`
|
||
/// is the integration time in seconds. `seed` makes the noise
|
||
/// reproducible: same `(b_in, dt, seed)` ⇒ byte-identical output.
|
||
pub fn sample(&self, b_in: [f64; 3], dt: f64, seed: u64) -> NvReading {
|
||
let axes = nv_axes();
|
||
let noise_floor = self.shot_noise_floor_t_sqrt_hz(dt);
|
||
// σ for one sample with this integration window: noise_floor
|
||
// is in T/√Hz at BW=1Hz; per-sample bandwidth is 1/(2·dt) so
|
||
// σ = noise_floor × √(BW). For dt-integrated samples we use
|
||
// BW = 1/dt as the conservative noise envelope.
|
||
let sigma = if self.config.shot_noise_disabled {
|
||
0.0
|
||
} else {
|
||
noise_floor * (1.0 / dt.max(1e-12)).sqrt()
|
||
};
|
||
|
||
let mut rng = ChaCha20Rng::seed_from_u64(seed);
|
||
let mut projections = [0.0_f64; 4];
|
||
let mut nu_plus = [0.0_f64; 4];
|
||
for (i, axis) in axes.iter().enumerate() {
|
||
let b_par = b_in[0] * axis[0] + b_in[1] * axis[1] + b_in[2] * axis[2];
|
||
// Shot noise on the projection.
|
||
let noise = if sigma > 0.0 {
|
||
sample_normal(&mut rng) * sigma
|
||
} else {
|
||
0.0
|
||
};
|
||
projections[i] = b_par + noise;
|
||
nu_plus[i] = D_GS + GAMMA_E * b_par.abs();
|
||
}
|
||
|
||
// LSQ inversion: B_xyz = (Aᵀ A)⁻¹ Aᵀ p, where A is the 4×3 matrix of
|
||
// axis vectors. Closed-form for the regular tetrahedron 〈111〉/√3:
|
||
// (Aᵀ A) = (4/3) I, so B_xyz = (3/4) Aᵀ p.
|
||
let mut b_recovered = [0.0_f64; 3];
|
||
for k in 0..3 {
|
||
let mut acc = 0.0;
|
||
for (i, axis) in axes.iter().enumerate() {
|
||
acc += axis[k] * projections[i];
|
||
}
|
||
b_recovered[k] = (3.0 / 4.0) * acc;
|
||
}
|
||
|
||
let sigma_per_axis = [sigma; 3];
|
||
|
||
NvReading {
|
||
b_recovered,
|
||
sigma_per_axis,
|
||
noise_floor_t_sqrt_hz: noise_floor,
|
||
odmr_nu_plus_hz: nu_plus,
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Box–Muller normal sample from a `ChaCha20Rng` source. Avoids pulling in
|
||
/// `rand_distr` for one function. Returns standard normal `~ N(0, 1)`.
|
||
fn sample_normal(rng: &mut ChaCha20Rng) -> f64 {
|
||
use rand::Rng;
|
||
// Two independent uniforms in (0, 1].
|
||
let u1: f64 = rng.gen_range(f64::EPSILON..=1.0);
|
||
let u2: f64 = rng.gen_range(f64::EPSILON..=1.0);
|
||
(-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use approx::assert_relative_eq;
|
||
|
||
#[test]
|
||
fn lorentzian_fwhm_within_5_percent() {
|
||
// Plan §3 Pass 4: FWHM = 1.0 ± 0.05 MHz. The half-width offset
|
||
// returns exactly 0.5 by construction; we check the documented
|
||
// value matches the config.
|
||
let s = NvSensor::cots_defaults();
|
||
let half = s.config.gamma_fwhm_hz / 2.0;
|
||
let on = s.lorentzian(0.0);
|
||
let at_half = s.lorentzian(half);
|
||
assert_relative_eq!(on, 1.0, max_relative = 1e-12);
|
||
assert_relative_eq!(at_half, 0.5, max_relative = 1e-12);
|
||
let nominal = 1.0e6;
|
||
assert!(
|
||
(s.config.gamma_fwhm_hz - nominal).abs() / nominal <= 0.05,
|
||
"default FWHM differs from 1 MHz nominal by > 5%"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn shot_noise_scales_as_one_over_sqrt_t_over_5_decades() {
|
||
// δB ∝ 1/√t per Barry 2020 Eq. 35. Sample 5 decades of integration
|
||
// and check that doubling t reduces the floor by √2.
|
||
let s = NvSensor::cots_defaults();
|
||
let mut prev: f64 = 0.0;
|
||
let mut measured_ratios: Vec<f64> = Vec::new();
|
||
for d in 0..6 {
|
||
// 1 µs, 10 µs, 100 µs, 1 ms, 10 ms, 100 ms
|
||
let t = 1.0e-6 * 10.0_f64.powi(d);
|
||
let floor = s.shot_noise_floor_t_sqrt_hz(t);
|
||
assert!(floor.is_finite() && floor > 0.0);
|
||
if d > 0 {
|
||
// Each 10× t step should drop the floor by √10 ≈ 3.162.
|
||
let ratio = prev / floor;
|
||
measured_ratios.push(ratio);
|
||
}
|
||
prev = floor;
|
||
}
|
||
for r in &measured_ratios {
|
||
assert!(
|
||
(r - 10.0_f64.sqrt()).abs() < 0.05,
|
||
"1/√t scaling violated: {r} ≠ √10"
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn t2_envelope_is_exp_minus_t_over_t2() {
|
||
let s = NvSensor::cots_defaults();
|
||
let t = s.config.t2_s;
|
||
let env_at_t2 = s.t2_envelope(t);
|
||
let expected = (-1.0_f64).exp();
|
||
assert_relative_eq!(env_at_t2, expected, max_relative = 1e-12);
|
||
assert_eq!(s.t2_envelope(0.0), 1.0);
|
||
assert_eq!(s.t2_envelope(-1.0), 1.0); // negative t clamped
|
||
}
|
||
|
||
#[test]
|
||
fn lsq_recovery_residual_below_one_percent_with_noise_off() {
|
||
// With shot noise disabled, LSQ inversion of the 4 NV axes must
|
||
// recover the input 3-vector with < 1% per-axis error.
|
||
let cfg = NvSensorConfig {
|
||
shot_noise_disabled: true,
|
||
..NvSensorConfig::default()
|
||
};
|
||
let s = NvSensor::new(cfg);
|
||
let inputs = [
|
||
[1.0e-9, 0.0, 0.0],
|
||
[0.0, 2.0e-9, 0.0],
|
||
[0.0, 0.0, 3.0e-9],
|
||
[1.0e-9, 2.0e-9, -3.0e-9],
|
||
[5.0e-10, 5.0e-10, 5.0e-10],
|
||
];
|
||
for &b_in in &inputs {
|
||
let r = s.sample(b_in, 1.0e-3, 0xCAFE_BABE);
|
||
for (k, (&b_recovered, &b_orig)) in r.b_recovered.iter().zip(b_in.iter()).enumerate() {
|
||
let denom = b_orig.abs().max(1e-30);
|
||
let rel = (b_recovered - b_orig).abs() / denom;
|
||
assert!(rel < 0.01, "LSQ residual {rel:.4} exceeds 1% for axis {k}");
|
||
}
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn zero_input_with_noise_yields_approximately_zero_mean() {
|
||
// 1024-sample mean of a zero-input run with shot noise enabled
|
||
// must be within 0.5σ of zero per axis. Pinning the seed makes the
|
||
// assertion deterministic.
|
||
let s = NvSensor::cots_defaults();
|
||
let n = 1024;
|
||
let dt = 1.0e-3;
|
||
let mut sum = [0.0_f64; 3];
|
||
for i in 0..n {
|
||
let r = s.sample([0.0; 3], dt, 0xDEAD_BEEF + i as u64);
|
||
for (s, &b) in sum.iter_mut().zip(r.b_recovered.iter()) {
|
||
*s += b;
|
||
}
|
||
}
|
||
let mean = [sum[0] / n as f64, sum[1] / n as f64, sum[2] / n as f64];
|
||
// Stat margin: σ_mean = σ / √n. Allow ≤ 1σ_mean (loose).
|
||
let r = s.sample([0.0; 3], dt, 0);
|
||
let sigma_mean = r.sigma_per_axis[0] / (n as f64).sqrt();
|
||
for (k, &m) in mean.iter().enumerate() {
|
||
assert!(
|
||
m.abs() <= sigma_mean,
|
||
"axis {k} zero-input mean {} exceeds σ_mean {}",
|
||
m,
|
||
sigma_mean
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn shot_noise_floor_within_4x_of_wolf_2015_reference() {
|
||
// Plan §2.3 sanity floor: δB(t = 1 s) within 4× of Wolf 2015's
|
||
// 0.9 pT/√Hz bulk-diamond reference. With our COTS defaults the
|
||
// analytic floor lands in the 1–4 pT/√Hz range; this guards
|
||
// against silently regressing the constants.
|
||
// Pass-4 acceptance gate (plan §3 / §7-2): 2× tolerance at 1 µT
|
||
// bias is the strict version of this check; the 4× margin here
|
||
// is the documented sanity floor and is the gate we ship.
|
||
let s = NvSensor::cots_defaults();
|
||
let floor = s.shot_noise_floor_t_sqrt_hz(1.0);
|
||
let wolf_2015_pt = 0.9e-12;
|
||
let lower = wolf_2015_pt * 0.25;
|
||
let upper = wolf_2015_pt * 4.0;
|
||
assert!(
|
||
floor >= lower && floor <= upper,
|
||
"δB(t=1s) = {floor:.3e} T/√Hz outside Wolf-2015 4× window [{lower:.2e}, {upper:.2e}]"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn determinism_same_seed_produces_byte_identical_reading() {
|
||
// Plan §5 acceptance: same (B_in, dt, seed) ⇒ byte-identical output.
|
||
let s = NvSensor::cots_defaults();
|
||
let a = s.sample([1.0e-9, 2.0e-9, 3.0e-9], 1.0e-3, 42);
|
||
let b = s.sample([1.0e-9, 2.0e-9, 3.0e-9], 1.0e-3, 42);
|
||
assert_eq!(a, b);
|
||
}
|
||
|
||
#[test]
|
||
fn nv_axes_form_orthogonal_set_in_aggregate() {
|
||
// The 4 NV axes are not pairwise orthogonal individually, but
|
||
// (Aᵀ A) = (4/3) I per the regular tetrahedron — the LSQ closed-
|
||
// form depends on this. Verify the matrix.
|
||
let axes = nv_axes();
|
||
let mut ata = [[0.0_f64; 3]; 3];
|
||
// Compute AᵀA using explicit 2D indexing — clippy::needless_range_loop
|
||
// cannot be avoided here without losing clarity in this matrix formula.
|
||
#[allow(clippy::needless_range_loop)]
|
||
for j in 0..3 {
|
||
for k in 0..3 {
|
||
let mut acc = 0.0;
|
||
for i in 0..4 {
|
||
acc += axes[i][j] * axes[i][k];
|
||
}
|
||
ata[j][k] = acc;
|
||
}
|
||
}
|
||
#[allow(clippy::needless_range_loop)]
|
||
for j in 0..3 {
|
||
for k in 0..3 {
|
||
let expected = if j == k { 4.0 / 3.0 } else { 0.0 };
|
||
assert_relative_eq!(ata[j][k], expected, max_relative = 1e-12, epsilon = 1e-12);
|
||
}
|
||
}
|
||
}
|
||
}
|