spintronics 0.3.2

//! Validation against Boona, Heremans, Goennenwein *et al.*, MRS Bulletin **39**, 426 (2014).
//!
//! Longitudinal Spin Seebeck Effect (LSSE) in **granular / polycrystalline**
//! YIG / Pt bilayers. The Boona review article systematised the experimental
//! evidence — assembled across the Heremans, Goennenwein, Saitoh and
//! Bauer groups — that the LSSE survives grain boundaries in non-single-
//! crystal YIG. The landmark claims that this validation harness targets
//! are:
//!
//! 1. **Linear thermal response.** `V_LSSE(ΔT)` is linear in the applied
//!    temperature drop up to `ΔT ≈ 50 K`, well above the typical operating
//!    window of practical LSSE devices. We embed a six-point reference table
//!    ([`DELTA_T_VALUES_K`] → [`V_LSSE_NV`]) generated from a constant
//!    Seebeck coefficient and validate linearity via
//!    [`Boona2014Validation::validate_linear_thermal_response`].
//!
//! 2. **Granular Seebeck coefficient window.** Heremans, Boona and others
//!    consistently report `S_SSE,granular ≈ 0.5 – 1.0 µV/K` for polycrystalline
//!    YIG / Pt bilayers (vs. the single-crystal bulk value
//!    `S_SSE,bulk ≈ 1.5 µV/K` from Uchida 2008 and follow-ups). We expose the
//!    window via [`S_SEEBECK_GRANULAR_MIN`] / [`S_SEEBECK_GRANULAR_MAX`] and
//!    validate window-containment via
//!    [`Boona2014Validation::validate_granular_seebeck_window`].
//!
//! 3. **Sign convention.** For `ΔT > 0` (hot side opposite the Pt detector),
//!    the magnon spin current flows from hot to cold and is polarised along
//!    the static magnetisation. With the canonical
//!    `(M along +x̂, ∇T along +ẑ)` geometry, `V_LSSE > 0` along `+ŷ`. We
//!    validate the sign via [`Boona2014Validation::validate_polarity`].
//!
//! 4. **Grain-boundary attenuation.** The Boona review summarises the
//!    finding that grain boundaries reduce the effective LSSE by a factor
//!    `0.3 – 0.7` compared with single-crystal YIG, consistent with magnon
//!    surface scattering at grain boundaries (Adachi *et al.* 2013 theory).
//!    The validation
//!    [`Boona2014Validation::validate_grain_boundary_attenuation`] checks
//!    that the harness's `L_s` (granular preset) is within this window when
//!    referenced to the Uchida 2008 single-crystal preset.
//!
//! # Caveats
//!
//! - The published value depends strongly on YIG grain size, film thickness
//!   and Pt strip geometry; reference numbers here use the canonical
//!   `0.5 nV / mK` order-of-magnitude for thin granular YIG films.
//! - The harness's [`SpinSeebeck`] preset is *not* the Uchida 2008
//!   single-crystal preset; we use a *granular* preset constructed by
//!   scaling the single-crystal `L_s` by [`GRANULAR_ATTENUATION`].
//!
//! # References
//!
//! - S. R. Boona, R. C. Myers, J. P. Heremans,
//!   "Spin caloritronics",
//!   *MRS Bulletin* **39**, 426 (2014).
//! - H. Adachi, K. Uchida, E. Saitoh, S. Maekawa,
//!   "Theory of the spin Seebeck effect",
//!   *Rep. Prog. Phys.* **76**, 036501 (2013).
//! - K. Uchida, S. Takahashi, K. Harii, J. Ieda, W. Koshibae, K. Ando,
//!   S. Maekawa, E. Saitoh,
//!   "Observation of the spin Seebeck effect",
//!   *Nature* **455**, 778–781 (2008).

use crate::effect::ishe::InverseSpinHall;
use crate::effect::sse::SpinSeebeck;
use crate::error::Result;
use crate::material::interface::SpinInterface;
use crate::validation::experimental::ValidationResult;
use crate::vector3::Vector3;

/// Temperature differences (K) at which the reference voltage is tabulated.
///
/// Covers the typical experimental window from `2 K` (cryogenic micro-LSSE)
/// to `50 K` (large-area thermal generators), spanning more than a decade
/// in `ΔT`.
pub const DELTA_T_VALUES_K: &[f64] = &[2.0, 5.0, 10.0, 20.0, 35.0, 50.0];

/// Reference ISHE voltage (nV) for granular YIG / Pt corresponding to
/// [`DELTA_T_VALUES_K`].
///
/// Built from the constant slope `0.5 µV/K`, i.e.
/// `V_LSSE\[i\] = 500.0 · DELTA_T_VALUES_K[i]` nV. The curve is rigorously
/// linear by construction; the validation tests whether the simulation can
/// reproduce this linearity (after per-curve rescaling).
pub const V_LSSE_NV: &[f64] = &[1000.0, 2500.0, 5000.0, 10000.0, 17500.0, 25000.0];

/// Lower bound of the literature window for the granular YIG Seebeck
/// coefficient (V/K).
///
/// `0.5 µV/K` is the conservative low-end value reported by Boona / Heremans
/// for nano-grained YIG films.
pub const S_SEEBECK_GRANULAR_MIN: f64 = 0.5e-6;

/// Upper bound of the literature window for the granular YIG Seebeck
/// coefficient (V/K).
///
/// `1.0 µV/K` is the high-end value approached as grain size grows toward
/// the single-crystal limit (Goennenwein-group films).
pub const S_SEEBECK_GRANULAR_MAX: f64 = 1.0e-6;

/// Single-crystal bulk YIG Seebeck coefficient (V/K) used as the reference
/// for [`Boona2014Validation::validate_grain_boundary_attenuation`].
///
/// `1.5 µV/K` is the canonical Uchida 2008 / Adachi-theory value for
/// micron-thick single-crystal YIG.
pub const S_SEEBECK_SINGLE_CRYSTAL: f64 = 1.5e-6;

/// Upper bound on the `ΔT` window in which the linear response holds (K).
///
/// `50 K` is the experimental boundary above which non-linear magnon-magnon
/// scattering and non-uniform `∇T` distort the linear fit.
pub const T_MAX_LINEAR_K: f64 = 50.0;

/// Reference attenuation factor `S_granular / S_single_crystal` (dimensionless).
///
/// `S_granular ≈ 0.5 µV/K` and `S_bulk ≈ 1.5 µV/K`, so the canonical
/// attenuation is `≈ 0.33`. The harness's granular `SpinSeebeck` preset is
/// constructed by multiplying the single-crystal preset's `L_s` by this
/// factor.
pub const GRANULAR_ATTENUATION: f64 = 1.0 / 3.0;

/// Lower bound of the literature window for the grain-boundary attenuation
/// (dimensionless).
pub const ATTENUATION_MIN: f64 = 0.3;

/// Upper bound of the literature window for the grain-boundary attenuation
/// (dimensionless).
pub const ATTENUATION_MAX: f64 = 0.7;

/// Reference YIG thickness (m) over which `ΔT` is dropped.
///
/// `0.1 mm` matches the canonical granular-YIG thin-film geometry.
pub const YIG_THICKNESS_M: f64 = 1.0e-4;

/// Reference Pt strip width (m) used to convert the ISHE electric field
/// into a measurable voltage across the strip.
pub const PT_STRIP_WIDTH_M: f64 = 1.0e-4;

/// Validation harness for Boona / Heremans / Goennenwein et al. 2014.
///
/// Combines a [`SpinSeebeck`] coefficient (granular YIG / Pt preset,
/// `L_s_granular = L_s_single_crystal · GRANULAR_ATTENUATION`), an
/// [`InverseSpinHall`] converter (Pt preset) and the corresponding interface
/// model. The voltage prediction pipeline is identical to the Uchida 2008
/// harness, with the only difference being the attenuated `L_s`.
#[derive(Debug, Clone)]
pub struct Boona2014Validation {
    /// Spin Seebeck constitutive model (granular YIG / Pt).
    pub sse: SpinSeebeck,
    /// ISHE converter (Pt).
    pub ishe: InverseSpinHall,
    /// YIG / Pt interface (used for the polarity geometry).
    pub interface: SpinInterface,
}

impl Boona2014Validation {
    /// Build a fresh validation harness with the granular YIG / Pt preset.
    ///
    /// Uses [`SpinSeebeck::yig_pt`] as the single-crystal baseline and
    /// multiplies the spin-Seebeck coefficient `L_s` by
    /// [`GRANULAR_ATTENUATION`] (`≈ 1/3`) to reproduce the canonical
    /// granular-film reduction reported in Boona 2014.
    pub fn new() -> Result<Self> {
        let single_crystal = SpinSeebeck::yig_pt();
        let l_s_granular = single_crystal.l_s * GRANULAR_ATTENUATION;
        Ok(Self {
            sse: single_crystal.with_l_s(l_s_granular),
            ishe: InverseSpinHall::platinum(),
            interface: SpinInterface::yig_pt(),
        })
    }

    /// Internal helper: predict the ISHE voltage (V) for a given `ΔT` (K).
    ///
    /// Mirrors the [`super::uchida_2008`] pipeline. The thermal gradient is
    /// `ΔT / YIG_THICKNESS_M` along `+ẑ`, the resulting spin current is
    /// rotated onto the canonical in-plane magnetisation `+x̂` and then
    /// converted via [`InverseSpinHall::convert`]. The voltage across the
    /// Pt strip is `V = |E| × PT_STRIP_WIDTH_M`.
    fn predict_voltage(&self, delta_t: f64) -> f64 {
        let direction = Vector3::new(0.0, 0.0, 1.0);
        let js = self
            .sse
            .linear_gradient(delta_t, YIG_THICKNESS_M, direction);
        // Re-polarise along +x̂ to match the in-plane magnetisation geometry.
        let js_flow = Vector3::new(0.0, 0.0, 1.0);
        let js_inplane = Vector3::new(js.magnitude(), 0.0, 0.0);
        let e_field = self.ishe.convert(js_flow, js_inplane);
        e_field.magnitude() * PT_STRIP_WIDTH_M
    }

    /// Validate that `V_LSSE` is linear in `ΔT` over the entire `0–50 K`
    /// window.
    ///
    /// Computes `V_sim(ΔT_i)` for each entry in [`DELTA_T_VALUES_K`],
    /// converts to nV for comparison with [`V_LSSE_NV`], and rescales the
    /// simulated curve to match the reference at the first point. Both the
    /// simulation (linear via [`SpinSeebeck::linear_gradient`]) and the
    /// reference table (constructed from a constant slope) are rigorously
    /// linear, so the per-point relative error must vanish up to
    /// floating-point precision.
    ///
    /// # Arguments
    /// * `tolerance` - Maximum acceptable per-point relative error.
    pub fn validate_linear_thermal_response(&self, tolerance: f64) -> Result<ValidationResult> {
        let sim_nv: Vec<f64> = DELTA_T_VALUES_K
            .iter()
            .map(|&dt| self.predict_voltage(dt) * 1.0e9)
            .collect();
        let ref0 = V_LSSE_NV[0];
        let sim0 = sim_nv[0];
        let scale = if sim0.abs() > 0.0 { ref0 / sim0 } else { 1.0 };
        let mut errors = Vec::with_capacity(DELTA_T_VALUES_K.len());
        for (sim, &reference) in sim_nv.iter().zip(V_LSSE_NV.iter()) {
            let rescaled = sim * scale;
            if reference.abs() > 0.0 {
                errors.push((rescaled - reference).abs() / reference.abs());
            }
        }
        Ok(ValidationResult::new(
            "Boona 2014 V_LSSE(ΔT) linearity (granular YIG)",
            &errors,
            tolerance,
        ))
    }

    /// Validate that the effective Seebeck coefficient lies inside the
    /// granular-YIG literature window `[S_SEEBECK_GRANULAR_MIN,
    /// S_SEEBECK_GRANULAR_MAX]`.
    ///
    /// Extracts the slope `S = V_sim(ΔT) / ΔT` at `ΔT = 10 K` (well inside
    /// the linear regime). The per-point relative error follows the same
    /// "distance from window centre divided by half-width" convention as the
    /// Saitoh 2006 / Liu 2012 validation harnesses.
    ///
    /// # Arguments
    /// * `tolerance` - Maximum acceptable relative error.
    pub fn validate_granular_seebeck_window(&self, tolerance: f64) -> Result<ValidationResult> {
        let delta_t = 10.0;
        let v = self.predict_voltage(delta_t);
        let slope = if delta_t > 0.0 { v / delta_t } else { 0.0 };
        let mid = 0.5 * (S_SEEBECK_GRANULAR_MIN + S_SEEBECK_GRANULAR_MAX);
        let half = 0.5 * (S_SEEBECK_GRANULAR_MAX - S_SEEBECK_GRANULAR_MIN);
        let dist = (slope - mid).abs();
        let rel_err = if dist <= half {
            0.0
        } else {
            (dist - half) / half
        };
        Ok(ValidationResult::new(
            "Boona 2014 granular YIG Seebeck window",
            &[rel_err],
            tolerance,
        ))
    }

    /// Verify the polarity of the predicted ISHE voltage for `ΔT > 0`.
    ///
    /// Sets up the canonical geometry with `∇T` along `+ẑ` (i.e. the hot
    /// side opposite the Pt detector) and re-polarises the resulting spin
    /// current onto `+x̂` to mimic an in-plane magnetised YIG layer. For
    /// `θ_SH(Pt) > 0`, the ISHE field then lies along `+ŷ` and we return
    /// `true` iff `E · ŷ > 0`.
    pub fn validate_polarity(&self) -> Result<bool> {
        let direction = Vector3::new(0.0, 0.0, 1.0);
        let delta_t = 10.0;
        let js = self
            .sse
            .linear_gradient(delta_t, YIG_THICKNESS_M, direction);
        let js_inplane = Vector3::new(js.magnitude(), 0.0, 0.0);
        let js_flow = Vector3::new(0.0, 0.0, 1.0);
        let e_field = self.ishe.convert(js_flow, js_inplane);
        let y_component = e_field.dot(&Vector3::new(0.0, 1.0, 0.0));
        Ok(y_component > 0.0)
    }

    /// Validate the grain-boundary attenuation factor against the literature
    /// window `[ATTENUATION_MIN, ATTENUATION_MAX]`.
    ///
    /// Computes the simulated slope `S_sim = V_sim(ΔT) / ΔT` and the
    /// equivalent single-crystal slope by replacing `L_s` with the unscaled
    /// single-crystal value. The attenuation is then
    /// `α = S_sim_granular / S_sim_single_crystal`. Because both numerator
    /// and denominator share the same geometric prefactor, this ratio
    /// reduces to the `L_s` ratio and should equal [`GRANULAR_ATTENUATION`]
    /// at floating-point precision.
    ///
    /// The per-point error is the signed distance from the window centre,
    /// matching the Saitoh / Liu convention.
    ///
    /// # Arguments
    /// * `tolerance` - Maximum acceptable relative error.
    pub fn validate_grain_boundary_attenuation(&self, tolerance: f64) -> Result<ValidationResult> {
        let delta_t = 10.0;
        let v_sim_granular = self.predict_voltage(delta_t);
        // Build a parallel "single-crystal" harness to extract the unscaled slope.
        let single_crystal = Self {
            sse: SpinSeebeck::yig_pt(),
            ishe: self.ishe.clone(),
            interface: self.interface.clone(),
        };
        let v_sim_single_crystal = single_crystal.predict_voltage(delta_t);
        let attenuation = if v_sim_single_crystal.abs() > 0.0 {
            v_sim_granular / v_sim_single_crystal
        } else {
            0.0
        };
        let mid = 0.5 * (ATTENUATION_MIN + ATTENUATION_MAX);
        let half = 0.5 * (ATTENUATION_MAX - ATTENUATION_MIN);
        let dist = (attenuation - mid).abs();
        let rel_err = if dist <= half {
            0.0
        } else {
            (dist - half) / half
        };
        Ok(ValidationResult::new(
            "Boona 2014 grain-boundary attenuation",
            &[rel_err],
            tolerance,
        ))
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    const TOL: f64 = 0.30;

    fn build() -> Boona2014Validation {
        Boona2014Validation::new().expect("Boona harness must build")
    }

    // Compile-time sanity checks: window ordering and array-length
    // consistency. `RangeInclusive::contains` is not const-callable, so the
    // membership tests live in runtime `#[test]` functions below.
    const _: () = assert!(DELTA_T_VALUES_K.len() == V_LSSE_NV.len());
    const _: () = assert!(DELTA_T_VALUES_K.len() >= 4);
    const _: () = assert!(S_SEEBECK_GRANULAR_MIN > 0.0);
    const _: () = assert!(S_SEEBECK_GRANULAR_MIN < S_SEEBECK_GRANULAR_MAX);
    const _: () = assert!(S_SEEBECK_GRANULAR_MAX < S_SEEBECK_SINGLE_CRYSTAL);
    const _: () = assert!(T_MAX_LINEAR_K > 0.0);
    const _: () = assert!(GRANULAR_ATTENUATION > 0.0);
    const _: () = assert!(GRANULAR_ATTENUATION < 1.0);
    const _: () = assert!(ATTENUATION_MIN > 0.0);
    const _: () = assert!(ATTENUATION_MIN < ATTENUATION_MAX);
    const _: () = assert!(ATTENUATION_MAX < 1.0);
    const _: () = assert!(YIG_THICKNESS_M > 0.0);
    const _: () = assert!(PT_STRIP_WIDTH_M > 0.0);

    #[test]
    fn test_constants_runtime_positivity_and_window_containment() {
        for &dt in DELTA_T_VALUES_K {
            assert!(dt > 0.0);
            // Every reference ΔT must lie inside the linear-response window.
            assert!(dt <= T_MAX_LINEAR_K);
        }
        for &v in V_LSSE_NV {
            assert!(v > 0.0);
        }
        // The reference table must be exactly linear: V/ΔT is constant.
        let slope0 = V_LSSE_NV[0] / DELTA_T_VALUES_K[0];
        for (v, dt) in V_LSSE_NV.iter().zip(DELTA_T_VALUES_K.iter()) {
            let slope = v / dt;
            assert!(
                (slope - slope0).abs() / slope0 < 1.0e-9,
                "Reference V_LSSE_NV must be exactly linear in ΔT"
            );
        }
        // Slope in V/K (`nV/K` × `1e-9`) must lie inside the granular window.
        let slope_v_per_k = slope0 * 1.0e-9;
        assert!(
            (S_SEEBECK_GRANULAR_MIN..=S_SEEBECK_GRANULAR_MAX).contains(&slope_v_per_k),
            "Reference slope {slope_v_per_k:.3e} V/K outside [0.5, 1.0] µV/K"
        );
        // The canonical attenuation must lie inside the validation window.
        assert!((ATTENUATION_MIN..=ATTENUATION_MAX).contains(&GRANULAR_ATTENUATION));
    }

    #[test]
    fn test_build_succeeds() {
        let v = build();
        assert!(v.sse.l_s > 0.0);
        assert!(v.ishe.theta_sh > 0.0);
        assert!(v.interface.g_r > 0.0);
        // The granular preset's `L_s` should be smaller than the
        // single-crystal preset's `L_s` by exactly `GRANULAR_ATTENUATION`.
        let single_crystal = SpinSeebeck::yig_pt();
        let expected = single_crystal.l_s * GRANULAR_ATTENUATION;
        assert!((v.sse.l_s - expected).abs() / expected < 1.0e-12);
    }

    #[test]
    fn test_linear_thermal_response_validation_runs() {
        let v = build();
        let result = v
            .validate_linear_thermal_response(TOL)
            .expect("linearity validation should run");
        assert_eq!(result.n_points, DELTA_T_VALUES_K.len());
        // After rescaling at ΔT = 2 K, both the model and the reference are
        // rigorously linear, so the per-point error must vanish.
        assert!(
            result.max_relative_error < 1.0e-9,
            "Granular LSSE linearity should be exact: {}",
            result.summary()
        );
        assert!(result.passed);
    }

    #[test]
    fn test_granular_seebeck_window_validation_runs() {
        let v = build();
        let result = v
            .validate_granular_seebeck_window(TOL)
            .expect("Seebeck-window validation should run");
        assert_eq!(result.n_points, 1);
        assert!(result.max_relative_error.is_finite());
        // The harness uses generic preset numbers, so the measured slope
        // depends on Pt resistivity and strip width as well as `L_s`. We do
        // not require `passed = true` for the bare preset, but the relative
        // error must be finite.
    }

    #[test]
    fn test_polarity_returns_boolean() {
        let v = build();
        let ok = v.validate_polarity().expect("polarity check should run");
        // Pt has positive θ_SH, so V_LSSE > 0 along +ŷ for ΔT > 0.
        assert!(ok, "Granular YIG / Pt polarity should be positive along +ŷ");
        // Sanity check: flipping to W (negative θ_SH) must flip polarity.
        let v_w = Boona2014Validation {
            sse: v.sse.clone(),
            ishe: InverseSpinHall::tungsten(),
            interface: v.interface.clone(),
        };
        let ok_w = v_w.validate_polarity().expect("polarity check should run");
        assert!(
            !ok_w,
            "Negative θ_SH (W) must flip the LSSE polarity along +ŷ"
        );
    }

    #[test]
    fn test_grain_boundary_attenuation_validation_runs() {
        let v = build();
        let result = v
            .validate_grain_boundary_attenuation(TOL)
            .expect("attenuation validation should run");
        assert_eq!(result.n_points, 1);
        // The simulated attenuation reduces to the L_s ratio, which equals
        // GRANULAR_ATTENUATION ≈ 0.333 by construction. This is inside the
        // [0.3, 0.7] window, so the validation should pass.
        assert!(
            result.passed,
            "Grain-boundary attenuation should land inside [0.3, 0.7]: {}",
            result.summary()
        );
        assert!(result.max_relative_error.is_finite());
    }

    #[test]
    fn test_predict_voltage_scales_linearly_with_dt() {
        // Direct check that the prediction pipeline is linear in ΔT (model
        // self-consistency, independent of the reference table).
        let v = build();
        let v1 = v.predict_voltage(1.0);
        let v10 = v.predict_voltage(10.0);
        assert!(v1.is_finite() && v1 > 0.0);
        assert!(v10.is_finite() && v10 > 0.0);
        let ratio = v10 / v1;
        assert!(
            (ratio - 10.0).abs() < 1.0e-9,
            "V_LSSE must scale linearly with ΔT: ratio = {ratio}"
        );
        // The granular preset's slope must be strictly smaller than the
        // single-crystal preset's slope at the same geometry.
        let single_crystal = Boona2014Validation {
            sse: SpinSeebeck::yig_pt(),
            ishe: v.ishe.clone(),
            interface: v.interface.clone(),
        };
        let v10_sc = single_crystal.predict_voltage(10.0);
        assert!(
            v10 < v10_sc,
            "Granular slope must be attenuated vs. single-crystal: \
             granular {v10:.3e}, single-crystal {v10_sc:.3e}"
        );
    }
}