From 80edf057ec2a0d3813d477dd556700fd95c6274e Mon Sep 17 00:00:00 2001 From: Felix Hekhorn Date: Sat, 21 Sep 2024 16:01:01 +0300 Subject: [PATCH] rust: Init ad.p.s.as2 --- .../polarized/spacelike.rs | 30 +-- .../polarized/spacelike/as2.rs | 206 ++++++++++++++++++ 2 files changed, 221 insertions(+), 15 deletions(-) create mode 100644 crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs index 8526b53b1..bb394deda 100644 --- a/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs +++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs @@ -5,23 +5,23 @@ use crate::harmonics::cache::Cache; use num::complex::Complex; use num::Zero; pub mod as1; -// pub mod as2; +pub mod as2; // pub mod as3; /// Compute the tower of the non-singlet anomalous dimensions. pub fn gamma_ns_qcd(order_qcd: usize, mode: u16, c: &mut Cache, nf: u8) -> Vec> { let mut gamma_ns = vec![Complex::::zero(); order_qcd]; gamma_ns[0] = as1::gamma_ns(c, nf); - // // NLO and beyond - // if order_qcd >= 2 { - // let gamma_ns_1 = match mode { - // PID_NSP => as2::gamma_nsp(c, nf), - // // To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here - // PID_NSM | PID_NSV => as2::gamma_nsm(c, nf), - // _ => panic!("Unkown non-singlet sector element"), - // }; - // gamma_ns[1] = gamma_ns_1 - // } + // NLO and beyond + if order_qcd >= 2 { + let gamma_ns_1 = match mode { + PID_NSP => as2::gamma_nsp(c, nf), + // To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here + PID_NSM | PID_NSV => as2::gamma_nsm(c, nf), + _ => panic!("Unkown non-singlet sector element"), + }; + gamma_ns[1] = gamma_ns_1 + } // // NNLO and beyond // if order_qcd >= 3 { // let gamma_ns_2 = match mode { @@ -45,10 +45,10 @@ pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Compl order_qcd ]; gamma_S[0] = as1::gamma_singlet(c, nf); - // // NLO and beyond - // if order_qcd >= 2 { - // gamma_S[1] = as2::gamma_singlet(c, nf); - // } + // NLO and beyond + if order_qcd >= 2 { + gamma_S[1] = as2::gamma_singlet(c, nf); + } // // NNLO and beyond // if order_qcd >= 3 { // gamma_S[2] = as3::gamma_singlet(c, nf); diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs new file mode 100644 index 000000000..3b3779b0f --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs @@ -0,0 +1,206 @@ +//! |NLO| |QCD|. + +use num::complex::Complex; +use num::Zero; + +use super::super::super::unpolarized::spacelike::as2::gamma_nsm as unpol_nsm; +use super::super::super::unpolarized::spacelike::as2::gamma_nsp as unpol_nsp; +use crate::constants::{CA, CF, TR, ZETA2, ZETA3}; +use crate::harmonics::cache::{Cache, K}; + +pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex { + unpol_nsp(c, nf) +} + +pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex { + unpol_nsm(c, nf) +} + +/// Compute the pure-singlet quark-quark anomalous dimension. +/// +/// Implements Eq. (A.3) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { + let n = c.n(); + let gqqps1_nfcf = (2. * (n + 2.) * (1. + 2. * n + n.powu(3))) / ((1. + n).powu(3) * n.powu(3)); + 4.0 * TR * (nf as f64) * CF * gqqps1_nfcf +} + +/// Compute the quark-gluon singlet anomalous dimension. +/// +/// Implements Eq. (A.4) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { + let n = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let Sp2m = c.get(K::S2mh); + + #[rustfmt::skip] + let gqg1_nfca = ( + (S1.powu(2) - S2 + Sp2m) * (n - 1.) / (n * (n + 1.)) + - 4. * S1 / (n * (1. + n).powu(2)) + - (-2. - 7. * n + 3. * n.powu(2) - 4. * n.powu(3) + n.powu(4) + n.powu(5)) / (n.powu(3) * (1. + n).powu(3)) + ) * 2.0; + #[rustfmt::skip] + let gqg1_nfcf = ( + (-(S1.powu(2)) + S2 + 2. * S1 / n) * (n - 1.) / (n * (n + 1.)) + - (n - 1.) + * (1. + 3.5 * n + 4. * n.powu(2) + 5. * n.powu(3) + 2.5 * n.powu(4)) + / (n.powu(3) * (1. + n).powu(3)) + + 4. * (n - 1.) / (n.powu(2) * (1. + n).powu(2)) + ) * 2.; + 4.0 * TR * (nf as f64) * (CA * gqg1_nfca + CF * gqg1_nfcf) +} + +/// Compute the gluon-quark singlet anomalous dimension. +/// +/// Implements Eq. (A.5) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { + Complex::zero() + // let N = c.n(); + // let S1 = c.get(K::S1); + // let S2 = c.get(K::S2); + // let Sp2m = c.get(K::S2mh); + // #[rustfmt::skip] + // let ggq1_cfcf = ( + // (2 * (S1**2 + S2) * (n + 2)) / (n * (n + 1)) + // - (2 * S1 * (n + 2) * (1 + 3 * n)) / (n * (1 + n) ** 2) + // - ((n + 2) * (2 + 15 * n + 8 * n**2 - 12.0 * n**3 - 9.0 * n**4)) + // / (n**3 * (1 + n) ** 3) + // + 8 * (n + 2) / (n**2 * (1 + n) ** 2) + // ) * 0.5; + // #[rustfmt::skip] + // let ggq1_cfca = -(-(-(S1 * *2) - S2 + Sp2m) * (n + 2) / (n * (n + 1)) + // - S1 * (12 + 22 * n + 11 * n * *2) / (3 * n * *2 * (n + 1)) + // + (36 + 72 * n + 41 * n * *2 + 254 * n * *3 + 271 * n * *4 + 76 * n * *5) + // / (9 * n * *3 * (1 + n) * *3)); + // #[rustfmt::skip] + // let ggq1_cfnf = + // (-S1 * (n + 2)) / (3 * n * (n + 1)) + ((n + 2) * (2 + 5 * n)) / (9 * n * (1 + n) * *2); + // 4 * CF * (CA * ggq1_cfca + CF * ggq1_cfcf + 4.0 * TR * nf * ggq1_cfnf) +} + +/// Compute the gluon-gluon singlet anomalous dimension. +/// +/// Implements Eq. (A.6) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex { + Complex::zero() + // let n = c.n(); + // let S1 = c.get(K::S1); + // let Sp1m = c.get(K::S1mh); + // let Sp2m = c.get(K::S2mh); + // let Sp3m = c.get(K::S3mh); + // let S1h = c.get(K::S1h); + // let g3n = c.get(K::G3); + // let SSCHLM = ZETA2 / 2 * (Sp1m - S1h + 2. / n) - S1 / n * *2 - g3 - 5 * ZETA3 / 8.; + // #[rustfmt::skip] + // let ggg1_caca = ( + // -4 * S1 * Sp2m + // - Sp3m + // + 8 * SSCHLM + // + 8 * Sp2m / (n * (n + 1)) + // + 2.0 + // * S1 + // * (72 + 144 * n + 67 * n**2 + 134 * n**3 + 67 * n**4) + // / (9 * n**2 * (n + 1) ** 2) + // - (144 + 258 * n + 7 * n**2 + 698 * n**3 + 469 * n**4 + 144 * n**5 + 48 * n**6) + // / (9 * n**3 * (1 + n) ** 3) + // ) * 0.5; + // #[rustfmt::skip] + // let ggg1_canf = ( + // -5 * S1 / 9 + // + (-3 + 13 * n + 16 * n**2 + 6 * n**3 + 3 * n**4) / (9 * n**2 * (1 + n) ** 2) + // ) * 4; + // #[rustfmt::skip] + // let ggg1_cfnf = (4 + 2 * n - 8 * n**2 + n**3 + 5 * n**4 + 3 * n**5 + n**6) / ( + // n**3 * (1 + n) ** 3 + // ); + // 4 * (CA * *2 * ggg1_caca + TR * nf * (CA * ggg1_canf + CF * ggg1_cfnf)) +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let gamma_qq = gamma_nsp(c, nf) + gamma_ps(c, nf); + [ + [gamma_qq, gamma_qg(c, nf)], + [gamma_gq(c, nf), gamma_gg(c, nf)], + ] +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::harmonics::cache::Cache; + use crate::{assert_approx_eq_cmplx, cmplx}; + use num::complex::Complex; + use num::traits::Pow; + use num::Zero; + use std::f64::consts::PI; + + const NF: u8 = 5; + + #[test] + fn physical_constraints() { + // qg_helicity_conservation + let mut c = Cache::new(cmplx!(1., 0.)); + assert_approx_eq_cmplx!(f64, gamma_qg(&mut c, NF), Complex::zero(), epsilon = 2e-6); + + // // qg momentum + // let mut c = Cache::new(cmplx!(1., 0.)); + // let gS1 = gamma_singlet(&mut c, NF); + // assert_approx_eq_cmplx!(f64, gS1[0][0], cmplx!(12. * TR * (NF as f64) * CF, 0.), epsilon = 1e-6); + } + + // #[test] + // fn N2() { + // // reference values are obtained from MMa + // let mut c = Cache::new(cmplx!(2., 0.)); + + // // ns+ + // assert_approx_eq_cmplx!( + // f64, + // gamma_nsp(&mut c, NF), + // cmplx!( + // (-112.0 * CF + 376.0 * CA - 64.0 * (NF as f64)) * CF / 27.0, + // 0. + // ), + // epsilon = 2e-6 + // ); + + // // ns- + // let check = ((34.0 / 27.0 * (-47.0 + 6. * PI.pow(2)) - 16.0 * ZETA3) * CF + // + (373.0 / 9.0 - 34.0 * PI.pow(2) / 9.0 + 8.0 * ZETA3) * CA + // - 64.0 * (NF as f64) / 27.0) + // * CF; + // assert_approx_eq_cmplx!( + // f64, + // gamma_nsm(&mut c, NF), + // cmplx!(check, 0.), + // epsilon = 2e-6 + // ); + + // // singlet sector + // let gS1 = gamma_singlet(&mut c, NF); + // // ps + // assert_approx_eq_cmplx!( + // f64, + // gamma_ps(&mut c, NF), + // cmplx!(-40.0 * CF * (NF as f64) / 27.0, 0.) + // ); + // // qg + // assert_approx_eq_cmplx!( + // f64, + // gS1[0][1], + // cmplx!((-74.0 * CF - 35.0 * CA) * (NF as f64) / 27.0, 0.) + // ); + // // gq + // assert_approx_eq_cmplx!( + // f64, + // gS1[1][0], + // cmplx!( + // (112.0 * CF - 376.0 * CA + 104.0 * (NF as f64)) * CF / 27.0, + // 0. + // ), + // epsilon = 1e-13 + // ); + // } +}