MSSMatNLO: MoGRe_v1.1.m

File MoGRe_v1.1.m, 29.7 KB (added by BenjF, 4 weeks ago)

MoGRe version 1.1

Line 
1(* ::Package:: *)
2
3(* ::Section:: *)
4(*Renormalization scheme definitions*)
5
6
7(* ::Subsection::Closed:: *)
8(*Declaring unrenormalized parameters*)
9
10
11DefineUnrenormalizedParameters[List[args__]]:=DefineUnrenormalizeOne/@{args};
12
13
14DefineUnrenormalizedParameters[args__]:=DefineUnrenormalizeOne/@{args};
15
16
17DefineUnrenormalizeOne[exp_]:=Block[{},SchemeRules=Append[SchemeRules,FR$delta[{exp},{}]->0];];
18
19
20(* ::Subsection::Closed:: *)
21(*Defining a dependent renormalisation constant*)
22
23
24AddRenormalizationCondition[obj_, rule_]:= Block[{}, SchemeRules=Append[SchemeRules,Rule[obj,rule]]];
25
26
27(* ::Subsection::Closed:: *)
28(*Declaring unrenormalized fields*)
29
30
31DeclareUnrenormalizedFields[List[args__]]:=DeclareUnrenormalizedOneField/@{args};
32
33
34DeclareUnrenormalizedFields[args__]:=DeclareUnrenormalizedOneField/@{args};
35
36
37DeclareUnrenormalizedOneField[exp_]:=Block[{},NoRenoFields=Append[NoRenoFields,exp];];
38
39
40(* ::Subsection::Closed:: *)
41(*Removing vanishing internal parameters*)
42
43
44EnforceZeroParameters[]:=PrePutIndices/@Cases[(Rule[#[[1]],#[[2]]]&/@IParamList)/.ParameterSwitchesAndSwaps[][[3]]/.Rule->myr,myr[_,0]]/.myr->Rule;
45
46
47(* ::Subsection::Closed:: *)
48(*Declaring real renormalization constants*)
49
50
51RealFieldRenormalization[{fields__}]:= RealFieldRenormalization/@{fields};
52
53
54RealFieldRenormalization[field_]:=Block[{}, RealRenoFields=Append[RealRenoFields,field]; CnumQ[FR$deltaZ[{field, field},__]]:=False;];
55
56
57RealFieldRenormalization[fields___]:= RealFieldRenormalization/@{fields};
58
59
60RealFieldRenormalization[]:= RealFieldRenormalization/@Select[Flatten[ClassMemberList/@PartList[[All,1,2]]],!(GhostFieldQ[#]===True||GoldstoneQ[#]===True)&];
61
62
63(* ::Subsection::Closed:: *)
64(*Removing an internal parameter renormalization constants*)
65
66
67RemovingInternalCst[param_]:=Block[{},If[FreeQ[InternalRemovalList,param], InternalRemovalList = Append[InternalRemovalList,param]];];
68
69
70ExecuteRemovingInternalCst[]:= DeleteDuplicates[Flatten[Cases[Global`RenormalizationRules["Internals"][[1]]/.Rule->myr,myr[FR$delta[{PrePutIndices[#]},{}],_]] &/@ InternalRemovalList,1]/.myr->Rule];
71
72
73(* ::Subsection::Closed:: *)
74(*Clearing the rules*)
75
76
77ClearRenormalizationScheme[]:=Block[{},
78  If[Length[FR$RmDblExt]===0,FR$RmDblExt={}];
79  Block[{}, CnumQ[FR$deltaZ[{#, _},__]]:=True;]&/@RealRenoFields;
80  SchemeRules = {}; InternalRemovalList = {}; Global`EnforceZeros = True; RealRenoFields={}; NoRenoFields = {};
81  Global`RenormalizationRules["Externals"]={}; Global`RenormalizationRules["Masses"]={}; Global`RenormalizationRules["Internals"]={{},{}};
82  Global`RenormalizationRules["Fields"]={}; Global`RenormalizationRules["Tadpoles"]={}; Global`RenormalizationRules["MassShifts"] = {};
83];
84
85
86ClearRenormalizationScheme[];
87
88
89(* ::Subsection::Closed:: *)
90(*Enforcing the scheme*)
91
92
93SchemeDependence[rules_]:=Block[{output=rules/.SchemeRules},
94  output=DeleteCases[output,Rule[a_,a_]];
95  Return[output];
96];
97
98
99SchemeDependenceInternal[{intdefs_, intrules_}]:=Block[{newintdefs, newintrules},
100  (* Simplifying the definitions and the rules according to the scheme *)
101  newintdefs = If[FreeQ[#[[2]],FR$delta],Rule[#[[1]],Simplify[#[[2]]]],#]&/@SchemeDependence[intdefs];
102  newintrules=intrules/.newintdefs;
103  newintrules=If[FreeQ[Coefficient[#[[2]],FR$CT],FR$delta],Rule[#[[1]],(#[[2]]/.FR$CT->0)+FR$CT*Simplify[Coefficient[#[[2]],FR$CT]]],#]&/@newintrules;
104  newintrules=SchemeDependence[newintrules];
105  newintrules=Rule[#[[1]], Collect[#[[2]],{FR$delta[__], FR$deltaZ[__]}]]&/@newintrules;
106  newintrules=Rule[#[[1]],#[[1]]+FR$CT*FR$delta[{#[[1]]},{}]]&/@newintrules;
107  newintrules = newintrules/.FR$delta[{Conjugate[aaa_]},{}]:>Conjugate[FR$delta[{aaa},{}]];
108  newintdefs=DeleteCases[newintdefs,Rule[_,0]];
109  (* output *)
110  Return[{PrePutIndices/@newintdefs,PrePutIndices/@newintrules}];
111];
112
113
114(* ::Section:: *)
115(*Misc*)
116
117
118(* ::Subsection:: *)
119(*Version information*)
120
121
122MoGRe`Version = "1.1";
123MoGRe`Date    = "25.10.2019"
124
125
126(* ::Subsection:: *)
127(*Welcome*)
128
129
130Welcome[]:=Block[{},
131  Print["  -- MoGRe: More General Renormalization in FeynRules --"];
132  Print["      > Version: " <> MoGRe`Version];
133  Print["      > Date: " <> MoGRe`Date];
134  Print["  "];
135  Print["      > Author: Benjamin Fuks (LPTHE / Sorbonne U.)"];
136  Print["      > Contact: fuks@lpthe.jussieu.fr"];
137  Print["  "];
138  Print["      > Reference: arXiv:1907.04898 [hep-ph]"];
139  Print["  "];
140];
141
142
143(* ::Subsection::Closed:: *)
144(*Definitions*)
145
146
147numQ[FR$CT]=True;  CnumQ[FR$CT]=False;
148
149
150(* ::Subsection:: *)
151(*Printing*)
152
153
154(* ::Text:: *)
155(*Set of methods for print statements*)
156
157
158TimeStamp[method_, msg_,t0_]:= Block[{ti=SessionTime[]},Print[Style[method,Orange,Bold],": ",msg ," done in ", ti-t0," seconds"];Return[ti];];
159SubMethodStamp[method_, msg_]:= Print[Style[method,Darker[Green],Bold],": starting ", msg];
160Error[ msg_]:=Print[Style["  ** Error: ",Darker[Red],Bold], msg];
161Warning[ msg_]:=Print[Style["  ** Warning: ",Darker[Purple],Bold], msg];
162DebugPrint[msg__]:=If[Global`MoGRe$Debug,Print[msg]];
163
164
165(* ::Subsection::Closed:: *)
166(*Parallelization*)
167
168
169(* ::Text:: *)
170(*Condition for parallelization (expressions longer than 40 terms and kernel availability)*)
171
172
173DoPara[expr_]:=Global`FR$Parallelize && Length[expr]>40 && $KernelCount>1;
174
175
176(* ::Text:: *)
177(*Options of the method*)
178
179
180Options[ParaExec]={Opts->MR$Null};
181
182
183(* ::Text:: *)
184(*Main method (two ways, depending whether the method to parallelize has options*)
185
186
187ParaExec[func_,arg_,OptionsPattern[]]:=Block[{myoptions=OptionValue[Opts],tmpres,inter},
188
189  tmpres=If[myoptions=!=MR$Null,
190    (* with options *)
191    DistributeDefinitions[myoptions];SetSharedFunction[func];
192    Table[inter=arg[[ii]]; ParallelSubmit[{ii,inter},$Output={}; tmpres=func[inter,Sequence@@myoptions]; $Output={OutputStream["stdout",1]}; tmpres],{ii,Length[arg]}],
193
194    (*no options *)
195    SetSharedFunction[func];
196    Table[inter=arg[[ii]]; ParallelSubmit[{ii,inter},$Output={}; tmpres=func[inter]; $Output={OutputStream["stdout",1]}; tmpres],{ii,Length[arg]}]
197  ];
198
199  (* output and exit *)
200  tmpres=WaitAll[tmpres];
201  Return[tmpres];
202];
203
204
205(* ::Subsection::Closed:: *)
206(*Get the class name of a particle*)
207
208
209ParticleClassName[particle_]:=If[MemberQ[PartList[[All,1,2]],particle], Cases[PartList[[All,1]],{_,particle}][[1,1]],Error["The particle does not exist"];Abort[]];
210
211
212(* ::Subsection::Closed:: *)
213(*Get the index type of the flavor index*)
214
215
216ParticleClassFlavorType[class_]:=FlavorIndex/.MR$ClassesRules[class];
217
218
219(* ::Section:: *)
220(*Checks*)
221
222
223(* ::Subsection::Closed:: *)
224(*Mixing*)
225
226
227(* ::Text:: *)
228(*Checks whether the mixing pattern is either a boolean, or a list of 2-tuples*)
229
230
231CheckRenormalizeOptions["Mixing", mixpattern_]:= Block[{},
232  If[BooleanQ[mixpattern],Return[]];
233  If[Not[ListQ[mixpattern]], Error["The mixing fields must be provided as a list of n-tuples"]; Abort[]];
234  If[Not[And@@(ListQ/@mixpattern)], Error["The mixing fields must be provided as a list of n-tuples"]; Abort[]];
235  If[Not[And@@((#===1)&/@(Length[Union[#]]&/@Map[Through,Map[MR$QuantumNumbers,#]&/@mixpattern,{2}]))],
236     Error["Imposing the one-loop mixing of fields with different quantum numbers"]; Abort[];]
237];
238
239
240(* ::Subsection::Closed:: *)
241(*Coupling orders*)
242
243
244(* ::Text:: *)
245(*Checks whether the selected coupling orders are consistent with the QCDOnly option*)
246
247
248CheckRenormalizeOptions["Couplings",couplingorders_]:=Block[{},
249  If[Not[ListQ[couplingorders]],Error["The CouplingOrders option must be a list of coupling orders"]; Abort[]];
250  If[Complement[couplingorders,M$InteractionOrderHierarchy[[All,1]]]=!={},Error["At least one unknown CouplingOrder"]; Abort[]];
251];
252
253
254(* ::Subsection::Closed:: *)
255(*Vevs*)
256
257
258VevDeclarations[]:= Block[{},
259  If[Not[And@@(!FreeQ[#[[1]]//.MR$Definitions,#[[2]]]&/@M$vevs)], Error["vev declarations (M$vevs) incompatible with (unphysical) field definitions."]; Abort[]];
260];
261
262
263(* ::Section:: *)
264(*Formatting the Lagrangian so that it could be processes by the renormalization method*)
265
266
267(* ::Text:: *)
268(*This method flavor-expands the Lagrangian (necessary as some fields of a given class can be massive or massless so that the renormalization is different.*)
269(*Next, the 4-scalar interactions are factored out in case they do not need to be renormalized (option to be provided by the user)*)
270
271
272FormattingLagrangian[lag_,no4S_,cano_]:=Block[{tmplag,lag4S={},l2pt,irules},
273  (* init *)
274  tmplag=PRIVATE`Listize[lag];
275
276  (* 1. Flavor expansion *)
277  tmplag=If[DoPara[tmplag],
278    DebugPrint["  ** Expanding indices over ", Global`FR$KernelNumber," cores"];
279    ParaExec[ExpandIndices,tmplag,Opts->{FlavorExpand->True}],
280    DebugPrint["  ** Expanding indices"];
281    Plus@@(ExpandIndices[#,FlavorExpand->True]&/@tmplag)
282  ];
283  tmplag=PRIVATE`Listize[Plus@@tmplag];
284 
285  (* 2. Removing 4-scalar interactions *)
286  If[no4S,
287    DebugPrint["  ** Removing four scalar interactions"];
288    lag4S=Select[tmplag,(GetFieldContent[#]/._?ScalarFieldQ -> S)==={S,S,S,S}&];
289    tmplag=Complement[tmplag,lag4S];
290  ];
291
292  (* Cleaning the two-point/one-point interactions *)
293  If[cano,
294    DebugPrint["  ** Inserting canonical two point terms..."];
295    tmplag= Select[tmplag,Length[GetFieldContent[#]]>2&];
296    l2pt = Block[{sym=If[SelfConjugateQ[#],2,1], idx = $IndList[#]/.Index[type_]:>Index[type,Unique[dum]]},
297     Which[
298      FermionQ[#]===True,I/sym anti[#].Ga[mu].del[#,mu]-Mass[#]/sym anti[#].# ,
299      VectorFieldQ[#]===True,GetQuadraticTerms[-1/(2 sym) If[Length[idx]==1,FS[#,mu,nu]FS[anti[#],mu,nu],FS[#,mu,nu,aa]FS[anti[#],mu,nu,aa]] +1/sym Mass[#]^2 anti[#][Sequence@@idx]#[Sequence@@idx]],
300      ScalarFieldQ[#]===True,1/sym del[#[Sequence@@idx],mu]del[anti[#][Sequence@@idx],mu]-1/sym Mass[#]^2 anti[#][Sequence@@idx]#[Sequence@@idx]/.fld_[]:>fld
301     ]]&/@Flatten[ClassMemberList/@Select[PartList[[All,1,2]],(GhostFieldQ[#]=!=True)&&(GoldstoneQ[#]=!=True)&]];
302    l2pt = ExpandIndices[#,FlavorExpand->True]&/@l2pt,
303    (* else *)
304    irules = DeleteCases[(Rule[#[[1]],#[[2]]]&/@IParamList)/.ParameterSwitchesAndSwaps[][[3]],Rule[_,0]];
305    irules = Join[ParameterSwitchesAndSwaps[][[2]],Flatten[irules/.ConditionalExpression[a_,_]->a]/.ParameterSwitchesAndSwaps[][[2]]];
306    irules=(If[MemberQ[ParameterSwitchesAndSwaps[][[1]][[All,2]],#1],Solve[#1==(#2//.irules),Cases[ParameterSwitchesAndSwaps[][[1]],{_,#1}][[1,1]]][[1,1]],#1->#2]&)@@@irules;
307    irules = Rule[#[[1]],#[[2]]//.MR$Definitions]&/@irules;
308    irules = Rule[#[[1]],#[[2]]//.irules]&/@irules;
309    DebugPrint["  ** Getting the two point terms..."];
310    l2pt=Join[DeleteCases[GetQuadraticTerms/@tmplag,0], Select[tmplag,Length[GetFieldContent[#]]===1&]];
311    tmplag=Complement[tmplag,PRIVATE`Listize[Plus@@l2pt]];
312    l2pt=Simplify/@ (l2pt//.Join[MR$Definitions,Rule[PrePutIndices[#[[1]]],#[[2]]]&/@irules]);
313  ];
314 
315  (* output *)
316  Return[{Join[tmplag,l2pt],lag4S}];
317];
318
319
320(* ::Section:: *)
321(*Field renormalization*)
322
323
324(* ::Subsection::Closed:: *)
325(*Shortcut to get a field renormalization constant*)
326
327
328ConstructDeltaZ[fld1_,fld2_,idx_,newidx_,tp___]:= FR$CT*FR$deltaZ[{fld1,fld2},{DeleteCases[Append[Cases[idx,Index[FindClassFlavor[fld1],_]],Index[FindClassFlavor[fld2],newidx]],Index[MR$Null,_]]}, tp];
329
330
331(* ::Subsection::Closed:: *)
332(*Main method to get the replacement rule 'bare field => renormalized field'*)
333
334
335(* ::Text:: *)
336(*The method starts by defining a unique sequence of indices that the given field carries.*)
337(*Then, all fields with which the given field can mix, at the loop-level, are determined.*)
338(*The renormalized field can then be computed, with off-diagonal terms corresponding to the potential mixing pattern.*)
339(*The results for the antifeld are then derived, and the replacement rule finally generated.*)
340
341
342RenormalizeField[field_]:= Block[{FreeIndices,SummedIndices,MixList,rfield=field,result,antiresult,newidx=Unique["idx"],spidx=Index[Spin,Unique["sp"]]},
343  (* 1. Indices *)
344  FreeIndices=List@@PrePutIndices[field[Sequence@@Table[Unique["idx"],{Length[$IndList[field]]}]]];
345
346  (* 2. Potiential mixing generated at the loop level *)
347  If[AntiFieldQ[field]===True,rfield=anti[field]];
348  MixList=DeleteCases[PRIVATE`MR$ClassNameList,_?(UnphysicalQ[#]===True&)];
349  MixList=DeleteCases[ MixList,_?(DeleteCases[$IndList[#],Index[FindClassFlavor[#]]]=!=DeleteCases[$IndList[rfield],Index[FindClassFlavor[rfield]]]&)];
350  MixList=(DeleteCases[ MixList,_?( Through[MR$QuantumNumbers[#]]=!=Through[MR$QuantumNumbers[rfield]]&)]);
351
352  (* 3. Main routine *)
353  result=(If[FermionQ[rfield]===True && GhostFieldQ[rfield]=!=True,
354    ConstructDeltaZ[rfield,#, FreeIndices,newidx,"L"]*ProjM[FreeIndices[[1]],spidx] + ConstructDeltaZ[rfield,#,FreeIndices,newidx,"R"] * ProjP[FreeIndices[[1]],spidx],
355    ConstructDeltaZ[rfield,#, FreeIndices,newidx]]*
356    (#[Sequence@@DeleteCases[$IndList[#]/.(FreeIndices/.Index[a_,b_]:>Rule[Index[a],Index[a,b]])/.{Index[Spin,_]->spidx,Index[FindClassFlavor[#],___]->Index[FindClassFlavor[#],newidx]},Index[MR$Null,_]]]/.#[]->#)
357  )&/@MixList;
358  result=Expand[(rfield[Sequence@@FreeIndices]/.rfield[]->rfield )+Plus@@result/2];
359
360  (* Majorana *)
361  If[FermionQ[rfield]===True && SelfConjugateQ[field], result = Expand[result/.{FR$deltaZ[{a_,b_}, {{}}, "R"]:>FR$deltaZ[{a,b}, {{}}, "L"]}]];
362
363  (* 4. Derivation of results for the antifield case *)
364  antiresult=result;
365  If[AntiFieldQ[field]===True,
366    result=result/.{fld_?(FieldQ[#]===True&)->anti[fld],FR$deltaZ[args__]->Conjugate[FR$deltaZ[args]],ProjP[a_,b_]->ProjM[b,a],ProjM[a_,b_]->ProjP[b,a]};
367    result=Refine[result/.FR$deltaZ[{fi_?(AntiFieldQ[#]===True&),gi_?(AntiFieldQ[#]===True&)},bla___]:>FR$deltaZ[{anti[fi],anti[gi]},bla],Assumptions->Element[FR$CT,Reals]],
368    (* else *)
369    antiresult=result/.{fld_?(FieldQ[#]===True&)->anti[fld],FR$deltaZ[args__]->Conjugate[FR$deltaZ[args]],ProjP[a_,b_]->ProjM[b,a],ProjM[a_,b_]->ProjP[b,a]};
370    antiresult=Refine[antiresult/.FR$deltaZ[{fi_?(AntiFieldQ[#]===True&),gi_?(AntiFieldQ[#]===True&)},bla___]:>FR$deltaZ[{anti[fi],anti[gi]},bla],Assumptions->Element[FR$CT,Reals]]
371  ];
372
373  (* 5. Creating elements for the replacement rule *)
374  SummedIndices=Union[Flatten[Complement[PRIVATE`ToIndexList[#]/.Index[_,a_]->a, FreeIndices/.Index[_,a_]->a]&/@List@@result]];
375  FreeIndices=#/.Index[a_,b_]:>Index[a,MyPattern[b,Blank[]]]&/@FreeIndices;
376  rfield=field[Sequence@@FreeIndices]/.a_?(FieldQ[#]===True&)[]->a/.MyPattern->Pattern;
377
378  (* 6. Outputting the replacement rule *)
379  result=If[field===anti[field], {MyRuleDelayed[rfield,MyModule[SummedIndices,result]]}, {MyRuleDelayed[rfield,MyModule[SummedIndices,result]],MyRuleDelayed[anti[rfield],MyModule[SummedIndices,antiresult]]}];
380  Return[result/.{MyRuleDelayed->RuleDelayed,MyModule[{},a_]->a}/.MyModule->Module];
381];
382
383
384(* ::Subsection::Closed:: *)
385(*Expansion over flavors*)
386
387
388(* ::Text:: *)
389(*This perform a full flavor expansion in the rules (to get one rule for each flavor) and in the expression (to have no implicit summed flavor index in the rules)*)
390(*First, this method checks what is the flavor index of the bare field and makes this expansion.*)
391(*Next, it gets the summed flavor indices in the rules, and makes this flavor expasion.*)
392
393
394FlavorExpFieldReno[field_, rule_]:=Block[{myfield=field,flavtype,flavindex,nflavors,newfields={field},newrules={rule},kill},
395  (* we need the field name and not the antifield name *)
396  If[AntiFieldQ[Head[field]], myfield=anti[field]];
397
398  (* treatment of the flavor of the field *)
399  If[PRIVATE`FlavoredQ[field],
400    (* We get the type of flavor index, its range and the name of the index as it apepars in the replacement rule*)
401    flavtype=ParticleClassFlavorType[ParticleClassName[Head[myfield]]];
402    flavindex=(Cases[List@@field, Index[flavtype,_]]/.Pattern->mypattern/.mypattern[a_,_]:>a)[[1]];
403    nflavors = IndexRange[Index[flavtype]][[-1]];
404 
405   (* from class to classmembers *)
406    newfields=Table[ApplyDefinitions[field/.Index[flavtype,_]->Index[flavtype,ii]],{ii,nflavors}]/.fld_?(FieldQ[#]===True&)[]->fld;
407    newrules=Table[
408      ApplyDefinitions[Refine[rule,Assumptions->Element[FR$CT,Reals]]/.flavindex->Index[flavtype,ii]/.FR$deltaZ[{a_,b_},{{Index[flavtype,ii],c___}},tp___]:>FR$deltaZ[{ClassMemberList[Head[myfield]][[ii]],b},{{c}},tp]],
409    {ii,nflavors}];
410    newrules=newrules/.fld_?(FieldQ[#]===True&)[]->fld
411  ];
412
413  (* Treating the implicit sums *)
414  newrules=newrules/.{
415    FR$deltaZ[{a_, cls_}, {{ix_}},tp___]*fld_[c___, ix_,d___]:>Plus@@Table[ApplyDefinitions[FR$deltaZ[{a, ClassMemberList[cls][[jj]]}, {{}},tp]*PrePutIndices[fld[c,jj,d]]],{jj,Length[ ClassMemberList[cls]]}]*kill[ix],
416    Conjugate[FR$deltaZ[{a_, cl_}, {{ix_}},tp___]]*fld_[c___, ix_,d___]:>Plus@@Table[ApplyDefinitions[Conjugate[FR$deltaZ[{a,ClassMemberList[cl][[jj]]},{{}},tp]]*PrePutIndices[fld[c,jj,d]]],{jj,Length[ClassMemberList[cl]]}]*kill[ix]
417  }/.fld_?(FieldQ[#]===True&)[]->fld/.kill[Index[_,a_]]->kill[a];
418  newrules=newrules/.{mymod[{a1___,ix_,a2___},bla_]:>mymod[{a1,a2},bla]/; Not[FreeQ[bla,kill[ix]]]}//.{kill[_]->1, mymod[{},bla_]->bla};
419
420  (* output *)
421  Return[Inner[RuleDelayed,newfields,newrules,List]/.mymod->Module];
422];
423
424
425(* ::Subsection::Closed:: *)
426(*Field mixing*)
427
428
429(* ::Text:: *)
430(*This method removes potential loop-induced mixings if not needed.*)
431(*By default, all potential mixings are kept, except of FieldMixing is set to False or contains a list of 2-tuples representing all fields that can mix*)
432
433
434TreatFieldMixing[renorule_,mix_]:=Block[{right = (renorule/.Module->mod)[[2]]},
435
436  (* Field mixing is not allowed *)
437  If[mix===False,Return[myrd[renorule[[1]],right/.{FR$deltaZ[{a_,b_},__]:>0/;a=!=b}]/.myrd->RuleDelayed/.mod->Module]];
438
439  (* Field mixing only allowed for specific fields *)
440  Return[myrd[renorule[[1]],right/.{FR$deltaZ[{flds__},__]:>0/;(Not[MatchQ[{flds},{a_,a_}]]&& Not[MemberQ[Sort/@Flatten[Subsets[#,{2}]&/@mix,1],Sort[{flds}]]] )}]/.myrd->RuleDelayed/.mod->Module];
441
442];
443
444
445(* ::Subsection::Closed:: *)
446(*Tadpole renormalization*)
447
448
449TadpoleShift[vev_,field_]:=Block[{ShiftedField},
450  (* Sanity check *)
451  If[Not[SelfConjugateQ[field]],Error["All fields getting a vev must be selfconjugate"];Abort[]];
452
453  (* Getting the shift *)
454  ShiftedField=Select[Expand[If[Im[Coefficient[field,vev]]===0,(field+HC[field])/.{vev->0},(field-HC[field])/.{vev->0}]],FieldQ[#]===True&];
455  CnumQ[FR$deltat[ShiftedField]]=False;
456  Return[Rule[ShiftedField,ShiftedField-FR$CT*FR$deltat[ShiftedField]/Mass[ShiftedField]^2]];
457];
458
459
460TadpoleShifts[fields_]:=Block[{RelevantFields = Select[fields,ScalarFieldQ[#]===True&],rules},
461  (*Sanity checks *)
462  VevDeclarations[];
463
464  (* Getting the field connected with the vevs and creating the rule; filtering all irrelevant fields*)
465  rules=TadpoleShift[Sequence@@#]&/@({#[[2]],#[[1]]/.MR$Definitions}&/@M$vevs);
466  Return[Select[rules,MemberQ[RelevantFields,#[[1]]]&]];
467];
468
469
470DeltaMShifts[lag_]:=Block[{massterms,shifts,FR$delta2,shiftedlag},
471  (* Extracting the mass shift induced by the normalization procedure *)
472  massterms =Select[lag,!FreeQ[#,FR$CT]&&FreeQ[#,FR$deltaZ]&& FreeQ[#,_?(GhostFieldQ[#]===True&)]&& FreeQ[#,_?(GoldstoneQ[#]===True&)]&& Length[GetFieldContent[#]]===2&];
473  If[massterms==={},massterms={{}}, $Output=OpenWrite[];massterms=(GetMassSpectrum[Plus@@massterms])[[1,2;;,{1,2}]];$Output={OutputStream["stdout",1]}];
474  shifts=Flatten[If[FreeQ[#,FR$delta],{},Solve[Simplify[#[[2]]/.{FR$deltat[_]->0,FR$delta->FR$delta2}]==#[[2]],{Cases[#[[2]],_FR$delta,\[Infinity]][[1]]}]/.FR$delta2->FR$delta]&/@massterms];
475 
476  (* Preparing the replacement rules and applying them *)
477  shifts=Join[shifts, (Expand[#/.shifts]&/@Global`RenormalizationRules["Internals"][[1]])];
478  shifts=Select[shifts,!MemberQ[Global`RenormalizationRules["Internals"][[1,All,2]],#[[2]]]&];
479  shiftedlag =Expand[#/.shifts]&/@lag;
480
481  (* output *)
482  Return[{shifts, PRIVATE`Listize[Expand[Plus@@shiftedlag]]}];
483];
484
485
486(* ::Subsection:: *)
487(*Wrapper*)
488
489
490(* ::Text:: *)
491(*Loop over all relevant fields; the lagrangian must be properly formatted*)
492
493
494Renormalization["Fields", lag_, couplingorders_]:= Block[{fields, orderrules,myrule},
495  (* Relevant orders *)
496  orderrules =myrule[#,1]&/@ Select[M$InteractionOrderHierarchy[[All,1]],Not[MemberQ[couplingorders,#]]&];
497
498  (* Getting the fields involved in the interactions related to the coupling orders under consideration *)
499  fields=Select[lag, Union[
500    Map[PRIVATE`GetIntOrder,If[ListQ[#],#,{#}]//.{ del[a_,_]->a,_?(FieldQ[#]===True&)[__]->1,Dot->Times,Index[_,i_]->i}]/.PRIVATE`GetIntOrder->IntOrder/.(orderrules/.myrule->Rule)]=!={1}&];
501  fields=Flatten[GetFieldContent/@fields];
502  fields=Union[If[AntiFieldQ[#],anti[#],#]&/@ DeleteCases[fields,_?(UnphysicalQ[#] || GhostFieldQ[#]===True &)]];
503  fields = Select[fields,!MemberQ[NoRenoFields,#]&];
504
505  (* Renormalization *)
506  fields = Flatten[RenormalizeField /@ fields];
507
508  (* exiting the function *)
509  Return[fields];
510];
511
512
513(* ::Section:: *)
514(*Parameter renormalization*)
515
516
517(* ::Subsection::Closed:: *)
518(*Get replacement rules*)
519
520
521(* ::Text:: *)
522(*This method check whether some external parameters and internal parameters must be traded.  It then gets a formatted list for the switches.*)
523(*It also format the parameter names (matrix elements with things as "1x1" for the "1,1" element must be replaced by things like M[1,1])*)
524
525
526ParameterSwitchesAndSwaps[]:=Block[{switches1,switches2,renames},
527
528  (* Getting and formatting the switches *)
529  switches1=DeleteCases[FR$LoopSwitches/.Index[_,b_]->b/.MR$Restrictions,Rule[_,0]];
530  switches2=DeleteCases[FR$RmDblExt/.Index[_,b_]->b/.MR$Restrictions,Rule[_,0]];
531  renames = Reverse/@DeleteCases[Union[Rule[#[[1]],#[[2]]]&/@PRIVATE`$ParamListtemp],Rule[a_,a_]]/.Index[_,b_]->b;
532
533  (* output *)
534  Return[{switches1,switches2,renames}];
535];
536
537
538(* ::Subsection::Closed:: *)
539(*Main method to get the replacement rules 'bare parameter' to 'renormalized parameter'*)
540
541
542(* ::Text:: *)
543(*The method starts by getting all necessary switches (trade-of of external/internal parameters, etc...).*)
544(*Then we start the real business:*)
545(*   - we get the list of External parameters.*)
546(*   - we focus on the internal parameters and deduce the list of those that can appear in the bare Lagrangian*)
547(*   - we then get the list of masses to renormalize (removing those of the non QCD fields in the QCD case)*)
548(* Finally, we renormalize the external parameters (g -> g + delta g) and deduce the renormalization relations for the internal parameters.*)
549
550
551Renormalization["Parameters", fields_]:=Block[{switches, Externals,Internals, Masses, IntRules, ZeroRules,simplifs,myRule},
552  (* Getting the switches *)
553  switches = ParameterSwitchesAndSwaps[];
554
555  (* External parameters (without the masses and after applying the switches *)
556  Externals=Sort[Flatten[EParamList[[All,2,All,2,1]]]]/.switches[[3]]/.Rule@@@switches[[1]]/.switches[[2]];
557  Externals=Select[Externals,Not[MemberQ[MassList[[2,All,2]],#]]&];
558  (CnumQ[FR$delta[{#},{}]]=CnumQ[#])&/@Externals;
559
560  (* Internal parameters (applying the switches where relevant) *)
561  Internals = If[Global`EnforceZeros,
562    DeleteCases[(Rule[#[[1]],#[[2]]]&/@IParamList)/.switches[[3]],Rule[_,0]],
563    (Rule[#[[1]],#[[2]]]&/@IParamList)/.switches[[3]]
564  ];
565  Internals=If[FreeQ[switches[[1,All,2]],#[[1]]],#,Solve[#[[1]]==(#[[2]]//.Internals),Cases[switches[[1]],{_,#[[1]]}][[1,1]] ][[1]]]&/@Internals;
566  Internals=Join[switches[[2]],Flatten[Internals/.ConditionalExpression[a_,_]->a]/.switches[[2]]];
567  Internals = Rule[#[[1]],#[[2]]//.MR$Definitions]&/@Internals;
568  Internals = Rule[#[[1]],#[[2]]//.Internals]&/@Internals;
569
570  (* Masses *)
571  Masses=DeleteCases[Union[Mass/@fields],0];
572  (CnumQ[FR$delta[{#},{}]]=False)&/@Masses;
573
574  (* From bare parameters to renormalized parameters *)
575  simplifs=Cases[SchemeRules/.Rule->myRule,myRule[_,0]]/.myRule->Rule;
576  Externals=Rule[#,#+FR$CT*FR$delta[{#},{}]]&/@Externals;
577  Masses=Rule[#,#+FR$CT*FR$delta[{#},{}]]&/@Masses;
578  IntRules =Select[ Internals/.Join[Externals, Masses],Not[FreeQ[#,FR$CT]]&];
579  ZeroRules =Rule[FR$delta[{#},{}],0]&/@Select[ Internals/.Join[Externals, Masses],FreeQ[#,FR$CT]&][[All,1]];
580  Internals=Flatten[If[CnumQ[#], List[Rule[#,#+FR$CT*FR$delta[{#},{}]],Rule[Conjugate[#],Conjugate[#]+FR$CT*Conjugate[FR$delta[{#},{}]]]],Rule[#,#+FR$CT*FR$delta[{#},{}]]]&/@Internals[[All,1]]];
581  Internals=DeleteCases[Internals/.simplifs/.ZeroRules,Rule[a_,a_]];
582  Internals = Internals/.FR$delta[{Conjugate[aaa_]},{}]:>Conjugate[FR$delta[{aaa},{}]];
583  IntRules=If[FreeQ[#/.simplifs/.ZeroRules,FR$CT],
584    Rule[#[[1]],#[[1]]],
585    Rule[#[[1]],Expand[Normal[Series[#[[2]]/.simplifs,{FR$CT,0,1}]]/(Normal[Series[#[[2]],{FR$CT,0,0}]]/#[[1]])]]]&/@IntRules;
586  IntRules = Rule[FR$delta[{#[[1]]},{}],Collect [#[[2]]-#[[1]],{ FR$delta[__]},Simplify]]/.FR$CT->1&/@IntRules;
587 
588  (* Output *)
589  Return[{(Externals/.simplifs),(Masses/.simplifs), Internals, IntRules}];
590];
591
592
593(* ::Section:: *)
594(*Tools*)
595
596
597(* ::Subsection::Closed:: *)
598(*Get all renormalization constants*)
599
600
601(* ::Text:: *)
602(*Allow to extract all the possible renormalization constants of the *)
603
604
605ParameterRenormalizationConstants[]:= Block[{RenoPrms,RenoFields,dummy,tmplag},
606  {RenoPrms,dummy,dummy,dummy}=Renormalization["Parameters",{}]; 
607  RenoPrms = DeleteCases[Coefficient[#,FR$CT]&/@RenoPrms[[All,2]],Conjugate[__]];
608  RenoPrms = Join[FR$delta[{#},{}]&/@Union[MassList[[2,All,2]]],RenoPrms];
609  Return[Union[RenoPrms]];
610];
611
612
613FieldRenormalizationConstants[lag_]:=Block[{dummy = Global`MoGRe$Debug, tmplag, RenoFields},
614  (* Verbose = off *)
615  Global`MoGRe$Debug=False;
616
617  (* Lagrangian formatting *)
618  tmplag = Expand[lag-(lag/.{(_?FieldQ)[__]->0,_?FieldQ->0})];
619  tmplag= FormattingLagrangian[tmplag,False][[1]];
620
621  (* Constants *)
622  RenoFields=Flatten[FlavorExpFieldReno[Sequence@@(#/.{Module->mymod})]&/@Renormalization["Fields",tmplag,OptionValue[Global`MoGRe$Renormalize,CouplingOrders]]];
623  If[OptionValue[Global`MoGRe$Renormalize,FlavorMixing]=!=True,RenoFields = TreatFieldMixing[#,OptionValue[Global`MoGRe$Renormalize,FlavorMixing]]&/@RenoFields];
624  RenoFields=Expand[Coefficient[#,FR$CT]]&/@RenoFields[[All,2]];
625  RenoFields=Union[DeleteCases[Select[Expand[#],!FreeQ[#,FR$deltaZ]&]&/@(RenoFields/.Plus->Sequence),Conjugate[_]]];
626
627  (* Output and verbose set to default again *)
628  Global`MoGRe$Debug = dummy;
629  Return[RenoFields];
630];
631
632
633RenormalizationConstants[lag_]:= InputForm[Join[FieldRenormalizationConstants[lag], ParameterRenormalizationConstants[]]];
634
635
636(* ::Section:: *)
637(*Main method*)
638
639
640ShiftRules[]:= Block[{irules,rules},
641  irules = ExecuteRemovingInternalCst[];
642  rules = Join[Global`RenormalizationRules["Externals"],Global`RenormalizationRules["Masses"],Global`RenormalizationRules["Internals"][[2]],Global`RenormalizationRules["Fields"]];
643  Return[rules/.irules];
644];
645
646
647Options[Global`MoGRe$Renormalize]={Exclude4Scalars->False,FlavorMixing->True, Global`CouplingOrders->M$InteractionOrderHierarchy[[All,1]], Global`CanonicalTwoPoints -> False};
648
649
650Global`MoGRe$Renormalize[lag_, OptionsPattern[]]:=Block[{time=SessionTime[],lasttime,tmplag,lag4S, dum1, dum2,ishift,tmpCC},
651
652  (* Initialization *)
653  Welcome[];
654  If[Global`MoGRe$Debug,SubMethodStamp["Renormalize[]","initialization"]];
655  FR$Loop=True;                (*Keep two point in the vertex list*)
656  CheckRenormalizeOptions["Mixing",OptionValue[FlavorMixing]];
657  CheckRenormalizeOptions["Couplings",OptionValue[CouplingOrders]];
658  tmplag = Expand[lag-(lag/.{(_?FieldQ)[__]->0,_?FieldQ->0})];
659  If[Global`EnforceZeros, tmplag = tmplag/.EnforceZeroParameters[]];
660  If[Global`MoGRe$Debug,lasttime=TimeStamp["Renormalize[]","initialization" ,time]];
661 
662  (* Formatting the Lagrangian *)
663  If[Global`MoGRe$Debug,SubMethodStamp["Renormalize[]","Lagrangian formatting"]];
664  {tmplag,lag4S} = FormattingLagrangian[tmplag,OptionValue[Exclude4Scalars],OptionValue[CanonicalTwoPoints]];
665  If[Global`MoGRe$Debug,lasttime=TimeStamp["Renormalize[]","Lagrangian formatting" ,lasttime]];
666 
667  (* Field renormalization *)
668  If[Global`MoGRe$Debug,SubMethodStamp["MoGRe$Renormalize[]","field renormalization"]];
669  DebugPrint["  ** Getting the replacement rules (bare into renormalized)"];
670  Global`RenormalizationRules["Fields"]=Flatten[FlavorExpFieldReno[Sequence@@(#/.{Module->mymod})]&/@Renormalization["Fields",tmplag, OptionValue[CouplingOrders]]];
671  If[OptionValue[FlavorMixing]=!=True,
672    DebugPrint["  ** Removing unnecessary field mixing"];
673    Global`RenormalizationRules["Fields"] = TreatFieldMixing[#,OptionValue[FlavorMixing]]&/@Global`RenormalizationRules["Fields"];
674  ];
675  DebugPrint["  ** tadpole renormalization"];
676  Global`RenormalizationRules["Tadpoles"] = TadpoleShifts[If[Head[#]=!=Symbol, Head[#], #]&/@Global`RenormalizationRules["Fields"][[All,1]]];
677  If[Global`MoGRe$Debug,lasttime=TimeStamp["MoGRe$Renormalize[]","field renormalization" ,lasttime]];
678
679  (* Parameter renormalization *)
680  If[Global`MoGRe$Debug,SubMethodStamp["MoGRe$Renormalize[]","parameter renormalization"]];
681  DebugPrint["  ** Getting the replacement rules (bare into renormalized)"];
682  {Global`RenormalizationRules["Externals"],Global`RenormalizationRules["Masses"],dum2,dum1}=Renormalization["Parameters",If[Head[#]=!=Symbol, Head[#], #]&/@Global`RenormalizationRules["Fields"][[All,1]]];
683  Global`RenormalizationRules["Internals"]={dum1,dum2};
684  DebugPrint["  ** Enforcing the renormalization scheme"];
685  Global`RenormalizationRules["Externals"]= SchemeDependence[Global`RenormalizationRules["Externals"]];
686  Global`RenormalizationRules["Masses"]   = SchemeDependence[Global`RenormalizationRules["Masses"]];
687  Global`RenormalizationRules["Internals"]= SchemeDependenceInternal[Global`RenormalizationRules["Internals"]];
688  If[Global`MoGRe$Debug,lasttime=TimeStamp["MoGRe$Renormalize[]","parameter renormalization" ,lasttime]];
689   
690  (* Lagrangian (with a special treatment for charge-conjugate fields) *)
691  If[Global`MoGRe$Debug,SubMethodStamp["MoGRe$Renormalize[]","Lagrangian renormalization"]];
692  DebugPrint["  ** Shifting parameters and fields"];
693  ishift = ShiftRules[];
694  tmplag = (#/.CC[fld_][inds__]:>tmpCC[fld[inds]]/.ishift/.Global`RenormalizationRules["Tadpoles"])&/@tmplag;
695  tmplag =Normal[Series[(Expand[#/.tmpCC->CC/.Dot->FR$Dot/.FR$Dot->Dot]),{FR$CT,0,1}]]&/@tmplag;
696  tmplag=PRIVATE`Listize[Expand[Plus@@tmplag]];
697  DebugPrint["  ** Absorbing all mass shifts in the renormalization constant"];
698  {Global`RenormalizationRules["MassShifts"],tmplag}=DeltaMShifts[tmplag]; 
699  tmplag = Collect[#,{FR$CT},Simplify]&/@tmplag;
700  If[Global`MoGRe$Debug,lasttime=TimeStamp["MoGRe$Renormalize[]","Lagrangian renormalization" ,lasttime]];
701
702  (* output *)
703  Return[Plus@@Join[tmplag,lag4S]];
704];