More work on the approx O(as^3) massive coefficient functions (thresh…

…old approximation implemented)
vbertone · Jul 28, 2024 · 23dff5b · 23dff5b
1 parent 1665ad9
commit 23dff5b
Show file tree

Hide file tree

Showing 2 changed files with 147 additions and 1 deletion.
diff --git a/inc/apfel/massivecoefficientfunctionsunp_sl.h b/inc/apfel/massivecoefficientfunctionsunp_sl.h
@@ -226,6 +226,44 @@ namespace apfel
   };
   ///@}
 
+  /**
+   * @defgroup NNLOthr Approximated NNLO massive coefficient functions near threshold
+   * @ingroup NCMassive
+   */
+  ///@{
+  /**
+   * @brief O(&alpha;<SUB>s</SUB><SUP>3</SUP>) gluon coefficient
+   * function for F2 near
+   * threshold. Reference:https://arxiv.org/pdf/1205.5727.
+   */
+  class CmTh23gNC: public Expression
+  {
+  public:
+    CmTh23gNC(int const& nf, double const& eta);
+    double Regular(double const& x) const;
+  private:
+    int     const _nf;
+    double  const _eta;
+    Cm21gNC const _c21g;
+  };
+
+  /**
+   * @brief O(&alpha;<SUB>s</SUB><SUP>3</SUP>) gluon coefficient
+   * function for FL near
+   * threshold. Reference:https://arxiv.org/pdf/1205.5727.
+   */
+  class CmThL3gNC: public Expression
+  {
+  public:
+    CmThL3gNC(int const& nf, double const& eta);
+    double Regular(double const& x) const;
+  private:
+    int     const _nf;
+    double  const _eta;
+    CmL1gNC const _cL1g;
+  };
+  ///@}
+
   /**
    * @defgroup NLOhq NLO massive heavy-quark-initiated coefficient functions
    * Collection of the massive coefficient functions for processes

diff --git a/src/structurefunctions/massivecoefficientfunctionsunp_sl.cc b/src/structurefunctions/massivecoefficientfunctionsunp_sl.cc
@@ -301,7 +301,115 @@ namespace apfel
     const double z   = eta * x;
     double xi = 4 * eta / ( 1 - eta );
     double wr = xi * ( 1 / z - 1 ) / 4 - 1;
-    return _eta * 16 * M_PI * xi * clnlobarq_(&wr, & xi) / z;
+    return _eta * 16 * M_PI * xi * clnlobarq_(&wr, &xi) / z;
+  }
+
+  //_________________________________________________________________________________
+  // Helper functions for the approximated O(as^3) coefficient
+  // functions. The correspond to Eqs. (3.10) and (3.11) of
+  // https://arxiv.org/pdf/1205.5727.
+  //_________________________________________________________________________________
+  double c0Th(double const& xi)
+  {
+    const double y  = sqrt(1 + 4 / xi);
+    const double L1 = log(1 + xi / 2);
+    const double L2 = log(2 + xi / 2);
+    const double L3 = log(sqrt(xi) * ( y - 1 ) / 2);
+    const double xip2 = 2 + xi;
+    const double xip4 = 4 + xi;
+    const double Li2 = dilog(-2 / xip2);
+    return CA * ( 50 - Pi2 + 12 * L3 / y + 4 * L3 * L3 + L1 * L1 + 6 * L2
+                  - 4 * L2 * L2 + 2 * Li2 + 48 / xip2 - 4 * L2 / xip2
+                  + 64 * L2 / xip2 / xip2 - 128 * L2 / ( xip2 * xip2 * xip4 )
+                  - 160 / xip2 / xip4 - 64 * L2 / xip2 / xip4
+                  + 128 / ( xip2 * xip4 * xip4 ) - 12 * ( 4 + zeta2 ) / xip4
+                  - 8 * L3 * L3 / xip4 + 64 / xip4 / xip4 )
+           + CF * ( -18. - 2. / 3. * Pi2 - 24. * L3 / y - 8. * L3 * L3
+                    + 2. * L1 * L1 - 6. * L2 + 4. * Li2 - 48. / xip2
+                    + 8. * L2 / xip2 + 360. / xip2 / xip4
+                    + 128. * L2 / xip2 / xip4 - 544. / (xip2 * xip4 * xip4)
+                    + 48. * L3 * L3 / xip4 - 8. * L1 * L1 / xip4
+                    + (44. + 40. * zeta2) / xip4 - 120. * L2 / xip2 / xip2
+                    + 256. * L2 / (xip2 * xip2 * xip4) - 16. * Li2 / xip4
+                    - 272. / xip4 / xip4 );
+  }
+  double c0Thbar(double const& xi)
+  {
+    return 4 * CA * ( 2 + log(1 + xi / 4) ) - 4 * TR / 3;
+  }
+  double ThExpansiong3(int const& nf, double const& xi, double const& z)
+  {
+    const double beta = sqrt(1 - 4 * z / xi / ( 1 - z ));
+    const double l  = log(beta);
+    const double l2 = l * l;
+    const double l3 = l * l2;
+    const double l4 = l * l3;
+    const double Lm  = - log(xi);
+    const double Lm2 = Lm * Lm;
+    const double ln2  = log(2);
+    const double ln22 = ln2 * ln2;
+    const double ln23 = ln2 * ln22;
+    const double c_log4 = 128 * CA * CA;
+    const double c_log3 = ( 768 * ln2 - 6464 / 9. ) * CA * CA + 128. / 9. * CA * nf + 128 * CA * CA * Lm;
+    const double c_log2 = ( 1728 * ln22 - 3232 * ln2 - 208 / 3. * Pi2 + 15520 / 9. ) * CA * CA
+                          + ( 64. * ln2 - 640 / 9. ) * CA * nf + 16 * CA * c0Th(xi) + 32 * CA * ( CF - CA / 2 ) * Pi2 / beta
+                          - ( ( - 512 * ln2 + 1136 / 3. ) * CA * CA - 32 / 3. * CA * nf + 16 * CA * c0Thbar(xi) ) * Lm
+                          + 32 * CA * CA * Lm2;
+    const double c_log_const = ( 1728 * ln23 - 4848 * ln22 + 15520 / 3. * ln2 - 208 * Pi2 * ln2
+                                 + 936 * zeta3 + 608. / 3. * Pi2 - 88856. / 27.) * CA * CA
+                               + ( 96 * ln22 - 640 / 3. * ln2 - 16 / 3. * Pi2 + 4592 / 27. ) * CA * nf
+                               - 32 * CF * ( CF - CA / 2 ) * Pi2 + ( 48 * ln2 - 40 ) * CA * c0Th(xi);
+    const double c_log_fracbeta = ( ( - 92 / 3. + 32 * ln2 ) * CA + 8 / 3. * nf) * ( CF - CA / 2 ) * Pi2;
+    const double c_log_Lm = - ( ( - 672 * ln22 + 976 * ln2 + 104 / 3. * Pi2 - 4160 / 9. ) * CA * CA
+                                + ( - 32 * ln2 + 320 / 9. ) * CA * nf + ( 48  * ln2 - 40 ) * CA * c0Thbar(xi)
+                                - 8 * CA * c0Th(xi) - 16 * CA * ( CF - CA / 2 ) * Pi2 / beta );
+    const double c_log_Lm2 = ( 64 * ln2 - 44 / 3. ) * CA * CA + 8 / 3. * CA * nf - 8 * CA * c0Thbar(xi);
+    const double c_log = c_log_const + c_log_fracbeta / beta + c_log_Lm * Lm + c_log_Lm2 * Lm2;
+    const double c_fracbeta = ( ( 8 * ln22 - 68 / 3. * ln2 + 8 / 3. * Pi2 - 658 / 9. ) * CA
+                                + ( 8 / 3. * ln2 - 20 / 9.) * nf + 2 * c0Th(xi)
+                                + ( 26 / 3. * CA + 4 / 3. * nf - 2 * c0Thbar(xi) ) * Lm ) * ( CF - CA / 2 ) * Pi2;
+    const double c_fracbeta2 = 4 / 3. * ( CF - CA / 2 ) * ( CF - CA / 2 ) * Pi2 * Pi2;
+    return c_log4 * l4 + c_log3 * l3 + c_log2 * l2 + c_log * l + c_fracbeta / beta + c_fracbeta2 / beta / beta;
+  }
+  double ThExpansiong3const(double const& xi)
+  {
+    return pow(c0Th(xi) + 36 * CA * pow(log(2), 2) - 60 * CA * log(2) - log(xi) * ( 8 * CA * log(2) - c0Thbar(xi) ), 2);
+  }
+
+  //_________________________________________________________________________________
+  CmTh23gNC::CmTh23gNC(int const& nf, double const& eta):
+    Expression(),
+    _nf(nf),
+    _eta(eta),
+    _c21g(Cm21gNC{eta})
+  {
+  }
+  double CmTh23gNC::Regular(double const& x) const
+  {
+    if (x >= 1)
+      return 0;
+    const double eta = this->_eta;
+    const double z   = eta * x;
+    const double xi  = 4 * eta / ( 1 - eta );
+    return _c21g.Regular(x) * ( ThExpansiong3(_nf, xi, z) + ThExpansiong3const(xi) );
+  }
+
+  //_________________________________________________________________________________
+  CmThL3gNC::CmThL3gNC(int const& nf, double const& eta):
+    Expression(),
+    _nf(nf),
+    _eta(eta),
+    _cL1g(CmL1gNC{eta})
+  {
+  }
+  double CmThL3gNC::Regular(double const& x) const
+  {
+    if (x >= 1)
+      return 0;
+    const double eta = this->_eta;
+    const double z   = eta * x;
+    const double xi  = 4 * eta / ( 1 - eta );
+    return _cL1g.Regular(x) * ( ThExpansiong3(_nf, xi, z) + ThExpansiong3const(xi) );
   }
 
   //_________________________________________________________________________________