Models with Thresholds in SPheno
General
In the SPheno Output of SARAH it is possible to include thresholds in the RGE evaluation at which heavy degrees of freedoms are removed from the spectrum. The implementation depends if the thresholds are also associated with the breaking of an additional gauge symmetry or now.
Thresholds without gauge symmetry breaking
Assumptions
If all scales have the same gauge structure, it is possible for SARAH to derive the RGEs for all scales from the RGEs for the highest scale by performing the following steps:
 For those fields which should be integrated out during the run new variables n_{gen}(Φ_{i}) are introduced, which define the number of generation of the heavy field Φ_{i}. All gauge group constants like the the Dynkin index summed over chiral superfields, S(R), are expressed as function of n_{gen}(Φ_{i}). These n_{gen}(Φ_{i}) are dynamically adjusted, when the energy scale crosses a threshold.
 It is also necessary to set the couplings which involve heavy fields to zero when a threshold is crossed. For example, the Yukawa type coupling of the formY^{i**j}Φ_{i}ϕ_{j}H involves three generations of the heavy fieldΦ. At the threshold ofΦ_{k}, thekth row ofY is set to zero. That happens similarly for all other superpotential and softbreaking parameters.
 It is assumed, that the masses of the scalar and fermionic component of a heavy superfield are the same, i.e. the masses are much larger than the softbreaking masses. Furthermore, it is assumed that the masses are given by a bilinear superpotential term.
SARAH included automatically the threshold corrections to gauge couplings and gaugino masses in the SPheno output.
Implementation
To include thresholds without gauge symmetry breaking, the following steps have to be performed:

The heavy fields must be deleted in the SARAH model definition: DeleteFields = {...};
This ensures, that the these particle are not take into account for the calculation of decays or loop corrections at the SUSY scale.

The thresholds have to be defined in
SPheno.m
: Thresholds = {{Scale1, {HeavyFields1}}, {Scale2, ... }};For all scales an entry in the array
Thresholds
has to be added. Each entry defines first the threshold scale, at second position a list with the heavy superfields is given. Also specific generations for a superfield can be given.
It is possible to define boundary conditions at each threshold scale for running up and down separately:
BoundaryConditionsUp[[/xx]] = { ...};
BoundaryConditionsDown[[/xx]] = { ...};
Example
The simplest example for a SUSY model with thresholds is the seesaw I version of the MSSM. The changes compared to the MSSM are as follows:

Note the presence of the Weinberg operator
WOp

The righthanded fields are deleted from the SUSY scale Lagrangian DeleteParticles={v};

The thresholds in SPheno.m are defined Thresholds={ {Abs[MvIN[1,1]],{v[1]}}, {Abs[MvIN[2,2]],{v[2]}}, {Abs[MvIN[3,3]],{v[3]}} };

When running down, boundary conditions to calculate the shift in the Weinberg operator are added for each threshold scale BoundaryConditionsUp=Table[{},{Length[Thresholds]}]; BoundaryConditionsDown=Table[{},{Length[Thresholds]}];
BoundaryConditionsDown[[/11]]={ {WOp[index1,index2],WOp[index1,index2]  Yv[1,index1] Yv[1,index2]/MassOfv[1]} }; BoundaryConditionsDown[[/22]]={ {WOp[index1,index2],WOp[index1,index2]  Yv[2,index1] Yv[2,index2]/MassOfv[2]} }; BoundaryConditionsDown[[/33]]={ {WOp[index1,index2],  Yv[3,index1] Yv[3,index2]/MassOfv[3]} };
Thresholds with gauge symmetry breaking
If the gauge structure at the different scales are different, each set of RGEs is calculated separately and this information is then combined into one consistent version of SPheno . This code includes routines for calculating finite shifts in the gauge couplings and gaugino mass parameters. As an example or for instanceS**O(10) motivated leftright models.
Implementation
The implementation of models with a threshold scale that involved a gauge symmetry breaking is more involved. In general, the following steps are necessary:

For each regime, a separate model file for SARAH has to be created. These model file have to be saved in the subdirectories
Regime1
,Regime2
, …of [SARAH Directory]/Models/[Model]/They are numbered beginning with the highest scale.

The SPheno input file for the higher scales must provide the following information:
 IntermediateScale = True
 RegimeNr = X
 A list of the heavy fields, which should be integrated out, the gauge sector below the threshold as well as the corresponding quantum numbers of the fields integrated. That’s needed to calculate the finite shifts at the threshold scale. For instance, the entries might read HeavyFields = {Field_1, Field_2,..}; NextGauge = {U[1], SU[2], SU[2], SU[3]}; NextQN = { {Field_1, 0, 2, 1, 1}, {Field_2, 1/3, 1, 2, 4}, ... };

All necessary information for combining the regimes to one SPheno is given in
SPheno.m
of the lowest scale.
IntermediateScale = False

RegimeNr = X

The threshold scales: ThresholdScales = ...

In the boundary conditions
index1
,index2
, …can be used for defining sums over indices. 
The usual information for SPheno , defined in the sec. [sec:sphenoinputfile].

When starting the SPheno output of the lowest scale, automatically all other scales are evolved. Note, to calculate the RGEs of the different regimes requires SARAH to start one additional Mathematica kernel. For passing the information between the different Mathematica kernels a directory Meta
in the model directory is created by SARAH . Also the screen output of Mathematica during the evaluation of the higher regimes is written to that directory (OutputRegimeX.m
). So, the user can supervise the progress and see potential error messages. The necessary information of each regime for writing the combined source code for SPheno at the end is saved by SARAH in the files RegimeX.m
.
Example
The entire discussion in the following is based on the model OmegaShort
delivered with SARAH. This model has a threshold scale SU(2)_L \times SU(2)R \times U(1){BL} \rightarrow SU(2)_L \times U(1)_Y

'''Boundary conditions: ''' We don't further specify here all incredients of the model but just listen some boundary conditions which are important when matching the model to the MSSM. \begin{aligned} \label{eq:boundary_LR_1} Y_d = Y_Q^1 \cos \theta_1  Y_Q^2 \sin \theta_1 \thickspace, &\qquad& Y_u =  Y_Q^1 \cos \theta_2 + Y_Q^2 \sin \theta_2 \thickspace, \\ Y_e = Y_L^1 \cos \theta_1  Y_L^2 \sin \theta_1 \thickspace, &\qquad& Y_\nu =  Y_L^1 \cos \theta_2 + Y_L^2 \sin \theta_2 \thickspace,\end{aligned} <br: whereR = sin(θ_{1} − θ_{2}). For the softtrilinear couplings,Y is replaced withT in the expressions. For the sfermionic soft masses, we have \begin{aligned} m_{q}^2 =m_{u^c}^2 = m_{d^c}^2 &=& m_{Q^c}^2 \thickspace,\\ m_{l}^2=m_{e^c}^2 &=& m_{L^c}^2 \thickspace, \\ M_L = M_R &=& M_2 \thickspace.\end{aligned} while in the Higgs sector we need the relations \begin{aligned} m_{H_d}^2 &=& \cos^2 \theta_1 (m_\Phi^2){11} + \sin^2 \theta_1 (m\Phi^2){22}  \sin \theta_1 \cos \theta_1 \left[ (m\Phi^2){12} + (m\Phi^2){21} \right] \thickspace, \\ m{H_u}^2 &=& \cos^2 \theta_2 (m_\Phi^2){11} + \sin^2 \theta_2 (m\Phi^2){22}  \sin \theta_2 \cos \theta_2 \left[ (m\Phi^2){12} + (m\Phi^2){21} \right] \thickspace, \end{aligned} In the gauge sector, we have to express the hypercharge coupling and the corresponding gaugino through \begin{aligned} g_1 & = & \frac{\sqrt{5} g_2 g{BL}}{\sqrt{2 g_2^2 + 3 g_{BL}^2}} \thickspace , \\ \label{eq:boundary_LR_2} M_1 & = & \frac{2 g_2^2 M_{BL} + 3 g_{BL}^2 M_R}{2 g_2^2 + 3 g_{BL}^2} \thickspace.\end{aligned}

Defining the SARAH models: one needs two model files
 Above the breaking scale the full model
 Below the breaking scale, only the particles and gauge groups of the MSSM survive as dynamical degrees of freedom. Therefore, the model file is almost identical that for the MSSM with only tiny differences:

The superpotential contains the Weinberg operator SuperPotential = ... + WOp/2 l.Hu.l.Hu;

The neutrinos are massive and mix among each other DEFINITION[EWSB][MatterSector]= { ...,{{FvL}, {FV, UV}}};

Input file for SPheno

AboveS**U(2)_{R} × U(1)_{B − L} breaking scale: First, the regime must be flagged as an intermediate scale, and SARAHmust be told its position in the scale hierarchy (counted from GUT to low scale).
RegimeNr = 1; IntermediateScale = True;
Afterwards, we give a list with all particles which are integrated out at the threshold scale after gauge symmetry breaking.
HeavyFields = {Hpm1R1, ChiR1, Cha1r1, hhR1, AhR1, FvR1, SVRr1, SH0r1[3], SHCr1[3], SO1r1, SDLpR1, SDLppR1, SDL0r1, SDRmmR1, DR3r1, DL1r1, DL2r1, DL3r1, H0r1, HCr1};
The numbers in square brackets indicate that only the third generation and above is integrated out.
In order to calculate the finite shifts of the gauge couplings and gaugino masses, it is necessary to define the gauge sector of the next scale
NextGauge
as well as the quantum number of the fields which are integrated out with respect to those gauge groups.NextGauge= {U[1], SU[2], SU[3]}; NextQN = { {Hpm1R1, 1, 1, 1}, {ChiR1, 0, 1, 1}, {Cha1r1, 1, 1, 1}, {hhR1, 0, 1, 1}, {AhR1, 0, 1, 1}, {FvR1, 0, 1, 1}, {SVRr1, 0, 1, 1}, {SH0r1, 1/2, 1, 1}, {SHCr1, 1/2, 2, 1}, {SO1r1, 0, 1, 1}, {SDLpR1, 1, 1, 1}, {SDLppR1, 2, 1, 1}, {SDL0r1, 1, 3, 1}, {SDRmmR1, 2, 1, 1}, {DR3r1, 2, 1, 1}, {DL1r1, 1, 1, 1}, {DL2r1, 1, 1, 1}, {DL3r1, 1, 3, 1}, {H0r1, 1/2, 1, 1}, {HCr1, 1/2, 2, 1} };
Finally, SARAHneeds information on the vacuum conditions. There are two different VEVs, and therefore, we need to choose two parameters which are fixed by the tadpole equations. As there is no closed analytical solution for them, we give an approximation obtained by neglecting the softbreaking terms.
ParametersToSolveTadpoles = {Mdelta, Momega}; UseGivenTadpoleSolution = True; SubSolutionsTadpolesTree={ Mdelta >  SignumMdelta ac1 vR/Sqrt[2], Momega >  SignumMomega ac1 vBL^2/(2 Sqrt[2] vR) }; SubSolutionsTadpolesLoop={};

BelowS**U(2)_{R} × U(1)_{B − L} breaking scale: The second scale is not an intermediate scale, and hence RegimeNr = 2; IntermediateScale = False;
We make the following choice of free parameters of that model: to the set of standard mSugra parameters (m_0,M_{1/2},A_0,\tan\beta,\text{sign}\mu), we addB_{0}, the superpotential parametera, the signs ofM_{Ω} andM_{Δ}, the two VEVsv_{R} andv_{B**L}, and the threshold scale. These are defined in the blocks
MINPAR
andEXTPAR
.MINPAR= { {1, m0}, {2, m12}, {3, TanBeta}, {4, SignumMu}, {5, Azero}, {6, Bzero}, {7, SignumMomega}, {8, SignumMdelta}, {9, aInput}}; EXTPAR = { {100, vRinput}, {101, vBLinput}, {200, TScale}};
As in the, MSSM we fix the numerical values ofμ andB_{μ} by the solving the two tadpole equations,
ParametersToSolveTadpoles = {\[Mu],B[\[Mu]]};
Furthermore, we use also the common definitions for the SUSY scale, see also sec. [sec:Example_{S}Pheno_{M}SSM].
RenormalizationScaleFirstGuess = m0^2 + 4 m12^2; RenormalizationScale = MSu[1]*MSu[6];
We use as condition for the GUT scale
ConditionGUTscale = {gBL==g2, g1 == g2};
The second entry is necessary because it might be that the masses of the fields breaking the leftright symmetry are above the unification scale which has to be handled seperatly. At the GUT scale we use boundary conditions motivated by minimal supergravity: all scalar softbreaking masses are proportional tom_{0}, the gaugino masses are proportional toM_{1/2}, the trilinear softbreaking couplings are given by the corresponding superpotential parameter timesA_{0}, and the bilinear softbreaking couplings are set toB_{0} times the superpotential parameter. The values for the coupling matricesf,α as well as forμ are read from the LesHouches input file.
BoundaryHighScale={ {T[YQ], Azero*YQ}, {T[YL], Azero*YL}, {f, LHInput[f]}, {T[f], Azero*LHInput[fm]}, {AlphaOm, LHInput[AlphaOm]}, {T[AlphaOm], Azero*LHInput[AlphaOm]}, {T[a], Azero*aInput}, {B[Mdelta], Bzero*Mdelta}, {B[Momega], Bzero*Momega}, {B[Mu3], Bzero*LHInput[Mu3]}, {mqL2, DIAGONAL m0^2}, {mqR2, DIAGONAL m0^2}, {mlL2, DIAGONAL m0^2}, {mlR2, DIAGONAL m0^2}, {mPhi2, DIAGONAL m0^2}, {mdeltaL2, m0^2}, {mdeltaLbar2, m0^2}, {mdeltaR2, m0^2}, {mdeltaRbar2, m0^2}, {momegaL2, m0^2}, {momegaR2, m0^2}, {MassB, m12}, {MassWL, m12}, {MassG, m12} };
To glue the both regimes, we need to define the appropriate boundary conditions. First, we initialize the arrays
ThresholdScales = {TSCALE}; BoundaryConditionsUp = Table[{},{Length[ThresholdScales]}]; BoundaryConditionsDown = Table[{},{Length[ThresholdScales]}];
and then encode the equations eqs. ([eq:boundary_{L}R_{1}])([eq:boundary_{L}R_{2}]). In order to keep the code short, we define
ST1 = Sin[Theta1]; CT1 = Cos[Theta1]; ST2 = Sin[Theta2]; CT2 = Cos[Theta2]; ST21 = Sin[Theta2Theta1]; CT21 = Cos[Theta2Theta1];
Using these abbreviations, the boundary conditions can be written as
BoundaryConditionsUp[[/11]] = { {YQ[index1,index2,1], (Yu[index1,index2] ST1 + Yd[index1,index2]ST2)/ST21 }, {YQ[index1,index2,2], (Yu[index1,index2] CT1 + Yd[index1,index2]CT2)/ST21 }, {YL[index1,index2,1], (Yv[index1,index2] ST1 + Ye[index1,index2]ST2)/ST21 }, {YL[index1,index2,2], (Yv[index1,index2] CT1 + Ye[index1,index2]CT2)/ST21 }, {gBL, Sqrt[2] g1 g2 /Sqrt[5 g2^2 3 g1^2]}, {Yv, LHInput[Yv]} }; BoundaryConditionsDown[[/11]] = { {vR, vRinput}, {vBL, vBLinput}, {a, aInput}, {Theta1, ArcTan[RealPart[((vR*AlphaOm[1,2])/2 + Mu3[1,2])/Mu3[2,2]]]}, {Theta2, ArcTan[RealPart[((vR*AlphaOm[1,2])/2 + Mu3[1,2])/Mu3[2,2]]]}, {g1, Sqrt[5] g2 gBL/Sqrt[2 g2^2 + 3 gBL^2]}, {MassB, (2 g2^2 MassB + 3 gBL^2 MassWL)/(2 g2^2 + 3 gBL^2)}, {MassWB, MassWL}, {Yd[index1,index2], YQ[index1,index2,1] CT1  YQ[index1,index2,2] ST1}, {Yu[index1,index2],  YQ[index1,index2,1] CT2 + YQ[index1,index2,2] ST2}, {Ye[index1,index2], YL[index1,index2,1] CT1  YL[index1,index2,2] ST1}, {Yv[index1,index2],  YL[index1,index2,1] CT2 + YL[index1,index2,2] ST2}, {T[Yd][index1,index2], T[YQ][index1,index2,1] CT1 T[YQ][index1,index2,2]ST1}, {T[Yu][index1,index2], T[YQ][index1,index2,1] CT2+T[YQ][index1,index2,2] ST2}, {T[Ye][index1,index2], T[YL][index1,index2,1] CT1 T[YL][index1,index2,2] ST1}, {mu2, mqR2}, {md2, mqR2}, {mq2, mqR2}, {me2, mlR2}, {ml2, mlR2}, {mHd2, CT1^2 mPhi2[1,1] + ST1^2 mPhi2[2,2]  ST1 CT1(mPhi2[1,2] + mPhi2[2,1])}, {mHu2, CT2^2 mPhi2[1,1] + ST2^2 mPhi2[2,2]  ST2 CT2(mPhi2[1,2] + mPhi2[2,1])}, {WOp, MatMul2[MatMul2[Yv,InverseMatrix[f],FortranFalse],Transpose[Yv], FortranFalse]/vR} };
Note thatY_{ν} does not appear in any of the model files and it is therefore necessary to fix the dimension of that matrix by hand
AdditionalVariablesSPheno={Yv[3,3]};
Several parameters are restricted to be real. In addition, it is helpful to choose initialization values for some parameters to stabilize the numerics in the first iteration
RealParameters = {TanBeta, vRinput,vBLinput,Theta1,Theta2,TScale}; InitializationValues = { {Mu3IN[1,1], (Mu3IN[1,2]^2  AlphaOmIN[1,2]^2 vRInput^2/4)/Mu3IN[2,2]}, {Theta1, ArcTan[RealPart[(Mu3IN[1,2]+AlphaOmIN[1,2]vRInput/2)/Mu3IN[2,2]]]}, {Theta2, ArcTan[RealPart[(Mu3IN[1,2]AlphaOmIN[1,2]vRInput/2)/Mu3IN[2,2]]]}, {Mdelta, aInput*SignumMdelta*vRinput/2 }, {Momega, SignumMomega*(aInput^2*vBLinput^2)/(8 Mdelta)} };


Getting the SPheno code: One needs to start the MSSM part of the model and run
MakeSPheno
. SARAH will start automatically a second Mathematica kernel to generate all necessary information for the leftright symmetric phase.