Rustify QED

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

@@ -1,10 +1,14 @@
 //! The unpolarized, space-like anomalous dimensions at various couplings power.
 
-use crate::constants::{PID_NSM, PID_NSP, PID_NSV};
+use crate::constants::{
+    ed2, eu2, PID_NSM, PID_NSM_ED2, PID_NSM_EU2, PID_NSP, PID_NSP_ED2, PID_NSP_EU2, PID_NSV,
+};
 use crate::harmonics::cache::Cache;
 use num::complex::Complex;
 use num::Zero;
+pub mod aem1;
 pub mod as1;
+pub mod as1aem1;
 pub mod as2;
 pub mod as3;
 
@@ -55,3 +59,26 @@ pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Compl
     }
     gamma_S
 }
+
+/// Compute the tower of the QED non-singlet anomalous dimensions.
+pub fn gamma_ns_qed(
+    order_qcd: usize,
+    order_qed: usize,
+    mode: u16,
+    c: &mut Cache,
+    nf: u8,
+) -> Vec<Vec<Complex<f64>>> {
+    let row = vec![Complex::<f64>::zero(); order_qcd + 1];
+    let mut gamma_ns = vec![row; order_qed + 1];
+    gamma_ns[1][0] = as1::gamma_ns(c, nf);
+    gamma_ns[0][1] = choose_ns_as_aem1(mode, c, nf);
+    gamma_ns
+}
+
+pub fn choose_ns_as_aem1(mode: u16, c: &mut Cache, nf: u8) -> Complex<f64> {
+    match mode {
+        PID_NSP_EU2 | PID_NSM_EU2 => eu2 * aem1::gamma_ns(c, nf),
+        PID_NSP_ED2 | PID_NSM_ED2 => ed2 * aem1::gamma_ns(c, nf),
+        _ => panic!("Unkown non-singlet sector element"),
+    }
+}

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/anomalous_dimensions/unpolarized/spacelike/as1aem1.rs

@@ -0,0 +1,142 @@
+//! The $O(a_s^1a_{em}^1)$ Altarelli-Parisi splitting kernels.
+use crate::cmplx;
+use num::complex::Complex;
+
+use crate::constants::{ed2, eu2, uplike_flavors, CA, CF, NC, TR, ZETA2, ZETA3};
+use crate::harmonics::cache::{Cache, K};
+
+/// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension.
+///
+/// Implements Eq. (36) of
+pub fn gamma_phq(c: &mut Cache, _nf: u8) -> Complex<f64> {
+    let N = c.n();
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+
+    #[rustfmt::skip]
+    let tmp_const =
+        2.0
+        * (
+            -4.0
+            - 12.0 * N
+            - N.powu(2)
+            + 28.0 * N.powu(3)
+            + 43.0 * N.powu(4)
+            + 30.0 * N.powu(5)
+            + 12.0 * N.powu(6)
+        ) / ((-1.0 + N) * N.powu(3) * (1.0 + N).powu(3));
+
+    #[rustfmt::skip]
+    let tmp_S1 = -4.0
+        * (10.0 + 27.0 * N + 25.0 * N.powu(2) + 13.0 * N.powu(3) + 5.0 * N.powu(4))
+        / ((-1.0 + N) * N * (1.0 + N).powu(3));
+
+    let tmp_S12 = 4.0 * (2.0 + N + N.powu(2)) / ((-1.0 + N) * N * (1.0 + N));
+    let tmp_S2 = 4.0 * (2.0 + N + N.powu(2)) / ((-1.0 + N) * N * (1.0 + N));
+
+    CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2)
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ quark-photon anomalous dimension.
+///
+/// Implements Eq. (26) of
+pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
+    let N = c.n();
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+
+    #[rustfmt::skip]
+    let tmp_const = -2.0
+        * (4.0
+            + 8.0 * N
+            + 25.0 * N.powu(2)
+            + 51.0 * N.powu(3)
+            + 36.0 * N.powu(4)
+            + 15.0 * N.powu(5)
+            + 5.0 * N.powu(6))
+        / (N.powu(3) * (1.0 + N).powu(3) * (2.0 + N));
+
+    let tmp_S1 = 8.0 / N.powu(2);
+    let tmp_S12 = -4.0 * (2.0 + N + N.powu(2)) / (N * (1.0 + N) * (2.0 + N));
+    let tmp_S2 = 4.0 * (2.0 + N + N.powu(2)) / (N * (1.0 + N) * (2.0 + N));
+
+    2.0 * (nf as f64) * CA * CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2)
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ gluon-photon anomalous dimension.
+///
+/// Implements Eq. (27) of
+pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex<f64> {
+    let N = c.n();
+    CF * CA
+        * (8.0 * (-4.0 + N * (-4.0 + N * (-5.0 + N * (-10.0 + N + 2.0 * N.powu(2) * (2.0 + N))))))
+        / (N.powu(3) * (1.0 + N).powu(3) * (-2.0 + N + N.powu(2)))
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ photon-gluon anomalous dimension.
+///
+/// Implements Eq. (30) of
+pub fn gamma_phg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    TR / CF / CA * (NC as f64) * gamma_gph(c, nf)
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ quark-gluon singlet anomalous dimension.
+///
+/// Implements Eq. (29) of
+pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
+    TR / CF / CA * (NC as f64) * gamma_qph(c, nf)
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ gluon-quark singlet anomalous dimension.
+///
+/// Implements Eq. (35) of
+pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
+    gamma_phq(c, nf)
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ photon-photon singlet anomalous dimension.
+///
+/// Implements Eq. (28) of
+pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
+    let nu = uplike_flavors(nf);
+    let nd = nf - nu;
+    cmplx!(4.0 * CF * CA * ((nu as f64) * eu2 + (nd as f64) * ed2), 0.)
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ gluon-gluon singlet anomalous dimension.
+///
+/// Implements Eq. (31) of
+pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex<f64> {
+    cmplx!(4.0 * TR * (NC as f64), 0.)
+}
+
+/// Compute the $O(a_s^1a_{em}^1)$ singlet-like non singlet anomalous dimension.
+///
+/// Implements Eqs. (33-34) of
+pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex<f64> {
+    let N = c.n();
+    let S1 = c.get(K::S1);
+    let S2 = c.get(K::S2);
+    let S3 = c.get(K::S3);
+    let S1h = c.get(K::S1h);
+    let S2h = c.get(K::S2h);
+    let S3h = c.get(K::S3h);
+    let S1p1h = c.get(K::S1ph);
+    let S2p1h = c.get(K::S2ph);
+    let S3p1h = c.get(K::S3ph);
+    let g3N = c.get(K::G3);
+    let g3Np2 = c.get(K::G3p2);
+
+    #[rustfmt::skip]
+    let result = 32.0 * ZETA2 * S1h - 32.0 * ZETA2 * S1p1h
+        + 8.0 / (N + N.powu(2)) * S2h
+        - 4.0 * S3h + (24.0 + 16.0 / (N + N.powu(2))) * S2
+        - 32.0 * S3 - 8.0 / (N + N.powu(2)) * S2p1h
+        + S1 * (16.0 * (3.0 / N.powu(2) - 3.0 / (1.0 + N).powu(2) + 2.0 * ZETA2) - 16.0 * S2h
+            - 32.0 * S2 + 16.0 * S2p1h )
+        + (-8.0 + N * (-32.0 + N * ( -8.0 - 3.0 * N * (3.0 + N) * (3.0 + N.powu(2)) - 48.0 * (1.0 + N).powu(2) * ZETA2)))
+        / (N.powu(3) * (1.0 + N).powu(3))
+        + 32.0 * (g3N + g3Np2) + 4.0 * S3p1h - 16.0 * ZETA3;
+
+    CF * result
+}

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,43 @@ pub const PID_NSM: u16 = 10201;
 
 /// non-singlet all-valence |PID|.
 pub const PID_NSV: u16 = 10200;
+
+/// QED |PID|. Need to give sensible names
+pub const PID_NSP_EU2: u16 = 10102;
+
+pub const PID_NSP_ED2: u16 = 10103;
+
+pub const PID_NSM_EU2: u16 = 10202;
+
+pub const PID_NSM_ED2: u16 = 10203;
+
+/// 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,
+    }
+}