QQbar_threshold Documentation ## Table of Contents## IntroductionQQbar_threshold is a C++ library that provides functions to compute observables near the production threshold of heavy quark pairs in \(e^+ e^-\) collisions. It includes N The QQbar_threshold library can also be accessed from Wolfram Mathematica through the QQbarThreshold package. ## Installation## LinuxThe easiest way to install QQbar_threshold is via the included installation script. The following software has to be available on the system: - a compiler complying with the C++11 standard (e.g. g++ 4.8 or higher);
- cmake (http://cmake.org/);
- GSL, the GNU scientific library (http://www.gnu.org/software/gsl);
- the QCD library (included in the installation script);
- the QCD library requires odeint from the boost libraries (http://www.boost.org/).
It is recommended to run the installation script in a separate build directory, e.g. wget https://www.hepforge.org/downloads/qqbarthreshold/install.tar.gz tar xzf install.tar.gz cd QQbar_threshold_source ./install.sh /my/path/ If During installation there is the opportunity to change some predefined physical constants like the W and Z masses and the default settings for some of the options discussed in the Options section. The defaults are defined in the header constants.hpp. Changing these values Before using the C++ library part, it may also be necessary to adjust certain environment variables. After installation to the base directory LIBRARY_PATH="/my/path/lib:$LIBRARY_PATH" LD_LIBRARY_PATH="/my/path/lib:$LD_LIBRARY_PATH" CPLUS_INCLUDE_PATH="/my/path/include:$CPLUS_INCLUDE_PATH" ## OS XUnder OS X, QQbar_threshold can be installed like under Linux, apart from two exceptions. First, the environment variable #!/bin/sh MATH_PATH=/Applications/Mathematica.app/Contents/MacOS/ DYLD_LIBRARY_PATH="$DYLD_LIBRARY_PATH:$MATH_PATH" $MATH_PATH/WolframKernel "$@" For the installation of the Mathematica package to work, the script file has to be in the executable path and must be named ## Installing more than one versionIt is possible to install more than one version of QQbar_threshold. This requires QQbar_threshold version 1.1 or later with the exception of version 2.0. The latest version has to be installed last. When using a specific version in C++ code only the header files for the version in question should be included. To prevent these header files from being overwritten by later installations it is therefore recommended to use different installation locations for different versions. The details of linking C++ code to a specific library version depend on the system and compiler in question. For example, with g++ under Linux one can use the flag Within Mathematica ## Basic usage and examplesIn this section, we give a brief overview over QQbar_threshold's main functionality and show several code examples. For the sake of a more accessible presentation we postpone the discussion of most details to later sections. The main observables that can be computed with QQbar_threshold are the total cross section and the energy levels and residues of quarkonium bound states. The cross section is calculated in picobarn, whereas all other dimensionful quantities are given in (powers of) GeV. While the examples below demonstrate the usage of specialised C++ header files, we also provide a header QQbar_threshold.hpp} which exposes all functionality offered by the QQbar_threshold library. Functions related to \(t\bar{t}\) production start with the prefix The C++ examples below have to be compiled with a reasonably recent compiler (complying with the C++11 standard) and linked to the QQbar_threshold library. For example, one could compile the first code snippet g++ -o resonance -std=c++11 resonance.cpp -lQQbar_threshold and run it with ./resonance The code for all examples will also be installed alongside with the library. Assuming again ## Mathematica usageWhile the following main text generally describes the usage in C++ programs, the C++ code examples are also followed by equivalent Mathematica code. After loading the package with Needs["QQbarThreshold`"] an overview over the available symbols can be obtained with Names["QQbarThreshold`*"] ## ResonancesQQbar_threshold can calculate the binding energy of the S-wave \(Q\bar{Q}\) resonances. As a first example, we calculate the binding energy of the \(\Upsilon(1S)\) resonance at leading order. The binding energy \(E^{\text{RS}}_1\) depends on the chosen renormalisation scheme and is defined by \(E^{\text{RS}}_1 = M_{\Upsilon(1S)} - 2m^{\text{RS}}_b\), where \(m^{\text{RS}}_b\) is the bottom-quark mass in the given scheme RS, and \(M_{\Upsilon(1S)}\) refers to the theoretically computed mass of the \(\Upsilon(1S)\). By default, the quark masses are defined in the PS-shift scheme (cf. the Mass schemes section). All masses and energies are in GeV, so the output implies that at leading order \(E^{\text{PS}}_1 = 0.273372\,\)GeV. #include <iostream> #include "QQbar_threshold/energy_levels.hpp" int main(){ namespace QQt = QQbar_threshold; const double mb_PS = 4.5; // bottom mass in PS scheme double E_1 = QQt::bbbar_energy_level( 1, // principal quantum number mb_PS, // renormalisation scale mb_PS, // quark mass QQt::LO // perturbative order ); std::cout << E_1 << '\n'; } Needs["QQbarThreshold`"]; With[ {mbPS = 4.5, scale = 4.5}, E1 = BBbarEnergyLevel[1, scale, mbPS, "LO"] ]; Print[E1]; For an unstable (e.g. top) quark, width corrections can be included by replacing the mass argument by Higher-order corrections can be included by changing the last argument to ## Cross sectionsOne of the most interesting phenomenological applications is the threshold scan of the \(t\bar{t}\) production cross section. The functions \begin{align} \label{eq:rt_def} \sigma_t(s) ={}& \sigma(e^+ e^- \to W^+W^-b\bar{b}X)\,,\\ \label{eq:rb_def} \sigma_b(s) ={}& \sigma(e^+ e^- \to b\bar{b})\,, \end{align} in picobarn. The name #include <iostream> #include "QQbar_threshold/xsection.hpp" int main(){ namespace QQt = QQbar_threshold; std::cout // QQt::ttbar_xsection(sqrt_s, {scales}, {mass, width}, order) << QQt::ttbar_xsection(340., {80., 350.}, {168., 1.4}, QQt::NLO) << '\n'; } Needs["QQbarThreshold`"]; Print[TTbarXSection[340., {80., 350.}, {168., 1.4}, "NLO"]]; Note that in this example two scales appear. As before, the first scale is the overall renormalisation scale. The second scale is due to the separation of resonant and nonresonant contributions to the cross section (cf. the Nonresonant cross section section). It exclusively appears in the top cross section, all other functions only take a single scale as their argument. ## GridsFor N #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" #include <iostream> int main(){ namespace QQt = QQbar_threshold; const double mu = 50.; const double mu_width = 350.; const double mt_PS = 168.; const double width = 1.4; for (double sqrt_s = 330.; sqrt_s < 345.; sqrt_s += 1.0) { std::cout << sqrt_s << '\t' sqrt_s ,{mu, mu_width} ,{mt_PS , width}, QQt::N3LO ) << '\n'; } } Needs["QQbarThreshold`"]; LoadGrid[GridDirectory <> "ttbar_grid.tsv"]; With[ { mu = 50., mtPS = 168., width = 1.4, muWidth = 350., order = "N3LO" }, Do[ Print[ sqrts, "\t", TTbarXSection[sqrts, {mu, muWidth}, {mtPS, width}, order] ], {sqrts, 330., 345., 1.} ] ]; Grids cover a specific range in the rescaled energy and width coordinates \begin{align} \label{eq:grid_coordinates} \tilde{E} ={}& - 4/(\alpha_s^2C_F^2m_Q)E\,,\\ \tilde{\Gamma} ={}& - 4/(\alpha_s^2C_F^2m_Q)\Gamma\,, \end{align} where \(m_Q\) is the pole quark mass, \(E = \sqrt{s} - 2m_Q\) the kinetic energy, and \(C_F = 4/3\). If the arguments of the cross section functions lead to rescaled coordinates outside the range covered by the grid an exception is thrown. In the Mathematica package an error message is displayed and the cross section function will return a symbolic After loading a grid, its coordinate range can be identified as shown in the following example. For the default top grid, the range is \(-19.6793 \le \tilde{E} \le 20.9913\) and \(-2.29592 \le \tilde{\Gamma} \le -0.655977\). For the default bottom grid we would find \(-37.1901 \le \tilde{E} \le -9.29752 \times 10^{-6}\) and \(-9.29752 \times 10^{-8} \le \tilde{\Gamma} \le -4.55432 \times 10^{-15}\) #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" #include <iostream> int main(){ namespace QQt = QQbar_threshold; auto range = QQt::grid_range(); std::cout << "Et range: (" << range.Et_min << ", " << range.Et_max << ")\n"; std::cout << "Gammat range: (" << range.Gammat_min << ", " << range.Gammat_max << ")\n"; } Needs["QQbarThreshold`"]; LoadGrid[GridDirectory <> "ttbar_grid.tsv"]; Print["Et range: ", {GridEtMin[], GridEtMax[]}]; Print["Gammat range: ", {GridGammatMin[], GridGammatMax[]}]; If no grid is loaded 0 will be returned for all coordinate limits. Note that at most one grid can be used at any given time; if a second grid is loaded it will replace the first one. Loading a grid is not thread safe, i.e. one should not try to load more than one grid in parallel. See the Parallelisation section for more details on parallelisation. ## ThresholdsFor convenience, there are also the functions #include <iostream> #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" int main(){ namespace QQt = QQbar_threshold; const double mu = 50.; const double mu_width = 350.; const double mt_PS = 168.; const double width = 1.4; const auto order = QQt::N3LO; for(double E = -3.; E < 5.; E += 1.0){ double sqrt_s = thr + E; std::cout << sqrt_s << '\t' sqrt_s, {mu, mu_width}, {mt_PS, width}, order ) << '\n'; } } Needs["QQbarThreshold`"]; LoadGrid[GridDirectory <> "ttbar_grid.tsv"]; With[ { mu = 50., muWidth = 350., mtPS = 168., width = 1.4, order = "N3LO" }, With[ {thr = TTbarThreshold[mu, mtPS, order]}, Do[ Print[ thr + energy, "\t", TTbarXSection[ thr + energy, {mu, muWidth}, {mtPS, width}, order ] ], {energy, -3., 5., 1.} ] ] ]; For the case that the centre-of-mass energy is not required anywhere else in the program, the C++ #include <iostream> #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" #include "QQbar_threshold/threshold.hpp" int main(){ namespace QQt = QQbar_threshold; std::cout << bbbar_xsection(QQt::threshold()+1e-3, 4.5, 4.5, QQt::N3LO) // equivalent to // << QQt::bbbar_xsection( // QQt::bbbar_threshold(4.5, 4.5, QQt::N3LO) + 1e-3, // 4.5, 4.5, QQt::N3LO // ) << '\n'; } In Mathematica, we can use the symbol Needs["QQbarThreshold`"]; LoadGrid[GridDirectory <> "bbbar_grid.tsv"]; Print[BBbarXSection[QQbarThreshold + 10^-3, 4.5, 4.5, "N3LO"]]; Note that for the bottom quark the width is considered to be zero and thus the continuum cross section is discontinuous at the production threshold. In fact, the expression for the cross section directly at threshold is currently not known to N \(^3\)LO and it is necessary to add a small positive offset to the energy, which was somewhat arbitrarily set to 1 MeV in the above examples. ## Basic optionsThe functions provided by QQbar_threshold support a plethora of optional settings to control their behaviour. In this section we only discuss a small selection; a short summary of all options followed by a comprehensive discussion is given in the Options section. The default settings for the options are given by As an example let us have another look at the \(t\bar{t}\) threshold scan. Here, we change the renormalisation scheme for the mass from the default PS-shift scheme to the MS scheme and discard all Standard Model corrections beyond QCD. #include <iostream> #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" int main(){ namespace QQt = QQbar_threshold; const double mu = 50.; const double mu_width = 350.; const double mt_MS = 160.; // mass mt(mt) in the MSbar scheme const double width = 1.4; QQt::options opt = QQt::top_options(); // MSbar scheme with mu_MSbar = mt_MS opt.mass_scheme = {QQt::MSshift, mt_MS}; // turn off all non-QCD corrections opt.beyond_QCD = QQt::SM::none; for(double sqrt_s = 330.; sqrt_s < 345.; sqrt_s += 1.0){ std::cout << sqrt_s << '\t' sqrt_s, {mu, mu_width}, {mt_MS, width}, QQt::N3LO, opt ) << '\n'; } } Needs["QQbarThreshold`"]; LoadGrid[GridDirectory <> "ttbar_grid.tsv"]; With[ { mu = 50., muWidth = 350., mtMS = 160., (* mass mt(mt) in the MSbar scheme *) width = 1.4, order = "N3LO" }, Do[ Print[ sqrts, "\t", TTbarXSection[ sqrts, {mu, muWidth}, {mtMS, width}, order, MassScheme -> {"MSshift", mtMS}, BeyondQCD -> None ] ], {sqrts, 330., 345., 1.} ] ]; Another useful option allows to modify the reference value for the strong coupling at the scale of the Z boson mass: #include <iostream> #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" int main(){ namespace QQt = QQbar_threshold; const double mu = 50.; const double mu_width = 350.; const double mt_PS = 168.; const double width = 1.4; QQt::options opt = QQt::top_options(); opt.alpha_s_mZ = 0.1174; for( double sqrt_s = 330.; sqrt_s < 345.; sqrt_s += 1.0){ std::cout << sqrt_s << '\t' sqrt_s, {mu, mu_width}, {mt_PS, width}, QQt::N3LO, opt ) << '\n'; } } Needs["QQbarThreshold`"]; LoadGrid[GridDirectory <> "ttbar_grid.tsv"]; With[ { mu = 50., muWidth = 350., mtPS = 168., width = 1.4, order = "N3LO" }, Do[ Print[ sqrts, "\t", TTbarXSection[ sqrts, {mu, muWidth}, {mtPS, width}, order, alphaSmZ -> 0.1174 ] ], {sqrts, 330., 345., 1.} ] ]; For debugging purposes it can be quite helpful to print the current option settings. The following example shows the default options for top-related functions. Note that for some options the default is signified by a physically meaningless sentinel value. #include <iostream> #include "QQbar_threshold/parameters.hpp" int main(){ std::cout << QQbar_threshold::top_options(); } In the Mathematica package the current option settings if different from the default are always visible through their specification in the function call. The default settings for a given function can be inspected with It is also possible to inspect the settings used internally for the actual calculation: #include <iostream> #include "QQbar_threshold/parameters.hpp" int main(){ namespace QQt = QQbar_threshold; const double sqrt_s = 340.; const double mu = 50.; const double mu_width = 350.; const double mt_PS = 168.; const double width = 1.4; const auto order = QQt::N3LO; // internal settings for top cross section std::cout << QQt::top_internal_settings( sqrt_s, {mu, mu_width}, {mt_PS, width}, order, ); // internal settings for energy levels and residues std::cout << QQt::top_internal_settings( mu, mt_PS, order, ); } Needs["QQbarThreshold`"]; With[ { sqrts = 340., mu = 50., muWidth = 350., mtPS = 168., width = 1.4, order = "N3LO" }, Print[ TopInternalSettings[sqrts, {mu, muWidth}, {mtPS, width}, order] ]; Print[TopInternalSettings[mu, mtPS, order]]; ]; Correspondingly, for bottom quarks the function ## Scheme conversionWhile the PS scheme is appropriate for the description of threshold observables, in other kinematic regions schemes like MS may be more suitable. For converting masses to the pole scheme, QQbar_threshold provides the functions #include <iostream> #include "QQbar_threshold/scheme_conversion.hpp" int main(){ namespace QQt = QQbar_threshold; const double mt_PS = 168.; const double mu = 50.; //convert to pole scheme //convert to MSbar scheme QQt::options opt = QQt::top_options(); double mt_MS = mt_PS; double delta_M; do{ opt.mass_scheme = {QQt::MSshift, mt_MS}; delta_M = QQt::top_pole_mass(mu, mt_MS, QQt::N3LO, opt) - mt_Pole; mt_MS -= delta_M; } while(std::abs(delta_M) > precision); std::cout << "pole mass: " << mt_Pole << '\n' << "MSbar mass: " << mt_MS << '\n'; } Needs["QQbarThreshold`"]; With[ {mtPS = 168., mu = 50., precision = 10^-4}, (* convert to pole scheme *) mtPole = TopPoleMass[mu, mtPS, "N3LO"]; (* convert to MSbar scheme *) mtMS = mtPS; For[ deltaM = TopPoleMass[ mu, mtMS, "N3LO", MassScheme -> {"MSshift", mtMS} ] - mtPole, Abs[deltaM] > precision, deltaM = TopPoleMass[ mu, mtMS, "N3LO", MassScheme -> {"MSshift", mtMS} ] - mtPole, mtMS -= deltaM ]; Print["pole mass: ", mtPole]; Print["MSbar mass: ", mtMS]; ]; ## Top quark widthThe width of the top quark is an external parameter of the #include <iostream> #include "QQbar_threshold/width.hpp" int main(){ namespace QQt = QQbar_threshold; const double mt_PS = 171.5; const double mu = 50.; std::cout << QQt::top_width(mu, mt_PS, QQt::N2LO) << '\n'; } Needs["QQbarThreshold`"]; With[ { mtPS = 171.5, mu = 50. }, Print[TopWidth[mu, mtPS, "N2LO"]]; ]; It should be noted that ## Structure of the cross sectionBefore discussing the optional settings in detail, we first give an overview over the structure of the cross section up to N ## Power countingIn PNRQCD, an expansion in \(\alpha_s \sim v \ll 1\) is performed, where \(v = [\sqrt{s}/m_Q - 2]^{1/2}\) is the non-relativistic velocity of the quarks and \(m_Q\) their pole mass. The Coulomb interaction leads to terms scaling with powers of \(\alpha_s/v \sim 1\), which are resummed to all orders. Concerning the electroweak interactions including the Higgs boson, we choose the power counting \(\alpha \sim y_t^2 \sim \alpha_s^2\) for the QED coupling constant and the top Yukawa coupling. We include pure QCD corrections up to N ## Resummation of QED effectsA large correction to the cross section arises from logarithmically enhanced initial-state radiation. I can be accounted for by performing a convolution with structure functions, which are currently known at leading-logarithmic (LL) accuracy. The total production cross section \(\sigma_Q^{\text{ISR}}\) for a \(Q\bar{Q}\) pair including initial-state radiation is then given by \begin{equation} \label{eq:sigma_ISR} \sigma_Q^{\text{ISR}}(s) = \int_0^1 dx_1 \int_0^1 dx_2 \Gamma_{ee}^{\text{LL}}(x_1)\Gamma_{ee}^{\text{LL}}(x_2) \sigma_Q(x_1 x_2 s)\,. \end{equation} This computationally expensive convolution is not done by default in QQbar_threshold. However, there is an example in the section Initial-state radiation. It is also customary to absorb large logarithmic corrections due to vacuum polarisation into a running QED coupling constant \(\alpha(\mu_\alpha)\), which coincides with the fine structure constant \(\alpha \equiv \alpha(0)\) in the Thomson limit. The cross section \(\sigma_Q\) is then proportional to \(\alpha(\mu_\alpha)^2\). We therefore factorise the cross sections \(\sigma_Q\) with \(Q=b,t\) as follows. \begin{equation} \label{eq:sigma_noQED} \sigma_Q = \frac{4\pi\alpha(\mu_\alpha)^2}{3s}R_Q\,. \end{equation} \(R_Q\) then depends on the QED coupling only through higher-order corrections. A further source of large logarithms is given by photon initial state radiation off the electron–positron pair. ## Nonresonant cross sectionSince for top quarks the width is non-negligible, it is necessary to consider the full process \(e^+ e^- \to W^+ W^- b\bar{b}\) instead of just the production of an on-shell \(t\bar{t}\) pair. A systematic analysis in the framework of unstable particle effective theory [5], [6] shows that the cross section can then be written as the sum of resonant and non-resonant production: \begin{equation} \label{eq:sigma_nonres} R_Q(s) = R_{\text{res}}(s) + R_{\text{non-res}}(s)\,. \end{equation} While the resonant part by construction only contains the contributions from top quarks near their mass shell, the invariant mass of the final state \(W\,b\) pair in the nonresonant part can be quite different from the top quark mass. In order to reduce such background contributions, it is possible to specify a cut on the invariant mass (see the Options section). The current implementation in QQbar_threshold includes the nonresonant cross section at N \begin{equation} \label{eq:sigma_nonres_dec} R_{\text{non-res}}(s) = R_{\text{non-res}}^{\text{NLO}}(s) + R_{\text{non-res}}^{\text{N}^2\text{LO}}(s)\,. \end{equation} Both the resonant and the non-resonant part are separately divergent. We remove the poles using MS subtraction and associate the remaining logarithms with a new scale \(\mu_w\) (cf. the Cross sections section). While these logarithms cancel order by order in the sum, a dependence on \(\mu_w\) remains in the present implementation at and N During the evaluation of the nonresonant cross section, interpolation on a precomputed grid is performed. While physical values of the W and the top quark mass are covered by a built-in grid, exotic parameter settings may require the generation of a custom nonresonant grid. This is covered in the Grid generation section. Custom grids can be loaded with ## Production channelsWhile the resonant quark pair is mostly produced in an S wave, there is also a subleading P-wave contribution starting at N \begin{equation} \label{eq:sigma_S_P} R_{\text{res}}(s) = R_S(s) + R_P(s)\,. \end{equation} \(R_S\) and \(R_P\) can be expressed in terms of the imaginary parts of the vector and axialvector polarisation functions, respectively. One obtains \begin{align} \label{eq:sigma_S} R_S(s) ={}& R_{S, \text{QCD}}(s) + R_{S, \text{EW}}(s)\,,\\ \label{eq:sigma_S_QCD} R_{S, \text{QCD}}(s) ={}& \big[{C^{(v)}_0}^2 + {C^{(a)}_0}^2\big]12\pi\,\text{Im}[\Pi_{\text{PR}}^{(v)}(s)] + \delta_{\Gamma_\epsilon/\epsilon} R(s) \,,\\ \label{eq:sigma_P} R_P(s) ={}& a_Q^2\big[a_e^2 + v_e^2\big]\frac{s^2}{(s-m_Z^2)^2}12\pi\,\text{Im}[\Pi_{\text{PR}}^{(a)}(s)]\,,\\ \label{eq:Cp_v} C^{(v)}_0 ={}& e_ee_Q + v_Q\,v_e\frac{s}{s - m_Z^2}\,,\\ \label{eq:Cp_a} C^{(a)}_0 ={}& -v_Q\,a_e\frac{s}{s - m_Z^2}\,, \end{align} where \(v_f\) is the vector coupling of a fermion to the Z boson and \(a_f\) the corresponding axialvector coupling given by \begin{equation} \label{eq:v_f_a_f} v_f = \frac{T_3^f - 2e_fs_w^2}{2s_wc_w}\,, \qquad a_f = \frac{T_3^f}{2s_wc_w}\,. \end{equation} \(e_f\) is the fermion charge in units of the positron charge, \(T_3^f\) its third isospin component, \(c_w = m_W/m_Z\) the cosine of the Weinberg angle, and \(s_w = (1 - c_w^2)^{1/2}\). \(R_{S, \text{EW}}(s)\) is the electroweak correction to S-wave production [19], [20], [21], [22], [3]; more details are given in the Electroweak section. The correction \(\delta_{\Gamma_\epsilon/\epsilon} R(s)\) arises in dimensional regularisation from the product of the \(d=4-2\epsilon\)-dimensional leading-order width with the \(1/\epsilon\) pole in the Coulomb Green function. Its explicit form is given in [3]. ## Pole resummationThe polarisation functions exhibit poles at \(E = E_N - i\Gamma\), where \(E = \sqrt{s} - 2m_Q\) is the kinetic energy, \(\Gamma\) the quark width, and \(E_N\) the (real) binding energy of the \(N\)th bound state. For reasons detailed in [4] the pole contributions to the polarisation functions should be resummed by subtracting the contribution expanded around the leading-order pole position and adding back the unexpanded contributions, i.e. \begin{equation} \label{eq:pole_resum} \Pi_{\text{PR}}^{(v)}(s) = \Pi^{(v)}(s) + \frac{N_C}{2m_Q^2} \sum_{N=1}^\infty\bigg\{ \frac{\big[Z_N\big]_{\text{expanded}}}{\big[E_N - E - i\Gamma\big]_{\text{unexpanded}}} - \bigg[\frac{Z_N}{E_N - E - i\Gamma}\bigg]_{\text{expanded}}\bigg\}\,, \end{equation} and similar for the axialvector polarisation function \(\Pi_{\text{PR}}^{(a)}(s)\) with the P-wave energy levels \(E_N^P\) and residues \(Z_N^P\). A more precise definition of \(Z_N\) is given in the Energy levels and residues section. It should be emphasised that in the limit of a vanishing width, i.e. for bottom quarks, pole resummation has no effect on the (continuum) cross section. In the actual implementation, it is of course not possible to evaluate the sum over the bound states up to infinity. The number of resummed poles is instead set via an option of the cross section functions and defaults to 6. From the scaling of the residues with \(N\) the resulting error on the cross section can be estimated to be comparable to the difference between resumming 4 and 6 poles and is typically at most about 2 per mille. ## Hard matching## QCD and HiggsThe polarisation functions without pole resummation are a product of hard current matching coefficients and the non-relativistic Green functions \(G, G^P\), \begin{align} \label{eq:Pi_v} \Pi^{(v)}(s) ={}& \frac{2N_c}{s}c_v\bigg[c_v - \frac{E + i \Gamma}{m_Q}\frac{d_v}{3}\bigg]G(E) + \dots\,, \\ \label{eq:Pi_a} \Pi^{(a)}(s) ={}& \frac{2N_c}{m_Q^2s}\frac{d-2}{d-1}c_a^2G^P(E) + \dots\,. \end{align} Here \(G(E)\) denotes the Green function at complex energy \(E + i\Gamma\). It is also understood that the products are consistently expanded to N \begin{align} \label{eq:hard_matching_exp} c_v ={}& 1 + \frac{\alpha_s(\mu)}{4\pi} c_v^{(1)} + \bigg(\frac{\alpha_s(\mu)}{4\pi}\bigg)^2 c_v^{(2)} + \bigg(\frac{\alpha_s(\mu)}{4\pi}\bigg)^3 c_v^{(3)}\notag\\ &+ \frac{y_Q^2}{2}\bigg[c_{vH}^{(2)} + \frac{\alpha_s(\mu)}{4\pi}c_{vH}^{(3)}\bigg]\,,\\ \label{eq:dv_exp} d_v ={}& d_v^{(0)} + \frac{\alpha_s(\mu)}{4\pi} d_v^{(1)}\,,\\ \label{eq:ca_exp} c_a ={}& 1 + \frac{\alpha_s(\mu)}{4\pi} c_a^{(1)}\,, \end{align} where numerically \(d_v^{(0)} = 1\). The index \(H\) indicates corrections where a Higgs boson couples exclusively to the heavy quark. \(\alpha_s(\mu)\) denotes the strong coupling constant in the MS scheme at the overall renormalisation scale \(\mu\). Explicit formulas for the coefficients can be found in [17], [9], [26]. ## ElectroweakThe electroweak correction to the cross section can be written (up to the N \begin{equation} \label{eq:sigma_EW} R_{S, \text{EW}}(s) = \frac{12N_c}{s}\alpha(\mu_\alpha)\Big\{c_v\text{Im}\Big[(C^{(v)}_0C_{\text{EW}}^{(v)} + C^{(a)}_0C_{\text{EW}}^{(a)})\,G_{\text{PR}}(E)\Big] + R_{\text{IS}}^{\text{conv}}(s) \Big\}\,. \end{equation} \(R_{\text{IS}}^{\text{conv}}(s)\) denotes the part of the QED initial-state corrections that is not captured by the convolution with the stucture functions \(\Gamma_{ee}^{\text{LL}}\) introduced in the Resummation of QED effects section. For details, see [3]. In contrast to the QCD and Higgs hard matching coefficients discussed in the QCD and Higgs section the electroweak Wilson coefficients \(C_{\text{EW}}^{(v)}\) and \(C_{\text{EW}}^{(a)}\) have a non-vanishing imaginary part that contributes to the cross section. They can be decomposed further into a pure QED contribution and corrections involving at least one \(W\), \(Z\), or Goldstone boson. We also include the Higgs-loop correction to the s-channel \(Z\) propagator in \(C^{(v, a)}_{\text{WZ}}\). \begin{equation} \label{eq:Cva_EW} C^{(v, a)}_{\text{EW}} = C^{(v, a)}_{\text{QED}} + C^{(v, a)}_{\text{WZ}}\,. \end{equation} Note that \(C_{\text{WZ}}^{(v)}\) and \(C_{\text{WZ}}^{(a)}\) do not contain corrections from Higgs bosons coupling exclusively to heavy quarks; these are instead absorbed into \(c_v\). As already mentioned in the Resummation of QED effects section, \(C_{\text{QED}}^{(v)}\) and \(C_{\text{QED}}^{(a)}\) do not include purely photonic corrections that couple only to the initial-state leptons, yet. \(G_{\text{PR}}(E)\) is the pole-resummed Green function (see also the Pole resummation section) given by \begin{equation} \label{eq:PR_GF} G_{\text{PR}}(E) = G(E) + \sum_{N=1}^\infty\bigg\{ \frac{\big[|\psi_N(0)|^2\big]_{\text{expanded}}}{\big[E_N - E - i\Gamma\big]_{\text{unexpanded}}} - \bigg[\frac{|\psi_N(0)|^2}{E_N - E - i\Gamma}\bigg]_{\text{expanded}}\bigg\}\,, \end{equation} where \(\psi_N(0)\) is the quarkonium wave function at the origin. Like in the QCD pole resummation, \(E_N\) denotes the binding energy to the same order as the total cross section, i.e. N According to the power counting outlined in the Power counting section, the electroweak correction first contributes at N Since the mass of the bottom quark lies significantly below the electroweak scale, corrections due to W, Z, and Higgs bosons should be considered in an effective (Fermi) theory, if at all. For this reason we discard all corrections contributing to \(C^{(v, a)}_{\text{WZ}}\) when computing the bottom production cross section. Technically, these corrections are excluded whenever \(m_Q < m_Z\), so for unphysical values of the top mass also the top production cross section would be affected. ## Green functionsThe expansion of the S-wave Green function to third order can be written as \begin{equation} \begin{split} \label{eq:GF_exp} G(E) ={}& \langle \mathbf{0} | \hat{G}_0(E) |\mathbf{0}\rangle + \langle \mathbf{0} | \hat{G}_0(E)\,i\,\delta V\, i\,\hat{G}_0(E) |\mathbf{0}\rangle \\ &+ \langle \mathbf{0} | \hat{G}_0(E)\,i\,\delta V\, i\,\hat{G}_0(E)\,i\,\delta V\, i\,\hat{G}_0(E) |\mathbf{0}\rangle\\ &+ \langle \mathbf{0} | \hat{G}_0(E)\,i\,\delta V\, i\,\hat{G}_0(E)\,i\,\delta V\, i\,\hat{G}_0(E)\,i\,\delta V\, i\,\hat{G}_0(E) |\mathbf{0}\rangle\\ &+\delta^{us}G(E) + \dots\,. \end{split} \end{equation} Here, \(\hat{G}_0(E)\) is the Green function operator of unperturbed PNRQCD and \(\delta V\) denotes a correction to the leading-order potential. \(\delta^{us}G(E)\) stands for the ultrasoft correction contributing only at third order. Again, this equation is to be understood as consistent expansions to N For the P-wave Green function an analogous formula holds. In this case, there is an additional insertion of \(p \cdot p'\), where \(p\) and \(p'\) are the momenta of the initial and final state. Only the first two terms (no perturbation and a single insertion) are needed at N ## PotentialsThe corrections to the potential up to third order can be classified in the following way: \begin{equation} \label{eq:delta_V} \delta V = \underbrace{\delta_C V + \delta_{\text{QED}} V}_{\text{NLO}} + \underbrace{\delta_{1/r^2} V + \delta_{\delta} V + \delta_p V + \delta \text{kin} + \delta \text{width} + \delta \text{width-kin}}_{\text{N}^2\text{LO}} + \underbrace{\delta_H V}_{\text{N}^3\text{LO}}\,, \end{equation} - \(\delta_C V\): Corrections to the colour Coulomb potential.
- \(\delta_{\text{QED}} V\): QED Coulomb potential.
- \(\delta_{1/r^2} V\): Potential proportional to \(1/r^2\) (equivalently \(1/m\)).
- \(\delta_{\delta} V\): Potential proportional to \(\delta(r)\) (equivalently \(1/m^2\)).
- \(\delta_{p} V\): Momentum-dependent potential.
- \(\delta_H V\): Potential due to Higgs exchange.
- \(\delta \text{kin}\): Kinetic energy correction.
- \(\delta \text{width}\): Width correction.
- \(\delta \text{width-kin}\): Mixed width-kinetic energy correction.
The braces indicate the order at which these potential corrections first appear. While all QCD corrections to the potentials are implemented up to N In this vein, we exclude the QED corrections to both the colour Coulomb potential and the delta potential \(\delta_{\delta} V\). For \(\delta_{\text{QED}} V\), the one-loop (order \(\alpha^2\)) QED contribution, corresponding to a N
For the case of a non-zero light quark (e.g. charm) mass, the potentials receive further contributions, which are known to N \begin{equation} \label{eq:delta_ml} \delta_C V = \delta_{C,0} V + \delta_{C,m_l} V\,, \end{equation} where \(\delta_{C,0} V\) corresponds to the contribution for a vanishing light quark mass \(m_l\). Because of the strong mass hierarchy, corrections to top-related observables due to a non-zero light-quark mass are negligible. However, in the case of the bottom quark the charm-quark mass is of the same order as the heavy-quark momentum. According to our power counting (the Power counting section) the charm contributions to the colour Coulomb potential are therefore formally of the same order as the contributions from massless quarks. In practice, charm-mass effects are found to be numerically small but computationally rather expensive, typically requiring numerical Mellin-Barnes integrations. Therefore, we liberally discard sub-leading effects during calculations. For the continuum cross section, at most the single insertion of the potential \(\delta_{C, m_l} V\) at NLO into the S-wave Green function is taken into account. The bound state energies and residues also contain the single insertion of the N ## Energy levels and residuesAs noted in the Pole resummation section, the S-wave energy levels \(E_N\) are given by the position of the poles in the vector polarisation function, or, equivalently, the S-wave Green function \(G(E)\). The residues of the Green function then correspond to the modulus squared of the wave function at the origin: \begin{equation} \label{eq:psi_N_def} G(E) \xrightarrow{E+i\Gamma \to E_N} \frac{|\psi_N(0)|^2}{E_N - E - i \Gamma}\,. \end{equation} Since \(\psi_N(0)\) is a factorisation scheme dependent quantity, the \begin{equation} \label{eq:ZN_def} Z_N = \frac{4m_Q^2}{s_N}c_v\bigg[c_v - \frac{E_N}{m_Q}\frac{d_v}{3}\bigg]|\psi_N(0)|^2 \end{equation} with \(s_N = (2m_Q + E_N)^2\). ## Mass schemesSo far, all formulas have been expressed in terms of the pole mass of the heavy quark. Mass values in other schemes RS can be converted via relations of the form \begin{equation} \label{eq:conv_pole} m_Q = m_Q^{\text{RS}} + \sum_{i=0}^{o} \delta m_i^{\text{RS}}\,, \end{equation} where \(0 \leq o \leq 3\) is the considered order according to the PNRQCD power counting summarised in the Power counting section. The pole mass can now be substituted in two different ways [2] : - First compute the numerical value of the pole mass from the mass value in the scheme RS by using the above relation. The value obtained for the pole mass in this way will strongly depend on the order \(o\). Then evaluate the expressions for the cross section, residues, and energy levels in the
*pole*scheme with the pole mass \(m_Q\) as determined above. This the*shift*prescription. Symbolically replace the pole mass in all expressions by the mass in the scheme RS via the above relation. Then perform a systematic expansion wherever \(\delta m_i^{\text{RS}}\) constitutes a small correction. To ensure both a consistent expansion and order-by-order renormalon cancellation, \(\delta m_0^{\text{RS}}\) has to be of order \(E\sim v^2\). Finally, insert the numeric value of the mass in the scheme RS. This corresponds to the *insertion*prescription.Deviating from this general rule, in the present version of the code the residue \(Z_N\) in the insertion scheme is computed as \((m_Q^{\text RS})^2 \times [Z_N/m_Q^2]_{\text{RS}}\) where the quantity in square brackets is transformed from the pole scheme to RS according to the general rule. That is, we perform the naive replacement \(m_Q \to m_Q^{\text{RS}}\) *without*the correction terms \(\delta m_i^{\text{RS}}\) in the factors \(4m_Q^2/s_N\) in the definition of the residue and also in \(N_C/(2m_Q^2)\) in the pole resummation This has no effect on the cross section, since \(m_Q\) cancels in the product of these factors. However, \(Z_N\) as defined above differs slightly from the value \([Z_N]_{\text{RS}}\) it would attain, if the general rule were applied to \(Z_N\) directly.The insertion prescription leads to unphysical oscillations of the cross section near threshold [2]. What is more, for the bottom quark the cross section is only defined in the sense of a distribution. This is due to the expansion of the threshold step function \(\theta(s - 4m_b^2)\) in \(\delta m_i^{\text{RS}} < \sqrt{s}\) for \(i > 0\). It is therefore recommended to use the shift prescription instead.
In both prescriptions, the energy variable \(E^\text{RS}\) is defined as \(E^\text{RS} = \sqrt{s} - 2m_Q^\text{RS}\), and similarly the binding energies are defined by the bound state masses minus \(2m_Q^\text{RS}\) so that \(s_N = (2m_Q^\text{RS} + E_N^\text{RS})^2\) is scheme independent (see also the Resonances section). So far, apart from the pole scheme, the following schemes are implemented in QQbar_threshold: - The potential-subtracted (PS) scheme [13] up to N
^{3}LO [7]. Corrections from a non-zero light-quark mass are only contained up to N^{2}LO [11]. We define the subtraction potential to not include any electroweak corrections. Because of this, the first-order QED corrections leads to a visible shift of the \(t\bar{t}\) cross section peak for fixed input PS mass, but contrary to QCD, higher-order QED and electroweak corrections are rapidly convergent. - The 1S scheme [23]. Up to N
^{3}LO the conversion formula to the pole scheme is given by\begin{equation} \label{eq:1S_to_pole} m_Q = m_Q^{\text{1S}} - \frac{E_1(m_Q)}{2} = m_Q^{\text{1S}} - \frac{E_1(m_Q^{\text{1S}})}{2} + \frac{E_1(m_Q^{\text{1S}})}{4} \frac{\partial E_1(m_Q^{\text{1S}})}{\partial m_Q^{\text{1S}}} + \dots\,. \end{equation} Since the 1S scheme is closely connected to the bound state energy levels, corrections are implemented to the same order, i.e. N^{3}LO for the QCD and Higgs corrections and N^{2}LO for the electroweak corrections. - The MS scheme in QCD. For this scheme, we keep \(\delta m_0^{\overline{\text{MS}}}\) at LO in (28), \(\delta m_1^{\overline{\text{MS}}}\) at NLO, and so on. Corrections are available at order \(\alpha_s^4\) [27], which corresponds to N
^{3}LO as required at the present highest accuracy. Since \(\delta m_0^{\overline{\text{MS}}}\) is of the same order as \(v\) (rather than \(v^2\) as in the PS and 1S scheme), only the shift prescription is self-consistent with the convention adopted here. We also include corrections from a non-zero light-quark mass to N^{2}LO [1]. We define this scheme via the pure QCD relation to the pole mass and therefore do not include any electroweak corrections to the mass conversion.
## OptionsOptions are set by passing an
In the C++ case, it is recommended to modify an object initialised with the helper functions Note that in many cases there is more than one option that disables certain parts. In case of conflicting settings, a contribution is discarded. For example, if all QED contributions are switched off through the The `contributions` : Specifies multiplicative factors for the potentials and current matching coefficients. For example, settingimplies that corrections due to the leading-order delta potential are discarded, but corrections from the next-to-leading delta potential are kept, i.e. multiplied by 1. The following table lists the relations to the definitions in the Potentials section.options opt;opt.contributions.v_delta = {{0., 1.}};C++ name Mathematica name Corrections `v_Coulomb` `vCoulomb` \(\{\delta_CV^{(1)}, \delta_CV^{(2)}, \delta_CV^{(3)}\}\) `v_delta` `vdelta` \(\{\delta_\delta V^{(0)}, \delta_\delta V^{(1)}\}\) `v_r2inv` `vr2inv` \(\{\delta_{1/r^2} V^{(1)}, \delta_{1/r^2} V^{(2)}\}\) `v_p2` `vp2` \(\{\delta_{p} V^{(0)}, \delta_{p} V^{(1)}\}\) `v_kinetic` `vkinetic` \(\{\delta\text{kin}\}\) `v_width2` `vwidth2` \(\{\delta\text{width}\}\) `v_width_kinetic` `vwidthkinetic` \(\{\delta\text{width-kin}\}\) `width_ep` `widthep` \(\{\delta_{\Gamma_\epsilon/\epsilon} R\}\) `ultrasoft` `ultrasoft` \(\{\delta^{us}G(E)\}\) `v_Higgs` `vHiggs` \(\{\delta_H V^{(0)}\}\) `v_QED_Coulomb` `vQEDCoulomb` \(\{\delta_{\text{QED}}V^{(0)}\}\) `cv` `cv` \(\{c_v^{(1)}, c_v^{(2)}, c_v^{(3)}\}\) `cv_Higgs` `cvHiggs` \(\{c_{vH}^{(2)}, c_{vH}^{(3)}\}\) `Cv_QED` `CvQED` \(\{C^{(v)}_{\text{QED}}\}\) `Ca_QED` `CaQED` \(\{C^{(a)}_{\text{QED}}\}\) `Cv_WZ` `CvWZ` \(\{C^{(v)}_{\text{WZ}}\}\) `Ca_WZ` `CaWZ` \(\{C^{(a)}_{\text{WZ}}\}\) `dv` `dv` \(\{d_v^{(0)},d_v^{(1)}\}\) `ca` `ca` \(\{c_a^{(1)}\}\) `nonresonant` `nonresonant` \(\{R_{\text{non-res}}^{\text{NLO}}, R_{\text{non-res}}^{\text{N}^2\text{LO}}\}\) `Contributions` option expects a list of all contributions with their coefficients:To facilitate the usage, theContributions -> {vCoulomb -> {1., 1., 1.},vdelta -> {0., 1.},...}`QQbarThreshold` package provides the auxiliary functions`ExceptContributions` and`OnlyContributions` which set the factors for all contributions that are not listed explicitly to 1 or 0, respectively. For example, to discard only the leading-order delta potential one could useContributions -> ExceptContributions[vdelta -> {0., 1.}]`alpha_s_mZ` or`alpha_s_mu0` specifies the input value for the strong coupling constant. If the option`alpha_s_mZ` is used, it is assumed that the given value corresponds to \(\alpha_s(m_Z)\).`alpha_s_mu0` specifies both a reference scale and the value of \(\alpha_s\) at that scale. For exampleoptions opt;opt.alpha_s_mu0 = {10., 0.22};sets \(\alpha_s(10\,\text{GeV}) = 0.22\). If both options are set, the value for `alpha_s_mZ` is ignored.The input value for the strong coupling is evolved automatically to the overall renormalisation scale using four-loop evolution. For bottom-related functions decoupling to the four-flavour theory is performed only if the input scale is above the decoupling scale `mu_thr` defined in the`constants.hpp` header. With the current default settings, decoupling is performed at twice the scale-invariant mass \(m_b^{\overline{\text{MS}}}(m_b^{\overline{\text{MS}}}) = 4.203\,\)GeV. Note that for top-related functions the input value for the strong coupling is always assumed to refer to the five-flavour theory and no decoupling is performed.The final values used for the actual calculations can be inspected with the `alpha_s_bottom` and`alpha_s_top` functions from the header alpha_s.hpp (or`alphaSBottom` ,`alphaSTop` in Mathematica), which take the renormalisation scale as their first argument and the value of either`alpha_s_mZ` or`alpha_s_mu0` as their second argument.`m_Higgs` : Specifies the value of the Higgs boson mass.`Yukawa_factor` : Specifies a multiplier for the top-quark Yukawa coupling. This can be used to parametrise a possible deviation from the Standard Model relation between the top-quark mass and the coupling to the Higgs boson.We assume that this deviation is caused by the dimension-6 operator \begin{equation} \label{eq:yukawa_mod_op} \Delta {\cal L} = - \frac{c_{\text{NP}}}{\Lambda^2}(\phi^\dagger\phi)(\bar{Q}_3i\sigma^2\phi^* t_R) + \text{h.c.}\,, \end{equation} which implies the relation [12] \begin{equation} \texttt{Yukawa_factor} = 1+\frac{c_\text{NP}}{\Lambda^2}\frac{v^3}{\sqrt{2} m_t}\,. \end{equation} While this operator modifies the coupling to the physical Higgs boson, the couplings to the Goldstone bosons remain unchanged, provided they are expressed in terms of the top-quark mass. The operator also generates four- and five-point vertices. Since we count the coupling \(c_\text{NP} v^2/\Lambda^2\) as N \(^2\)LO, similar to \(\alpha\) and \(y_t^2\) (c.f. the Power counting section), these vertices contribute only through a Higgs tadpole diagram to the top self-energy, which has no effect in the top mass renormalisation schemes adopted here. Hence in the present approximation, the only effect of the dimension-6 operator is a rescaling of the Yukawa coupling. `resonant_only` : If set to`true` , the nonresonant contribution to the cross section is discarded.`invariant_mass_cut` : Specifies an invariant mass cut for the nonresonant contribution to the cross section. The invariant mass of each \(W\,b\) pair in the final state is restricted to the region between \(m_t - \texttt{invariant_mass_cut}\) and \(m_t + \texttt{invariant_mass_cut}\). By default, the loosest possible cut \(\texttt{invariant_mass_cut} = m_t - m_W\) is taken, which yields the total cross section.`ml` : Specifies the value of the light (e.g. charm or bottom) quark mass. The input value should be the MS quark mass at the overall renormalisation scale \(\mu\). This option only affects the mass of light quarks in virtual corrections, the bottom quarks in the final state of the process \(e^+ e^- \to W^+ W^- b\bar{b}\) are always assumed to be massless.`r4` : Specifies the four-loop coefficient for the conversion between the pole and the MS scheme. More precisely,`r4` is defined by the relation\begin{equation} \label{eq:r4_def} m_Q = m_Q^{\overline{\text{MS}}}\big(m_Q^{\overline{\text{MS}}})\big)\big[1 + r_1a_s + r_2a_s^2 + r_3a_s^3 + `r4`\,a_s^4 + {\cal O}\big(a_s^5\big)\big]\,, \end{equation} where \(a_s = \alpha_s^{(n_l)}\big(m_Q^{\overline{\text{MS}}}\big)/\pi\). Note that only the \(n_l\) massless quark flavours contribute to the running of the strong coupling; in the notation of [28] `r4` is given by`r4` \(=c_m^{(4)} + \delta c_m^{(4)}\) with \(\mu = m_Q^{\overline{\text{MS}}}\big(m_Q^{\overline{\text{MS}}}\big)\).This option is only relevant if the MS scheme was chosen and only affects the conversion at N ^{3}LO.`alpha` : Specifies the value of the QED coupling constant at the scale`mu_alpha` . This mainly affects the overall normalisation factor of the cross section but also all electroweak corrections.`mu_alpha` : Specifies the scale for the QED coupling constant.
Let us comment on the usage of the two previous options by adopting two examples. (I) We want to compute the top cross section using a different value for \(\alpha(m_Z)\), say \(\alpha(m_Z) = 1/130\). (II) We think that the default scale \(\mu_\alpha = m_Z\) is too low and want to use \(\alpha(m_t)\) as input. In scenario (I) we set `resum_poles` : Specifies the number of bound states that are resummed into the vector polarisation function and the electroweak contribution to the cross section. At N^{3}LO the current maximum value is 6; at lower orders there is no such upper limit. This option does not affect the pole resummation for the axialvector polarisation function, where in the current version always the three leading poles are resummed.`beyond_QCD` : Specifies the Standard Model corrections beyond QCD that should be taken into account. More precisely, each setting defines the Lagrangian of the underlying full theory, from which the higher-order corrections in the effective nonrelativistic theory are then derived. For example, with the setting`beyond_QCD = SM::Higgs` corrections involving Higgs bosons coupling to the heavy quarks are added to the usual PNRQCD corrections. Note that this option does not affect the leading-order production process, i.e. s-channel production via a virtual \(Z\) or photon is still taken into account with the above setting, although higher-order corrections due to photons or \(Z\) bosons are disabled. Nonresonant production is also unaffected by this option.Here is a chart with the possible settings: \(\mathcal{L}_{\text{QCD}}\) and \(\mathcal{L}_{\text{SM}}\) denote the usual QCD and Standard Model Lagrangians. Furthermore, we use \begin{align} \label{eq:lagrangians} \mathcal{L}_{\text{QED}} ={}& -\frac{1}{4} F_{\mu\nu}F^{\mu\nu} + \sum_{l \in \text{leptons}} \overline{\psi}_l i \partial_\mu \gamma^\mu \psi_l - \sum_{f \in \text{fermions}} e_f \overline{\psi}_f A_\mu \gamma^\mu \psi_f \,,\\ \mathcal{L}_{\text{Higgs}} ={}& \frac{1}{2} (\partial_\mu H)^2 - \frac{1}{2}m_H^2H^2 - \sqrt{\frac{\lambda}{2}}m_HH^3 - \frac{\lambda}{4}H^4 - \frac{y_t}{\sqrt{2}} \overline{t} t H\,, \end{align} where \(\lambda = \frac{\pi\alpha m_H^2}{2m_W^2s_w^2}\). `mass_scheme` : Specifies the renormalisation scheme. Depending on the scheme it is also mandatory to specify one or more associated scales.Here is a list of the available schemes:
The pole and 1S schemes do not have any associated scales and can be set with options opt; opt.mass_scheme = {Pole}; and similar code for MassScheme -> "Pole" The options opt; opt.mass_scheme = {MSshift, 160.}; in C++ or MassScheme -> {"MSshift", 160.} in Mathematica. Finally, the potential-subtracted schemes have two associated scales \(\mu_f\) and \(\nu\), where the first scale \(\mu_f\) specifies the factorisation scale defining the cut-off in the subtraction term and the second scale \(\nu\) is the scale appearing in logarithms that originate from the infrared singularity in the static potential. If only a single scale is given, \(\nu\) is set to \(\mu_f\): options opt; opt.mass_scheme = {PSshift, 20., 20.}; opt.mass_scheme = {PSshift, 20.}; // equivalent to the line before `production` : Specifies which production channels are taken into account. The possible settings are:`photon_only` : Production only via a virtual photon. This effectively discards the P-wave contribution, the second term in the vector production operator, and the axialvector production operator In addition, box corrections and corrections to the production via a virtual Z boson are discarded in the electroweak contribution to the cross section.`S_wave_only` : Production only via an S-wave photon or Z, i.e. the P-wave contribution is discarded.`all` : All possible production channels.
The corresponding Mathematica settings are `"PhotonOnly"` ,`"SWaveOnly"` , and`"All"` (or`All` ).`expand_s` : Specifies the treatment of the overall factor \(1/s\) in the polarisation functions defined in the Production channels section and the electroweak correction. If set to`true` , \(s = (2m_Q + E)^2\) is expanded in \(E/m_Q \sim v^2 \ll 1\) to the appropriate order. This also affects the prefactor \(1/s_N\) in the residue \(Z_N\).`double_light_insertion` : Specifies whether double insertions of the light-quark potential correction \(\delta_{C, m_l}V\) are taken into account. This option only affects the calculation of the energy levels and residues; in the continuum cross section double insertions are always neglected.The impact on observables that can be computed reliably within perturbation theory is typically small. For example, in the determination of the bottom quark mass from the 10th moment of the cross section the end result is changed by about \(0.1\) per mille, compared to an overall change of around \(0.5\) per mille when combined with the dominant single insertions [11]. In infrared-sensitive quantities like binding energies and residues of bound states with principal quantum number \(N > 5\) the effects can become significant. The calculation of the double insertions is computationally very expensive, requiring the evaluation of an infinite sum over integrals (cf. [11]). Therefore, setting this option to `true` leads to a considerable slowdown by a factor between \(100\) and \(1000\). In order to avoid an even more severe slowdown as well as numerical instabilities, the current implementation only computes the first few terms of the sum, so that the result is not very precise. Furthermore, the implementation is not thread safe (cf. the Parallelisation section).`ISR_const` : Specifies whether the non-logarithmic initial-state correction \(R_{\text{IS}}^{\text{conv}}\) should be included. This should be set to`true` if and only if also the logarithmically enhanced contribution are included by computing the convolution of the cross section with the luminosity function.
Finally, here are the default values for the option settings:
## Advanced usageIn the following we discuss several more complicated examples that require some knowledge of the options discussed in the Options section and the structure of the cross section outlined in the Structure of the cross section section. ## Initial-state radiationInitial-state radiation requires the computationally expensive convolution with structure functions. Therefore, this correction is not included automatically. This example shows how initial-state radiation can be added manually. After defining the luminosity function \begin{equation} \label{eq:lumi_def} {\cal L}(x) = \int_y^1 \frac{dy}{y} \, \Gamma_{ee}^\text{LL}(y)\Gamma_{ee}^\text{LL}(x/y) \end{equation} with the electron structure functions \(\Gamma_{ee}^\text{LL}\) the cross section after initial-state radiation is given by \begin{equation} \label{eq:ISR_lumi} \sigma_Q^{\rm ISR}(s) = \int_0^1dx \,{\cal L}(x)\,\sigma_Q(xs)\,. \end{equation} Here, \(\sigma_Q\) is the partonic cross section including the non-logarithmic initial-state radiation correction. The non-logarithmic correction can be included with the option setting options opt; opt.ISR_const = true; in C++ and In principle, the convolution integral covers the whole energy range from zero to the nominal center-of-mass energy. However, our prediction for the cross section is only valid in the vicinity of the threshold. Sufficiently below the threshold the actual cross section becomes negligible. This implies that we can introduce a lower cut-off \(x_\text{min}\) in the integral. In the following we choose \(x_\text{min} = (330\,\text{GeV})^2/s\). A further, purely numerical problem arises from the integrable divergence of the luminosity function for \(x \to 1\). In order to eliminate this divergence, we can change the integration variable to \(t = (1-x)^\beta\) and write the cross section as \begin{equation} \label{eq:ISR_lumi_mod} \sigma_{\rm ISR}(s) = \int_0^{t_\text{max}}\!\!dt \,\bar{{\cal L}}(t)\,\sigma^{\text{conv}}\big(x(t)\, s\big)\,,\qquad \bar{{\cal L}}(t) = \frac{(1-x)^{1-\beta}{\cal L}(x)}{\beta}\,, \qquad x(t) = 1 - t^\frac{1}{\beta}\,, \end{equation} with the modified luminosity function \(\bar{{\cal L}}(t)\) and a cut-off \(t_\text{max} = (1 - x_\text{min})^\beta\). The function \(\beta = -2\alpha(\mu_\alpha)/\pi[\log(m_e^2/s) + 1]\) is available in QQbar_threshold as Finally, version 2 of QQbar_threshold provides an The following C++ code prints the cross section \(\sigma = 0.591737\,\text{pb}\) after initial state radiation for a center-of-mass energy of \(\sqrt{s} = 344\,\text{GeV}\), including all known perturbative corrections: #include <iostream> #include <cmath> #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" #include "QQbar_threshold/structure_function.hpp" #include "QQbar_threshold/integrate.hpp" #include "QQbar_threshold/constants.hpp" int main(){ namespace QQt = QQbar_threshold; constexpr double sqrt_s = 344.; constexpr double mu = 80.; constexpr double mu_width = 350.; constexpr double mt_PS = 171.5; constexpr double width = 1.33; QQt::options opt = QQt::top_options(); opt.ISR_const = true; const auto integrand = [=](double t){ const double x = 1 - std::pow(t, 1/beta); std::sqrt(x)*sqrt_s, {mu, mu_width}, {mt_PS, width}, QQt::N3LO, opt ); return L*sigma; }; constexpr double x_min = 330.*330./(sqrt_s*sqrt_s); const double t_max = std::pow(1 - x_min, beta); std::cout << QQt::integrate(integrand, 0, t_max) << '\n'; } The corresponding Mathematica code is Needs["QQbarThreshold`"]; LoadGrid[GridDirectory <> "ttbar_grid.tsv"]; sqrts = 344.; mu = 80.; muWidth = 350.; mtPS = 171.5; width = 1.33; order = "N3LO"; beta = ISRLog[sqrts, alphamZ]; xmin = (330./sqrts)^2; tmax = (1-xmin)^beta; Print[ NIntegrate[ ModifiedLuminosityFunction[t, beta]* TTbarXSection[ Sqrt[1-t^(1/beta)] * sqrts, {mu, muWidth}, {mtPS, width}, order, ISRConst -> True ], {t, 0, tmax} ] ]; Numerically, the structure function under the replacement \(\beta \to 2\beta\) is very close to the luminosity function. The same holds for the modified versions of both functions. Indeed, substituting ## Wave function at the originWhile the quarkonium wave function at the origin is not a physical observable, it serves as a good example to illustrate some of the more advanced options. Our starting point is the definition of the residue \(Z_N\). The following assumes that we wish to determine the wave function at the origin in the PS-shift scheme and that the #include <iostream> #include "QQbar_threshold/scheme_conversion.hpp" #include "QQbar_threshold/energy_levels.hpp" #include "QQbar_threshold/residues.hpp" int main(){ namespace QQt = QQbar_threshold; const double mb_PS = 4.5; const double mu = 4.5; QQt::options opt = QQt::bottom_options(); opt.contributions.cv = {0., 0., 0.}; opt.contributions.dv = {0., 0.}; 1, mu, mb_PS, QQt::N3LO, opt ); const double s_1 = (2*mb_PS + E_1_PS)*(2*mb_PS + E_1_PS); mu, mb_PS, QQt::N3LO, opt ); std::cout << s_1/(4*mb_pole*mb_pole)*QQt::bbbar_residue( 1, mu, mb_PS, QQt::N3LO, opt ) << '\n'; } Needs["QQbarThreshold`"]; With[ { mbPS = 4.5, mu = 4.5, opt = Contributions -> ExceptContributions[ cv -> {0, 0, 0}, dv -> {0, 0} ] }, E1PS = BBbarEnergyLevel[1, mu, mbPS, "N3LO", opt]; s1 = (2*mbPS + E1PS)^2; mbPole = BottomPoleMass[mu, mbPS, "N3LO", opt]; Print[s1/(4*mbPole^2)*BBbarResidue[1, mu, mbPS, "N3LO", opt]]; ]; Since \(s_N=(2m_Q + E_N)^2 = (2m_Q^{\text{PS}} + E_N^{\text{PS}})^2\) is a scheme-independent quantity we could also have computed the binding energy in the pole scheme and combined it with the pole mass. This method works analogously in the other shift schemes. In an insertion scheme, we proceed the same way. However, in this case one needs to multiply by \(s_N/(4(m_Q^{\text{RS}})^2)\). ## ParallelisationWhile the observables for a single given set of parameters can be calculated rather quickly, computing e.g. the cross section over a range of centre-of-mass energies and parameter settings can become somewhat time consuming. In such a situation parallelisation can lead to significant speed-ups. As an example, we show how to perform a threshold scan similar to the one discussed in the Cross sections section, but including scale variation. Our strategy is to generate a list containing the centre-of-mass energy, the cross section obtained with a default scale of \(80\,\)GeV, the minimum and the maximum cross section obtained through scale variation for each point in parallel. Since the mechanisms typically used in C++ vs. Mathematica programs are rather different, we discuss these languages separately in the following sections. Since we use threads for parallelisation it may be necessary to add additional flags for compilation. The g++ -o parallel -std=c++11 parallel.cpp -pthread -lQQbar_threshold As in the previous examples we use an abbreviation for the somewhat unwieldy QQbar_threshold namespace: namespace QQt = QQbar_threshold; First, we set up a struct comprising the minimum, maximum, and default cross sections obtained for a given centre-of-mass energy: struct xs_point{ double sqrt_s; double xs_default; double xs_min; double xs_max; }; Since we call the cross section function many times with mostly the same arguments, it is also useful to define an auxiliary function that has only the energy and the scale as remaining arguments: double xsection(double sqrt_s, double mu){ static constexpr double mu_width = 350.; static constexpr double mt = 171.5; static constexpr double width = 1.33; return QQt::ttbar_xsection( sqrt_s, {mu, mu_width}, {mt, width}, QQt::N3LO ); } Using this function, we can then define the scale variation for a single centre-of-mass energy. For the sake of simplicity, we use a rather crude sampling of the cross section to estimate the extrema. xs_point xsection_scale_variation(double sqrt_s){ static constexpr double mu_default = 80.; static constexpr double mu_min = 50.; static constexpr double mu_max = 350.; static constexpr double mu_step = 5.; xs_point result; result.sqrt_s = sqrt_s; result.xs_default = xsection(sqrt_s, mu_default); result.xs_min = result.xs_default; result.xs_max = result.xs_default; for(double mu = mu_min; mu < mu_max; mu += mu_step){ double current_xsection = xsection(sqrt_s, mu); if(current_xsection < result.xs_min){ result.xs_min = current_xsection; } else if(current_xsection > result.xs_max){ result.xs_max = current_xsection; } } return result; } The code to actually calculate the scale variation for various energies in parallel is then rather short: std::vector<std::future<xs_point>> results; for(double sqrt_s = 340.; sqrt_s < 349.; sqrt_s += 0.2){ results.emplace_back( std::async( std::launch::async, xsection_scale_variation, sqrt_s ) ); } For each energy, a new thread is launched for computing the scale variation. While this can be rather inefficient in practice, it still ilustrates how parallelisation can be achieved in principle. To complete the example, we should also include the appropriate headers, load a grid, and produce some output. While all these steps are straightforward, some care should be taken when loading the grid. As already stated in the Cross sections section, the Finally, here is the full code for the parallelised threshold scan with scale variation: #include <iostream> #include <iomanip> #include <vector> #include <thread> #include <future> #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" namespace QQt = QQbar_threshold; struct xs_point{ double sqrt_s; double xs_default; double xs_min; double xs_max; }; double xsection(double sqrt_s, double mu){ static constexpr double mu_width = 350.; static constexpr double mt = 171.5; static constexpr double width = 1.33; return QQt::ttbar_xsection( sqrt_s, {mu, mu_width}, {mt, width}, QQt::N3LO ); } xs_point xsection_scale_variation(double sqrt_s){ static constexpr double mu_default = 80.; static constexpr double mu_min = 50.; static constexpr double mu_max = 350.; static constexpr double mu_step = 5.; xs_point result; result.sqrt_s = sqrt_s; result.xs_default = xsection(sqrt_s, mu_default); result.xs_min = result.xs_default; result.xs_max = result.xs_default; for(double mu = mu_min; mu < mu_max; mu += mu_step){ double current_xsection = xsection(sqrt_s, mu); if(current_xsection < result.xs_min){ result.xs_min = current_xsection; } else if(current_xsection > result.xs_max){ result.xs_max = current_xsection; } } return result; } int main(){ std::vector<std::future<xs_point>> results; for(double sqrt_s = 340.; sqrt_s < 349.; sqrt_s += 0.2){ results.emplace_back( std::async( std::launch::async, xsection_scale_variation, sqrt_s ) ); } std::cout << std::fixed; std::cout << "sqrt_s \tcentral \tmin \tmax\n"; for(auto & res: results){ xs_point current_point = res.get(); std::cout << current_point.sqrt_s << '\t' << current_point.xs_default << '\t' << current_point.xs_min << '\t' << current_point.xs_max << '\n'; } } ## MathematicaFor parallelisation in Mathematica, we have to ensure that each kernel knows all relevant definitions and has its own precomputed grid: Needs["QQbarThreshold`"]; LaunchKernels[]; ParallelEvaluate[Needs["QQbarThreshold`"]]; ParallelEvaluate[LoadGrid[GridDirectory <> "ttbar_grid.tsv"]]; This is very different from the parallelisation at thread level we used in the C++ case. In the present example, the kernels are independent processes and therefore unable to share a common grid. It is not only safe, but even necessary to load multiple copies of a grid at the same time. For the actual threshold scan, we first define an auxiliary function for the cross section with fixed top quark properties and width scale: XSection[sqrts_, mu_] := With[ {muWidth = 350, mt = 171.5, width = 1.33}, TTbarXSection[sqrts, {mu, muWidth}, {mt, width}, "N3LO"] ]; To perform the scale variation, we use Mathematica's built-in XSectionScaleVariation[sqrts_] := Module[ { muDefault = 80, muMin = 50, muMax = 350, xsDefault, xsMin, xsMax, mu }, xsDefault = XSection[sqrts, muDefault]; xsMin = NMinValue[ {XSection[sqrts, mu], muMin <= mu <= muMax}, mu, Method -> "SimulatedAnnealing" ]; xsMax = NMaxValue[ {XSection[sqrts, mu], muMin <= mu <= muMax}, mu, Method -> "SimulatedAnnealing" ]; Return[{sqrts, xsDefault, xsMin, xsMax}]; ]; After this, the code for the actual parallelised scan is again rather compact: results = ParallelTable[ XSectionScaleVariation[sqrts], {sqrts, 340, 349, 0.2} ]; We complete the example by adding some code for the output: LaunchKernels[]; Needs["QQbarThreshold`"]; ParallelEvaluate[Needs["QQbarThreshold`"]]; ParallelEvaluate[LoadGrid[GridDirectory <> "ttbar_grid.tsv"]]; XSection[sqrts_, mu_] := With[ {muWidth = 350, mt = 171.5, width = 1.33}, TTbarXSection[sqrts, {mu, muWidth}, {mt, width}, "N3LO"] ]; XSectionScaleVariation[sqrts_] := Module[ { muDefault = 80, muMin = 50, muMax = 350, xsDefault, xsMin, xsMax, mu }, xsDefault = XSection[sqrts, muDefault]; xsMin = NMinValue[ {XSection[sqrts, mu], muMin <= mu <= muMax}, mu, Method -> "SimulatedAnnealing" ]; xsMax = NMaxValue[ {XSection[sqrts, mu], muMin <= mu <= muMax}, mu, Method -> "SimulatedAnnealing" ]; Return[{sqrts, xsDefault, xsMin, xsMax}]; ]; results = ParallelTable[ XSectionScaleVariation[sqrts], {sqrts, 340, 349, 0.2} ]; PrependTo[results, {"sqrt_s", "central", "min", "max"}]; Print[TableForm[results]]; ## Moments for nonrelativistic sum rulesNext to threshold scans for \(t\bar{t}\) production, the calculation of moments for \(\Upsilon\) sum rules is one of the key applications for QQbar_threshold. It is conventional to consider the normalised cross section \(R_b = \sigma_b/\sigma_{pt}\) with \(\sigma_{pt} = 4\pi\alpha(\mu_\alpha)^2/(3s)\), which can be calculated with the \begin{equation} \label{eq:M_def} {\cal M}_n = \int_0^\infty \frac{R_b(s)}{s^{n+1}}\,. \end{equation} Splitting the moments into the contribution from the narrow \(\Upsilon\) resonances and the remaining continuum contribution we obtain \begin{equation} \label{eq:M_res_cont} {\cal M}_n = \frac{12 \pi^2 N_c e_b^2}{m_b^2} \sum_{N=1}^\infty \frac{Z_N}{s_N^{2n+1}} + \int_{4m_b^2}^\infty ds\,\frac{R_b(s)}{s^{n+1}}\,. \end{equation} We now show how this formula can be evaluated with QQbar_threshold. Since discussing numerical integration is clearly outside the scope of this work, we assume the existence of a C++ header integral.hpp that provides a suitable g++ -o moments -std=c++11 moments.cpp -lQQbar_threshold -lgsl -lgslcblas For the sake of simplicity, we use the standard bottom options and only keep the bottom quark mass and \(n\) as free parameters in our example. We can then define auxiliary functions for the energy levels, residues, and the continuum cross section: static constexpr double pi = 3.14159265358979323846264338328; double Z(int N, double mb_PS){ return QQt::bbbar_residue(N, mb_PS, mb_PS, QQt::N3LO); } double E(int N, double mb_PS){ return QQt::bbbar_energy_level(N, mb_PS, mb_PS, QQt::N3LO); } double Rb(double s, double mb_PS){ try{ return QQt::bbbar_R_ratio(std::sqrt(s), mb_PS, mb_PS, QQt::N3LO); } catch(std::out_of_range){ return std::numeric_limits<double>::quiet_NaN(); } } Note that for the bottom cross section we have to handle the case that the centre-of-mass energy is outside the region covered by the precomputed grid. In fact, we assume that the value obtained for the continuum integral is reliable in spite of the limited grid size. For a more careful analysis, this assumption should of course be checked, for example by varying the upper bound of the integral. Let us first consider the resonance contribution. For the prefactor, we have to convert the input mass from the PS scheme to the pole scheme, which can be done with the function double M_resonances(int n, double mb_PS){ static constexpr int N_c = 3; double sum_N = 0.0; for(int N = 1; N <= 6; ++ N){ sum_N += Z(N, mb_PS)*std::pow(2*mb_PS + E(N, mb_PS), -2*n-1); } return 12*pi*pi*N_c*e_b*e_b/(mb_pole*mb_pole)*sum_N; } Since at N For the continuum contribution, we encounter the problem that typical C++ integration routines can only evaluate integrals over a finite interval. To deal with this we perform a substitution, e.g. \(s = s(x) = 4m_b^2 + \frac{x}{1-x}\), so we have \begin{equation} \label{eq:M_cont_subst} \int_{4m_b^2}^\infty ds\, \frac{R_b(s)}{s^{n+1}} = \int_0^1 \frac{dx}{(1 - x)^2} \frac{R_b\big(s(x)\big)}{s(x)^{n+1}}\,. \end{equation} The continuum moments can then be computed with the following code: double M_continuum(int n, double mb_PS){ double const s0 = pow(QQt::bbbar_threshold(mb_PS, mb_PS, QQt::N3LO), 2); auto s = [=](double x){ return s0 + x/(1-x); }; auto integrand = [=](double x){ return Rb(s(x), mb_PS)*std::pow(s(x), -n-1)*std::pow(1 - x, -2); }; return integral(0, 1, integrand); } All that remains is then to add up both contributions and add the standard boilerplate code for including headers, loading the grid, etc. It is also convenient to rescale the moments to be of order one. For this, we multiply them by a factor of \((10\,\text{GeV})^{2n}\). Finally, here is the complete code for our example: #include <cmath> #include <limits> #include <iostream> #include "QQbar_threshold/load_grid.hpp" #include "QQbar_threshold/xsection.hpp" #include "QQbar_threshold/threshold.hpp" #include "QQbar_threshold/residues.hpp" #include "QQbar_threshold/energy_levels.hpp" #include "QQbar_threshold/scheme_conversion.hpp" #include "QQbar_threshold/integrate.hpp" namespace QQt = QQbar_threshold; static constexpr double pi = 3.14159265358979323846264338328; double Z(int N, double mb_PS){ return QQt::bbbar_residue(N, mb_PS, mb_PS, QQt::N3LO); } double E(int N, double mb_PS){ return QQt::bbbar_energy_level(N, mb_PS, mb_PS, QQt::N3LO); } double Rb(double s, double mb_PS){ try{ return QQt::bbbar_R_ratio(std::sqrt(s), mb_PS, mb_PS, QQt::N3LO); } catch(std::out_of_range){ return std::numeric_limits<double>::quiet_NaN(); } } double M_resonances(int n, double mb_PS){ static constexpr int N_c = 3; double sum_N = 0.0; for(int N = 1; N <= 6; ++ N){ sum_N += Z(N, mb_PS)*std::pow(2*mb_PS + E(N, mb_PS), -2*n-1); } return 12*pi*pi*N_c*e_b*e_b/(mb_pole*mb_pole)*sum_N; } double M_continuum(int n, double mb_PS){ double const s0 = pow(QQt::bbbar_threshold(mb_PS, mb_PS, QQt::N3LO), 2); auto s = [=](double x){ return s0 + x/(1-x); }; auto integrand = [=](double x){ return Rb(s(x), mb_PS)*std::pow(s(x), -n-1)*std::pow(1 - x, -2); }; return QQt::integrate(integrand, 0, 1); } double M(int n, double mb_PS){ return pow(10, 2*n)*(M_resonances(n, mb_PS) + M_continuum(n, mb_PS)); } int main(){ std::cout << M(10, 4.5322) << '\n'; } In units of \((10\,\text{GeV})^{-20}\), we find \({\cal M}_{10} = 0.264607\), which is in good agreement with the experimental value \({\cal M}_{10} = 0.2648(36)\). In other words, our determination of the PS mass agrees very well with the central value \(m_b^{\text{PS}} = 4.532\,\)GeV found in [11]. There, QED effects were treated differently and a non-zero mass was used for the light quark. The numerical effect on the extracted PS mass is negligible. For Mathematica, the structure is quite similar. To avoid clashes with built-in functions, we have renamed some of the variables. Needs["QQbarThreshold`"]; ZN[i_, mbPS_] := BBbarResidue[i, mbPS, mbPS, "N3LO"]; EN[i_, mbPS_] := BBbarEnergyLevel[i, mbPS, mbPS, "N3LO"]; Rb[s_, mbPS_] := BBbarRRatio[Sqrt[s], mbPS, mbPS, "N3LO"]; Mresonances[n_, mbPS_] := With[ { Nc = 3, eb = eD, mbPole = BottomPoleMass[mbPS, mbPS, "N3LO"] }, 12*Pi^2*Nc*eb^2/mbPole^2*Sum[ ZN[i, mbPS]/(2*mbPS + EN[i, mbPS])^(2*n+1), {i, 1, 6} ] ]; Mcontinuum[n_, mbPS_] := Module[ {s0, s, x}, s0 = BBbarThreshold[mbPS, mbPS, "N3LO"]^2; s[x_] := s0 + x/(1 - x); Return[ NIntegrate[ Rb[s[x], mbPS]/(s[x]^(n+1)*(1 - x)^2), {x, 0, 1}, AccuracyGoal -> 5 ] ]; ]; M[n_, mbPS_] := 10^(2*n)*(Mresonances[n, mbPS] + Mcontinuum[n, mbPS]); LoadGrid[GridDirectory <> "bbbar_grid_large.tsv"]; Print[M[10, 4.5322]]; Here, we obtain a slightly different value of \({\cal M}_{10} = 0.264706\), which stems from a different integration algorithm. Since this change can be offset by changing the input mass by a small amount of about \(0.1\,\)MeV, we can conclude that integration errors are under control. ## Grid generationDepending on the chosen input parameters, the precomputed grids distributed with QQbar_threshold may not be sufficient. In these cases custom grids can be generated with the ## Top and bottom gridsThe main function provided by this package is QQbarCalcGrid[ Energy -> {MinEnergy, MaxEnergy, EnergyStep}, Width -> {MinWidth, MaxWidth, WidthStep}, "GridFileName" ]; <<QQbarGridCalc.m; LaunchKernels[]; QQbarCalcGrid[ Energy -> {-1, 1, 1}, Width -> {1.5, 1.6, 0.1}, "top_grid_example.tsv" ]; Note that loading the package typically takes several minutes. The calculation of the grids themselves is even more time-consuming, so we restrict the examples to very small and coarse grids and suggest to rely on parallelisation as much as possible. The energy and width ranges always refer to reference values for the quark mass and the strong coupling, specified with the In some cases it is desirable to have grids that are relatively coarse in one region, e.g. at high energies, and much finer in another region. To this end it is possible to directly specify the energy and width points when calling QQbarCalcGrid[ Energy -> {{EnergiesPts ... }}, Width -> {{WidthPts ... }}, "GridFileName" ]; The following example shows how a bottom grid with a higher resolution close to the threshold can be generated: <<QQbarGridCalc.m; LaunchKernels[]; emin = 10^-6; emax = 10; n = 3; epoints = Module[ {stepfact}, stepfact = (emax/emin)^(1/(n - 1)); Table[N[emin*stepfact^i], {i, 0, n - 1}] ]; QQbarCalcGrid[ Energy -> {epoints}, Width -> {{10^-10, 10^-8}}, "bottom_grid_example.tsv", QuarkMass -> 5, AlphaS -> 0.25 ]; Note that the numerical evaluation requires at least a small non-vanishing width, which is internally set to \(10^{-9}\) for bottom quarks. For such a small width it is not possible (and not very useful) to calculate grid points with negative energies. Finally, QQbarCalcGrid[ Energy -> {...}, Width -> {...}, "GridFileName", Comments -> { "Comment in the first line of the grid file", "Comment in the second line of the grid file" } ]; The default setting ## Nonresonant grids## Grids at NLOThe second function in QQbarCalcNonresonantGrid[ MassRatio -> {...}, Cut -> {...}, "GridFileName" ]; As with <<QQbarGridCalc.m; LaunchKernels[]; QQbarCalcNonresonantGrid[ MassRatio -> {0.15, 0.30, 0.01}, Cut -> {0, 1, 0.01}, "non-resonant_grid.tsv" ]; ## Grids at NNLOQQbar_threshold comes with a built-in N The code for producing the N docker load -i https://www.hepforge.org/downloads/qqbarthreshold/nnlo_nonres_grid_entry.tar.gz and run with docker run -it amaier/nnlo_nonres_grid_entry <x> <yw> with For convenience there is also a Mathematica interface to the calculation of the N QQbarCalcNNLONonresonantGridEntry[x, yw, Verbose -> True] which returns a list of the grid entries (without the coordinates Especially for large values \(y_w \sim 1\) the calculation of the automated contribution can fail, in which case a slight change in the input parameters may help. In practice, this is not a severe problem as the automated contribution becomes essentially constant in this region. The precision of the automated calculation is not very high, and the values obtained can easily deviate from the ones in the reference grid by around 10%. Since the N Generated on Fri Apr 27 2018 08:24:08 for QQbar_threshold by 1.8.13 |