deep_causality_physics 0.7.0

/*
 * SPDX-License-Identifier: MIT
 * Copyright (c) 2023 - 2026. The DeepCausality Authors and Contributors. All Rights Reserved.
 */

//! Quantum Chromodynamics (QCD) kernels for SU(3) gauge theory.
//!
//! Provides fundamental QCD operations including:
//! - Gell-Mann matrices (SU(3) generators)
//! - Structure constants (f^{abc})
//! - Covariant derivative
//! - Wilson loop (confinement order parameter)
//! - Confinement potential (string tension)
//! - Running coupling (asymptotic freedom)
//!
//! All kernels are generic over a real field `R: RealField`. The `f64`-defined
//! numeric constants (½, √3, 4π, …) are cast into `R` via [`crate::real_from_f64`].

use crate::PhysicsError;
use crate::real_from_f64;
use deep_causality_algebra::RealField;
use deep_causality_num::FromPrimitive;

// ============================================================================
// Constants: Gell-Mann Matrices
// ============================================================================

/// The 8 Gell-Mann matrices λ^a (a = 1..8), generators of SU(3).
///
/// Each matrix is 3x3 complex, stored as [Re(00), Im(00), Re(01), Im(01), ...]
/// for row-major order with 9 complex entries = 18 f64 values.
///
/// For simplicity, we store only the real parts since most components are real.
/// Non-zero imaginary parts are handled via the structure constants.
///
/// Returns: 8 matrices, each as `[3][3]` = 9 real values (imaginary handled separately).
pub fn gell_mann_matrices<R: RealField + FromPrimitive>() -> [[R; 9]; 8] {
    let zero = R::zero();
    let one = R::one();

    // λ_1: off-diagonal in (1,2) positions
    let lambda1 = [zero, one, zero, one, zero, zero, zero, zero, zero];

    // λ_2: anti-symmetric imaginary (1,2) - stored as pure real with sign convention
    // Actually has imaginary parts: (0,-i,0), (i,0,0), (0,0,0)
    // We store the magnitude; use structure constants for phases
    let lambda2 = [zero, zero, zero, zero, zero, zero, zero, zero, zero]; // Imaginary - use structure constants

    // λ_3: diagonal (1, -1, 0)
    let lambda3 = [one, zero, zero, zero, -one, zero, zero, zero, zero];

    // λ_4: off-diagonal in (1,3) positions
    let lambda4 = [zero, zero, one, zero, zero, zero, one, zero, zero];

    // λ_5: anti-symmetric imaginary (1,3)
    let lambda5 = [zero, zero, zero, zero, zero, zero, zero, zero, zero]; // Imaginary

    // λ_6: off-diagonal in (2,3) positions
    let lambda6 = [zero, zero, zero, zero, zero, one, zero, one, zero];

    // λ_7: anti-symmetric imaginary (2,3)
    let lambda7 = [zero, zero, zero, zero, zero, zero, zero, zero, zero]; // Imaginary

    // λ_8: diagonal (1, 1, -2) / sqrt(3)
    let inv_sqrt3 = real_from_f64::<R>(1.0 / 3.0_f64.sqrt());
    let lambda8 = [
        inv_sqrt3,
        zero,
        zero,
        zero,
        inv_sqrt3,
        zero,
        zero,
        zero,
        real_from_f64::<R>(-2.0) * inv_sqrt3,
    ];

    [
        lambda1, lambda2, lambda3, lambda4, lambda5, lambda6, lambda7, lambda8,
    ]
}

// ============================================================================
// Constants: SU(3) Structure Constants
// ============================================================================

/// Returns the totally antisymmetric SU(3) structure constants f^{abc}.
///
/// Non-zero values (and their antisymmetric permutations):
/// - f^{123} = 1
/// - f^{147} = f^{165} = f^{246} = f^{257} = f^{345} = f^{376} = 1/2
/// - f^{458} = f^{678} = √3/2
///
/// Returns: Function that computes f^{abc} for given indices (1-indexed).
pub fn structure_constant<R: RealField + FromPrimitive>(a: usize, b: usize, c: usize) -> R {
    // Canonical non-zero structure constants (1-indexed in physics convention)
    // We convert to 0-indexed internally
    let half = real_from_f64::<R>(0.5);
    let sqrt3_half = real_from_f64::<R>(3.0_f64.sqrt() * 0.5);

    // Normalize to sorted order with sign tracking for antisymmetry
    let mut indices = [a, b, c];
    let mut sign = R::one();

    // Bubble sort with sign tracking
    for i in 0..2 {
        for j in 0..2 - i {
            if indices[j] > indices[j + 1] {
                indices.swap(j, j + 1);
                sign = -sign;
            }
        }
    }

    let [i, j, k] = indices;

    // Non-zero structure constants (1-indexed)
    let value = match (i, j, k) {
        (1, 2, 3) => R::one(),
        (1, 4, 7) => half,
        (1, 5, 6) => -half, // Note sign
        (2, 4, 6) => half,
        (2, 5, 7) => half,
        (3, 4, 5) => half,
        (3, 6, 7) => -half, // Note sign
        (4, 5, 8) => sqrt3_half,
        (6, 7, 8) => sqrt3_half,
        _ => R::zero(),
    };

    sign * value
}

/// Returns all non-zero structure constants as a list of (a, b, c, f^abc).
pub fn all_structure_constants<R: RealField + FromPrimitive>() -> Vec<(usize, usize, usize, R)> {
    let half = real_from_f64::<R>(0.5);
    let sqrt3_half = real_from_f64::<R>(3.0_f64.sqrt() * 0.5);

    vec![
        (1, 2, 3, R::one()),
        (1, 4, 7, half),
        (1, 6, 5, half), // Permutation of (1,5,6)
        (2, 4, 6, half),
        (2, 5, 7, half),
        (3, 4, 5, half),
        (3, 7, 6, half), // Permutation of (3,6,7)
        (4, 5, 8, sqrt3_half),
        (6, 7, 8, sqrt3_half),
    ]
}

// ============================================================================
// Kernels
// ============================================================================

/// Computes the gauge covariant derivative: $D_\mu \psi = \partial_\mu \psi + i g A_\mu \psi$.
///
/// For SU(3), $A_\mu = A_\mu^a T^a$ where $T^a = \lambda^a / 2$.
///
/// # Arguments
/// * `psi` - Field value (color triplet, 3 complex components = 6 real values).
/// * `psi_gradient` - Ordinary derivative $\partial_\mu \psi$ (4 spacetime × 6 color = 24 values).
/// * `gluon_field` - Gluon potential $A_\mu^a$ (4 spacetime × 8 color = 32 values).
/// * `coupling` - QCD coupling constant $g$.
///
/// # Returns
/// * `Ok(Vec<R>)` - Covariant derivative $D_\mu \psi$ (4 × 6 = 24 values).
pub fn covariant_derivative_kernel<R>(
    psi: &[R],          // 6 values: (Re, Im) for each color
    psi_gradient: &[R], // 24 values: 4 spacetime × 6 color
    gluon_field: &[R],  // 32 values: 4 spacetime × 8 color adjoint
    coupling: R,
) -> Result<Vec<R>, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    // Validate input dimensions
    if psi.len() != 6 {
        return Err(PhysicsError::DimensionMismatch(format!(
            "Psi must have 6 components (3 complex), got {}",
            psi.len()
        )));
    }
    if psi_gradient.len() != 24 {
        return Err(PhysicsError::DimensionMismatch(format!(
            "Psi gradient must have 24 components, got {}",
            psi_gradient.len()
        )));
    }
    if gluon_field.len() != 32 {
        return Err(PhysicsError::DimensionMismatch(format!(
            "Gluon field must have 32 components, got {}",
            gluon_field.len()
        )));
    }

    let zero = R::zero();
    let half = real_from_f64::<R>(0.5);
    let matrices = gell_mann_matrices::<R>();
    let mut result = vec![zero; 24];

    // For each spacetime direction mu
    for mu in 0..4 {
        // Start with ordinary derivative
        for c in 0..6 {
            result[mu * 6 + c] = psi_gradient[mu * 6 + c];
        }

        // Add gauge term: i * g * A_mu^a * (lambda^a / 2) * psi
        // Simplified: we add the diagonal contributions from λ_3 and λ_8
        for a in 0..8 {
            let a_mu_a = gluon_field[mu * 8 + a];
            let lambda = &matrices[a];

            // Matrix-vector multiply (lambda/2) * psi
            // For real matrices acting on complex vector
            for i in 0..3 {
                for j in 0..3 {
                    let m_ij = lambda[i * 3 + j] * half; // T^a = λ^a / 2
                    if m_ij != zero {
                        // Multiply complex: (0 + i * g * A * m_ij) * psi[j]
                        // = i * g * A * m_ij * (psi_re + i * psi_im)
                        // = i * g * A * m_ij * psi_re - g * A * m_ij * psi_im
                        let psi_re = psi[j * 2];
                        let psi_im = psi[j * 2 + 1];
                        let factor = coupling * a_mu_a * m_ij;

                        // Re part: -g * A * m * psi_im
                        result[mu * 6 + i * 2] -= factor * psi_im;
                        // Im part: +g * A * m * psi_re
                        result[mu * 6 + i * 2 + 1] += factor * psi_re;
                    }
                }
            }
        }
    }

    // Check for numerical issues
    if result.iter().any(|v| !v.is_finite()) {
        return Err(PhysicsError::NumericalInstability(
            "Non-finite result in covariant derivative".into(),
        ));
    }

    Ok(result)
}

/// Computes a simplified Wilson loop trace for confinement analysis.
///
/// The Wilson loop is $W(C) = \text{Tr}[\mathcal{P} \exp(i g \oint_C A_\mu dx^\mu)]$.
///
/// For a rectangular loop of area $A$, in confining phase: $\langle W \rangle \sim \exp(-\sigma A)$.
///
/// This simplified version computes: $W \approx \text{Tr}[\exp(i g \sum_i A_i \cdot \Delta x_i)]$
///
/// # Arguments
/// * `gluon_values` - Gluon field values at each path segment (N segments × 8 color).
/// * `path_lengths` - Length of each path segment (N values).
/// * `coupling` - QCD coupling constant $g$.
///
/// # Returns
/// * `Ok(R)` - Wilson loop trace (should be ~3 for trivial loop, decays for confinement).
pub fn wilson_loop_kernel<R>(
    gluon_values: &[R],
    path_lengths: &[R],
    coupling: R,
) -> Result<R, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    let three = real_from_f64::<R>(3.0);
    let num_segments = path_lengths.len();
    if gluon_values.len() != num_segments * 8 {
        return Err(PhysicsError::DimensionMismatch(format!(
            "Expected {} gluon values for {} segments, got {}",
            num_segments * 8,
            num_segments,
            gluon_values.len()
        )));
    }

    if num_segments == 0 {
        return Ok(three); // Tr(I) = 3 for SU(3)
    }

    // For small loops, approximate: W ≈ 1 - (g²/2) * sum(A² * dl²)
    // This captures the area law behavior for confinement
    let mut phase_sum = R::zero();

    for i in 0..num_segments {
        let dl = path_lengths[i];

        // Sum of |A|² at this segment
        let mut a_squared = R::zero();
        for a in 0..8 {
            let a_a = gluon_values[i * 8 + a];
            a_squared += a_a * a_a;
        }

        phase_sum += a_squared * dl * dl;
    }

    // Wilson loop with quadratic approximation
    let neg_half = real_from_f64::<R>(-0.5);
    let w = three * (neg_half * coupling * coupling * phase_sum).exp();

    if !w.is_finite() {
        return Err(PhysicsError::NumericalInstability(
            "Non-finite Wilson loop result".into(),
        ));
    }

    Ok(w)
}

/// Computes the linear confining potential: $V(r) = \sigma r + V_0$.
///
/// The string tension $\sigma \approx 0.18 \text{ GeV}^2$ in physical units.
///
/// # Arguments
/// * `distance` - Separation distance $r$ in fm (or natural units).
/// * `string_tension` - String tension $\sigma$ in GeV² (typically ~0.18).
/// * `coulomb_term` - Optional Coulomb coefficient for short-range: $-\alpha/r$.
///
/// # Returns
/// * `Ok(R)` - Potential energy $V(r)$ in GeV.
pub fn confinement_potential_kernel<R>(
    distance: R,
    string_tension: R,
    coulomb_term: Option<R>,
) -> Result<R, PhysicsError>
where
    R: RealField,
{
    let zero = R::zero();
    if distance <= zero {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "Distance must be positive".into(),
        ));
    }

    if !distance.is_finite() || !string_tension.is_finite() {
        return Err(PhysicsError::NumericalInstability(
            "Non-finite input to confinement potential".into(),
        ));
    }

    // Linear confining term
    let mut v = string_tension * distance;

    // Optional Coulomb term for short-range behavior
    if let Some(alpha) = coulomb_term {
        if !alpha.is_finite() {
            return Err(PhysicsError::NumericalInstability(
                "Non-finite Coulomb coefficient".into(),
            ));
        }
        v -= alpha / distance;
    }

    Ok(v)
}

/// Computes the running QCD coupling constant $\alpha_s(Q^2)$ (asymptotic freedom).
///
/// One-loop beta function: $\alpha_s(Q^2) = \frac{4\pi}{(11 - 2n_f/3) \ln(Q^2/\Lambda_{QCD}^2)}$
///
/// # Arguments
/// * `q_squared` - Momentum transfer squared $Q^2$ in GeV².
/// * `lambda_qcd` - QCD scale $\Lambda_{QCD}$ in GeV (typically ~0.2).
/// * `n_flavors` - Number of active quark flavors (typically 3-6).
///
/// # Returns
/// * `Ok(R)` - Running coupling $\alpha_s(Q^2)$.
pub fn running_coupling_kernel<R>(
    q_squared: R,
    lambda_qcd: R,
    n_flavors: u32,
) -> Result<R, PhysicsError>
where
    R: RealField + FromPrimitive,
{
    let zero = R::zero();
    if q_squared <= zero {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "Q² must be positive".into(),
        ));
    }

    if lambda_qcd <= zero {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "Lambda_QCD must be positive".into(),
        ));
    }

    let lambda_squared = lambda_qcd * lambda_qcd;
    if q_squared <= lambda_squared {
        return Err(PhysicsError::PhysicalInvariantBroken(
            "Q² must be greater than Λ_QCD² for perturbative regime".into(),
        ));
    }

    // One-loop beta function coefficient (computed in f64 since n_flavors is integral).
    let b0_f64 = 11.0 - (2.0 * n_flavors as f64) / 3.0;
    if b0_f64 <= 0.0 {
        return Err(PhysicsError::PhysicalInvariantBroken(format!(
            "Too many flavors ({}): b0 = {} <= 0",
            n_flavors, b0_f64
        )));
    }
    let b0 = real_from_f64::<R>(b0_f64);

    let log_ratio = (q_squared / lambda_squared).ln();
    let four_pi = real_from_f64::<R>(4.0 * std::f64::consts::PI);
    let alpha_s = four_pi / (b0 * log_ratio);

    if !alpha_s.is_finite() || alpha_s <= zero {
        return Err(PhysicsError::NumericalInstability(
            "Invalid running coupling (non-finite or non-positive)".into(),
        ));
    }

    Ok(alpha_s)
}