diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib
index a8a049ac1..0035b7a39 100644
--- a/crates/ekore/refs.bib
+++ b/crates/ekore/refs.bib
@@ -94,3 +94,16 @@ @article{Buza1996wv
     pages = "301--320",
     year = "1998"
 }
+@article{Bierenbaum2009zt,
+    author = "Bierenbaum, Isabella and Blumlein, Johannes and Klein, Sebastian",
+    title = "{The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production}",
+    eprint = "0901.0669",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "DESY-08-187, SFB-CPP-08-107, IFIC-08-68",
+    doi = "10.1016/j.physletb.2009.01.057",
+    journal = "Phys. Lett. B",
+    volume = "672",
+    pages = "401--406",
+    year = "2009"
+}
diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs
index 01b7025f8..3e8a1eb8f 100644
--- a/crates/ekore/src/bib.rs
+++ b/crates/ekore/src/bib.rs
@@ -1,4 +1,4 @@
-//! List of References (autogenerated on 2024-07-18T17:42:47.650964).
+//! List of References (autogenerated on 2024-09-10T11:22:16.645053).
 
 #[allow(non_snake_case)]
 /// The Three loop splitting functions in QCD: The Nonsinglet case
@@ -95,3 +95,15 @@ pub fn Ball2015tna() {}
 ///
 /// DOI: [10.1007/BF01245820](https:dx.doi.org/10.1007/BF01245820)
 pub fn Buza1996wv() {}
+
+#[allow(non_snake_case)]
+/// The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production
+///
+/// Bierenbaum, Isabella and Blumlein, Johannes and Klein, Sebastian
+///
+/// Published in: Phys. Lett. B 672 (2009), 401--406
+///
+/// e-Print: [0901.0669](https://arxiv.org/abs/0901.0669)
+///
+/// DOI: [10.1016/j.physletb.2009.01.057](https:dx.doi.org/10.1016/j.physletb.2009.01.057)
+pub fn Bierenbaum2009zt() {}
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
index 79b866da8..c9d624e84 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
@@ -1,4 +1,4 @@
-//! |NLO| |QCD|
+//! |NLO| |QCD|.
 
 use num::complex::Complex;
 use num::Zero;
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
index 686260be8..c01f54932 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
@@ -1,17 +1,17 @@
-//! |NNLO| |QCD|
+//! |NNLO| |QCD|.
 
 use num::complex::Complex;
 use num::traits::Pow;
 use num::Zero;
 
-use crate::cmplx;
 use crate::constants::{CA, CF, TR, ZETA2, ZETA3};
 use crate::harmonics::cache::{Cache, K};
 
 use crate::operator_matrix_elements::unpolarized::spacelike::as1;
 
 /// |NNLO| light-light non-singlet |OME|.
-/// It is given in Eq.() of
+///
+/// Implements Eq. (B.4) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -41,8 +41,9 @@ pub fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     CF * TR * (a_qq_l2 * L.pow(2) + a_qq_l1 * L + a_qq_l0)
 }
 
-/// |NNLO| heavy-light pure-singlet |OME|
-/// It is given in Eq.() of
+/// |NNLO| heavy-light pure-singlet |OME|.
+///
+/// Implements Eq. (B.1) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S2 = c.get(K::S1);
@@ -80,8 +81,10 @@ pub fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     CF * TR * (a_hq_l2 * L.pow(2) + a_hq_l1 * L + a_hq_l0)
 }
 
-/// |NNLO| heavy-gluon |OME|
-/// It is given in Eq.() of
+/// |NNLO| heavy-gluon |OME|.
+///
+/// Implements Eq. (B.3) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
+/// The expession for ``A_Hg_l0`` comes form [\[Bierenbaum:2009zt\]](crate::bib::Bierenbaum2009zt).
 pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -179,8 +182,9 @@ pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     a_hg_l2 * L.pow(2) + a_hg_l1 * L + a_hg_l0
 }
 
-/// |NNLO| gluon-quark |OME|
-/// It is given in Eq.() of
+/// |NNLO| gluon-quark |OME|.
+///
+/// Implements Eq. (B.5) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -205,8 +209,9 @@ pub fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     CF * TR * (a_gq_l2 * L.pow(2) + a_gq_l1 * L + a_gq_l0)
 }
 
-/// |NNLO| gluon-gluon |OME|
-/// It is given in Eq.() of
+/// |NNLO| gluon-gluon |OME|.
+///
+/// Implements Eq. (B.7) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -254,7 +259,7 @@ pub fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     a_gg_l2 * L.pow(2) + a_gg_l1 * L + a_gg_l0
 }
 
-/// |NNLO| singlet |OME|
+/// |NNLO| singlet |OME|.
 pub fn A_singlet(c: &mut Cache, nf: u8, L: f64, is_msbar: bool) -> [[Complex<f64>; 3]; 3] {
     let A_hq_2 = A_hq_ps(c, nf, L);
     let A_qq_2 = A_qq_ns(c, nf, L);
@@ -274,7 +279,7 @@ pub fn A_singlet(c: &mut Cache, nf: u8, L: f64, is_msbar: bool) -> [[Complex<f64
     ]
 }
 
-/// |NNLO| non-singlet |OME|
+/// |NNLO| non-singlet |OME|.
 pub fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
     [
         [A_qq_ns(c, nf, L), Complex::<f64>::zero()],