Rustify QED

crates/ekore/src/anomalous_dimensions/unpolarized/spacelike.rs

@@ -4,6 +4,7 @@ use crate::constants::{PID_NSM, PID_NSP, PID_NSV};
 use crate::harmonics::cache::Cache;
 use num::complex::Complex;
 use num::Zero;
+pub mod aem1;
 pub mod as1;
 pub mod as2;
 pub mod as3;

crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/aem1.rs

@@ -0,0 +1,130 @@
+//! |LO| |QED|.
+use num::complex::Complex;
+use num::Zero;
+
+use crate::constants::{charge_combinations, ed2, eu2, uplike_flavors, CF, NC, TR};
+use crate::harmonics::cache::Cache;
+
+use crate::anomalous_dimensions::unpolarized::spacelike::as1;
+
+/// Compute the leading-order photon-quark anomalous dimension.
+///
+/// Implements Eq. (2.5) of
+pub fn gamma_phq(c: &mut Cache, nf: u8) -> Complex<f64> {
+    as1::gamma_gq(c, nf) / CF
+}
+
+/// Compute the leading-order quark-photon anomalous dimension.
+///
+/// Implements Eq. (2.5) of
+pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
+    as1::gamma_qg(c, nf) / TR * (NC as f64)
+}
+
+/// Compute the leading-order photon-photon anomalous dimension.
+///
+/// Implements Eq. (2.5) of
+pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
+    let nu = uplike_flavors(nf);
+    let nd = nf - nu;
+    (4.0 / 3.0 * (NC as f64) * ((nu as f64) * eu2 + (nd as f64) * ed2)).into()
+}
+
+/// Compute the leading-order non-singlet QED anomalous dimension
+///
+/// Implements Eq. (2.5) of
+pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex<f64> {
+    as1::gamma_ns(c, nf) / CF
+}
+
+/// Compute the leading-order singlet QED anomalous dimension matrix
+///
+/// Implements Eq. (2.5) of
+pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 4]; 4] {
+    let cc = charge_combinations(nf);
+
+    let gamma_ph_q = gamma_phq(c, nf);
+    let gamma_q_ph = gamma_qph(c, nf);
+    let gamma_nonsinglet = gamma_ns(c, nf);
+
+    [
+        [
+            Complex::<f64>::zero(),
+            Complex::<f64>::zero(),
+            Complex::<f64>::zero(),
+            Complex::<f64>::zero(),
+        ],
+        [
+            Complex::<f64>::zero(),
+            gamma_phph(c, nf),
+            cc.e2avg * gamma_ph_q,
+            cc.vue2m * gamma_ph_q,
+        ],
+        [
+            Complex::<f64>::zero(),
+            cc.e2avg * gamma_q_ph,
+            cc.e2avg * gamma_nonsinglet,
+            cc.vue2m * gamma_nonsinglet,
+        ],
+        [
+            Complex::<f64>::zero(),
+            cc.vde2m * gamma_q_ph,
+            cc.vde2m * gamma_nonsinglet,
+            cc.e2delta * gamma_nonsinglet,
+        ],
+    ]
+}
+
+/// Compute the leading-order valence QED anomalous dimension matrix
+///
+/// Implements Eq. (2.5) of
+pub fn gamma_valence(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
+    let cc = charge_combinations(nf);
+
+    [
+        [cc.e2avg * gamma_ns(c, nf), cc.vue2m * gamma_ns(c, nf)],
+        [cc.vde2m * gamma_ns(c, nf), cc.e2delta * gamma_ns(c, nf)],
+    ]
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::{assert_approx_eq_cmplx, cmplx};
+    use num::complex::Complex;
+    use num::Zero;
+    const NF: u8 = 5;
+
+    #[test]
+    fn number_conservation() {
+        const N: Complex<f64> = cmplx!(1., 0.);
+        let mut c = Cache::new(N);
+        let me = gamma_ns(&mut c, NF);
+        assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
+    }
+
+    #[test]
+    fn quark_momentum_conservation() {
+        const N: Complex<f64> = cmplx!(2., 0.);
+        let mut c = Cache::new(N);
+        let me = gamma_ns(&mut c, NF) + gamma_phq(&mut c, NF);
+        assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
+    }
+
+    #[test]
+    fn photon_momentum_conservation() {
+        const N: Complex<f64> = cmplx!(2., 0.);
+        let mut c = Cache::new(N);
+
+        for nf in 2..7 {
+            let nu = uplike_flavors(nf);
+            let nd = nf - nu;
+            assert_approx_eq_cmplx!(
+                f64,
+                eu2 * gamma_qph(&mut c, nu) + ed2 * gamma_qph(&mut c, nd) + gamma_phph(&mut c, nf),
+                cmplx!(0., 0.),
+                epsilon = 2e-6
+            );
+        }
+    }
+}

crates/ekore/src/constants.rs

@@ -1,4 +1,5 @@
 //! Global constants.
+use std::unimplemented;
 
 /// The number of colors.
 ///
@@ -20,6 +21,16 @@ pub const CA: f64 = NC as f64;
 /// Defaults to $C_F = \frac{N_C^2-1}{2N_C} = 4/3$.
 pub const CF: f64 = ((NC * NC - 1) as f64) / ((2 * NC) as f64);
 
+/// Up quark charge square.
+///
+/// Defaults to $e_u^2 = 4./9$
+pub const eu2: f64 = 4. / 9.;
+
+/// Down quark charge square.
+///
+/// Defaults to $e_d^2 = 1./9$
+pub const ed2: f64 = 1. / 9.;
+
 /// Riemann zeta function at z = 2.
 ///
 /// $\zeta(2) = \pi^2 / 6$.
@@ -41,3 +52,34 @@ pub const PID_NSM: u16 = 10201;
 
 /// non-singlet all-valence |PID|.
 pub const PID_NSV: u16 = 10200;
+
+/// compute the number of up flavors
+pub fn uplike_flavors(nf: u8) -> u8 {
+    if nf > 6 {
+        unimplemented!("Selected nf is not implemented")
+    }
+    nf / 2
+}
+
+pub struct ChargeCombinations {
+    pub e2avg: f64,
+    pub vue2m: f64,
+    pub vde2m: f64,
+    pub e2delta: f64,
+}
+
+pub fn charge_combinations(nf: u8) -> ChargeCombinations {
+    let nu = uplike_flavors(nf) as f64;
+    let nd = (nf as f64) - nu;
+    let e2avg = (nu * eu2 + nd * ed2) / (nf as f64);
+    let vue2m = nu / (nf as f64) * (eu2 - ed2);
+    let vde2m = nd / (nf as f64) * (eu2 - ed2);
+    let e2delta = vde2m - vue2m + e2avg;
+
+    ChargeCombinations {
+        e2avg,
+        vue2m,
+        vde2m,
+        e2delta,
+    }
+}