RGESolver: a C++ library to perform renormalization group evolution in the Standard Model Effective Theory

Renormalization group evolution above the electroweak scale is a crucial ingredient in the phenomenology of the Standard Model Effective Theory. The RGESolver open-source C++ library performs the evolution at leading order for dimension-six operators in the most general flavour scenario (assuming lepton and baryon number conservation). Given its efficiency, RGESolver can be used to include the effects of renormalization group evolution in extensive phenomenological analyses in the framework of the Standard Model Effective Theory.


Introduction
The Standard Model (SM) is one of the biggest scientific successes of our time: it describes three (weak, strong and electromagnetic) of the four currently known elementary interactions in nature, closing a long path started between the 19th and the 20th century.The particle content of the SM was completed in 2012 with the discovery of a Higgs-like boson with properties consistent with the SM within current experimental errors.In spite of its success, however, the SM leaves several phenomena unexplained.For example, neutrino masses, the origin of the baryon asymmetry in the universe and dark matter suggest the presence of physics beyond the SM.The absence of direct evidence of new particles at energies O (TeV) allows us to parametrize the effects of possible heavy New Physics (NP), lying beyond the reach of the LHC for direct production, with an effective field theory, known as the Standard Model Effective Field Theory (SMEFT) [1,2].The SMEFT Lagrangian contains the SM Lagrangian plus a complete set of independent higher-dimensional operators.Working in the framework of the SMEFT makes it possible to search for NP through its virtual effects in a general, model-independent way.
Going from the SM to the SMEFT entails dramatic phenomenological consequences, since the SM enjoys several accidental symmetries that are potentially broken by higherdimensional operators, such as Baryon (B) and Lepton (L) number conservation, or the absence of tree-level Flavour Changing Neutral Currents (FCNC).Phenomenology requires the coefficients of those higher-dimensional operators that violate accidental symmetries of the SM to be tiny, implying in turn that any NP not too far from the electroweak (EW) scale must be at least approximately invariant under the accidental symmetries of the SM.However, while SM interactions will not generate B or L violation perturbatively if the corresponding operators in the SMEFT have vanishing coefficients, FCNC's will always be generated, even if NP is invariant under the full U(3) 5 flavour symmetry group of SM gauge interactions, due to the SM Yukawa couplings.Therefore, the flavour properties of the SMEFT Wilson coefficients must be specified at the NP scale, and then the coefficients must be evolved using Renormalization Group Equations (RGE's) down to the scale relevant for the processes of interest in order to compute the NP contributions.Conversely, a phenomenological bound on low-scale Wilson coefficients can be turned into a bound on the coefficients at the NP scale, allowing to extract information on the viable values of SMEFT coefficients and, hopefully, on the symmetries of the NP models of interest.
Assuming B and L conservation but a general flavour structure, the SMEFT has 2499 independent operators, so that the full system of RGE's involves more than 2500 parameters.At Leading Order (LO), the renormalization group (RG) evolution is dictated by the Anomalous Dimension Matrices (ADM's) of the SMEFT operators, which in general receive contributions from scale-dependent gauge and Yukawa couplings, as well as from the Higgs self-coupling.Since the ADM contributions proportional to different couplings are in general non-commuting, an analytic resummation of logarithmic contributions cannot be achieved and the numerical solution of the full system of equations is the only possibility, in particular when the NP scale is much heavier than the EW one.RGESolver is an open-source C++ library that performs the RG evolution of the SMEFT Wilson coefficients in a fast and easy-to-use manner, as detailed below.RGESolver will be also integrated in HEPfit [3], a flexible open-source tool which, given the Standard Model or any of its extensions, allows to fit the model parameters to a given set of experimental observables and to obtain predictions for observables.
This paper is organized as follows: we briefly introduce the theoretical framework of the SMEFT in Section 2, presenting the notation used in RGESolver.In Section 3 we describe in more detail the structure of the library, with a specific focus on the handling of the flavour structure and on the implementation of the numerical solution of the RGE's.Section 4 is devoted to describing the usage of RGESolver: we discuss the installation procedure and we describe how to use the basic functionalities of the library, together with a list of its most important methods.We discuss the efficiency of the library and the comparison with DSixTools ( [4,5]) in Section 5. We present our conclusions in Section 6.

Theoretical framework: the SMEFT
As discussed above, the absence of new particles at energies at or above the TeV scale allows us to parametrize the effects of physics beyond the SM with a tower of dimension D > 4, Lorentz and gauge-invariant operators [1,2].The resulting effective field theory is the SMEFT.The SM gauge group is (2.1) The three factors represent color, weak isospin and hypercharge gauge group.The quantum numbers of SM fields under G SM are listed in Table 1.The Lagrangian of the Standard Model is (following the notation used in Refs.[6][7][8]) With this notation, when spontaneous symmetry breaking occurs, the physical Higgs boson acquires a mass given by m2 H = 2λv 2 , with v = ( √ 2G F ) −1/2 ∼ 246 GeV.The sum over ψ is extended to all fermion fields in the SM, both left-handed (q, l) and right-handed (u, d, e).The gauge covariant derivative is D µ = ∂ µ +ig 1 yB µ +ig 2 t I W I µ +ig 3 T A G A µ , with T A and t I = τ I /2 1 the generators in the fundamental representation of SU(3) and SU(2) respectively: Finally, H and X (with X = B, W I , G A ) are defined via the totally antisymmetric Levi-Civita symbol, with ε 12 = 1 and ε 0123 = 1: Fermions appear in n g = 3 different generations, so every fermion field carries, together with the gauge indices, a generation (or flavour) index that runs from 1 to 3. We did not write explicitly neither of these indices in (2.2) for the sake of simplicity.Yukawa couplings Y ψ are thus complex matrices in flavour space.
Table 2: The 59 independent dimension-six operators built from SM fields which conserve baryon and lepton number.At dimension six a complete basis of independent and gauge invariant operators that conserve B and L is given by the so-called Warsaw basis, defined in [2] and displayed in Table 2. Their independence means that no linear combination of them and their Hermitian conjugates vanishes due to the equations of motion (up to total derivatives).The equations of motion are used at O 1/Λ 2 level, so they can be derived from L SM alone.The Warsaw basis consists of 59 operators, for a total of 2499 independent parameters in the generic flavour scenario (see Appendix A of Ref. [8]).The operators are divided in 8 classes, depending on their field content, which schematically are where X stands for a gauge field-strength tensor (or its dual), ψ for a fermion field, D for a derivative and H is the Higgs field.
To conclude the discussion about the SMEFT we note that higher dimensional operators contribute to observables not only through their matrix element, but also through their contribution to the running of SM couplings.An example is given by the diagrams in Figure 1.Coherently with the power counting, their effects in the running of λ are suppressed by a factor of m 2 H /Λ 2 .The contributions to the running of the SM parameters due to SMEFT operators are available in Ref. [6].

Flavour basis and back-rotation
In the SM a unitary transformation for each fermion field in flavour space, does not alter the structure of the Lagrangian, provided the Yukawa matrices are redefined as (2.7) This freedom allows to diagonalize simultaneously two out of the three matrices (Y e and one between Y u and Y d ). 2 Denoting with Ŷ a diagonal Yukawa matrix we can define two reference bases: the down basis where Y e = Ŷe , Y d = Ŷd , Y u = Ŷu V and the up basis where Y e = Ŷe , Y u = Ŷu , Y d = Ŷd V † , with V a unitary matrix.In the SM, V is the Cabibbo-Kobayashi-Maskawa matrix.In the SMEFT this is not true, since the mass matrices are not simply proportional to the Yukawa couplings due to the contributions from dimension-six operators.With a slight abuse of notation, in the following we still refer to the matrix V as the CKM matrix.It is worth noticing that the matrix V is not the one appearing in the interaction vertices between quarks and W ± bosons, again due to effects from higher dimensional operators.A more detailed discussion is available in Section 2 of Ref. [9].
Notice that the up and down bases are not stable under renormalization group evolution, even in the SM.For example, the SM β function for Y d (available in Ref. [10]) , which is always non-diagonal in flavour space.Starting from a diagonal Y d at the scale µ I (down basis) leads to a non-diagonal Y d at a different scale µ F .To go back into the down basis, a further flavour rotation is thus required.Clearly, the same holds for the up basis.In the SMEFT the flavour transformations in Equation (2.6) imply not only a redefinition of the Yukawa couplings, but also a rotation of the Wilson coefficients.For example, the coefficient It should be noticed that the flavour rotations are not uniquely determined by the singular value decomposition.The three Yukawa matrices are diagonalized via six unitary matrices: If U ψ , V ψ satisfy the relation (2.8), also U ψ φ ψ , V ψ φ ψ with φ ψ = diag e α ψ 1 , e α ψ 2 , e α ψ 3 satisfy it.The phase matrices φ u , φ d are determined unambiguosly once the phase convention for the CKM matrix is chosen, while, in the absence of right-handed neutrinos, the phase of the lepton matrices cannot be determined.A possibility to fix the phases would be using one of the SMEFT coefficients that involve leptons, but this would not be a general solution.In fact, the user may want to set that coefficient to 0, frustrating the choice of phase convention.To overcome this problem, we choose φ e = diag e −i arg[(Ue) 11 ] , e −i arg[(Ue) 22 ] , e −i arg[(Ue) 33 ] .This choice ensures that an evolution between two energy scales close to each other produces a small change in the coefficients.The user is then free to perform the rotation in the lepton sector choosing any other convention.To go into the up basis, the rotation matrix R q in (2.6) must be set equal to V u φ u .To go into the down basis one must instead set R q = V d φ d .Obviously, in both bases one sets

Solution of the renormalization group equations
The RGE for a parameter x i , that in this context is a SMEFT coefficient C i or a SM parameter, is that defines the β-function for the parameter.In general, when the theory has multiple parameters, the β-function for x i can contain terms proportional to integer powers of x i , but also terms proportional to other parameters x j , i = j.In this article we always consider the β-functions at one-loop level.
The β-functions of SMEFT coefficients (as well as SMEFT contributions to the running of SM parameters) are available in Refs.[6][7][8], the SM β-functions of gauge couplings can be found for example in Ref. [11] and the β-functions of Yukawa couplings, m 2 H and λ in Ref. [10]. 3he β-function for SMEFT coefficients must be linear in the coefficients themselves by power counting.In fact, the product of two dimension-six SMEFT coefficients would be suppressed by a factor of 1/Λ 4 , negligible at O m 2 H /Λ 2 .This allows to write the RGE's for the C i s introducing the anomalous dimension matrix Γ (ADM), (2.10) The indices i, j in the previous equation run on the whole set of independent coefficients: Γ is a square 2499 × 2499 matrix.Clearly, the β-function for the i-th coefficient is A simple approximation consists in neglecting the µ dependence of the right-hand side of Equation (2.10) and taking only linear order terms in ln (µ F /µ I ), yielding (2.11) While this drastic simplification may be appropriate for cases in which the two energy scales are not too different, so that the second term in Equation (2.11) is a small correction, in general the log on the right-hand side will become large, calling for a resummation of log-enhanced terms.
To provide a more accurate solution of the system of first-order differential equations in Equation (2.10) is a non-trivial task since the anomalous dimension matrix depends on µ through the SM couplings, the running of which is in turn influenced by SMEFT operators.The anomalous dimension matrix can be written as where the matrices on the right-hand side do not depend on the renormalization scale.
In general the Γ (g i ) matrices cannot be diagonalized simultaneously: an exact analytic solution of the system is then impossible.One could proceed with a hybrid scheme as commonly done in the Effective Theory for Weak interactions (WET) below the EW scale, where QCD and QED corrections are implemented by resumming the logs generated by the strong coupling only.However, in the SMEFT the top Yukawa coupling is also large, calling for the resummation of both strong and Yukawa contributions.Thus, a numerical approach is required to obtain precise results in the general case.4More details about the implementation of this method are given in Section 3.2.

Description of the library
RGESolver is implemented in C++.We use the GNU Scientific Library for the integration of the RGE's (more information about the numerical methods used are given below).The whole code (including input & output) handles separately real and imaginary parts of the parameters.

Flavour structure
Operators (or equivalently their coefficients) in the SMEFT can be divided in different categories according to their symmetry properties as in Table 3, 5 where we partially modify the notation used in Table 14  We work splitting real and imaginary part for each C i .For example, the real part of the coefficient of an Hermitian 2F operator is a (real) symmetric matrix, while the imaginary part is a (real) antisymmetric matrix.Also the real and imaginary part of symmetric 4F operators are in different categories.This is not the case for non-Hermitian operators: the real and imaginary parts of their coefficient do not have any constraints, and they belong to categories 1 (2F) and 5 (4F).Clearly it is convenient to store and evolve only the independent parameters, but for ease of use it is convenient to be able to access the coefficients with any combination of indices, not only the independent ones.We achieve this flexibility by defining specific symmetry-aware getters and setters for each category.For example, setting Re[C 1,2  Hl(1) ] = 0.5 via S.SetCoefficient("CHl1R",0.5,1,2) (where S is an instance of the RGESolver class) will lead to S.GetCoefficient("CHl1R",2,1) returning the value 0.5 due to Hermiticity.

Numerical implementation of renormalization group evolution
As already mentioned, we use the routines provided by the GNU Scientific Library (GSL , [13]) to perform the integration of the RGEs.Let y be the vector that contains all the independent parameters of the SMEFT.Introducing t = ln(µ/GeV), the system in Equation (2.9) can be rewritten as where the right-hand side of the equation depends on t only through y.When all the symmetries of the operators are considered, the coefficients of operators at dimension-six level are completely determined through 2499 independent real parameters.The three gauge couplings g 1 , g 2 , g 3 , the Higgs quartic coupling λ and the mass m H of the Higgs boson are real scalars.Yukawa couplings Y u , Y d and Y e are 3 × 3 complex matrices, each of which has 18 parameters. 6This yelds 59 further parameters, raising the total number to 2558: this is the dimension of y.The system in Equation (2.9) is thus a 2558-dimensional system of first order coupled differential equations.We use an adaptive stepsize routine to solve the system: the stepsize h is changed throughout the integration in order to match the estimated local error E i with a user-defined error level, to obtain the maximum efficiency.The desired error level D i for each component is determined through four parameters a y , a dydt , ǫ abs and ǫ rel with the expression The second term in brackets in the previous expression ensures a correct estimation also for situations where some components y i are very close to zero.The error is determined not only from the value y i but also from its increment.In particular, we use a y =a dydt = 1.
E i can be estimated with the step-doubling technique, that consists in advancing the solution from t to t + 2h in two different ways (performing two steps of length h or one step of length 2h) and taking the error as the difference between the two.A more sofisticated and efficient error estimate is possible for the embedded integration methods, where the same evaluations of the function f are used to compute two different values of the solution at t + h.The error is taken as the difference between the two values (see Section 16.2 of Ref. [14]).
If E i exceeds the desired error D i by more than 10% for any component, the stepsize is reduced according to where 0.9 is a safety factor (the error can only be estimated, not accurately determined), q is the consistency order of the method (i.e. the local error scales as h q+1 ) and E/D = max i (E i /D i ) is the maximum ratio of estimated and desired error among the components.If, instead, the estimated error is lower than the desired one (precisely, when E/D < 50%) the stepsize is increased according to To avoid uncontrolled changes in the stepsize, the overall scaling factor is limited to the range (1/5, 5).In RGESolver the explicit embedded Runge-Kutta-Fehlberg method is used (with initial stepsize h = ln(µ/Λ)/100), for which q = 4.
It should be stressed that Equation (3.2) refers to a local error: there is no simple relation that connects the four parameters to an estimate of the global error affecting the final result.We tested the accuracy of the adaptive stepsize integration comparing it with a fixed-step integration with an high number of steps.Using µ I = Λ = 1000 TeV and µ F = 250 GeV as energy scales, ǫ abs = 10 −16 and ǫ rel = 10 −4 for the adaptive integration and N FS = 1000 steps for the fixed stepsize integration we obtained percentual differences 10 −5 .The initial conditions for SM parameters at the high energy scale are obtained evolving them from µ ∼ v to Λ with the pure SM β-functions (namely neglecting SMEFT contributions).The initial conditions for SMEFT coefficients are O 1/Λ 2 , as prescribed by the power counting.

Installation
RGESolver is a free software released under the GNU General Public License.The download can be performed directly from the command line, typing in the terminal: git clone https://github.com/silvest/RGESolver--recursive More details can be found on the dedicated GitHub webpage.The extended documentation is also available here.

Class methods
Here we briefly describe the most important methods available in RGESolver.The detailed documentation is available on GitHub.

Evolution
• void Evolve(std::string method, double muI, double muF): performs the oneloop renormalization group evolution from µ = muI to µ = muF muI and muF are expressed in GeV).The current values of the SMEFT parameters are taken as initial condition at µ = muI.After the evolution, they are updated with the new values at µ = muF.The available methods for the solution of the RGE's are "Approximate" and "Numeric".The first method is faster than the latter, but less accurate, as explained in Section 2.2.
• void EvolveSMOnly(std::string method, double muI, double muF): same as Evolve , but only for the SM parameters.The user should use this method instead of Evolve when interested in pure SM running.Using this function is the same of using Evolve with all the SMEFT coefficients set to 0, but it is faster since it computes only the evolution for the SM parameters.
• void EvolveToBasis(std::string method, double muI, double muF, std::string basis): same as Evolve, but performs in addition the flavour back-rotation described in Section 2.1 to go into the selected basis ("UP" or "DOWN").After the evolution, the CKM matrix is computed.
• void GenerateSMInitialConditions(...): generates the initial conditions for Standard Model parameters (g 1 , g 2 , g 3 , λ, m 2 h , Y u , Y d , Y e ) at the scale muFin (in GeV), using one-loop pure SM beta functions.At such scale, the CKM matrix is computed.The generation can be done starting from user-defined low-energy input or using the default values summarized in Table 4 (the arguments to be passed to the function are different in the two cases, see the documentation for the details and Subsection 4.3 for an example).In case of user-defined input, this method should be used with usual fermion hierarchy (smallest mass for the 1st generation and greatest mass for the 3rd without mass degeneracy for all up and down quarks and for charged leptons).

Input & Output
Getters and setters take as first argument a string (name), that identifies the corresponding parameter.The names that must be used to correctly invoke these functions are reported in Tables 5, 6, 7 and 8.The functions for non-scalar coefficients take as additional arguments the flavour indices (2 or 4), that must be in the interval [0:2].
• void SetCoefficient(std::string name, double val): setter function for the scalar parameters.Sets the parameter name equal to val. Coeff.Name • double GetCoefficient(std::string name): getter function for the scalar parameters.Returns the parameter name.
• void SetCoefficient(std::string name, double val, int i, int j): setter function for the parameters with two flavour indices.Sets the parameter name[i,j] equal to val.
• double GetCoefficient(std::string name, int i, int j): getter function for the parameters with two flavour indices.Returns the parameter name[i,j].
• void SetCoefficient(std::string name, double val, int i, int j, int k, int l): setter function for the parameters with four flavour indices.Sets the parameter name[i,j,k,l] equal to val.
• double GetCoefficient(std::string name, int i, int j, int k, int l): getter function for the parameters with four flavour indices.Returns the parameter name[i,j,k,l].
• double GetCKMAngle(std::string name), double GetCKMPhase(): getter methods to access matrix parameters.These methods should be called only after methods that choose a specific flavour basis (as GenerateSMInitialConditions() or EvolveToBasis()), in which case the CKM matrix is updated.
• double GetCKMRealPart(int i, int j), double GetCKMImagPart(int i, int j): getter functions for the real and imaginary part of the (i, j) element of the CKM matrix.These methods should be called only after methods that choose a specific flavour basis (as GenerateSMInitialConditions() or EvolveToBasis()), in which case the CKM matrix is updated.The indices must be in the range [0:2].
• void SaveOutputFile(std::string filename, std::string format): saves the current values of the SMEFT parameters in the file filename.Currently, the only available format is Susy Les Houches Accord "SLHA" ( [15]).

Numerical integration parameters
Here we list the methods related to the parameters that control the numerical evolution, as described in Section 3.2: • void Setepsrel(double epsrel) : sets ǫ rel equal to epsrel (the default value is ǫ rel = 5 × 10 −3 ).
• Reset(): resets all the SMEFT coefficients to 0 and the SM parameters to their default value (in the up basis).ǫ abs and ǫ rel are reset to their default value.This function should be called after the evolution, if the same instance of the class is used for another run. S.EvolveToBasis("Numeric",Lambda,muLow,"DOWN"); Since this method computes also the CKM matrix, we can access its parameters or its elements using the dedicated getter functions.For example, the following lines print onscreen sin θ 12 at the scale muLow: std::cout<<"Sin(theta12)("<<muLow<< " GeV): "<< S.GetCKMAngle("s12")<< std::endl;

Execution times and comparison with DSixTools
We tested our code against DSixTools ([4, 5]), a Mathematica package for the handling of the SMEFT (and the Low-energy Effective Field Theory, but we focus on the SMEFT part of the package).DSixTools solves the SMEFT RGE's numerically and via the approximation described in Section 2.2, as RGESolver does.We performed an evolution from Λ to µ Low = 250 GeV with both codes to compare the values of the coefficients after the evolution.The initial conditions at the high energy scale Λ for the SM parameters are obtained using pure SM RGE's to run them up to such scale.We set O 1/Λ 2 initial conditions for every SMEFT coefficient at µ = Λ. 8The result of this comparison is displayed in Table 9.
To compare the two results, we introduce ∆ max = max|C R i (µ Low ) − C D i (µ Low )|/(1/Λ 2 ), with C R(D) i (µ Low ) the evolved coefficient computed by RGESolver (DSixTools).We observe that the two codes produce the same results up to O 10 −5 /Λ 2 (O 10 −6 /Λ 2 ) in the case of numeric solution (approximate solution).Clearly, the agreement is better for the simpler solution, namely the approximate one.Let us briefly discuss RGESolver's efficiency performances.We report in Table 10 the execution times that we measured (referring to a PC whose specifications are shown in Table 11).The execution times consider not only the evolution, but also the input/output (setting the initial conditions for the SMEFT parameters and accessing them after the evolution).These times should not be taken as exact, but as an order-of-magnitude estimate.This result should be compared to the execution times of DSixTools, that can be up to O (20 s) (O (10 s)) for the numeric (approximate) solution.This package performs a consistency check on the input that affects heavily the execution time for large inputs.Such impact is ∼ 8 s for the input used in the comparison between the two programs.Conversely, considering just 3 non vanishing independent operators at µ = Λ, the execution time is O (10 s) (O (2 s)) for the numeric (approximate) solution.

Conclusions
RGESolver is an open-source C++ library that performs the renormalization group evolution in the context of the SMEFT at dimension-six level, in the most general flavour scenario.Only operators that preserve baryon and lepton number are considered.RGESolver was designed to be easy to use, highly efficient and suitable to perform extensive phenomenological analysis.To this aim, it will be included in HEPfit ( [3]), a multipurpose and flexible analysis framework that can be used for fitting models to experimental and theoretical constraints.RGESolver outperforms DSixTools in execution time by two orders of magnitude, while retaining an excellent accuracy.Further details and the full documentation can be found on GitHub.
p, r, s, t are fermion flavour indices, j, k are indices in the fundamental representation of SU(2) W , I, J, K (A, B, C) are indices in the adjoint representation of SU(2) W (SU(3) C ) and greek letters (µ, ν . . . ) are Lorentz indices.Contraction of indices in the fundamental representation of SU(3) C is implicit.
SMEFT diagram contributing to the β function of λ with terms O m 2 H C H . β function of λ with terms O λm 2 H C HD ,O λm 2 H C H .

Figure 1 :
Figure 1: Two SMEFT diagrams that contribute to the running of the Higgs quartic coupling λ.The solid square denotes a SMEFT vertex and the dot denotes a SM vertex.

Table 9 :
Comparison between DSixTools and RGESolver: the table on the left refers to the approximate solution, the one on the right to the numeric solution.

Table 3 :
[5]Ref.[5].Symmetry categories for operators in the SMEFT.nF indicates the number of flavour indices for each category.

Table 4 :
SM parameters used by default to generate SM initial conditions at an arbitrary scale.The scale at which these parameters are given is µ = 173.6GeV.We follow UTfit for what concerns the conventions for the CKM matrix.

Table 5 :
Standard Model parameters inRGESolver.The labels in the left column must be used with the GetCoefficient method.The ones in the right column must be used with the GetCKMAngle method.

Table 10 :
Execution times of C++ programs that use RGESolver.

Table 11 :
Technical specifications of the computer used for the test.