RGEs of the MSSM - problems with Mathematica [closed]

Question

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.

Answer 1

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\)]\)]"}]