Division by 0 in MakeSPheno

[JPEllis]

In MakeSPheno, when it is "Processing couplings for 2-loop effective potential", the replacement

Code: Select all

temp = Select[temp /.subZeroGaugeLess ,#[[1,2,1]]=!=0&];

can cause divisions by zero error even if the limit is well defined.

Assuming that the terms are smooth in the neighbourhood of the couplings being equal to zero, then this can be fixed with the following patch:

Code: Select all

--- a/SARAH-4.8.5/Package/SPheno/SPhenoCoupling.m
+++ b/SARAH-4.8.5/Package/SPheno/SPhenoCoupling.m
@@ -99,7 +99,7 @@ CouplingUsedForEffPot=True; (* in order to apply the correct color sum *)
  listBrokenGaugeCouplings=DeleteCases[Transpose[BetaGauge][[1]],strongCoupling];
 (* listBrokenGaugeCouplings=Extract[Gauge,Join[#,{4}]&/@(Position[SGauge /. A_[{b__}]\[Rule]A,#][[1]]&/@Select[SGauge /. A_[{b__}]\[Rule]A,FreeQ[Particles[SA`CurrentStates],#]&])]; *)
 subZeroGaugeLess=Table[listBrokenGaugeCouplings[[i]]->0,{i,1,Length[listBrokenGaugeCouplings]}];
-temp = Select[temp /. subZeroGaugeLess,#[[1,2,1]]=!=0&];
+temp = Select[Fold[Limit[#1, #2] &, temp, subZeroGaugeLess] ,#[[1,2,1]]=!=0&];
 
 temp=SPhenoCouplingList[temp];
 CouplingUsedForEffPot=False;
@@ -110,7 +110,7 @@ WriteSPhenoAllCouplings[SPhenoCouplings4P,parametersAll4P,namesAll4P,"CouplingsF
 
 DynamicCouplingsEffpot="3 point vertices";
 AllRelevant=getAllRelevantCouplings[VertexListNonCC];
-AllRelevant =  Select[AllRelevant /. subZeroGaugeLess,If[Length[#[[1]]]===3,(#[[1,2,1]]=!=0) ||( #[[1,3,1]]=!=0),#[[1,2,1]]=!=0]&];
+AllRelevant =  Select[Fold[Limit[#1, #2] &, AllRelevant, subZeroGaugeLess] ,If[Length[#[[1]]]===3,(#[[1,2,1]]=!=0) ||( #[[1,3,1]]=!=0),#[[1,2,1]]=!=0]&];
 temp=SPhenoCouplingList[AllRelevant];
 SPhenoCouplings3P=temp[[1]];
 parametersAll3P=temp[[2]];

[FStaub]

I don't see how there can ever be terms ~1/g in the vertices. That looks to me that there is another problem with the model.

[JPEllis]

I didn't have 1/g, but I did have terms with "g1 g2 / sqrt(g1^2 + g2^2)" which, depending on the particular context, would either produce a division by 0 error, a 0/0 error, or result in a term with 0 * infinity when g1 and g2 are naively set to zero. The limits are perfectly fine though.

[FStaub]

Ok, these terms show up because you make some replacements for the Weinberg angle, it seems. However, I'm not sure if I can recommend that because it's a purely tree-level relation. In particular, for some precision observables it's better to use the on-shell value (\sin\Theta_W as function of the measured M_Z and M_W) instead of the DR/MS value, see for instance the discussion in 1208.0934 for the B-decays.

[FStaub]

PS: Thanks for the patch, of course. I'll just need to test that it works fine with all models and that it doesn't lead to a huge delay for more complicated models. In that case, one would need to add a check if the simple replacement is working or if one has really to take the limit.

[JPEllis]

I suspect that doing a check if the substitution works and then trying the limit will only end up doubling the evaluation time (if that is indeed ever a bottle neck). In my experience, the limit will work nearly just as fast as a direct substitution if the substitution is well defined. In cases where the substitution is not well defined, then the limit can take much longer time

I tested this with the following two functions (one numerical, the other symbolic) and the substitution is only very slightly slower:

Code: Select all

(* Numerical version *)
f[x_, y_] := 
  Tr[(x RandomReal[{-1, 1}, {100, 100}]) .
    (y RandomReal[{-1, 1}, {100, 100}]) .
    (x RandomReal[{-1, 1}, {100, 100}]) .
    (y RandomReal[{-1, 1}, {100, 100}])] /
   Sqrt[Tr[(x RandomReal[{-1, 1}, {100, 100}]).
     (y RandomReal[{-1, 1}, {100, 100}]) .
     (x RandomReal[{-1, 1}, {100, 100}]) .
     (y RandomReal[{-1, 1}, {100, 100}])]];
(* Symbolic version *)
f[x_, y_] := 
  Tr[(x RandomInteger[{-1, 1}, {100, 100}]) .
    (y RandomInteger[{-1, 1}, {100, 100}]) .
    (x RandomInteger[{-1, 1}, {100, 100}]) .
    (y RandomInteger[{-1, 1}, {100, 100}])] /
   Sqrt[Tr[(x RandomInteger[{-1, 1}, {100, 100}]).
     (y RandomInteger[{-1, 1}, {100, 100}]) .
     (x RandomInteger[{-1, 1}, {100, 100}]) .
     (y RandomInteger[{-1, 1}, {100, 100}])]];
     
f[x, y] /. {x -> 0, y -> 0} // Timing
g[x, y] /. {x -> 0, y -> 0} // Timing

Fold[Limit[#1, #2] &, f[x, y], {x -> 0, y -> 0}] // Timing
Fold[Limit[#1, #2] &, g[x, y], {x -> 0, y -> 0}] // Timing

and the output:

Code: Select all

{7.77,    Indeterminate}
{6.59333, Indeterminate}
{8.38333, 0.}
{7.05333, 0}

Still better to check this against all the models though.

[FStaub]

Ok, these terms show up because you make some replacements for the Weinberg angle, it seems. However, I'm not sure if I can recommend that because it's a purely tree-level relation. In particular, for some precision observables it's better to use the on-shell value (\sin\Theta_W as function of the measured M_Z and M_W) instead of the DR/MS value, see for instance the discussion in 1208.0934 for the B-decays.

To clarify that: my recommendation is actually not to use any Dependence by default for any parameter when generating the SPheno output. Thus, just keep the Weinberg angle, and it will be calculated/used appropriately by SPheno:

• - in most cases where the MS/DR values is needed from the diagonalisation of the vector boson mass matrix
  - in some cases like the calculation of quark flavour observables and \delta\rho the suitable, numerical values are used

Cheers,
Florian