Introducing RGBeta: a Mathematica package for the evaluation of renormalization group \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \beta $$\end{document}β-functions

In completely generic four-dimensional gauge-Yukawa theories, the renormalization group \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \beta $$\end{document}β-functions are known to the 3–2–2 loop order in gauge, Yukawa, and quartic couplings, respectively. It does, however, remain difficult to apply these results to realistic models without the use of dedicated computer tools. We describe a procedure for extracting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \beta $$\end{document}β-functions using the general results and introduce RGBeta, a dedicated Mathematica package for extracting the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\text {MS}} $$\end{document}MS¯ \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \beta $$\end{document}β-functions in broad classes of models. The package and example notebooks are available from the GitHub repository.


Introduction
The renormalization group (RG) functions are fundamental quantities in quantum field theories (QFTs), governing how the dynamics of models change with energy. The β-functions, in particular, determine the flow of couplings g with the renormalization scale μ. They are a staple of BSM physics used in grand unification to extrapolate low-energy couplings to the unification scale or to generate mass spectra from high-scale input in supersymmetric models. Other applications involve the study of the ultimate ultraviolet fates of models in the search of fundamental theories (see e.g. [1,2]). It is, therefore, not surprising that the β-functions were computed to the first few orders in perturbation theory several decades ago. In the general case of four-dimensional gauge-Yukawa theories, they were calculated to the 2-loop order in both gauge, Yukawa, and quartic couplings already in the 80's [3][4][5][6]. What is perhaps more surprising is that it took another 20 years for the 3-loop result for the gauge β-function [7], which represents the current state-of-the-art in generic theories.
One contributing cause for this measured pace might be that although the general results for the β-function have long been known, it is an often time-consuming endeavor to apply the results to specific models. Typically, model builders write, a e-mail: thomsen@itp.unibe.ch (corresponding author) and work with, their models in terms of various matter fields, each in their separate irreducible representation of the model gauge group. On the other hand, the generic Lagrangian employed in the general result is framed in terms of but one real scalar and one Weyl fermion multiplet in reducible representations of the symmetry, encompassing all the matter fields of the model. Matching the specific model onto the general theory is in principle straightforward; however, in practice, it almost always turns into an arduous exercise and is the primary obstacle in applying the general formulas. Typical computations in the generic framework would be done by explicit construction of the coupling tensors, which are mostly large and sparse. In the SM for instance, the generalized Yukawa coupling Y ai j is a sparse 4 × 45 × 45 tensor and the fermion gauge generator (T A ) i j an 11 × 45 × 45 tensor. Every term in the β-functions corresponds to a monomial in these tensors, contracting all internal indices, and computations quickly become computationally expensive and/or cumbersome to set up.
For all of the above reasons, a number of computer tools have been developed for extracting model-specific βfunctions, of which RGBeta is but the latest. There is the general-purpose Mathematica 1 package SARAH 4 [8,9], which includes a 2-loop implementation of the β-functions based on the original 2-loop formulas [4,10]. Then there is the dedicated python code PyR@TE 3 [11][12][13], which in its latest iteration, based on the completely general set of basis tensors [14,15] up to the 3-loop order in the gauge β-function, is blazingly fast. With the PyLie module based on SUSYNO [16], PyR@TE 3 is capable of dealing with arbitrary gauge group representations, but it does not have the flexibility of e.g. using SU(N ) with N left unspecified. Finally, during the completion of this manuscript, we became aware of the independent, dedicated Mathematica package ARGES [17] with a significant overlap with the scope of RGBeta (even as it pertains to dedicated computation tools for RG functions) and a good degree of flexibility in the implementation of the gauge groups. Internally ARGES is based on the original RG formulation [3][4][5]7,10,18], where RGBeta is based on Ref. [14] implemented with the structure deltas of Ref. [19]. More noticeable to the user, ARGES implements the running of vacuum expectation values [20,21], whereas RGBeta includes kinetic mixing between Abelian gauge fields. Also recently, the C++ library RGE++ was introduced to provide an easyto-use framework for numerically solving RG equations in an efficient manner [22]. Although most of the new features of RGBeta has been covered in the latest PyR@TE release and ARGES, we still believe it worthwhile to release this latest tool in the line. With multiple available tools, the user can ultimately choose whichever fits them and their project best. Furthermore, with multiple independent codes, it becomes possible to cross-check the outputs. This is sorely needed considering that in all cases reliability of the output is subject to the user not making any errors in the (often complicated) input. There is also the black-box factor to consider, that is, the implementation in the various programs cannot easily be tested by the users.
RGBeta was originally developed to test the general βfunction basis used and developed in Ref. [14] and is used in ongoing work on higher-order β-functions [23]. It is an implementation of the generally applicable MS β-function formulas for gauge-Yukawa theories presented in Ref. [14], which in turn is an extension of the general 2-loop results [3-6,10] and 3-loop gauge β-function [7] to allow for completely generic gauge groups. The new formulas use a revised basis, which easily generalizes to the case with multiple Abelian gauge group factors and has been checked for inconsistencies with the Weyl consistency conditions (cf. Refs. [24][25][26]). RGBeta uses an implicit construction of the general coupling tensors to avoid large sparse tensors in the evaluation, and all group index contractions are evaluated with a handful of simplification rules and Fierz identities. Although this limits the scope of the package to fields in some of the more common representations of the ordinary Lie groups (cf. Sect. 3.4), it does come with the benefit of a usually quick evaluation time and the possibility of leaving group indices unspecified, meaning that it can handle e.g. SU(N ) groups without fixing N . This makes RGBeta ideal for using in the Mathematica notebook environment, where it is easy for the user to manipulate the output and/or make changes to the model on the fly.
The rest of the paper is organized as follows: the following section details the matching procedure from a specific model to the general framework. It provides background and insight into the underlying procedure but can be skipped by the impatient reader, who wants to learn how to use RGBeta. Section 3 explains the basic objects of the package, installation, and implementation of group theory. Next, Sect. 4 outlines the routines of RGBeta and how to implement a model, using the Standard Model (SM) as a detailed example. Finally, we round off with a short conclusion while the appendix contains more full-fledged documentation of the routines.

Generic four-dimensional renormalizable theories
The derivation of the general formulas for β-functions presupposes a certain generic formulation of the QFT. To include all marginal couplings of any four-dimensional theory, the generic Lagrangian (GL) used in the formulation of the general β-function results is compactly put on the form [14] The Fermion field i is of the Weyl type, as any Dirac spinor can be decomposed into two 2-component spinors. All fermion fields are gathered into this one multiplet in a, generally speaking, reducible representation of the flavor and gauge group, and lower-case Latin indices starting from i, j, . . . are used for the fermion indices. The scalar field a is taken to be real: complex fields can be decomposed to real fields but not vice versa. a , labeled with small Latin letters beginning with a, b, . . ., is one multiplet containing all scalar fields and is also generically a reducible representation of the gauge and flavor groups. Finally, all gauge fields are gathered into a multiplet A A μ with upper-case Latin indices A, B, . . . running over all adjoint representations of the individual product groups. F A μν is the field-strength tensor associated with the corresponding multiplet.
As for the couplings, the Yukawa coupling Y ai j is symmetric in the two fermion indices i and j and the quartic coupling abcd is completely symmetric. The gauge couplings have been absorbed into the gauge kinetic term and placed in the coupling matrix a −1 AB , which is proportional to the identity in the case of a simple gauge group and otherwise block diagonal (except for kinetic mixing terms between Abelian groups). 2 Thus, the covariant derivatives of the matter fields read In contrast to the GL, the model builder's Lagrangian (MBL), as it is commonly used in model building, is written in terms of multiple fields in irreducible representations of the gauge and global symmetry groups of the model. It takes the form where n labels the various products of the gauge group G = × n G n , and all repeated indices are summed over. This notation is generalized and -unfortunately -rather incomprehensible at a glance. 3 The Greek subscripts on the fields label the various irreducible representations of the scalar and fermion fields. Each field has a set of indices collectively denoted with a calligraphic letter, which comprises all gauge and flavor indices of the respective fields. As an example, the SM left-handed quark doublet would be written as ψ I q q with the collective index I q = (i, c, α) labeling generation, color, and isospin. For the gauge fields, A U n μ n is the field of product group G n with corresponding adjoint index U n . To apply known results and extract the β-functions for the coupling in the MBL, one must construct a mapping between it and the GL.

Mapping onto the generic Lagrangian
Let us proceed with the construction of a mapping of the MBL onto the GL. One approach would be to do an explicit construction of the a and i multiplets with each entry being filled with a particular gauge and flavor component of the fields φ α and ψ ρ . 4 Instead, we will employ and expand on the use of the structure deltas introduced by Mølgaard [19] in an approach that allows for treating most sums (matrix multiplications) of the program in a mostly implicit manner.
The peculiar form of the gauge kinetic term in the GL (2.1) is a compact way of including all the product groups of a generic semi-simple gauge group G = × n G n . The gauge couplings have been placed in the kinetic term of the gauge fields by rescaling the gauge fields A U n nμ → g −1 n A U n nμ . Next, all the gauge fields are collected into a single multiplet A A μ . The structure deltas [A n ] are introduced as projection operators 3 No doubt this is what we get for recklessly constructing a generic specific theory: a contradiction in terms. 4 See e.g. App. D of Ref. [18] for how this would be done in the SM. The reader will no doubt appreciate the effort underlying such explicit constructions and why one might wish to avoid this. satisfying As projection operators they fulfill the completeness relations Using the structure delta, one may write the gauge kinetic term (2.8) for the gauge couplings and group structure constants, respectively. One may think of a −1 AB as a block diagonal matrix, where gauge invariance ensures that the value of the coupling is the same for all the gauge fields across the adjoint representation of any one of the product groups G n . The gauge self-interactions always depend on the generalized structure constant F ABC , which also vanish between gauge fields of different product groups.
Moving on to the interactions with matter fields, the structure deltas are implemented with the purpose of projecting out the specific fields (in the phenomenological sense) φ α or ψ ρ from the general multiplets or , respectively. We define them by where e.g.
2 accounts for the different normalization of the kinetic terms of the real and complex scalar fields. Since no degrees of freedom are lost or gained from packing the fields into multiplets, the structure deltas obey the usual summation rules for projection matrices: As one would expect, the complex scalars satisfy 5 Using the structure deltas, the interaction terms of the MBL can be recast as (2.13) Through direct comparison with the GL (2.1), one then finds The sums over {i, j} and {a, b, c, d} are taken to be over all 2 and 24 permutations of the indices, respectively. In a similar manner, the kinetic terms can be made to match the GL form by inserting an identity with a pair of structure deltas. For the fermions it follows that which allows for identifying the gauge generator (T A ) i j of the fermion field multiplet. For the scalars there is a small complication due to the matching from real to complex fields.
Observe first that the scalar kinetic term can be put on the form Only then do we match to the GL, with (2.18) From here we can identify the generators (T A ) ab of the generic scalar multiplet. They are antisymmetric in a, b as indeed they should be for any real representation. It should be mentioned that any MBL could, of course, contain real scalars in addition to (or instead of) the complex ones. The formalism presented here is easily extended to accommodate such cases by identifying [φ α ] and [φ † α ] for the appropriate fields and getting rid of the corresponding √ 2. One great advantage of using the structure deltas is that they never need to be constructed explicitly as matrices. It is sufficient to invoke the completeness relations whenever GL indices are contracted to obtain a sum of contractions with MBL indices.

β-functions
Having established a mapping procedure between the MBL and GL, the known results for the Beta functions [3][4][5][6][7]10] can be applied to determine the β-functions. Specifically, the implementation of the β-functions in RGBeta uses the form given in [14], as it supports completely generic gauge groups including the effects of kinetic mixing up to the 3-loop order. In the GL notation, the general formula for the β-functions produces in terms of the couplings a AB , Y ai j , and abcd and gauge group generators (T A ) ab , (T A ) i j , and F ABC . As we saw, all of these objects can be expressed with their MBL counterparts and structure deltas, leaving only the task of projecting the generic β-functions back onto the specific couplings of the MBL.
The structure deltas are applied to map back from GL βfunctions to the β-functions of the specific MBL couplings. For the gauge couplings we invert Eq. (2.7), giving For the Yukawa and quartic interactions, a similar comparison to Eqs. (2.14) and (2.15) shows that Here we must introduce symmetry factors S n to compensate for double counting when multiple fields coincide. 6 Contrary to the gauge β-function, both the Yukawa and quartic βfunctions still have free indices. Some of these are coupling indices, such as the generation indices on the SM Yukawa couplings, and should be retained. However, all free gauge (and flavor indices irrelevant to said coupling) will have to be projected out. One should keep in mind a potential problem in projecting out the remaining symmetry indices of Eqs. (2.21) and (2.22). If the gauge and global symmetry representations of the fields allow for multiple invariants, the naive projection of one coupling may have a non-trivial overlap with the others. In such cases, one must take care to define linear combinations that pick out particular couplings. One example where this becomes relevant is in the single-trace/doubletrace terms in quartic potentials.

Treatment of other renormalization group functions
We should like to make some additional comments regarding what RG functions can be extracted with the methods discussed in this section. First of all, as pointed out in Ref. [14], the formulation of the gauge Lagrangian with a matrix gauge coupling generalizes easily to include kinetic mixing terms between multiple Abelian gauge groups, making it possible to treat any gauge group. In contrast to previous formulations of the kinetic-mixing as something other than the gauge coupling [27,28], we can include it by allowing for a matrix gauge coupling between all Abelian field-strength tensors in the gauge potential: where r, s run over all Abelian gauge fields. This includes all information of the mixing between the groups, and h rs can be embedded as a block in the general gauge coupling a AB .
Thus far we have only discussed the β-functions of the marginal couplings. Relevant couplings are also allowed in renormalizable theories, and their running is also governed by β-functions. The most general Lagrangian term for the relevant couplings take the form in the GL or in the MBL. The structure delta technique can be applied to establish a mapping between the GL and MBL relevant terms in a similar manner to what was just discussed for the marginal couplings. The actual β-functions for the relevant terms can be recovered from the Yukawa and quartic β-functions using the dummy field method [10,29]. It entails introducing nondynamic faux scalar fields to the relevant couplings to match them with marginal couplings whose β-functions can then be used for the relevant couplings. The main pitfall is that as the faux field is not dynamical, it does not receive any field-strength renormalization; all such contributions have to be removed from the corresponding β-function. At the 2loop order, one can unambiguously identify what terms in the β-function basis corresponds to the (neutral) scalar field anomalous dimensions, and these can therefore be removed without any specific knowledge of the underlying loop calculation. This method was also used by Sartore [15] for the new β-function basis [14].
The anomalous dimension of the matter fields are also known up to the 2-loop order. For a matter field η, the bare field is related to the renormalized field by where Z η is the field-strength renormalization. The anomalous dimension of the field is then given as The general 2-loop expressions found in Luo et al. [10] are easily adapted to the notation of Ref. [14] used in RGBeta.
It should be pointed out that γ η is gauge dependent, and that they have been computed using the R ξ gauges (with ξ = 0 being the Lorenz/Landau gauge and ξ = 1 the Feynman gauge).

RGBeta in a nutshell
RGBeta is a Mathematica package aimed at allowing theorists to easily extract the MS β-functions of their favorite models, be it for BSM physics or field-theoretic applications.
The package implements the state-of-the-art β-functions for general four-dimensional renormalizable theories, which are known up to loop order 3-2-2, for gauge, Yukawa, and quartic couplings, respectively. The intent is to provide Mathematica users with an easy-to-use, one-stop implementation of the β-function formulas directly in Mathematica, which is already widely used in the community for many computeralgebra tasks. Specifying a model to the package can be done in a rather compact manner and can easily be done directly in a Mathematica notebook, in which the user can proceed to extract and manipulate the various β-functions. The βfunctions can be obtained modularly for the users who wish to experiment with them in the Mathematica notebook environment.
The core of RGBeta is set up to evaluate the tensor structures of the general β-function formulas with Einstein summation conventions and various group identities. This approach is similar to what a person might do and is effective because many of the index contractions are in fact Kronecker deltas or group generators, which can be contracted with the pattern matching functionality of Mathematica. The performance of this strategy is generally very good, evaluating e.g. the full set of SM β-functions to highest loop order in less than 10 seconds on a decent laptop. 7 , 8 On top of that, it allows the user to keep the number of colors or generations arbitrary during the evaluation. The price of the implementation strategy is that it does require hard-coding various representation-specific identities. For this reason, only a selection of common representations have presently been implemented; cf. Sect. 3.4. 7 With an Intel Core TM i7-7700HQ CPU @ 2.8 GHz. 8 Other models with matter in the fundamental representations are similarly fast, though needless to say that evaluation time increases with the complexity of the model. The main obstacle to a fast evaluation is when the scalar quartic sector becomes complicated, either because of there being a large number of quartic couplings in the model or because of scalars in non-fundamental representations.

Definition of the β-functions
It pays to be explicit about the definition of the β-functions used in RGBeta, as many different normalizations are used throughout the literature. To reiterate, the package exclusively provides the MS (or equivalently MS) β-functions. In all cases we associate a loop expansion to the β-functions, writing For all non-gauge couplings, both relevant and marginal, the β-functions are defined as the logarithmic derivative of the coupling wrt. the renormalization scale: and where the ordering of flavor indices corresponds to the ordering of fields in the associated coupling. In all cases, the coupling will be the reference name of the β-function used in the package.
The only deviation from this pattern of β-function definitions is the gauge couplings. We will assume that the kinetic terms of the gauge fields and their β-functions are normalized as for the fields of non-Abelian gauge groups or the Abelian gauge group in a model with at most one Abelian group. For a model with multiple Abelian gauge groups, which therefore features kinetic mixing, the kinetic term for all the Abelian gauge fields is written jointly as where h i j is the symmetric coupling matrix. 9 The β-function β h of the coupling matrix contains all information of the running of couplings and mixing. With the normalization of the gauge fields employed in this notation, the covariant derivative of a matter field can be written as where T A n are Hermitian generators of the non-Abelian representations and q i Abelian charges. In particular, using a coupling matrix for the Abelian kinetic term allows for keeping the coupling of the Abelian gauge fields to the matter fields "diagonal".
The matter field anomalous dimensions are loop expanded as well, parameterizing They are defined in Eq. (2.26), and the indices in RGBeta are all calculated in the R ξ gauges.

Installation
RGBeta can be installed directly to the /Applications folder in Mathematica's § ¤

¥
In all cases, RGBeta ought to be loaded into a freshly initialized Mathematica kernel to avoid clashing symbol definitions.

Group theory
RGBeta is set up to evaluate the tensor structures in the β-function formulas in a manner similar to what a person would do if asked to use the 1-loop formulas by hand for a particular model (something which is typically doable with a bit of patience). That is to say, it avoids all explicit summation/matrix multiplication, except over the various field types, in favor of relying on implicit summation with Einstein summation convention for gauge and flavor indices. This is possible, and fast, for several of the common representations for the ordinary Lie groups, where gauge indices are typically contracted with Kronecker deltas or generators of the fundamental representations, which can be dealt with using Fierz identities. This approach has the benefit that it allows for keeping e.g. the index of the gauge group unspecified throughout the evaluation. Using implicit summation, RGBeta is set up with a number of symbols with internal summation rules set as Mathematica up-values, 10 a list of which is presented in Table 1. The most basic symbol used in the package is the Kronecker delta δ ab , which is represented by del[rep, a, b] with the first argument specifying the type of the indices a, b. The summation conventions assigned to the Kronecker delta makes contractions evaluate as e.g. for any representation rep. Identically named indices of different types do not contract, so e.g.  Table 2.
The tensor structures used for the 3-2-2 β-functions have been chosen such as to minimize the occurrence of nontrivial group structures stemming from the gauge groups [14]. Despite this, there remains a couple of them that cannot be evaluated by identifying factors of quadratic Casimir operators or Dynkin indices. These tensors are dealt with by implementing the Fierz identities of the fundamental representations of the groups (Table 3). In the SM, these are necessary in the 3-loop gauge β-functions and the quartic β-function already from 1-loop order. Table 1 List of various symbols used in RGBeta and their meaning. For details of the various representation we refer to Table 2 Symbol Interpretation

Bar[x]
Complex conjugation, x * . It is used general purpose for both fields, couplings, and group representations.

Trans[x]
Symbolic transposition of matrix couplings.

del[rep, a, b]
The Kronecker delta δ ab with indices running in the group representation or over flavor indices as specified by rep.

delA2[group, i, a, b]
The Clebsch-Gordan coefficient between indices i in A 2 and a, b in N (or i in A 2 and  a, b in N) representations of group.

delIndex[rep, a, b]
The Kronecker delta δ ab , where one of a, b is an integer pointing to a definite index value, e.g. δ a2 . The Levi-Civita symbol with indices in the given representation.

tGen[rep, A, a, b].
The Hermitian group generator T Aa b for the corresponding representation. The first index, A, is always taken to be in the adjoint representation, while the other two belong to rep.

Dim[rep]
The dimension of a representation or flavor index. The gauge representations are mostly predefined, but the dimension of flavor indices can be set by the user.

Matrix[A, B,. . .][a, b]
The matrix product of multiple couplings with 2 or fewer indices, e.g. (AB) ab = A ac B cb . Couplings with 1 index are interpreted as column vectors.

Tensor[A][a,. . .]
The tensor coupling, A, with three or more indices. No coupling contractions involving tensors are supported.
The need for Fierz identities and systematic symbolic treatment of group invariants are the main obstacles to implementing arbitrary representations of the gauge groups. At present RGBeta is, thus, restricted to the representations and groups listed in Table 2. In addition to fundamental and adjoint representations, N and G, of the ordinary groups, it includes the two-index symmetric and antisymmetric, S 2 and A 2 , of SU(N ); the traceless two-index symmetric S 2 of SO(N ); and the two-index antisymmetric that vanishes upon contraction with ab , A 2 , of Sp(N ). For all these representations, the generators can be decomposed in terms of the generators of the fundamental representation after which the Fierz identities can be employed. The U n (1) groups support representations with any charge assignments. The limit to these few, if frequently used, representations of the Lie groups is the primary limitation of RGBeta. None of the exceptional Lie groups or their representations are presently implemented.
The 2-index representations A 2 and S 2 are implemented using a single index in RGBeta (when referring to them as 2-index representation, it is because they can be written in terms of 2 fundamental indices, similar to how an adjoint representation can be written with a fundamental and an antifundamental index). We need a way to contract e.g. the A 2 label with two fundamental indices to use as an invariant in various interactions. For this purpose RGBeta includes the Clebsch-Gordan coefficients delA2 and delS2 (Table 1). These are really just shortcuts and could easily be constructed using the regular del by accessing the "subindices" a [1] and a [2] of e.g. the A 2 index a. Thus, delA2[G, i, a, [1], b])/2 for any group G with an A 2 representation (similarly for S 2 ). It is always preferred to use delA2 and delS2 directly, as they have additional internal identities.

Validation of the program
Like any other computer tool, RGBeta is for all intents a black box to the user. It does not matter that the underlying theory is sound, if there is a mistake in the implementation, something which can easily go unnoticed. To avoid this kind of errors, we have validated RGBeta against several results from literature and PyR@TE 3. The RGBeta implementation of all the models used for validation are included with the package in /Documentation/Sample_Models.nb.
First, we compared directly with the SM and type-III 2HDM (two-Higgs-doublet model) computation [30] up to order 3-2-1 with matrix Yukawa couplings. This comparison allowed us to settle a discrepancy in the 2-loop Yukawa β-function between Refs. [4,6,10,14], settling in favor of Refs. [6,14] (recall that RGBeta is based on Ref. [14]). 11 The terms in question are degenerate in the SM but not in the 2HDM. Furthermore, a comparison with the 2HDM result [30] revealed a typo in the 3-loop gauge β-functions of their result, which has since been fixed. We also have agreement with the 2-loop SM quartic result of Ref. [31].
To test the handling of unspecified gauge and flavor groups, we reproduced the 3-2-1 β-functions of the Litim-Sannino model, keeping N f and N c as free variables, and found complete agreement with Ref. [2]. The implementation of gauge kinetic mixing was tested in the SM with a gauged U(1) B−L symmetry (see e.g. [32]) against the PyR@TE 3 results [11]. Likewise, we compared the results for the SU(5) grand unified theory (GUT). 12 In both cases complete agreement was found up to order 3-2-2.

Using RGBeta
This section describes how to use RGBeta in some detail. We heartily recommend prospective users to have a look at the tutorial notebook /Documentation/Tutorial. nb, which illustrates some practical use cases in the form of the SM and Litim-Sannino model [2] with comments.

Package routines
To use the RGBeta package, one should start by specifying a model. First, the user should define symmetries of the model with the routines This is all that is necessary for defining the model. Once the model has been loaded into the kernel, the RG functions can be extracted with a minimum of effort using the routines The Appendix contains more detailed documentation for the various routines. Here, we explain the basic use of the package in the next subsections with a detailed example.

Setting up a model
We detail the core use of RGBeta with the SM as a concrete example, as it will be familiar to most/all users. The package is set to work with Weyl spinors to be able to treat all fermions in the same manner. The SM matter fields have charge assignments q ∈ (3, 2, 1/6),ū ∈ (3, 1, 2/3),d ∈ (3, 1, −1/3), ∈ (1, 2, 1/2),ē ∈ (1, 1, −1), There are three Yukawa couplings given by where c is a color index, α, β SU(2) L indices, and i, j generation indices. The reason for specifying the Yukawa couplings with right-handed (Hermitian-conjugated) fermions is to match the conventional definition of the couplings in the Dirac notation. 13 The SM also contains a scalar potential with a quartic Higgs interaction and a Higgs mass term: The RGBeta package automates most of the process of getting the β-functions. It is, however, unavoidable that the user will have to input the model in a precise manner. Specifying gauge groups and scalar and fermion field content is fairly straightforward. Arguably, the most difficult part in using RGBeta is specifying the Yukawa and quartic interaction: the user must manually specify how the flavor and gauge indices are contracted, which is a potential source of errors.
The first thing to do when defining a model is to specify the gauge groups to RGBeta with AddGaugeGroup. For each product group, one needs simply to specify the coupling, group name, and Lie group. The choice of Lie group will then set up the group invariants and generator properties of the supported representations. The unique name of the group is used for referencing the group representations. For the SM, the SU(3) c × SU(2) L × U(1) Y gauge group is added with § ¤ Having set up the gauge group, the next step is adding the matter content with AddFermion and AddScalar. It is sufficient to know the charges and flavor indices of the fields to do so. As previously mentioned, all fermions must be added as left-handed chiral fields. If any of the fields in the model are given as right-handed fermions, simply add it as a lefthanded spinor in the conjugate representation under the symmetries. The fermions are specified with a field name, and all their non-trivial gauge charges are given as a list with the GaugeRep option (the default being no gauge charges). A full list of the available representations is given in Table 2. Given that the fields come in three generation, they are also given a single flavor index with the FlavorIndices option. 14 To add the fermions to the model, we call § ¤

¦ ¥
The last line may seem a bit curious. It is simply there to specify the dimension of the flavor index "gen." In the SM there are, of course, three generations, but it is often kept as a free parameter in β-function computations. Either can be used in RGBeta. We proceed to add the Higgs field in a similar manner. It has just one generation, so there is no need to give it any flavor indices: § ¤

¦ ¥
The package assumes scalar fields to be complex-valued by default, as is the SM Higgs field. If we wished to add a real scalar field, we should give the option SelfConjugate→ True . In models with exact global symmetries, the flavor group can be added with DefineLieGroup, and the field representations can be specified by using the group representation in theFlavorIndices option for the fields. The most involved part of defining the model is setting up the couplings. In particular, one must take great care when defining the contraction of all field indices and specify the indices of the couplings. The particulars of this should be passed to the package in terms of pure functions. It is important to always remember to surround pure functions given as options with parentheses. Otherwise, the function call will

¦ ¥
The first AddYukawa call specifies the up-type Yukawa coupling. The first two arguments gives the coupling name, yu, and the fields entering the interaction, H, q, u, the first of which is always taken to be the scalar. Do bear in mind that the order of the fields matter for the construction of the index contractions. The Yukawa couplings are defined with left-handed spinors by convention. It is, therefore, the conjugated terms of the Yukawa Lagrangian (4.2), which are used as input.
The GroupInvariant option should be given a function with three arguments that specify the contraction between the indices of the fields. As the name suggests, this should be an invariant of the gauge and global symmetries. It is a function because it will be called every time the coupling appears in the β-function formulas to give appropriate labels to the invariant. For the up-type Yukawa, the del[SU3c@fund, #2, #3] part of the group invariant establishes that the fundamental SU(3) c indices of the second and third fields are contracted with a Kronecker delta. Furthermore, the fundamental SU(2) L index of the Higgs (the first field) is contracted by an -invariant of SU(2) L with the left-handed quark (second field) as specified by eps [SU2L@fund, #1, #2]. 15 The remaining free indices of the fields are the generation indices of the quarks. The coupling itself carries these indices, which we spell out with the option . This option should be passed a function of 3 arguments (one for each field) that returns a list of the indices of the coupling. In practice we can almost think of the arguments #1, #2, and #3 as the indices of the fields, and it is not a problem to have e.g. multiple #2 indices when they belong to different representations: there is no ambiguity.
The last option passed to AddYukawa is more of a qualityof-life option than strictly necessary. The Yukawa couplings in RGBeta are always given in terms of left-handed spinors as in the MBL (2.3), but conventionally the SM defines the Yukawa couplings on the right-handed spinors, and the com-plex conjugate of the coupling appears with the left-handed fermions. To tell RGBeta that the un-conjugated coupling should be placed with the right-handed spinors, we need to pass it the option Chirality→ Right (Left by default). The other two Yukawa couplings follow in mostly the same manner as the up-type Yukawa, the main difference being that they involve the complex conjugated Higgs field. This is specified by using Bar@H for the scalar field in AddYukawa. With the conjugated Higgs field, the SU(2) L contraction is done with a Kronecker delta instead of the antisymmetric invariant.
The ¦ ¥ Again, the function takes the coupling and the fields (four scalars this time) as arguments. In contrast to the matrix Yukawa couplings of the SM, λ is a scalar coupling, so no coupling indices have to be passed to the function. The 1/2 normalization of the Higgs self-coupling in Eq. (4.3) is put in the group invariant. RGBeta never assumes any normalization factors even when, as is the case with the SM, there is a symmetry factor associated with the coupling.
In ¦ ¥ completely parallel to the implementation of the other couplings. This concludes the implementation of the SM.
Before proceeding to extract the SM β-functions, we would like comment on cases with other invariants in the couplings. In some cases, it is possible to construct new invariants using the group invariants already implemented in RGBeta. In such cases, the user can manually define the relevant contractions, but it is important to avoid repeating index names in the internal contractions. The best way to avoid this is to use a combination of SetDelayed and Module to ensure unique naming of the indices every time the invariant is called. An example is SU(5) GUT, where there is an adjoint scalar , which has two independent invariants with four fields, one of them being ( a a ) 2 . For the second, independent invariant, one often uses a trace over 4 fundamental generators, which can be implemented in the program as The full SU(5) implementation is included with the package in /Documentation/Sample_Models.nb.

Producing the β-functions
Once the SM model has been loaded into the Mathematica kernel, extracting the β-function is as simple as anything. To obtain, for instance, the 2-loop β-function of the downtype Yukawa coupling, one simply has to call the function BetaFunction[yd, 2]. This function returns the β-function as defined in Eqs. (3.2-3.5).
A typical use of BetaFunction (selected for minimality) will look like  $i is used to denote an index of the first fermion of the coupling and $j is an index of the second fermion (as they are given when specifying the coupling). 16 With how we implemented the ye coupling, these are the l and e generation indices, respectively. For further manipulation with the β-function, one will typically wish to remove such explicit indices to leave the matrix structure implicit. This can be done with the Finalize routine, which can also be used to substitute Matrix/vector couplings with a list of their entries. We should mention that BetaTerm can be used to single out the contribution to a β-function at a particular loop order, β ( ) g . In many models, there will be multiple singlets in the product of the same four scalars allowing for multiple couplings. The quintessential example of this behavior is the single-/double-trace coupling of scalars in the fundamental representation of two different groups. In such events, the naive operators used in RGBeta to project out specific couplings will mix the couplings in question and BetaTerm /BetaFunction cannot be relied on to produce the correct βfunctions for these couplings. This behavior can also be seen with CheckProjection, which will mix the couplings in question. In such cases, one should rely on QuarticBetaFunctions , which extracts all quartic β-functions simultaneously and accounts for the mixing by identifying the correct linear combination of the projection operators.
This concludes our discussion of the functionality of RGBeta. Examples of further manipulation of the βfunctions can be found in /Documentation/Tutorial. nb. Here one can also find an example of how matrix couplings can be parametrized to obtain e.g. the β-function of the charm Yukawa coupling.

Summary and conclusions
In recent years there has been a renewed interest in higherorder β-functions for model building and precision physics. The theory of the general β-function has also matured to the point where we now have a completely general formalism that treats the couplings in a unified way, correcting several mistakes along the way [7,14,15,29]. RGBeta is a Mathematica package that leverages these theory advances along with the structure delta approach to Lagrangian matching [19] in a minimal yet fairly general tool for the automatic computation of β-functions in four-dimensional renormalizable models.
The RGBeta package presented in this paper is an attempt at striking a balance between having an intuitive and minimal way of implementing models (the problematic part of using any RG tool) and maintaining a good degree of generality. It is fast enough that it can typically be evaluated directly in the notebook environment, allowing the user to experiment with the resulting β-functions. It also has the advantage that group indices can be kept arbitrary. It is our ambition to implement the full set of β-functions up to order 4-3-2 in RGBeta pending ongoing work with Davies and Herren [23]. ing out the subtleties in the dummy-field method for extracting the β-functions for the relevant couplings. Likewise, Florian Herren was very helpful in tracking down the issue in the 2HDM β-functions and carefully reading the manuscript. This work has received funding from the Swiss National Science Foundation (SNF) through the Eccellenza Professorial Fellowship "Flavor Physics at the High Energy Frontier" project number 186866.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The Mathematica package can be found in the external GitHub repository linked to in the article.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 . defines a Yukawa coupling in the model. coupling specifies the coupling constant of the Yukawa interaction.
{phi, psi1, psi2} is the list consisting of the scalar phi and two fermions psi1 and psi2 of the interaction. The scalar can be conjugated with Bar.

Options
GroupInvariant → (1 &) is a pure function of 3 arguments defining the group invariants of the coupling. CouplingIndices → (Null &) is a pure function of 3 arguments that specifies the tensor indices of the coupling. Chirality → Left sets whether the coupling appears with left-handed or right-handed fields in the Lagrangian.

AnomalousDimension[field, loop, Options]
gives the full anomalous dimension of a field up to the given loop order. field the field one wishes to find the anomalous dimensions of. loop is either 1 or 2 specifying to what loop order.

Options
RescaledCouplings → False determines whether the couplings should all be rescaled with g → 4πg, y → 4π y, and λ → (4π) 2 λ in the output.

AnomalousDimTerm[field, loop]
gives the term γ ( ) φ in the anomalous dimension. field the field one wishes to find the anomalous dimensions of. loop is either 1 or 2 specifying the loop order.

BetaFunction[coupling, loop, Options]
gives the full β-function of a coupling up to the given loop order. coupling the coupling one wishes to find the β-function of. loop is an integer specifying to what loop order.

Options
FourDimensions → True determines if the β-function is strictly four-dimensional or whether it should include the zeroth order O( ) term in (4 − ) dimensions. RescaledCouplings → False determines whether the couplings should all be rescaled with g → 4πg, y → 4π y, and λ → (4π) 2 λ.

BetaTerm[coupling, loop]
gives the term β ( ) g in the β-function. coupling is a symbol corresponding to the coupling of the relevant β-function. loop is an integer specifying the loop order.

CheckProjection[coupling]
returns the value of the internal projection operator applied to the general coupling. coupling specifies the coupling to be checked.

DefineLieGroup[groupName, lieGroup[n]]
sets up all group constants associated to the a Lie group of the specified kind. This can be used e.g. if the model contains a global symmetry. groupName specifies the name of the added gauge group. lieGroup[n] specifies the type Lie Group. It can either be SU [n], Sp[n], or SO[n] with n either an integer or symbol.

Finalize[expr, Options]
further manipulates the expression (typically a β-function), removing indices unnecessary indices and allowing for couplings substitutions. expr an expression on the form such as given by BetaTerm.

Options
Parametrizations → {} a set of substitution rules, to replace the coupling symbols with e.g. coupling matrices. BarToConjugate → False replaces the RGBeta head Bar with the standard Mathematica Conjugate, giving an expression more suitable for further numerical analysis.

QuarticBetaFunctions[loop, Options]
gives the full β-functions of all quartic couplings using fully diagonalized projectors. loop is an integer specifying to what loop order.

Options
FourDimensions → True determines if the β-function is strictly four-dimensional or whether it should include the zeroth order O( ) term in (4 − ) dimensions. RescaledCouplings → False determines whether the couplings should all be rescaled with g → 4πg, y → 4π y, and λ → (4π) 2 λ. instructs RGBeta to treat one or more symbols (typically couplings) as being real by setting Bar@symb = symb. symbol the symbol that will be defined to be real.