RGErun 2.2.3

Author: Kristjan Kannike. Created: 2004-08-07. Changed: 2017-12-13.
Available under the GNU General Public Licence <http://www.gnu.org/copyleft/gpl.html>.
General framework for running renormalisation group equations.

The RGErun Package

BeginPackage["RGErun`"];

Preamble

β::usage="β[x,model,thresholds] is used to enter the beta function for the variable x in the model; the beta function may depend the thresholds below the current energy scale if different effective field theories are used; if not, you can set it to _.";

run::usage="run[thresholds,variables,direction,model] runs the RGEs between thresholds in the direction (up or down), gives variables as functions of the energy scale and optionally matches variables at thresholds or prints out the values of the variables and their predefined functions at the thresholds together with a plot of running couplings.";

printOut::usage="printOut is an option for run which specifies whether the results and matched results at each scale should be printed out together with a plot of the running variables; the default value is printOut -> True";

matches::usage="matches is an option for run which specifies the matching conditions for the running variables at thresholds.";

functions::usage="functions is an option for run which specifies functions of the running variables (e.g. eigenvalues of a mass matrix) that should be printed out at thresholds.";

𝕀::usage="𝕀 must be used in lieu of IdentityMatrix[n] in the beta functions.";

CirclePlus::usage="Use ⊕ instead of + to add a constant vector or matrix in a beta function.";

CircleTimes::usage="Use ⊗ instead of * to multiply element-wise by constant vector or matrix in a beta function.";

getRGEs::usage="getRGEs[name] imports the β-functions defined in the output of SARAH (name is the directory of SARAH output) or PyR@TE (name is the file with PyR@TE Mathematica output).";

returnVariables::usage="returnVariables is an option for getRGEs that specifies whether to return a list of variables defined in a PyR@TE Mathematica output.";
replaceNames::usage="replaceNames is an option for getRGEs to replace variable names defined in the PyR@TE output.";

updateInstalled::usage="updateInstalled[directory] replaces the RGErun.m package in the user's Mathematica Applications folder with the one in the directory (the default is NotebookDirectory[]).";

Options[getRGEs]={returnVariables->False,replaceNames->{}};
getRGEs[name_:"",model_:"",OptionsPattern[]]:=Block[{variableSubs,AllRGEs,DimensionParameter,conj,Tp,Adj,trace,MatMul,pi},
Which[name==""&&And@@Map[Head[#]===List&,{BetaGauge,BetaLijkl,BetaYijk,BetaTijk,BetaMuij,BetaBij,BetaVEV}],(* SARAH definitions in the Global context *)
AllRGEs=
Join[BetaGauge,BetaLijkl,BetaYijk,BetaTijk,BetaMuij,BetaBij,BetaVEV],
DirectoryQ[name],(* Directory with SARAH output *)
AllRGEs=Flatten[Get[FileNameJoin[{name,#}]]&/@{"BetaGauge.m","BetaLijkl.m","BetaYijk.m","BetaTijk.m","BetaMuij.m","BetaBij.m","BetaVEV.m"},1],
True,(* PyR@TE output file *)
AllRGEs=Get[name]
];
If[Last[Dimensions[AllRGEs]]==3,(* Convert from SARAH format *){First[#],{(1/(16 (π^2))),((((1/(16 (π^2)))))^2)}.Rest[#]}&/@AllRGEs];
variableSubs=Map[#->Head[#]&,Select[AllRGEs[[All,1]],!AtomQ[#]&]];Scan[(β[First[#],model,_]:=Last[#])&,AllRGEs/.pi->π/.variableSubs/.OptionValue[replaceNames]//.{conj->Conjugate,Tp->Transpose,Adj->ConjugateTranspose,trace[ms__]:>Tr[Dot[ms]],MatMul[ms__][__]:>Dot[ms]}];
If[OptionValue[returnVariables]==True,AllRGEs[[All,1]]/.variableSubs/.OptionValue[replaceNames]]
];

updateInstalled[directory_:NotebookDirectory[]]:=Block[{newVersion=FileNameJoin[{directory,"RGErun.m"}],
installedVersion=FileNameJoin[{$UserBaseDirectory,"Applications","RGErun.m"}]},
If[FileExistsQ[newVersion],
If[FileExistsQ[installedVersion],DeleteFile[installedVersion]];
CopyFile[newVersion,installedVersion];
Print[StringForm["Copied `1` to `2`.",newVersion,installedVersion]],
Print["No RGErun.m in ",directory]];
];

Begin["`Private`"];

Auxiliary Constants & Functions

SetAttributes[{𝕀},Protected];
SetAttributes[{CirclePlus,CircleTimes},{HoldAll,Protected}];
chopδ=0 (10^(-20));
runδ=(10^(-30));

numericQ::usage="numericQ[expr] returns True if expr is a numeric quantity or a full array of numeric quantities and False otherwise.";
numericQ[x_]:=NumericQ[x]||ArrayQ[x,_,NumericQ];

flatten::usage="flatten[x] flattens x if x is a list, else returns {x}.";
flatten[xs_]:=If[Dimensions[N[xs]]=={},{xs},Flatten[xs]];

threadEqual::usage="threadEqual[xs,ys] returns equations formed from flatten'ed xs and ys";
threadEqual[xs_,ys_]:=Thread[flatten[xs]==flatten[ys]];

numericsNamesQ::usage="numericsNamesQ[variables,scale] returns the list of variables that are numeric at the scale.";
numericAtScaleQ[variables_,scale_]:=Select[variables,numericQ[#[scale]]&];

complexInterpolationFunctionQ::usage="complexInterpolationFunctionQ[if] returns True if the InterpolationFunction is complex anywhere and False otherwise.";
complexInterpolationFunctionQ[if_]:=Cases[Chop[Im[if[[4,3]]]],Except[0|0.]]≠{};

split::usage="Split a list into parts such that the length of each part is given by the second argument.";
split[list_,parts_]:=Map[Take[list,#1+{1,0}]&,Partition[FoldList[Plus,0,parts],2,1]];
(* Taken from Programming With Mathematica, An Introduction by Paul Wellin *)

printSeparator::usage="printSeparator[n] prints a row of n —'s";
printSeparator[n_:50,c_:"—"]:=Print[StringJoin[ConstantArray[c,n]]];

style::usage="style[text,name] sets the text to a predefined style.";
style[text_,"header"]:=Style[text,FontSize->16];

display::usage="display[variable] returns the variable in MatrixForm if the variable is a vector or matrix, otherwise it returns the variable";
display[variable_]:=If[ArrayQ[variable],MatrixForm[Chop[N[variable],chopδ]],Chop[N[variable],chopδ]];

printIntermediateResults::usage="printIntermediateResults[thresholds, scale, variables, functions,done] prints out intermediate results in runup if printOut->True; done specifies whether the variables are \"Run\" etc.";
printIntermediateResults[thresholds_,scale_,variables_,functions_,done_:"Run",unit_:"GeV"]:=Block[{},
Print["Results @ "<>StringJoin[Riffle[thresholds," = "]]<>" = ",N[scale]," ",unit];
Print[done<>" variables:"];
Print[Sequence@@Flatten[Riffle[{ToString[#]," = ",display[#[scale]]}&/@variables,
", "]]];
If[functions≠{},
Print[done<>" functions:"];
Print[Sequence@@Flatten[Riffle[{ToString[#]," = ",display[#[scale]]}&/@functions,
", "]]];
];
];

checkInitialConditions::usage="checkInitialConditions[variables] checks whether all the initial conditions to run are defined and numerical at the scale.";
checkInitialConditions[variables_,scale_]:=With[{numerics=numericAtScaleQ[variables,scale]},If[numerics≠variables,Message[run::nonnum,ToString/@Complement[variables,numerics]]];
Complement[variables,numerics]];

groupThresholds::usage="groupThresholds[thresholds] takes as arguments the list of \"threshold name\" -> scale rules and returns the {sorted names, sorted mass scales} list.";
groupThresholds[thresholds_]:=Block[{thrs=GatherBy[SortBy[thresholds,Last],Last]},
{Map[First,thrs,{2}],Map[Last@First[#]&,thrs]}
];

ts::usage="ts[massScales] returns logs of ratios of mass levels from a vector.";
ts[massScales_]:={Log[#1/#2],runδ}&@@@Partition[massScales,2,1];(* The value (10^(-30)) instead of 0 as an upper bound avoids an NDSolve bug in v. 10.1 *)

runBetweenThresholds::usage="runBetweenThresholds[scale,tMin,tMax,thresholds,variables], where the initial values for the variables must be given at scale, solves the equations in the interval from tMin to tMax.";
runBetweenThresholds[scale_,tMin_,tMax_,model_,thresholds_,variables_,opts:OptionsPattern[]]:=Block[{eqs,vars,sub,t,sol,event,end},
sub=
ReleaseHold[With[{name=Symbol[ToString[#]<>"$"]},Hold[Sequence[#'[t]->Array[name[##]'[t]&,Dimensions[N[#[scale]]]],#->
Array[name[##][t]&,Dimensions[N[#[scale]]]]]]]&/@variables];
eqs=Flatten[ReleaseHold[Hold@Sequence[
threadEqual[#'[t]/.sub/.𝕀:>IdentityMatrix[First[Dimensions[N[#[scale]]]]],β[#,model,thresholds]/.sub/.𝕀:>IdentityMatrix[First[Dimensions[N[#[scale]]]]]/.{CirclePlus->Plus,CircleTimes->Times}],threadEqual[#/.sub/.t->tMin,#[scale]]]&/@variables]];
vars=Flatten[With[{name=Symbol[ToString[#]<>"$"]},Array[name[##]&,Dimensions[N[#[scale]]]]]&/@variables];
event[x_]:=Apply[And,Abs[#[x]]<4 π&/@vars];
sol=First[NDSolve[eqs,vars,{t,tMin,tMax},Method->{"EventLocator","Event"->event[t]},FilterRules[{opts},Options[NDSolve]]]];
end =sol[[1,2,1,1,-1]];(* End of the domain of the interpolation function; running up! *)
{sol,end}
];

pieceTogether::usage="pieceTogether[variables,solutions,scales] pieces together each variable and its Derivatives as a function of the renormalisation scale from the solutions to the RGEs.";
pieceTogether[variables_,solutions_,scales_,model_]:=(
Scan[
(#[μ_?NumericQ,model]:=Block[{intervals=Interval/@Partition[scales,2,1]},
Piecewise[
With[{var=#,varS=Symbol[ToString[#]<>"$"]},
MapThread[{
 Hold[(* To avoid evaluation of functions in wrong intervals *)
With[{scale=#2[[1,2]]},
Array[
varS[##][Log[μ]-Log[scale]]&,
Dimensions[N[var[scales[[1]]]]]
]
]/.#1
],
IntervalMemberQ[#2,μ]
}&,
{solutions,intervals}]],
Which[μ<scales[[1]],#[scales[[1]]],μ>scales[[-1]],#[scales[[-1]]]]]//ReleaseHold
];
#[μ_?NumericQ]:=Evaluate[#][μ,model];
)&,
variables];
Scan[
(Derivative[n_][#][μ_?NumericQ,model]:=Block[{intervals=Interval/@Partition[scales,2,1]},
Piecewise[
With[{var=#,varS=Symbol[ToString[#]<>"$"]},
MapThread[{
 Hold[(* To avoid evaluation of functions in wrong intervals *)
With[{scale=#2[[1,2]]},
Array[
(1/(μ^n)) StirlingS1[n,Range[n]].With[{is=##},Derivative[#][varS[is]][Log[μ]-Log[scale]]&/@Range[n]]&,
Dimensions[N[var[scales[[1]]]]]
]
]/.#1
],
IntervalMemberQ[#2,μ]
}&,
{solutions,intervals}]],
Array[0&,Dimensions[N[#[scales[[1]]]]]]]//ReleaseHold
];
Derivative[n_][#][μ_?NumericQ]:=Derivative[n][Evaluate[#]][μ,model];
)&,
variables];
);

Run & Display

plotStyles={Black,ColorData["Legacy"]["LightCadmiumRed"],ColorData["Legacy"]["PermanentGreen"],ColorData["Crayola"]["Blue"],ColorData["Legacy"]["Apricot"],ColorData["Legacy"]["BlueViolet"],ColorData["Legacy"]["Gold"],Black};

(*optimalPosition::usage="optimalPosition[variables,minX,maxX] gives the rough position between minX and maxX where the variables are most distant from each other.";optimalPosition[variables_,areComplexes_,minX_,maxX_]:=Block[{dist},
dist[list_]:=Total[Abs[Subtract@@#]&/@Subsets[list,{2}]];
SortBy[Block[{x=#},{x,dist[Flatten@MapThread[If[And@@#2==False,#1[x],ReIm[#1[x]]]&,{variables,areComplexes}]]}]&/@Rest[Most[Exp[Subdivide[Log[minX],Log[maxX],10]]]],Last][[-1,1]]
]*)

SetAttributes[{printOut},Protected];
Options[run]={printOut->True,matches->{},functions->{}}~Join~Options[NDSolve];
run::nonnum="The initial condition(s) for `1` is/are not numerical and/or not full arrays.";
run[thresholds_,variables_,direction_:"up",model_:"",opts:OptionsPattern[]]:=Block[{dir=Which[direction=="up",1,direction=="down",-1],thrs,scales,intervals,t,tEnd,thr,runthr,matchedHere,runBetween,piecewiseVariables,areComplexes,normalThresholds=Normal[thresholds](*,position*)},
{thrs,scales}=groupThresholds[normalThresholds];
intervals=Partition[scales,2,1];
intervals=Which[direction=="up",intervals,direction=="down",Reverse[intervals]];
Evaluate[thr/@scales]=thrs;
Evaluate[runthr/@scales]=FoldList[Join,First[thrs],Rest[thrs]];
(* Since there is no running at the first scale, it is undefined *)
Evaluate[t/@Which[direction=="up",Rest[scales],direction=="down",Most[scales]]]=ts[scales];
(* Run in each interval and print out results *)
runBetween=Block[{tMin,tMax,int=#2,runVariables},
If[checkInitialConditions[variables,int[[dir]]]=={},
{tMin,tMax}={t[int[[-dir]]][[dir]],t[int[[-dir]]][[-dir]]} ;
{runVariables,tEnd}=runBetweenThresholds[int[[dir]],tMin,tMax,model,runthr[int[[dir]]],variables,FilterRules[{opts},Options[NDSolve]]];

(*******************)
Print["tEnd = ",tEnd];
If[tEnd≠runδ,Throw[tEnd]];


Scan[(#[int[[-dir]]]=With[{name=Symbol[ToString[#]<>"$"]},Array[name[##][tMax]&,Dimensions[N[#[scales[[dir]]]]]]]/.runVariables)&,variables];
If[OptionValue[printOut]==True,printSeparator[];
printIntermediateResults[thr[int[[-dir]]],int[[-dir]],variables,OptionValue[functions]]];
matchedHere={};

Scan[Module[{mtch=#[[2]][model,int[[-dir]],thr[int[[-dir]]],direction]},
If[numericQ[mtch],#[[1]][int[[-dir]]]=mtch;AppendTo[matchedHere,#[[1]]]]]&,OptionValue[matches]];


If[OptionValue[printOut]==True&&matchedHere≠{},printIntermediateResults[thr[int[[-dir]]],int[[-dir]],matchedHere,OptionValue[functions],"Matched"];
printSeparator[]];
Append[#1,runVariables]
]
]&;

(* Do not evaluate anything if some initial conditions are ill defined: *)
If[checkInitialConditions[variables,scales[[dir]]]=={},
If[OptionValue[printOut]==True,
printSeparator[];
Print["RUNNING "<>ToUpperCase[direction]<>" IN MODEL "<>model];
 ];
If[OptionValue[printOut]==True,printIntermediateResults[thr[scales[[dir]]],scales[[dir]],variables,OptionValue[functions],"Initial"]];

Catch[piecewiseVariables=Fold[runBetween,{},intervals]]//Print;

If[OptionValue[printOut]==True,
printSeparator[]
 ];
If[direction=="down",piecewiseVariables=Reverse[piecewiseVariables]];

pieceTogether[variables,piecewiseVariables,scales,model];

(* Which variables have complex values? *)
areComplexes=Map[Or@@Map[complexInterpolationFunctionQ[Last[#]]&,#]&,Transpose[piecewiseVariables]];
areComplexes=split[areComplexes,Times@@Dimensions[#[scales[[1]]]]&/@variables];


If[OptionValue[printOut]==True&&checkInitialConditions[variables,scales[[-dir]]]=={},
(* Show a plot with all the running variables *)

(* Calculate the optimal position(s) for the labels *)
(*positions=Log[optimalPosition[variables,areComplexes,scales[[1]],scales[[-1]]]];*)
(*positions=optimalPositions[fs_,{{xMin_,xMax_},{yMin_,yMax_}}];*)

Print[
Show[Reverse@Sort@Flatten@MapThread[
Block[{variable=#1,areComplex=#2,colour=#3,position=RandomReal[Log[scales[[{1,-1}]]]],dimensions=Dimensions[N[#1[scales[[1]]]]],indices},
indices=If[dimensions=={},{""},Flatten[Array[Style[StringJoin[Riffle[ToString/@{##},","]],Red]&,dimensions]]];
MapIndexed[
{
LogLinearPlot[
With[{vari=flatten[variable[μ]][[Sequence@@#2]]},If[#1,{Re[vari],Im[vari]},vari]],{μ,scales[[1]],scales[[-1]]},AxesLabel->{"μ/GeV",""},AxesOrigin->{scales[[1]],0},PlotStyle->colour,PlotRange->All,Method->{"AxesInFront"->False}],
Block[{index=#2},With[{vari=flatten[variable[(E^position)]][[Sequence@@index]],ind=indices[[Sequence@@index]]},
Graphics[Text[If[dimensions=={},First[#],Subscript[First[#],ind]],{position,Last[#]},Background->White]]&/@If[#1,{{StringForm["Re ``",variable],Re[vari]},{StringForm["Im ``",variable],Im[vari]}},{{ToString[variable],vari}}]]
]
}&,
areComplex]
]&,
{variables,areComplexes,Flatten[ConstantArray[plotStyles,Quotient[Length[variables],Length[plotStyles]]]]~Join~Take[plotStyles,Length[variables]~Mod~Length[plotStyles]]}
],
PlotRange->All]
];
printSeparator[]
 ];
];
];

End[];

EndPackage[];

Change Log

TO DO: Switching RGEs not only by switching scales, but also models (splice together running variables of different models). Or: switching models at scales: a model or an EFT is defined between two scales. Because usually the RGEs for different EFTs are calculated separately and splicing them together is extra effort/source of error.
Instead of thresholds, models??
Model can be valid from 0... m1, m1...m2, m2...∞. In case of ∞, run up to the Planck mass.

2017-12-13: 2.2.3: Bug fix: matching work with new “up” and “down”.
2017-07-31: 2.2.2: Added model as an optional argument to variables; if two models “1” and “2” have a λ, then λ[μ] will equal the last evaluated of λ[μ,"1"] and λ[μ,"2"].
2017-07-30: 2.2.1: 𝕀 does not need an explicit argument anymore: where it is used, the variable is assumed to be a square matrix so 𝕀 gets first of its dimensions.
2016-09-01: 2.2.0:  Changed up and down to strings “up” and “down”.
2016-04-12: 2.1.9: Below the lowest or the highest scale, any variable now has the constant value at that scale and its derivative is zero.
2015-06-04: 2.1.8: Updated usage messages for run and β.
2015-05-29: 2.1.7: Changed getRGEs to read in β-functions from the Mathematica output of the SARAH package.
2015-05-28: 2.1.6: Renamed getPyRATERGEs to getRGEs.
2015-05-19: 2.1.5: Added model as an argument to run, β, match and getRGEs.
2015-05-19: 2.1.4: Changed matches and functions from arguments to run to its options.
2015-05-19: 2.1.3: Changed the format of β-functions: instead of β[λ][t,thresholds] it is β[λ,model,thresholds].
2015-05-18: 2.1.2: Put the code that updates the package in the Mathematica Applications directory in the updateInstalled function.
2015-05-18: 2.1.1: Changed the name of getPyrateRGEs into getPyRATERGEs to conform with the official name of that package.
2015-05-18: 2.1.0: Avoid an NDSolve bug in Mathematica 10.1: set (Subscript[t, max])=(10^(-30)) instead of (Subscript[t, max])=0 in the ts function.
2015-05-15: 2.0.9: Thresholds can now be entered as an Association.
2015-04-08: 2.0.8: Updated to the look of Mathematica 10.
2014-07-27: 2.0.7: Bug fix: made plotting Re and Im parts more robust; added Method->{"AxesInFront"->False}.
2014-07-27: 2.0.6: Added getPyrateRGEs to read in and convert β-functions from the Mathematica output of the PyR@TE package.
2014-07-26: 2.0.5: Bug fix: converted the cells of the usage messages of 𝕀, CirclePlus and CircleTimes into initialization cells.
2014-06-18: 2.0.4: Plotting the real and imaginary part of complex couplings.
2014-06-18: 2.0.3: All order derivatives of running couplings.
2013-10-03: 2.0.2: Removed the the maxSteps and workPrecision options to run and made it accept the Options of NDSolve.
2013-10-02: 2.0.1:  Bug fix: changed everywhere Dimensions[…] to Dimensions[N[…]] to ensure that a scalar with an exact initial condition such as (π/2) is internally treated as a scalar.

Examples

You can copy the package RGErun.m from the current directory to your Mathematica Applications directory (evaluate ‘The RGErun Package’ section above first):

updateInstalled[]

Read in the package:

Needs["RGErun`"];

?run

Running of the Gauge Couplings, the Higgs λ and (Subscript[y, t]) in the Standard Model

The RGes of the gauge couplings and (Subscript[y, t]) are at 2-loop level. For simplicity, the RGE for λ is at 1-loop level, as it makes little numerical difference. Notice that run accepts all the options of NDSolve.

mt=173.34;mP=1.22 (10^19);mHi=(10^20) mP;
thresholds={"(Subscript[m, t])"->mt,"TeV"->1000,"(Subscript[m, P])"->mP};
variables={g1,g2,g3,yt,λ};
g1[mt]=(Sqrt[(5/3)]) 0.35940;
g2[mt]=0.64754;
g3[mt]=1.1666;
yt[mt]=0.95096;
λ[mt]=0.12879;

β[g1,"SM",_]:=(1/(16 (π^2)))(41/10)(g1^3)+((((1/(16 (π^2)))))^2) (g1^3) ((44/5) (g3^2)+(27/10) (g2^2)+(199/50) (g1^2)-(17/10) (yt^2));

β[g2,"SM",_]:=(1/(16 (π^2)))(-(19/6)) (g2^3)+((((1/(16 (π^2)))))^2) (g2^3)(12 (g3^2)+(35/6) (g2^2)+(9/10) (g1^2)-(3/2) (yt^2));

β[g3,"SM",_]:=(1/(16 (π^2)))(-7)(g3^3)+((((1/(16 (π^2)))))^2) (g3^3) (-26 (g3^2)+(9/2) (g2^2)+(11/10) (g1^2)-2 (yt^2));

β[yt,"SM",_]:=(1/(16 (π^2))) ((9/2) (yt^2)-(17/20) (g1^2)-(9/4) (g2^2)-8 (g3^2)) yt+((((1/(16 (π^2)))))^2)  yt ((yt^2) (-12 (yt^2)-12 λ+36 (g3^2)+(225/16) (g2^2)+(393/80) (g1^2))+6 (λ^2)-108(g3^4)-(23/4) (g2^4)+(1187/600) (g1^4)+9 (g2^2) (g3^2)+(19/15) (g1^2) (g3^2)-(9/20) (g1^2) (g2^2));

β[λ,"SM",_]:=(1/(16 (π^2)))  (((9 (g2^4))/8)+((3 (g2^2) (g1^2))/4)+((3 (g1^4))/8)-(9 (g2^2)+3 (g1^2)-12 (yt^2)) λ+24 (λ^2)-6 (yt^4));
run[thresholds,variables,"up","SM",printOut->True,MaxSteps->(10^4),WorkingPrecision->MachinePrecision];

Now variables have values at (Subscript[m, P]) and we can run down from high energy again.

run[thresholds,variables,"down","SM",printOut->True,MaxSteps->(10^4),WorkingPrecision->MachinePrecision];

After the running, the variables are available as functions of the renormalisation scale, e.g.

g1[1000]

λ[1000]

You can also take their derivatives with respect to the renormalisation scale:

λ'[1000]

λ''[1000]

Below the lowest scale and above the highest scale, the variables are set to the value at that scale and their derivatives are set to zero for definiteness:

λ[0]

λ[0]==λ[mt]

λ'[0]

The name of the model is an optional argument to the running variables:

λ[1000,"SM"]

g2'[2000,"SM"]

The running variables without the model name specified are given by the model that was evaluated last.

A variation of the running of the Standard Model couplings, with all the gauge couplings in one vector variable g. On the plot, vector and matrix indices are shown in red in labels. Note that for element-wise multiplication by constant vectors or matrices one has to use ⊗ (CircleTimes).

mt=173.34;mP=1.22 (10^19);
thresholds={"(Subscript[m, t])"->mt,"TeV"->1000,"(Subscript[m, P])"->mP};
variables={g,yt,λ};
g[mt]={(Sqrt[(5/3)]) 0.35940,0.64754,1.1666};
yt[mt]=0.95096;
λ[mt]=0.12879;

(* Use ⊗ for element-wise multiplication with constant vectors or matrices *)
β[g,"SM",_]:=(1/(16 (π^2))){(41/10), -(19/6), -7}⊗(g^3)+((((1/(16 (π^2)))))^2) (g^3)⊗((({{(199/50), (27/10), (44/5)}, {(9/10), (35/6), 12}, {(11/10), (9/2), -26}})).(g^2)-{(17/10),(3/2),2} (yt^2));

β[yt,"SM",_]:=(1/(16 (π^2))) ((9/2) (yt^2)-{(17/20) ,(9/4) ,8}.(g^2)) yt+((((1/(16 (π^2)))))^2)  yt ((yt^2) (-12 (yt^2)-12 λ+{(393/80),(225/16),36}.(g^2))+6 (λ^2)+(g^2).(({{(1187/600), -(9/40), (19/30)}, {-(9/40), (23/4), (9/2)}, {(19/30), (9/2), -108}})).(g^2));

β[λ,"SM",_]:=(1/(16 (π^2)))  ((g^2).(({{(3/8), (3/8), 0}, {(3/8), (9/8), 0}, {0, 0, 0}})).(g^2)-({3,9,0}.(g^2)-12 (yt^2)) λ+24 (λ^2)-6 (yt^4));

run[thresholds,variables,"up","SM"]

You can evaluate vector or matrix variables and their derivatives transparently as well:

g[(10^8)]

g[(10^8),"SM"]

g'[(10^8)]

Standard Model Extended by a Singlet: the Scalar Threshold Effect

Matching is demonstrated by the scalar threshold effect that stabilises the electroweak vacuum [arXiv:1203.0237]. The heavy singlet S is integrated out at its mass scale (Subscript[m, S]) and gives a contribution to the Higgs quartic coupling λH via matching. Note that the RGEs are different below and above the scale (Subscript[m, S]).

mt=173.34;mS=(10^8);mP=1.22 (10^19);
thresholds={"(Subscript[m, t])"->mt,"(Subscript[m, S])"->mS,"(Subscript[m, P])"->mP};
variables={g1,g2,g3,yt,λH,λHS,λS};
g1[mt]=(Sqrt[(5/3)]) 0.35940;
g2[mt]=0.64754;
g3[mt]=1.1666;
yt[mt]=0.95096;
λH[mt]=0.12879;
λHS[mt]=0.1;
λS[mt]=0.1;

β[g1,"SMS",_]:=(1/(16 (π^2)))(41/10)(g1^3)+((((1/(16 (π^2)))))^2) (g1^3) ((44/5) (g3^2)+(27/10) (g2^2)+(199/50) (g1^2)-(17/10) (yt^2));

β[g2,"SMS",_]:=(1/(16 (π^2)))(-(19/6)) (g2^3)+((((1/(16 (π^2)))))^2) (g2^3)(12 (g3^2)+(35/6) (g2^2)+(9/10) (g1^2)-(3/2) (yt^2));

β[g3,"SMS",_]:=(1/(16 (π^2)))(-7)(g3^3)+((((1/(16 (π^2)))))^2) (g3^3) (-26 (g3^2)+(9/2) (g2^2)+(11/10) (g1^2)-2 (yt^2));

β[yt,"SMS",_]:=(1/(16 (π^2))) ((9/2) (yt^2)-(17/20) (g1^2)-(9/4) (g2^2)-8 (g3^2)) yt+((((1/(16 (π^2)))))^2)  yt ((yt^2) (-12 (yt^2)-12 λH+36 (g3^2)+(225/16) (g2^2)+(393/80) (g1^2))+6 ((λH)^2)-108(g3^4)-(23/4) (g2^4)+(1187/600) (g1^4)+9 (g2^2) (g3^2)+(19/15) (g1^2) (g3^2)-(9/20) (g1^2) (g2^2));

β[λH,"SMS",thresholds_]:=(1/(16 (π^2)))  (((9 (g2^4))/8)+((3 (g2^2) (g1^2))/4)+((3 (g1^4))/8)-(9 (g2^2)+3 (g1^2)-12 (yt^2)) λH+24 ((λH)^2)-6 (yt^4)+If[MemberQ[thresholds,"(Subscript[m, S])"],1,0] 4 ((λHS)^2));

β[λHS,"SMS",thresholds_]:=(1/(16 (π^2))) If[MemberQ[thresholds,"(Subscript[m, S])"],1,0] ((1/2) (12 (yt^2)-(9/5) (g1^2)-9 (g2^2)) λHS+4 λHS (3 λH+2 λS)+8 ((λHS)^2));

β[λS,"SMS",thresholds_]:=(1/(16 (π^2)))  If[MemberQ[thresholds,"(Subscript[m, S])"],1,0] (8 ((λHS)^2)+20 ((λS)^2));

(* The matching functions take as arguments the scale of the threshold, the name(s) of the threshold and the direction of the running (up or down) *)
match[λH]["SMS",scale_,thresholds_,dir_]:=Which[
dir=="up"&&MemberQ[thresholds,"(Subscript[m, S])"],
λH[scale]+(((λHS[scale])^2)/(λS[scale])),
dir=="down"&&MemberQ[thresholds,"(Subscript[m, S])"],
λH[scale]-(((λHS[scale])^2)/(λS[scale]))
];

run[thresholds,variables,"up","SMS",matches->{λH->match[λH]}];

The Higgs quartic coupling λH receives a threshold correction that stabilises the electroweak vacuum. The couplings λS and λHS are defined and run only above the threshold (Subscript[m, S]).

Show[LogLinearPlot[{λH[μ],λS[μ],λHS[μ]},{μ,mt,mP},AxesOrigin->{mt,0},PlotStyle->{ColorData["Crayola"]["Red"],ColorData["Crayola"]["Green"],ColorData["Crayola"]["Blue"]},AxesLabel->{"μ/GeV"}],Graphics[Text[ToString[#],{Log[(10^11)],#[(10^11)]},Background->White]]&/@{λH,λS,λHS}]

Import β-Functions from SARAH and PyR@TE

SARAH is a Mathematica package for analysis of particle physics models that can calculate two-loop RGEs. PyR@TE is a Python package to calculate two-loop RGEs. With a few changes, the output of both packages can be used with run.

Import β-Functions From PyR@TE

For an example, download the RGEsOutput_numerics.m that contains the one-loop RGEs for the two Higgs doublet model as defined in PyR@TE into the same directory as the current notebook. Use the function getRGEs to import RGEs. The returnVariables option makes getRGEs to return the list of variables; the replaceNames option lets you to replace variable names.

Needs["RGErun`"]

SetDirectory[NotebookDirectory[]];

variables=getRGEs["RGEsOutput_numerics.m","2HDM",returnVariables->True,replaceNames->{gSU2L->g2,gSU3c->g3,lambda1->λ1,lambda2->λ2,lambda3->λ3,lambda4->λ4,lambda5->λ5,lambda6->λ6,lambda7->λ7}]

The RGE for λ6, for example, is

β[λ6,"2HDM",_]

Import β-Functions From SARAH

Since some variable names in SARAH clash with RGErun, you should calculate your RGEs first, exit the kernel, and then import them with getRGEs. For an example, download the RGEs for the singlet extended Standard Model or SSM. Unpack the ZIP file into the same directory as the current notebook. The returnVariables option makes getRGEs to return the list of variables; the replaceNames option lets you to replace variable names.

Needs["RGErun`"]

SetDirectory[NotebookDirectory[]];

getRGEs[FileNameJoin[{"SSM","RGEs"}],"SSM",returnVariables->True,replaceNames->{Lambda1->λ,LambdaS->λS,K1->μSH,K2->λSH}]

For example, the RGE for the singlet self-coupling λS is

β[λS,"SSM",_]