1 Extensions to CAS by Rubi
Computer Algebra Systems (CAS) such as Mathematica (Wolfram Research, Inc., 2018) and SymPy (Meurer et al., 2017) (the popular opensource alternative implemented in Python), have builtin symbolic integral routines. Rulebased Integration (Rubi) developed by Rich and Scheibe (2018)
is principally a package (designed for Mathematica) that provides a method of symbolic integration organized by decision tree pattern matching, which matches the form of the integral against known integral rules. Rubi comprises 6700+ rules, collated from familiar favourites
Abramowitz and Stegun (1964); Beyer (1991); Gradstejn and Ryzik (1994) and in doing so it offers not only a means of integrating, but a growing complete reference for integration rules. These rules are in humanreadable form with cross references to Rubi rule numbers, and to the source. Rubi can also print the rules applied at each stage of solving the integral – a useful technique for pedagogical and diagnostic purposes.Without proper consideration it may not be obvious why Rubi marks a significant improvement to effectively solving integrals. The effectiveness of these routines have been independently investigated by Abbasi (2018) with the results presented in Table 1. Comparing Rubi 4.15.2, Mathematica 11.3 and SymPy 1.1.1 Abbasi (2018) divides the quality of integral’s antiderivatives into four groups.Group A consists of integrals that were easily solved: where the antiderivative is optimal in quality and leafsize. Group B is the group of integrals which were solved, but the leafsize twice that of optimal. Group C’s integrals were solved, but the solution contains hypergeometric functions, special functions or imaginary units while the optimal antiderivative does not. Finally Group F are all integrals which cannot be solved by the CAS. See Abbasi (2018) for more details.
System  % A grade  % B grade  % C grade  % F grade 
Rubi 4.15.2  99.76  0.08  0.06  0.1 
Mathematica 11.3  75.37  8.46  15.81  2.67 
SymPy 1.1.1  30.29  0  0  69.71 
Adapted: Abbasi (2018) pg. 6 
Rulebased integration is the focus of much attention in development not only by Rich and Scheibe (2018), but by others. For example SymPy 1.1.1 currently fairs comparatively poorly in symbolic integration to other CAS (Table 1). However, the rules and implementation (pattern matching in a decision tree) behind Rubi are currently being developed into SymPy, see SymPy and Collaborators (2017) for details. This would clearly improve the quality of this open source alternative.
2 Motivating Example: Quark Mass Renormalization Group Equation
Having examined how powerful Rubi is as a symbolic integration tool, we now explore how it can be applied in computation, using the quark mass renormalization group equation as a motivating example. This example is purposefully chosen, as will become apparent later, because the integral of interest is challenging for CAS.
Along with the strong coupling, the quark masses are fundamental parameters of Quantum Chromodynamics and it is therefore important to accurately know their numerical values. Further, it is important to know the scale dependence of these values.
In QCD, as in Quantum Electrodynamics (QED), one removes the present divergences with a technique known as renormalization. A nonphysical renormalization scale parameter, , is introduced in the renormalization procedure to represent the point at which one performs the subtraction of the divergences to render the amplitudes finite. Both the renormalized coupling , and the quark masses, , depend on the renormalization scheme used to define the theory, and on the scale parameter, . If we set approximately equal to the scale of the momentum transfer Q in a particular interaction, becomes the effective strength of the strong coupling for that interaction (Tanabashi et al., 2018). Throughout this paper, we make use of the physical energy scale parameter, (where ). The scale dependence of and is governed by corresponding renormalization group equations (RG equations) which rely on QCD’s anomalous dimensions as input.
The strong coupling, , satisfies the differential RGE (Davier et al., 2006):
(1) 
where the function is known up to , and ( is the gauge coupling of QCD). Given the renormalization point, the function describes how the strong coupling depends on the momentum transfer.
The quark masses, , satisfy the differential RGE (Davier et al., 2006):
(2) 
where the function is an anomalous dimension and is known up to . The dependence of and in Eqs.(1)(2) is implicit i.e. and .
The coefficients of the function, which are now known to fiveloop order (Baikov et al., 2017; Luthe et al., 2017; Herzog et al., 2017), are given (for three active quark flavours) by: , , etc. While the function coefficients, also currently known to fiveloop order (Baikov et al., 2014; Chetyrkin, 1997; Vermaseren et al., 1997), are: , , etc., for three flavours. Throughout this paper we work in the modified minimal subtraction scheme () (Veltman et al., 1972; Bardeen et al., 1978). This renormalization scheme is the most commonly used scheme in QCD perturbation theory.
It is important to be aware that there are consequences to crossing flavour thresholds. We concentrate on deriving the series expansion of Eq.(2) in the light quark sector (the up, down and strangequark). If we proceed to higher energies (into the heavy quark region), the renormalization scale crosses quark mass flavour thresholds: finite threshold corrections appear (Chetyrkin et al., 1998) and the scale dependence of the mass then needs to be matched above and below the threshold. One has to specify a new initial condition for the running coupling constant at each threshold. In principle, Eq.(8) is valid for any number of quarks , between two thresholds provided the correct initial values are used. For the light quark sector we set .
The recent calculation of the function to fiveloop order by Baikov et al. (2017); Luthe et al. (2017); Herzog et al. (2017), has ensured that the series expansion of the running quark mass can now be calculated to fiveloop order. Previously this series expansion has been calculated by Chetyrkin et al. (1997) to fourloop order, which built on the threeloop order calculation done by Kniehl (1996). Chishtie et al. (2018) provide a recent exposition into the topic to fourloop order without the use of CAS, but stop short of explicitly providing the series expansion. The perturbative series solution to Eq.(2) involves performing a Taylor expansion of at some reference scale , in powers of . To the third and fourthloop this calculation is a fairly trivial exercise. At higher loop orders, however, this computation becomes more difficult and CAS such as Mathematica and SymPy struggle to intuitively solve the RG equation without the additional use of Rubi. We outline the method for using Rubi to find a perturbative solution to Eq.(2) in the following section. The derivation is purely symbolic.
3 The Perturbative Series Expansion of
The quark mass RG equation (Eq.(2)) can be identified as a linearly separable differential equation. As such, we are able to exactly solve for given the coefficients of the and functions to a certain order. The exact solutions to leading and nexttoleading order are given in Kniehl (1996). However, it is difficult to obtain the exact solution of at higher orders, and this becomes a numerical procedure. Therefore, it is more lucid to solve the renormalization group equations in terms of a power expansion; since this type of solution provides insight into the renormalization scheme dependence of the running quark mass on the energy scale parameter , at higher powers. This is important in accurately determining the light quark mass at a chosen scale. Hence, we proceed with determining a perturbative series expansion of Eq.(2).
This is achieved by dividing Eq.(2) by Eq.(1) and linearly separating the differentials to yield
(3) 
Integrating Eq.(3) leads to
(4) 
Which can be easily rearranged to find
(5) 
where is the initial condition.
Both Mathematica and Rubi can be used in attempts to solve the integral in Eq.(5). What is of interest is how each of these CAS approach solving the chosen problem. In terms of the integral classification we introduced in Section 1
, we can classify the integral in Eq.(
5) as a Group F integral, which means that Rubi and Mathematica are unable to solve the integral analytically. Naively using Mathematica’s inbuilt integration function immediately yields an answer in terms of a RootSum object. Mathematica then struggles to find the definite integral (and series expansion) due to infinities arising from the logarithmic terms in this RootSum object. Mathematica’s solution has a low interpretability and it’s not what part of the integration process eventually yields the RootSum object.Comparatively Rubi’s attempt at the integral is a partial solution involving lower order integrals. The key advantage Rubi offers here is in simplification and clarity in identifying the unevaluated sections of the problem. Rubi performs the integral stepwise while printing the integration rule that it employs at each stage – which is is worth emphasising. This allows the researcher to focus on what Rubi does not know. Should the researcher find an analytical solution for these unknown integrals it is easy to develop the appropriate rule and submit it to the Rubi GitHub project. The high interpretability of Rubi’s attempt means that these remaining integrals can then be suitably approximated. Finding the series expansion from this point is straightforward.
Rubi’s attempt at the indefinite version of the integral in Eq.(5) yields
(6) 
where
(7) 
The integrals , do not at present have an analytic solution in terms of algebraic functions – at least, they are unknown by Rubi. At this stage, however, Mathematica is able to rewrite these integrals in terms of RootSum objects (without logarithmic divergences) that can be suitably simplified when the series expansion is performed.
The definite integral of Eq.(6) is found simply by using the Fundamental Theorem of Calculus^{1}^{1}1An assumption of the Fundamental Theorem of Calculus is that the function to be integrated must be continuous. In the present case, the integrand is a rational function and therefore continuous up to isolated poles in the complex plane.. The upper bound of the definite integral, the scale dependent strong coupling , is rewritten as its perturbative solution in terms of some known (e.g. at the taulepton mass scale) up to (Davier et al., 2006). The resulting definite integral is quite lengthy, despite some simplification occurring between polynomial sums arising from the integrals in Eq.(6). It can be viewed in the supplementary Mathematica notebook.
Finally focusing on Eq.(5), we exponentiate the definite integral, and perform a series expansion at some reference scale .
Reordering the perturbative solution in terms of yields
(8) 
where .
This is the updated series expansion of the quark mass renormalization group equation to fiveloop order. Up to threeloop order Eq.(8) agrees exactly with Kniehl (1996), and up to fourloop order with Chetyrkin et al. (1997).
For three active quark flavours, substituting the known values of the and coefficients into Eq.(8) results in
(9) 
with the Riemann zetafunction, in the light quark sector and .
4 Evaluating the Accuracy of the Series Expansion of
Eq.(8) is the perturbative series expansion of the running quark mass in powers of , with the initial value to fiveloop order. We are now interested in the effect of the latest loop order (i.e. the term). To do this we compare Eq.(8) with the fourloop series determined by Chetyrkin et al. (1997). Alternatively, we could directly numerically integrate Eqs.(1)(2) to find the running coupling and running quark mass (noting that discontinuities arise at flavour thresholds). This is the method employed in RunDec, a Mathematica (and C) package used for the decoupling and running of the strong coupling constant and quark masses developed by Chetyrkin et al. (2000) and now in its third version.
Fig.(1) provides a local error analysis, by plotting the difference between the direct numerical integration of Eq.(2) for the running of and i). the perturbative series solution to fiveloop order (Eq.(8)), ii). the perturbative series solution to fourloop order (Chetyrkin et al., 1997). Varying the energy scale between and in increments of 0.001, describes 4001 points at which to evaluate . We set the initial quark mass condition to be (Dominguez et al., 2018), where is defined as
(10) 
We have also made use of the strong coupling constant which is found using the perturbative series expansion of the strong coupling RG equation (Davier et al., 2006) with the initial condition (Pich, 2017).
The direct numerical integration approach can be used as a reference from which we statistically compare how well it is approximated by the fourloop (Chetyrkin et al., 1997) and by the fiveloop (Eq.(8)) series expansion. Two common statistical evaluation criteria: Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) are used to provide a global error analysis. The Root Mean Squared Error is defined as , and the Mean Absolute Error is calculated as .
Where is our reference i.e. calculated by directly numerically integrating Eq.(2) at each point in the range described; and is the quark mass calculated using the perturbative series solution to either the four or fiveloop order at a particular point within the range.
The MAE can be interpreted as the average error rate, while the RMSE is more sensitive to a large deviation between the function and the reference function at a single point. The MAE and RMSE for the four and fiveloop perturbative series solution are given in Table (2). The largest absolute deviation is also given, in order to provide context for the MAE and RMSE values.
Statistic  Fiveloop series solution  Fourloop series solution 
Mean Abs. Error  0.0015  0.0079 
Root Mean Squared Error  0.0029  0.0146 
Largest Abs. Deviation  0.0116  0.0578 
Smallest Abs. Deviation  0  0 
From Fig.(1). and the low MAE and RMSE in Table (2), we conclude that the fiveloop perturbative series solution for the quark mass does not deviate significantly from the direct numerical integration of the mass RG equation. Hence the correction to the series solution of the quark mass RG equation is a valuable addition.
5 Validation of Results
It has been established in Section 3 we find the peturbative series solution of the running quark mass by focusing on the separated RG equation (Eq.(5)): integrating the given integral, using Rubi, exponentiating the result, followed by Taylor expanding around a reference point . Mathematica, in comparison, can not find the peturbative solution in this way, since it fails to evaluate the integral in Eq.(5) into a useful form.
This aside, we can validate the result through a different method, this time more agreeable to Mathematica. The method relies on: i). interchanging the limiting processes (performing the series expansion before integration)^{2}^{2}2To be able to replace the integrand with its Taylor expansion, and then integrate term by term; a sufficient criterion is that the series expansion converges uniformly. This is indeed the case here, with the higher order terms having a decreasing contribution., ii). calculating the indefinite integral, followed by using the FTC to find the definite integral, iii). performing a second Taylor expansion after exponentiating the resultant integral.
The first stage of this process is to series expand the integrand of Eq.(5) which yields
(11) 
where the higher order terms are given in the supplementary Mathematica notebook.
Eq.(10) is now easily integrated with respect to . After taking the definite integral and exponentiating, we then perform a second Taylor expansion in order to yield the resultant series solution. See the supplementary Mathematica notebook for further details.
While the double series expansion may seem nonintuitive, it is able to reproduce the quark mass series expansion to fiveloop order. Which, in turn, provides validation to the central method of this paper – using Rubi.
6 Concluding Remarks
The case for using Rubi as a tool in this situation, and in other Science, Technology, Engineering and Mathematics (STEM) research areas, is thus: it provides a lucid and intuitive approach to solving integrals, which other CAS systems are often unable to solve directly. We have shown this through the motivating example of series solution of quark mass renormalization group equation.
Acknowledgements:
The authors wish to thank Hubert Spiesberger for insightful discussions, and Giulia Zanderighi for a helpful comment about interchanging the process of integration and series expansion.
Notice: The Mathematica code used is attached as a supplementary resource.
References
 (1)

Abbasi (2018)
Abbasi, N. (2018), ‘Computer algebra
independent integration tests’.
https://www.12000.org/my_notes/CAS_integration_tests/reports/rubi_4_15_2/index.pdf  Abramowitz and Stegun (1964) Abramowitz, M. and Stegun, I. (1964), Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables, United States National Bureau of Standards.
 Baikov et al. (2014) Baikov, P., Chetyrkin, K. and Kühn, J. (2014), ‘Quark mass and field anomalous dimensions to ’, Journal of High Energy Physics 2014(10), 76.
 Baikov et al. (2017) Baikov, P., Chetyrkin, K. and Kühn, J. (2017), ‘Fiveloop running of the QCD coupling constant’, Physical review letters 118(8), 082002.
 Bardeen et al. (1978) Bardeen, W., Buras, A., Duke, D. and Muta, T. (1978), ‘Deepinelastic scattering beyond the leading order in asymptotically free gauge theories’, Physical Review D 18(11), 3998.
 Beyer (1991) Beyer, W. (1991), CRC standard mathematical tables and formulae, 29. ed. edn, CRC Press, Boca Raton [u.a.].
 Chetyrkin (1997) Chetyrkin, K. (1997), ‘Quark mass anomalous dimension to ’, Physics Letters B 404, 161–165.

Chetyrkin et al. (1997)
Chetyrkin, K., Kniehl, B. and Sirlin, A. (1997), ‘Estimations of order
and corrections to massdependent observables’, Physics Letters B 402(34), 359–366.  Chetyrkin et al. (1998) Chetyrkin, K., Kniehl, B. and Steinhauser, M. (1998), ‘Decoupling relations to and their connection to lowenergy theorems’, Nuclear Physics B 510(12), 61–87.
 Chetyrkin et al. (2000) Chetyrkin, K., Kühn, J. and Steinhauser, M. (2000), ‘Rundec: A mathematica package for running and decoupling of the strong coupling and quark masses’, arXiv preprint hepph/0004189 .
 Chishtie et al. (2018) Chishtie, F., McKeon, D. and Sherry, T. (2018), ‘A systematic expansion of running couplings and masses’, arXiv preprint arXiv:1806.02534 .
 Davier et al. (2006) Davier, M., Höcker, A. and Zhang, Z. (2006), ‘The physics of hadronic tau decays’, Reviews of modern physics 78(4), 1043.
 Dominguez et al. (2018) Dominguez, C., Mes, A. and Schilcher, K. (2018), ‘Upand downquark masses from QCD sum rules’, arXiv preprint arXiv:1809.07042 .
 Gradstejn and Ryzik (1994) Gradstejn, I. and Ryzik, I. (1994), Table of integrals, series, and products, 5th ed. edn, Acad. Press, New York [u.a.].
 Herzog et al. (2017) Herzog, F., Ruijl, B., Ueda, T., Vermaseren, J. and Vogt, A. (2017), ‘The fiveloop beta function of yangmills theory with fermions’, Journal of High Energy Physics 2017(2), 90.
 Kniehl (1996) Kniehl, B. (1996), ‘Dependence of electroweak parameters on the definition of the topquark mass’, Zeitschrift für Physik C: Particles and Fields 72(3), 437.
 Luthe et al. (2017) Luthe, T., Maier, A., Marquard, P. and Schröder, Y. (2017), ‘Complete renormalization of QCD at five loops’, Journal of High Energy Physics 2017(3), 20.
 Meurer et al. (2017) Meurer, A., Smith, C., Paprocki, M., Čertík, O., Kirpichev, S., Rocklin, M., Kumar, A., Ivanov, S., Moore, J. and Singh, S. (2017), ‘Sympy: Symbolic computing in python’, PeerJ Computer Science 3, e103.
 Pich (2017) Pich, A. (2017), Precision physics with QCD, in ‘EPJ Web of Conferences’, Vol. 137, EDP Sciences, p. 01016.

Rich and Scheibe (2018)
Rich, A. and Scheibe, P. (2018),
‘Rulebased integration (Rubi)’.
https://rulebasedintegration.org/ 
SymPy and Collaborators (2017)
SymPy and Collaborators (2017),
‘Rubi integrator’.
https://github.com/sympy/sympy/issues/12233  Tanabashi et al. (2018) Tanabashi, M. et al. (2018), ‘Particle data group review’, Phys. Rev. D 98, 030001.
 Veltman et al. (1972) Veltman, M. et al. (1972), ‘Regularization and renormalization of gauge fields’, Nuclear Physics B 44(1), 189–213.
 Vermaseren et al. (1997) Vermaseren, J., Larin, S. and Van Ritbergen, T. (1997), ‘The 4loop quark mass anomalous dimension and the invariant quark mass’, Physics Letters B 405(34), 327–333.
 Wolfram Research, Inc. (2018) Wolfram Research, Inc. (2018), ‘Mathematica, Version 11.3’. Champaign, IL.
Comments
There are no comments yet.