# 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
show 1 more comment