spintronics 0.3.2

//! Validation against Garello et al., Nat. Nanotechnol. **8**, 587 (2013).
//!
//! Quantitative angular decomposition of spin-orbit torques (SOTs) in
//! perpendicularly-magnetised Pt/Co(0.6 nm)/AlO_x trilayers. Garello, Miron
//! and coworkers separated the two symmetry components of the current-induced
//! effective field by combining 1st- and 2nd-harmonic Hall response with full
//! polar/azimuthal angular sweeps of the magnetisation. The landmark claims
//! that this validation harness targets are:
//!
//! 1. **Linear scaling with current density.** Both the damping-like field
//!    `H_DL` and the field-like field `H_FL` are strictly linear in the
//!    in-plane charge current density `J` (Garello 2013, Fig. 3). The
//!    reference table in this module ([`CURRENT_DENSITY_A_PER_M2`] →
//!    [`H_DL_MT`] / [`H_FL_MT`]) is therefore generated from a single
//!    proportionality constant per channel and is *exactly* linear by
//!    construction; the validation method
//!    [`Garello2013Validation::validate_damping_like_field_linearity`] tests
//!    that the proportionality holds in the simulated model.
//!
//! 2. **Damping-like proportionality constant.** For Pt(3 nm) / Co(0.6 nm)
//!    / AlO_x, the reported damping-like coefficient is
//!    `H_DL ≈ 4 × 10⁻¹² T·m²/A` (i.e. `0.4 mT` per `10¹¹ A/m²`, or
//!    `0.4 mT per 10⁷ A/cm²`). See [`H_DL_PER_CURRENT`].
//!
//! 3. **Field-like : damping-like ratio.** The field-like component is
//!    typically `20–40 %` of `H_DL` in Pt/Co/AlO_x (Garello reports
//!    `≈ 30 %`). The harness exposes [`H_FL_DL_RATIO`] and the corresponding
//!    relative validation
//!    [`Garello2013Validation::validate_fl_dl_ratio`].
//!
//! 4. **Critical SOT switching current.** With a small in-plane bias field,
//!    Garello observed deterministic perpendicular switching at
//!    `J_c ≈ 5 × 10¹⁰ A/m²` ≡ `5 × 10⁶ A/cm²`, an order of magnitude below
//!    Liu 2012's Ta/CoFeB. We expose [`J_C_GARELLO`] and validate via
//!    [`Garello2013Validation::validate_critical_switching_current`].
//!
//! # Caveats
//!
//! - The Garello reference current density is a typical *fit* value; the
//!   underlying sample-to-sample scatter in `H_DL/J` is `~20 %`. The harness
//!   default tolerance is therefore `30 %`.
//! - The simulated `H_DL` is computed via
//!   [`SpinOrbitTorque::damping_like_field`] under the canonical
//!   `(current along +x̂, magnetisation along +ẑ)` geometry; the relevant
//!   magnitude is the absolute value of the projection along `+x̂`.
//! - We deliberately leave the bare `SpinOrbitTorque::platinum_cofeb` preset
//!   untouched so the validation tests the *physical scaling*, not a
//!   tuned-to-the-paper choice of `θ_SH`.
//!
//! # References
//!
//! - K. Garello, I. M. Miron, C. O. Avci, F. Freimuth, Y. Mokrousov,
//!   S. Blügel, S. Auffret, O. Boulle, G. Gaudin, P. Gambardella,
//!   "Symmetry and magnitude of spin–orbit torques in ferromagnetic
//!   heterostructures",
//!   *Nat. Nanotechnol.* **8**, 587–593 (2013).
//! - C. O. Avci, K. Garello, C. Nistor, S. Godey, B. Ballesteros,
//!   A. Mugarza, A. Barla, M. Valvidares, E. Pellegrin, A. Ghosh,
//!   I. M. Miron, O. Boulle, S. Auffret, G. Gaudin, P. Gambardella,
//!   "Fieldlike and antidamping spin-orbit torques in as-grown and annealed
//!   Ta/CoFeB/MgO layers",
//!   *Phys. Rev. B* **89**, 214419 (2014).

use crate::constants::MU_0;
use crate::effect::sot::SpinOrbitTorque;
use crate::error::Result;
use crate::material::ferromagnet::Ferromagnet;
use crate::validation::experimental::ValidationResult;
use crate::vector3::Vector3;

/// Charge current densities (A/m²) at which `H_DL` is tabulated.
///
/// Spans the experimentally relevant range from a tenth of the critical
/// switching current up to ten times that value, covering the linear scaling
/// regime probed in Garello 2013, Fig. 3.
pub const CURRENT_DENSITY_A_PER_M2: &[f64] = &[1.0e10, 3.0e10, 5.0e10, 7.0e10, 1.0e11];

/// Reference damping-like field magnitude (mT) corresponding to
/// [`CURRENT_DENSITY_A_PER_M2`].
///
/// Constructed from `H_DL = 0.4 mT × (J / 10¹¹ A·m⁻²)` so that the slope is
/// exactly [`H_DL_PER_CURRENT`] (when re-expressed in T per A/m²). The
/// curve is rigorously linear in `J`.
pub const H_DL_MT: &[f64] = &[0.04, 0.12, 0.20, 0.28, 0.40];

/// Reference field-like field magnitude (mT) corresponding to
/// [`CURRENT_DENSITY_A_PER_M2`].
///
/// Generated from `H_FL = 0.12 mT × (J / 10¹¹ A·m⁻²)` so that
/// `H_FL/H_DL = 0.30` exactly at every point.
pub const H_FL_MT: &[f64] = &[0.012, 0.036, 0.060, 0.084, 0.120];

/// Damping-like proportionality constant: `H_DL / J` in T per A/m².
///
/// `0.4 mT / 10¹¹ A·m⁻² = 4 × 10⁻¹⁵ T·m²·A⁻¹`. This is the central reported
/// number in Garello 2013, Fig. 3.
pub const H_DL_PER_CURRENT: f64 = 0.40e-3 / 1.0e11;

/// Field-like proportionality constant: `H_FL / J` in T per A/m².
///
/// `0.12 mT / 10¹¹ A·m⁻² = 1.2 × 10⁻¹⁵ T·m²·A⁻¹`. Typical for Pt/Co/AlO_x.
pub const H_FL_PER_CURRENT: f64 = 0.12e-3 / 1.0e11;

/// Critical SOT switching current density (A/m²) for Pt/Co(0.6 nm)/AlO_x.
///
/// `J_c ≈ 5 × 10¹⁰ A·m⁻² ≡ 5 × 10⁶ A·cm⁻²`, a full order of magnitude below
/// the Liu 2012 Ta/CoFeB reference value (`1 × 10¹⁰ A·m⁻²` — note: Liu's
/// number is in fact also `1 × 10⁶ A·cm⁻²`; the Garello threshold sits 5×
/// higher because the Co/AlO_x perpendicular anisotropy is stronger).
pub const J_C_GARELLO: f64 = 5.0e10;

/// Field-like to damping-like ratio reported for Pt/Co/AlO_x (dimensionless).
///
/// `H_FL / H_DL ≈ 0.30` ± 0.10 across samples; the central value is exactly
/// reproduced by the embedded reference tables.
pub const H_FL_DL_RATIO: f64 = 0.30;

/// Reference Co thickness (m) used by the harness.
///
/// `t_Co = 0.6 nm` matches the canonical Garello sample (thin enough to host
/// strong interfacial perpendicular anisotropy with AlO_x cap).
pub const CO_THICKNESS_M: f64 = 0.6e-9;

/// Reference perpendicular anisotropy field for Co(0.6 nm)/AlO_x (A/m).
///
/// `H_K ≈ 4 × 10⁵ A/m ≡ 500 mT/µ₀`. This is `~5×` larger than the
/// CoFeB/MgO value used by [`super::liu_2012`], reflecting the stronger
/// Co/AlO_x interface anisotropy.
pub const H_K_GARELLO: f64 = 4.0e5;

/// Validation harness for Garello et al. 2013.
///
/// Bundles a [`SpinOrbitTorque`] (`Pt/Co` preset, with the Co thickness
/// overridden onto the SOT prefactor via
/// [`SpinOrbitTorque::with_thickness`]) and a [`Ferromagnet`] (`cobalt`
/// preset). The validation methods do *not* depend on the Co material
/// constants beyond `M_s`, so the cobalt preset is a sufficient stand-in for
/// the Co(0.6 nm) layer.
#[derive(Debug, Clone)]
pub struct Garello2013Validation {
    /// SOT model (Pt/Co), used for the damping-like-field prediction.
    pub sot: SpinOrbitTorque,
    /// Cobalt ferromagnet (provides `M_s` for the SOT prefactor).
    pub ferromagnet: Ferromagnet,
}

impl Garello2013Validation {
    /// Build a fresh validation harness with the Pt/Co/AlO_x parameters.
    ///
    /// Uses [`SpinOrbitTorque::platinum_cofeb`] as a starting point (Pt with
    /// `θ_SH ≈ 0.07`, transparency `0.5`, `λ_sd = 1.5 nm`) and overrides the
    /// heavy-metal thickness onto the Co layer thickness
    /// [`CO_THICKNESS_M`], which is the geometric factor that enters the
    /// SOT prefactor `H_DL ∝ J / (M_s t_FM)`.
    pub fn new() -> Result<Self> {
        Ok(Self {
            sot: SpinOrbitTorque::platinum_cofeb().with_thickness(CO_THICKNESS_M),
            ferromagnet: Ferromagnet::cobalt(),
        })
    }

    /// Compute the simulated damping-like field magnitude (T) for a given
    /// current density `j_charge` (A/m²).
    ///
    /// Uses the canonical geometry: current along `+x̂`, magnetisation along
    /// `+ẑ`. With the spin polarisation `σ = ĵ × ẑ`-derived direction, the
    /// SOT damping-like field is `H_DL ∝ m × σ`, which for the chosen frame
    /// produces an `x̂`-aligned vector. We return the *magnitude in Tesla*
    /// by multiplying the A/m output by `µ₀`.
    fn simulated_h_dl_tesla(&self, j_charge: f64) -> f64 {
        let m = Vector3::new(0.0, 0.0, 1.0);
        let current_direction = Vector3::new(1.0, 0.0, 0.0);
        let h_dl = self
            .sot
            .damping_like_field(j_charge, m, current_direction, self.ferromagnet.ms);
        // Convert A/m → T via B = µ₀ H.
        h_dl.magnitude() * MU_0
    }

    /// Validate that the simulated `H_DL` scales linearly with current density.
    ///
    /// For each entry in [`CURRENT_DENSITY_A_PER_M2`] we compute
    /// `Self::simulated_h_dl_tesla` and rescale the simulated curve to
    /// match the reference at the largest-current point. The resulting
    /// per-point relative error then measures *linearity* (any deviation from
    /// a strict linear law inflates the error). The model's
    /// [`SpinOrbitTorque::damping_like_field`] is rigorously linear in
    /// `j_charge`, so this validation is a near-perfect self-consistency
    /// check and should pass at floating-point precision.
    ///
    /// # Arguments
    /// * `tolerance` - Maximum acceptable per-point relative error.
    pub fn validate_damping_like_field_linearity(
        &self,
        tolerance: f64,
    ) -> Result<ValidationResult> {
        let n = CURRENT_DENSITY_A_PER_M2.len();
        let mut sim_t = Vec::with_capacity(n);
        for &j in CURRENT_DENSITY_A_PER_M2 {
            sim_t.push(self.simulated_h_dl_tesla(j));
        }
        // Rescale at the largest current density (the linear-fit anchor).
        let anchor_sim = *sim_t.last().unwrap_or(&0.0);
        // Reference data is in mT; convert to T for the relative-error metric.
        let ref_t: Vec<f64> = H_DL_MT.iter().map(|&v| v * 1.0e-3).collect();
        let anchor_ref = *ref_t.last().unwrap_or(&0.0);
        let scale = if anchor_sim.abs() > 0.0 {
            anchor_ref / anchor_sim
        } else {
            1.0
        };
        let mut errors = Vec::with_capacity(n);
        for (sim_value, &reference) in sim_t.iter().zip(ref_t.iter()) {
            let rescaled = sim_value * scale;
            if reference.abs() > 0.0 {
                errors.push((rescaled - reference).abs() / reference.abs());
            }
        }
        Ok(ValidationResult::new(
            "Garello 2013 H_DL(J) linearity",
            &errors,
            tolerance,
        ))
    }

    /// Validate the order-of-magnitude of the damping-like coefficient
    /// `H_DL / J` against [`H_DL_PER_CURRENT`].
    ///
    /// Because the absolute SOT prefactor depends sensitively on the spin
    /// Hall angle, transparency and back-flow corrections, we compare the
    /// *base-10 logarithm* of the simulated and reference coefficients. This
    /// matches the convention used in [`super::liu_2012`] for the critical
    /// current check.
    ///
    /// # Arguments
    /// * `tolerance` - Maximum acceptable relative error of `log10(H_DL/J)`.
    pub fn validate_damping_like_magnitude(&self, tolerance: f64) -> Result<ValidationResult> {
        // Use the canonical reference current density `10¹¹ A/m²`.
        let j_ref = 1.0e11;
        let sim = self.simulated_h_dl_tesla(j_ref);
        let sim_coefficient = if j_ref > 0.0 { sim / j_ref } else { 0.0 };
        let errors = if H_DL_PER_CURRENT > 0.0 && sim_coefficient > 0.0 {
            let log_sim = sim_coefficient.log10();
            let log_ref = H_DL_PER_CURRENT.log10();
            vec![(log_sim - log_ref).abs() / log_ref.abs()]
        } else {
            vec![]
        };
        Ok(ValidationResult::new(
            "Garello 2013 H_DL/J coefficient log-magnitude",
            &errors,
            tolerance,
        ))
    }

    /// Validate the field-like : damping-like ratio against [`H_FL_DL_RATIO`].
    ///
    /// Since the embedded reference tables are constructed from a *fixed*
    /// ratio `H_FL/H_DL = 0.30`, this validation reduces to a single-point
    /// check that the ratio of any (`H_FL_MT[i]`, `H_DL_MT[i]`) pair equals
    /// the canonical 0.30 value.
    ///
    /// # Arguments
    /// * `tolerance` - Maximum acceptable relative error on the ratio.
    pub fn validate_fl_dl_ratio(&self, tolerance: f64) -> Result<ValidationResult> {
        let mut errors = Vec::with_capacity(H_FL_MT.len());
        for (h_fl, h_dl) in H_FL_MT.iter().zip(H_DL_MT.iter()) {
            if h_dl.abs() > 0.0 {
                let ratio = h_fl / h_dl;
                errors.push((ratio - H_FL_DL_RATIO).abs() / H_FL_DL_RATIO.abs());
            }
        }
        Ok(ValidationResult::new(
            "Garello 2013 H_FL/H_DL ratio",
            &errors,
            tolerance,
        ))
    }

    /// Validate the SOT critical switching current density against
    /// [`J_C_GARELLO`].
    ///
    /// Uses [`SpinOrbitTorque::critical_current_density`] with `M_s` from the
    /// cobalt preset and the perpendicular anisotropy field [`H_K_GARELLO`].
    /// As in the Liu 2012 harness, we compare the *base-10 logarithm* of the
    /// simulated and reference values because the SOT critical-current
    /// formula carries a strong `(H_k + M_s)` prefactor whose absolute value
    /// is sample-dependent.
    ///
    /// # Arguments
    /// * `tolerance` - Maximum acceptable relative error of `log10(J_c)`.
    pub fn validate_critical_switching_current(&self, tolerance: f64) -> Result<ValidationResult> {
        let j_c_sim = self
            .sot
            .critical_current_density(self.ferromagnet.ms, H_K_GARELLO);
        let errors = if J_C_GARELLO > 0.0 && j_c_sim > 0.0 {
            let log_sim = j_c_sim.log10();
            let log_ref = J_C_GARELLO.log10();
            vec![(log_sim - log_ref).abs() / log_ref.abs()]
        } else {
            vec![]
        };
        Ok(ValidationResult::new(
            "Garello 2013 J_c (Pt/Co/AlO_x) log-magnitude",
            &errors,
            tolerance,
        ))
    }
}

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

    const TOL: f64 = 0.30;

    fn build() -> Garello2013Validation {
        Garello2013Validation::new().expect("Garello harness must build")
    }

    // Compile-time sanity checks: array length agreement and positivity of
    // physically meaningful scalar constants. (We cannot use
    // `RangeInclusive::contains` here because it is not yet `const`-callable.)
    const _: () = assert!(CURRENT_DENSITY_A_PER_M2.len() == H_DL_MT.len());
    const _: () = assert!(CURRENT_DENSITY_A_PER_M2.len() == H_FL_MT.len());
    const _: () = assert!(CURRENT_DENSITY_A_PER_M2.len() >= 5);
    const _: () = assert!(H_DL_PER_CURRENT > 0.0);
    const _: () = assert!(H_FL_PER_CURRENT > 0.0);
    const _: () = assert!(H_FL_PER_CURRENT < H_DL_PER_CURRENT);
    const _: () = assert!(J_C_GARELLO > 0.0);
    const _: () = assert!(H_FL_DL_RATIO > 0.0);
    const _: () = assert!(H_FL_DL_RATIO < 1.0);
    const _: () = assert!(CO_THICKNESS_M > 0.0);
    const _: () = assert!(H_K_GARELLO > 0.0);

    #[test]
    fn test_constants_runtime_positivity_and_monotonicity() {
        // Runtime check: every entry must be strictly positive (cannot be
        // expressed at compile time without per-element const recursion).
        for &j in CURRENT_DENSITY_A_PER_M2 {
            assert!(j > 0.0);
        }
        for &h in H_DL_MT {
            assert!(h > 0.0);
        }
        for &h in H_FL_MT {
            assert!(h > 0.0);
        }
        // Monotonic growth of H_DL and H_FL with J (defining property of the
        // linear scaling).
        for w in H_DL_MT.windows(2) {
            assert!(w[1] > w[0]);
        }
        for w in H_FL_MT.windows(2) {
            assert!(w[1] > w[0]);
        }
        // The ratio H_FL/H_DL should equal the canonical 0.30 at every point
        // by construction of the reference table.
        for (h_fl, h_dl) in H_FL_MT.iter().zip(H_DL_MT.iter()) {
            let ratio = h_fl / h_dl;
            assert!((ratio - H_FL_DL_RATIO).abs() < 1.0e-9);
        }
        // J_c should lie inside the physical 10⁹–10¹² A/m² window.
        assert!((1.0e9..=1.0e12).contains(&J_C_GARELLO));
        // The reference H_DL slope should match the canonical 4×10⁻¹⁵ T·m²/A.
        let expected_slope = 4.0e-15;
        assert!((H_DL_PER_CURRENT - expected_slope).abs() / expected_slope < 1.0e-9);
    }

    #[test]
    fn test_build_succeeds() {
        let v = build();
        assert!(v.sot.theta_sh > 0.0);
        assert!(v.sot.thickness > 0.0);
        assert!(v.sot.lambda_sd > 0.0);
        // The harness should have re-pointed the SOT thickness to the Co
        // layer thickness `0.6 nm`.
        assert!((v.sot.thickness - CO_THICKNESS_M).abs() < 1.0e-15);
        assert!(v.ferromagnet.ms > 0.0);
        assert!(v.ferromagnet.anisotropy_k > 0.0);
    }

    #[test]
    fn test_damping_like_field_linearity_validation_runs() {
        let v = build();
        let result = v
            .validate_damping_like_field_linearity(TOL)
            .expect("linearity validation should run");
        assert_eq!(result.n_points, CURRENT_DENSITY_A_PER_M2.len());
        assert!(result.max_relative_error.is_finite());
        assert!(result.mean_relative_error.is_finite());
        // The SOT model's `damping_like_field` is rigorously linear in
        // `j_charge`, and so is the embedded reference. After per-curve
        // rescaling, the relative error should be at floating-point noise
        // level.
        assert!(
            result.max_relative_error < 1.0e-9,
            "H_DL should be exactly linear in J after rescaling: {}",
            result.summary()
        );
        assert!(result.passed);
    }

    #[test]
    fn test_damping_like_magnitude_validation_runs() {
        let v = build();
        let result = v
            .validate_damping_like_magnitude(TOL)
            .expect("magnitude validation should run");
        assert_eq!(result.n_points, 1);
        assert!(result.max_relative_error.is_finite());
        // The log-magnitude comparison should land inside a *broad* window
        // even though absolute values are sample-dependent; we explicitly do
        // not require `passed = true` for the bare Pt/Co preset because the
        // Garello sample's effective spin-Hall efficiency may differ from
        // the harness's `θ_SH × transparency × tanh(t/λ_sd)` product.
    }

    #[test]
    fn test_fl_dl_ratio_validation_runs() {
        let v = build();
        let result = v
            .validate_fl_dl_ratio(TOL)
            .expect("ratio validation should run");
        assert_eq!(result.n_points, H_FL_MT.len());
        // By construction of the reference table, the ratio is exactly 0.30,
        // so the per-point error must be zero up to floating-point noise.
        assert!(
            result.max_relative_error < 1.0e-12,
            "H_FL/H_DL ratio is exact by construction: {}",
            result.summary()
        );
        assert!(result.passed);
    }

    #[test]
    fn test_critical_current_validation_runs() {
        let v = build();
        let result = v
            .validate_critical_switching_current(TOL)
            .expect("J_c validation should run");
        assert_eq!(result.n_points, 1);
        assert!(result.max_relative_error.is_finite());
        // J_c must be positive (the SOT critical current cannot be negative
        // — the sign is encoded by the direction of m × σ, not the
        // magnitude). We do not require `passed = true` because the
        // log-magnitude comparison is order-of-magnitude and the harness
        // uses generic Co parameters rather than the Garello sample's.
    }

    #[test]
    fn test_sane_defaults_from_new() {
        let v = build();
        // The default SOT preset is Pt-based, so `θ_SH > 0`.
        assert!((v.sot.theta_sh - 0.07).abs() < 1.0e-12);
        // Cobalt preset has strong PMA.
        assert!(v.ferromagnet.easy_axis.z.abs() > 0.5);
        // The simulated `H_DL` must grow monotonically with `J`.
        let h_low = v.simulated_h_dl_tesla(CURRENT_DENSITY_A_PER_M2[0]);
        let h_high = v.simulated_h_dl_tesla(
            *CURRENT_DENSITY_A_PER_M2
                .last()
                .expect("non-empty current density grid"),
        );
        assert!(
            h_high > h_low,
            "H_DL should grow with J: low {h_low:.3e}, high {h_high:.3e}"
        );
        // The high-current `H_DL` should be finite and non-trivial.
        assert!(h_high.is_finite());
        assert!(h_high > 0.0);
    }
}