# RGEs of the MSSM - problems with Mathematica

I'm having some troubles with the trilinear soft couplings of the MSSM RGEs. I've used the ones written in Martin's supersymmetry primer and I run them using mathematica, if I do so without taking into account the soft trilinear terms (by putting them to zero at all scales and not considering that they get generated even for a[M_SUSYbreaking]=0) everything works well. I get unification of gauge coupling constants at some high scale and the running of other terms seems to be OK.

When I consider the running of trilinear terms too things don't work at all. Not only the running of soft terms change in a strange way, also ِYukawa terms and gauge couplings change and I know this can't be true as there is no soft trilinear dependence in these adimensional couplings.

This obviously points towards a bug in the code or something related but I've checked it many times and I don't know where the mistake or the typo could be.

I've also tried to run the adimensional parameters first and then using the numerical solution I get solve the soft parameters, by doing this I think that, at least, I can force the adimensional parameters to have a good behaviour and then see what the result for the other things is. The problem is that I get a Mathematica error that I don't know how to solve,

"Computed derivatives do not have dimensionality consistent with the initial conditions."

I've seen that these can come when trying to solve a system of differential equations using solutions from the numerical resolution of another system of diff equations (solution of NSolve into another NSolve). I don't know how could I do this. Any thoughts?

I copy the code here in case there is something really obvious that I'm missing. I don't know, maybe is just a typo that I haven't seen yet but I really don't see where it could be. As I guess the problem is with Mathematica maybe someone has dealt with a similar situation in the past even if is not writing RGEs.

Qmax = 35.
tanbeta = 10.
htzero = 173/174 Sqrt[1 + tanbeta^2]/tanbeta
hbzero = 4/174 Sqrt[1 + tanbeta^2]
gaugino = 1.5
sfermions = 1.5
higgsmassparam = 1.5
Exp[Qmax] .174
Xu[t_] := 2 ht[t]^2 ( mhu2[t] + mq2[t] + mu2[t]) + 2 at[t]^2
Xd[t_] := 2 hb[t]^2 (mhd2[t] + mq2[t] + md2[t]) + 2 ab[t]^2
system = NDSolve[{ht'[t] == 1/(16 \[Pi]^2) ht[t] (6 ht[t]^2 + hb[t]^2 - 16/3 g3[t]^2 - 3 g2[t]^2 - 13/15 g1[t]^2), ht[0] == htzero, hb'[t] == 1/(16 \[Pi]^2) hb[t] (6 hb[t]^2 + ht[t]^2 - 16/3 g3[t]^2 - 3 g2[t]^2 - 7/15 g1[t]^2), hb[0] == hbzero, g3'[t] == -3 /(16 \[Pi]^2) g3[t]^3, g3[0] == Sqrt[4 \[Pi] 0.118], g2'[t] == 1 /(16 \[Pi]^2) g2[t]^3, g2[0] == Sqrt[4 \[Pi]/30], g1'[t] == 33/5 /(16 \[Pi]^2) g1[t]^3, g1[0] == Sqrt[4 \[Pi]/60], M3'[t] == -6 /(16 \[Pi]^2) g3[t]^2 M3[t], M3[Qmax] == gaugino, M2'[t] == 2 /(16 \[Pi]^2) g2[t]^2 M2[t], M2[Qmax] == gaugino, M1'[t] == 66/5 /(16 \[Pi]^2) g1[t]^2 M1[t], M1[Qmax] == gaugino, mhu2'[t] == 3/16/\[Pi]^2  Xu[t] - 1/16/\[Pi]^2 ( 6 g2bis[t]^2 M2bis[t]^2 + 6/5 g1bis[t]^2 M1bis[t]^2), mhu2[Qmax] == higgsmassparam, mhd2'[t] == 3/16/\[Pi]^2  Xd[t] - 1/16/\[Pi]^2 (6 g2bis[t]^2 M2bis[t]^2 + 6/5 g1bis[t]^2 M1bis[t]^2), mhd2[Qmax] == higgsmassparam, mq2'[t] == 1/16/\[Pi]^2 (Xu[t] + Xd[t]) - 1/16/\[Pi]^2 (32/3 g3bis[t]^2 M3bis[t]^2 + 6 g2bis[t]^2 M2bis[t]^2 + 2/15 g1bis[t]^2. M1bis[t]^2), mq2[Qmax] == sfermions, mu2'[t] == 2/16/\[Pi]^2 Xu[t] - 1/16/\[Pi]^2 (32/3) g3bis[t]^2 M3bis[t]^2 - 1/16/\[Pi]^2 (32/15) g1bis[t]^2 M1bis[t]^2, mu2[Qmax] == sfermions, md2'[t] == 2/16/\[Pi]^2 Xd[t] - 1/16/\[Pi]^2 (32/3) g3bis[t]^2 M3bis[t]^2 - 1/16/\[Pi]^2 (8/15) g1bis[t]^2 M1bis[t]^2,  md2[Qmax] == sfermions, at'[t]==(1/16Pi^2)((at[t](18htbis[t]^2+hbbis[t]^2-(16/3)g3bis[t]^2-3g2bis[t]^2-(13/15)g1bis[t]^2)+2ab[t]hbbis[t]htbis[t]+ htbis[t]((32/3)g3bis[t]^2M3[t]+6g2bis[t]^2M2[t]+(26/15)g1bis[t]^2M1[t]))),at[Qmax]==0, ab'[t]==(1/16Pi^2)((ab[t](18hbbis[t]^2+htbis[t]^2-(16/3)g3bis[t]^2-3g2bis[t]^2-(13/15)g1bis[t]^2)+2at[t]htbis[t]hbbis[t]+\hbbis[t]((32/3)g3bis[t]^2M3[t]+6g2bis[t]^2M2[t]+(26/15)g1bis[t]^2M1[t]))),ab[Qmax]==0}, {g3[t], g2[t], g1[t], ht[t], hb[t], M1[t], M2[t], M3[t],mhu2[t], mhd2[t], mq2[t], mu2[t], md2[t],at[t],ab[t]}, {t, 0, Qmax}]
g3bis[t_] = g3[t] /. system

g2bis[t_] = g2[t] /. system

g1bis[t_] = g1[t] /. system

htbis[t_] = ht[t] /. system

hbbis[t_] = hb[t] /. system

M3bis[t_] = M3[t] /. system

M2bis[t_] = M2[t] /. system

M1bis[t_] = M1[t] /. system

mhu2bis[t_] = mhu2[t] /. systema

mhd2bis[t_] = mhd2[t] /. systema

atbis[t_]=at[t]/.system

abbis[t_]=ab[t]/.system

Plot[{g1bis[z], g2bis[z], g3bis[z]}, {z, 0, Qmax}, Frame -> True]

Plot[{M1bis[z], M2bis[z], M3bis[z]}, {z, 0, Qmax}, Frame -> True]

Plot[{mhu2bis[z], mhd2bis[z]}, {z, 0, Qmax}, Frame -> True]


As you can see I've taken the third family approximation and neglected lepton terms too. These code works without putting at[t] and ab[t] terms and dependencies, when not doing so, nothing works.

-
Would scicomp.stackexchange.com or mathematica.stackexchange.com be a better home for this question? – Qmechanic Aug 1 '13 at 20:18
For mhu2bis[] and mhd2bis[] you use the replacement rule systema, which doesn't seem to be defined (might be overlooking something). – Vibert Aug 1 '13 at 21:03
ups! sorry this is from some tests I was making... I use system for everything. – user27764 Aug 1 '13 at 21:17
yes maybe I should post these on the other forums, I just thought that maybe had dealt with something similar trying to run the MSSM RGEs before. – user27764 Aug 1 '13 at 21:18
I haven't run into this in this context, but I've run into similar things in Mathematica before, and found that things worked much better when I rewrote them in Python (using SciPy). I know that's not a very satisfying answer-- presumably real Mathematica experts could figure out how to fix it. – Matt Reece Aug 1 '13 at 22:44

Typos stops your evaluation of trilinears at[t] and ab[t] - you wrote (1/16Pi^2) as loop factor for these terms by mistake, and blow up at[0] and ab[0] seriously.

There are also 2 other typos, one appears in the beta function of ab[t] \hbbis[t] should be hbbis[t]; the other one is systema should be system for the interpolation of mhu2 and mhd2.

I feel, it's not necessary to use -bis variables in the beta-function. I put my revised version below, it works in my computers (Mathematica 9 and Mathematica 10)

Also, the original code omitted contribution from $\tau$ lepton, e.g. $y_\tau,a_\tau,m_{L_3}^2,m_{e_3}^2$--They are not neglectable in large $\tan\beta$ regime. In this updated version I add these missing parts as well as the Higgs $\mu$ and $B_\mu$ terms to refine the result. To make the code more readable, I defined some new variables and added illustrations along with them.

This code is able to reproduce unification of gauge couplings and running of soft breaking terms, which agree with Martin's Supersymmetry Primer p61 Fig. 6.8 (LHS: Martin, RHS: Code):

and p105 Fig. 8.4 (LHS: Martin, RHS: Code):

Qmax = 35.;(* t=ln[Q/(top quark mass)]*)
0.174*Exp[Qmax]
vevhiggs = 0.246;(* in TeV*)
vevsusy =
vevhiggs/Sqrt[2];(*vevsusy=Sqrt[vevu^2+vevd^2]*)
tanbeta = 10.;
vevu = vevsusy*tanbeta/Sqrt[1 + tanbeta^2];
vevd = vevsusy/Sqrt[1 + tanbeta^2];
tquark = 0.174;
bquark = 0.00430;
taulep = 0.00178;
htzero = tquark/vevu;
hbzero =
bquark/vevd;
htauzero = taulep/vevd;

gaugino = 0.6; (* Subscript[m, 1/2] GUT universal gaugino mass*)
\
susymu = 0.83 ;(* \[Mu] GUT supersymmetric higgs \[Mu] term*)
msf = \
0.2;(* m0 GUT universal sfermion mass*)
higgsmassparam =
Sqrt[msf^2 + susymu^2];(* GUT universal soft higgs term*)
at0 = 0;
ab0 = 0;
atau0 = 0; (* au,ad,atau are trilinear scalar Yukawa couplings*)
b0 = \
0; (* Subscript[B, \[Mu]] soft higgs term*)
\[Kappa] =
1/(16 \[Pi]^2);
Xu[t_] := 2 ht[t]^2 (mhu2[t] + mq2[t] + mu2[t]) + 2 at[t]^2;
Xd[t_] := 2 hb[t]^2 (mhd2[t] + mq2[t] + md2[t]) + 2 ab[t]^2;
Xtau[t_] := 2 htau[t]^2 (mhd2[t] + mq2[t] + md2[t]) + 2 ab[t]^2;
strm2[t_] :=
mhu2[t] - mhd2[t] + (mq2[t] - mL2[t] - 2 mu2[t] + md2[t] + me2[t]);
system = NDSolve[{
ht'[t] == \[Kappa] ht[
t] (6 ht[t]^2 + hb[t]^2 - 16/3 g3[t]^2 - 3 g2[t]^2 -
13/15 g1[t]^2), ht[0] == htzero,
hb'[t] == \[Kappa] hb[
t] (6 hb[t]^2 +(*ht[t]^2+htau[t]^2*)-16/3 g3[t]^2 -
3 g2[t]^2 - 7/15 g1[t]^2), hb[0] == hbzero,
htau'[
t] == \[Kappa] htau[
t] (4 htau[t]^2 + 3 hb[t]^2 - 3 g2[t]^2 - 9/5 g1[t]^2),
htau[0] == htauzero,
smu'[t] == \[Kappa] smu[
t] (3 ht[t]^2 + 3 hb[t]^2 - 3 g2[t]^2 - 3/5 g1[t]^2),
smu[0] == susymu, g3'[t] == -3 \[Kappa] g3[t]^3,
g3[0] == Sqrt[4 \[Pi] 0.118],
g2'[t] == \[Kappa] g2[t]^3, g2[0] == Sqrt[4 \[Pi]/30],
g1'[t] == 33 \[Kappa]/5 g1[t]^3, g1[0] == Sqrt[4 \[Pi]/60],
M3'[t] == -6 \[Kappa] g3[t]^2 M3[t], M3[Qmax] == gaugino,
M2'[t] == 2 \[Kappa] g2[t]^2 M2[t], M2[Qmax] == gaugino,
M1'[t] == 66/5 \[Kappa] g1[t]^2 M1[t], M1[Qmax] == gaugino,
mhu2'[t] == \[Kappa] (3 Xu[t] - 6 g2[t]^2 M2[t]^2 -
6/5 g1[t]^2 M1[t]^2 + 3/5 g1[t]^2 strm2[t]),
mhu2[Qmax] == higgsmassparam^2,
mhd2'[
t] == \[Kappa] (3 Xd[t] + Xtau[t] - 6 g2[t]^2 M2[t]^2 -
6/5 g1[t]^2 M1[t]^2 - 3/5 g1[t]^2 strm2[t]),
mhd2[Qmax] == higgsmassparam^2,
mq2'[t] == \[Kappa] (Xu[t] + Xd[t] - 32/3 g3[t]^2 M3[t]^2 -
6 g2[t]^2 M2[t]^2 - 2/15 g1[t]^2*M1[t]^2 +
1/5 g1[t]^2 strm2[t]), mq2[Qmax] == msf^2,
mu2'[t] == \[Kappa] (2 Xu[
t] - (32/3) g3[t]^2 M3[t]^2 - (32/15) g1[t]^2 M1[t]^2 -
4/5 g1[t]^2 strm2[t]), mu2[Qmax] == msf^2,
md2'[t] ==
2 \[Kappa] (Xd[
t] - (32/3) g3[t]^2 M3[t]^2 - (8/15) g1[t]^2 M1[t]^2 +
2/5 g1[t]^2 strm2[t]), md2[Qmax] == msf^2,
mL2'[t] == \[Kappa] (Xtau[t] -
6 g2[t]^2 M2[t]^2 - (6/5) g1[t]^2 M1[t]^2 -
3/5 g1[t]^2 strm2[t]), mL2[Qmax] == msf^2,
me2'[t] == \[Kappa] (2 Xtau[t] - (24/5) g1[t]^2 M1[t]^2 +
6/5 g1[t]^2 strm2[t]), me2[Qmax] == msf^2,
at'[t] == \[Kappa] (at[
t] (18*ht[t]^2 + hb[t]^2 - (16/3) g3[t]^2 -
3*g2[t]^2 - (13/15)*g1[t]^2) + 2 ab[t] hb[t] ht[t] +
ht[t] ((32/3) g3[t]^2 M3[t] +
6 g2[t]^2 M2[t] + (26/15) g1[t]^2 M1[t])), at[Qmax] == at0,
ab'[t] == \[Kappa] (ab[
t] (18 hb[t]^2 + ht[t]^2 + htau[t]^2 - (16/3) g3[t]^2 -
3 g2[t]^2 - (7/15) g1[t]^2) + 2 at[t] ht[t] hb[t] +
2 atau[t] htau[t] hb[t] +
hb[t] ((32/3) g3[t]^2 M3[t] +
6 g2[t]^2 M2[t] + (14/15) g1[t]^2 M1[t])),
ab[Qmax] == ab0,
atau'[
t] == \[Kappa] (atau[
t] (12 htau[t]^2 + 3 hb[t]^2 - 3 g2[t]^2 - (9/5) g1[t]^2) +
6 ab[t] hb[t] htau[t] +
htau[t] (6 g2[t]^2 M2[t] + (18/15) g1[t]^2 M1[t])),
ab[Qmax] == atau0,
b'[t] == \[Kappa] (b[
t] (3 ht[t]^2 + 3 hb[t]^2 + htau[t]^2 - 3 g2[t]^2 -
3/5 g1[t]^2) +
smu[t] (6 at[t] ht[t] + 6 ab[t] hb[t] + 2 atau[t] htau[t] +
6 g2[t]^2 M2[t] + 6/5 g1[t]^2 M1[t])), b[Qmax] == b0},
{g3[t], g2[t], g1[t], ht[t], hb[t], htau[t], smu[t], M1[t], M2[t],
M3[t], mhu2[t], mhd2[t], mq2[t], mu2[t], md2[t], mL2[t], me2[t],
at[t], ab[t], atau[t], b[t]}, {t, 0, Qmax}];

g3bis[t_] = g3[t] /. system;
g2bis[t_] = g2[t] /. system;
g1bis[t_] = g1[t] /. system;

htbis[t_] = ht[t] /. system;
hbbis[t_] = hb[t] /. system;
htaubis[t_] = htau[t] /. system;
smubis[t_] = smu[t] /. system;

M3bis[t_] = M3[t] /. system;
M2bis[t_] = M2[t] /. system;
M1bis[t_] = M1[t] /. system;

mhubis[t_] = Sqrt[mhu2[t]] /. system;
mhdbis[t_] = Sqrt[mhd2[t]] /. system;

mQbis[t_] = Sqrt[mq2[t]] /. system;
mubis[t_] = Sqrt[mu2[t]] /. system;
mdbis[t_] = Sqrt[md2[t]] /. system;
mLbis[t_] = Sqrt[mL2[t]] /. system;
mebis[t_] = Sqrt[me2[t]] /. system;

atbis[t_] = at[t] /. system;
abbis[t_] = ab[t] /. system;
ataubis[t_] = atau[t] /. system;
bbis[t_] = Sqrt[b[t]] /. system;

Plot[{g1bis[z], g2bis[z], g3bis[z]}, {z, 6.35, Qmax},
AxesOrigin -> {6.35, 0.4},
PlotLabel ->
"runing gauge couplings \!$$\*SubscriptBox[\(g$$, $$s$$]\), \
\!$$\*SubscriptBox[\(g$$, $$w$$]\), g'",
AxesLabel -> {"ln[Q/\!$$\*SubscriptBox[\(m$$, $$t$$]\)]"}]
Plot[{htbis[z], hbbis[z], htaubis[z]}, {z, 0, Qmax},
AxesOrigin -> {0, 0.4},
PlotLabel ->
"runing Yukawa couplings \!$$\*SubscriptBox[\(y$$, $$t$$]\), \
\!$$\*SubscriptBox[\(y$$, $$\(b$$$$,$$\)]\) \!$$\*SubscriptBox[\(y$$, \
$$\[Tau]$$]\)",
AxesLabel -> {"ln[Q/\!$$\*SubscriptBox[\(m$$, $$t$$]\)]"}]
Plot[{M1bis[z], M2bis[z], M3bis[z], mhubis[z], mhdbis[z], mQbis[z],
mubis[z], mdbis[z], mLbis[z], mebis[z]}, {z, 6.35, Qmax},
AxesOrigin -> {6.35, 0},
PlotLabel ->
"runing soft mass parameters \!$$\*SubscriptBox[\(M$$, $$1, 2, 3$$]\
\), \!$$\*SubscriptBox[\(m$$, $$hu, hd$$]\), \
\!$$\*SubscriptBox[\(m$$, $$Q, t, b, L, \[Tau]$$]\)",
AxesLabel -> {"ln[Q/\!$$\*SubscriptBox[\(m$$, $$t$$]\)]"}]

-
This could be quite useful. It would be really good if you made it clearer how to use it (e.g. set the parameters at the GUT scale) and validated it by e.g. comparing against SOFTSUSY or by reproducing a plot in Martin's SUSY primer. – innisfree Jan 27 at 5:18
See my updated version. I refined the code and add some illustrations. Now it's able to reproduce plots in Martin's SUSY primer. – Di Liu Jan 29 at 0:22
Cool. I just checked the code myself. It worked fine and I put in the plots generated into the answer. – innisfree Jan 29 at 3:51