Manual 5.0
User Manual: Pdf
Open the PDF directly: View PDF 
.
Page Count: 56
| Download | |
| Open PDF In Browser | View PDF |
Date: May 15, 2018
The user’s manual, version 5.0
D. Barducci1 , G. Bélanger2 , J. Bernon3 , F. Boudjema2 , J. Da Silva2 ,
A. Goudelis4 , S. Kraml5 , U. Laa6 A. Pukhov7 , A. Semenov8 ,
B. Zaldivar2,9 .
1) SISSA and INFN, Sezione di Trieste, via Bonomea 265, 34136 Trieste, Italy
2) UGA, USMB, CNRS, LAPTh, F-74940 Annecy, France
3) Institute for Advanced Studies, The Hong Kong University of Science and
Technology, Clear Water Bay, Kowloon, Hong Kong S.A.R, China
4)LPTHE, Sorbonne Universités, UPMC, CNRS, F-75252 Paris Cedex, France
5) Laboratoire de Physique Subatomique et de Cosmologie, Université Grenoble-Alpes,
CNRS/IN2P3, 53 Avenue des Martyrs, F-38026 Grenoble, France
6) Monash University, Melbourne, Victoria 3800 Australia.
7) Skobeltsyn Inst. of Nuclear Physics, Moscow State Univ., Moscow 119992, Russia
8) Joint Institute for Nuclear Research (JINR) 141980, Dubna, Russia
9) Departamento de Fı́sica Teórica, UAM, 28049 Madrid, Spain
Abstract
We give an up-to-date description of the micrOMEGAs functions. Only the routines which are available for the users are described. Examples on how to use these
functions can be found in the sample main programs distributed with the code.
Contents
1 Introduction
3
2 Discrete symmetry in micrOMEGAs.
3
3 Downloading and compilation of micrOMEGAs.
3.1 File structure of micrOMEGAs. . . . . . . . . . . . .
3.2 Compilation of CalcHEP and micrOMEGAs routines. .
3.3 Module structure of main programs. . . . . . . . . .
3.4 Compilation of codes for specific models. . . . . . .
3.5 Command line parameters of main programs. . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
4 Global Parameters and constants
4
4
5
6
7
7
8
5 Setting of model parameters, spectrum calculation, parameter display. 10
6 Relic density calculation.
6.1 Switches and auxilary routines . . . . . . . . . . . . . . . . . . . . .
6.2 Calculation of relic density for one-component Dark Matter models.
6.3 Calculation of relic density for two-component Dark Matter models.
6.4 Calculation of relic density for freeze-in. . . . . . . . . . . . . . . .
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
12
12
13
16
17
7 Direct detection.
7.1 Amplitudes for elastic scattering . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Scattering on nuclei . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 Auxiliary routines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
19
20
21
8 Indirect detection
8.1 Interpolation and display of spectra . .
8.2 Annihilation spectra . . . . . . . . . .
8.3 Distribution of Dark Matter in Galaxy.
8.4 Photon signal . . . . . . . . . . . . . .
8.5 Propagation of charged particles. . . .
21
21
22
23
24
25
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
9 Neutrino signal from the Sun and the Earth
25
9.1 Comparison with IceCube results . . . . . . . . . . . . . . . . . . . . . . . 27
10 Cross sections and decays.
27
11 Tools for model independent analysis
29
12 Constraints from colliders
12.1 The Higgs sector . . . . . . . . . . . . . . . .
12.1.1 Lilith . . . . . . . . . . . . . . . . . . .
12.1.2 HiggsBounds and HiggsSignals . . . . .
12.1.3 Automatic generation of interface files
12.2 Searches for New particles . . . . . . . . . . .
12.2.1 SModelS . . . . . . . . . . . . . . . . .
12.2.2 Other limits . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
31
31
31
32
33
34
34
36
13 Additional routines for specific
13.1 MSSM . . . . . . . . . . . . .
13.2 The NMSSM . . . . . . . . .
13.3 The CPVMSSM . . . . . . . .
13.4 The UMSSM . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
36
37
39
39
40
.
.
.
.
.
.
.
41
41
42
42
43
43
46
46
models
. . . . .
. . . . .
. . . . .
. . . . .
.
.
.
.
14 Tools for new model implementation.
14.1 Main steps . . . . . . . . . . . . . . . . .
14.2 Automatic width calculation . . . . . . .
14.3 Using LanHEP for model file generation.
14.4 QCD functions . . . . . . . . . . . . . .
14.5 SLHA reader . . . . . . . . . . . . . . .
14.5.1 Writing an SLHA input file . . .
14.6 Routines for diagonalisation. . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
15 Mathematical tools.
47
A An updated routine for b → sγ in the MSSM
49
2
1
Introduction
micrOMEGAs is a code to calculate the properties of cold dark matter (CDM) in a generic
model of particle physics. First developed to compute the relic density of dark matter, the
code also computes the rates for dark matter direct and indirect detection. micrOMEGAs
computes CDM properties in the framework of a model of particle interactions presented in
CalcHEP format [1]. It is assumed that the model is invariant under a discrete symmetry
like R-parity (even for all standard particles and odd for some new particles including
the dark matter candidate) which ensures the stability of the lightest odd particle (LOP).
Similarly in two-component dark matter models, a discrete symmetrie that guarantees the
stability of the lightest particle in each of the two dark matter sectors is assumed. The
CalcHEP package is included in micrOMEGAs and used for matrix elements calculations.
All annihilation and coannihilation channels are included in the computation of the relic
density. This manual gives an up-to-date description of all micrOMEGAs functions. The
methods used to compute the different dark matter properties are described in references
[2–7]. These references also contain a more complete description of the code. In the
following the cold dark matter candidate also called LOP or weakly-interactive massive
particle (WIMP) will be denoted by χ. Starting with version 5.0, micrOMEGAs also allows
to compute the abundance of feebly interacting dark matter candidates (FIMP) through
the freeze-in mechanism [8]
micrOMEGAs contains both C and Fortran routines. Below we describe only the Cversion of the routines, in general we use the same names and the same types of argument
for both C and Fortran functions. We always use double(real*8) variables for float point
numbers and int(INTEGER) for integers. In this manual we use FD for file descriptor
variables, the file descriptors are FILE* in C and channel number in Fortran. For Cfunctions which return values different from integer and real number, the first parameter
in the corresponding Fortran subroutine presents the return value. The symbol & before
the names of variables in C-functions stands for the address of the variable. It is used
for output parameters. In Fortran calls there is no need for & since all parameters are
passed via addresses. In C programs one can substitute NULL for any output parameter
which the user chooses to ignore. In Fortran one can substitute cNull, iNull, r8Null
for unneeded parameters of character, integer and real*8 type respectively.
A few C-functions use pointer variables that specify an address in the computer memory. Because pointers do not exist in Fortran, we use an INTEGER*8 variable whose length
is sufficient to store a computer address.
The complete format for all functions can be found in include/.h (for C) or include/_f.h
(for Fortran). Examples on how to use these functions are provided in the MSSM/main.c[F]
file.
2
Discrete symmetry in micrOMEGAs.
micrOMEGAs exploits the fact that models of dark matter exhibit a discrete symmetry
and that the fields of the model transform as φ → ei2πXφ φ where the charge |Xφ | < 1.
The particles of the Standard Model are assumed to transform trivially under the discrete
symmetry, Xφ = 0. In the following all particles with charge Xφ 6= 0 will be called odd and
the lightest odd particle will be stable. If neutral, it can be considered as a DM candidate.
3
Typical examples of discrete symmetries used for constructing single DM models are Z2
and Z3 . Multi-component DM can arise in models with larger discrete groups. A simple
example is a model with Z2 ×Z2′ symmetry, the particles charged under Z2 (Z2′ ) will belong
to the first (second) dark sector. The lightest particle of each sector will be stable and
therefore a potential DM candidate. Another example is a model with a Z4 symmetry.
The two dark sectors contain particles with Xφ = ±1/4 and Xφ = 1/2 respectively. The
lightest particle with charge 1/4 is always stable while the lightest particle of charge 1/2
is stable only if its decay into two particles of charge 1/4 is kinematically forbidden.
micrOMEGAs assumes that all odd particles have names starting with ’~’, for example,
~o1 for the lightest neutralino. In versions 4.X, to distinguish the particles with different
transformation properties with respect to the discrete group, that is particles belonging
to different ’dark’ sectors, we use the convention that the names of particles in the second
’dark’ sector starts with ’~~’. Note that micrOMEGAs does not check the symmetry of the
Lagrangian, it assumes that the name convention correctly identifies all particles with the
same discrete symmetry quantum numbers. For models with FIMPs, new particles are
considered to be in thermal equilibrium with the SM bath (B) unless explicitly defined as
being feeble, ie belonging to F. Both B and F can contain odd or even particles.
3
Downloading and compilation of micrOMEGAs.
To download micrOMEGAs, go to
http://lapth.cnrs.fr/micromegas
and unpack the file received, micromegas_5.0.tgz, with the command
tar -xvzf micromegas_5.0.tgz
This should create the directory micromegas_5.0/ which occupies about 67Mb of disk
space. You will need more disk space after compilation of specific models and generation
of matrix elements. In case of problems and questions
email: micromegas@lapth.cnrs.fr
3.1
File structure of micrOMEGAs.
calc
calchep.ini
Makefile
README
CalcHEP_src/
Packages/
clean
man/
newProject
sources/
include/
lib/
MSSM model directory
MSSM/
calculator
specify the fonts for graphics in CalCHEP
to compile the kernel of the package
short description on how to run the code
generator of matrix elements for micrOMEGAs
external codes
to remove compiled files
contains the manual: description of micrOMEGAs routines
to create a new model directory structure
micrOMEGAs code
include files for micrOMEGAs routines or external codes
contains library micromegas.a when micrOMEGAs is compiled
4
Makefile
to compile the code and executable for this model
main.c[pp] main.F files with sample main programs
lib/
directory for routines specific to this model
Makefile
to compile the auxiliary code library lib/aLib.a
*.c *.f *.h *.inc
source codes of auxiliary functions
lanhep/
directory containing lanhep source model files
work/
CalcHEP working directory for the generation of
matrix elements
Makefile
to compile the library work/work aux.a
models/
directory for files which specifies the model
files *1.mdl are used in micrOMEGAs sessions. Other *.mdl files
are intended for CalcHEP interactive sessions
vars1.mdl free variables
func1.mdl constrained variables
prtcls1.mdl particles
lgrng1.mdl Feynman rules
tmp/
auxiliary directories for CalcHEP sessions
results/
so_generated/ storage of matrix elements generated by CalcHEP
Directories of other models which have the same structure as MSSM/
NMSSM/
Next-to-Minimal Supersymmetric Model [9, 10]
CPVMSSM/
MSSM with complex parameters [11, 12]
UMSSM/
U(1) extensions of the MSSM [13, 14]
IDM/
Inert Doublet Model [15]
LLL_singlet/
Simplified model with singlet charged lepton and real scalar DM
LHM/
Little Higgs Model [16]
SingletDM/
Singlet scalar DM model with Z2 symmetry [17]
Z3IDM/
Inert doublet model with Z3 discrete symmetry [18, 19]
Z4IDSM/
Inert doublet and singlet model with Z4 symmetry [18, 19]
ZpPortal/
Simplified model with a Z’ portal and fermion DM
mdlIndep/
For model independent computation of DM signals
Other models can be downloaded on the web, http://lapth.cnrs.fr/micromegas,
for example : RHNM, a right-handed Neutrino Model [20], SM4, a toy model with a 4th
generation of leptons and neutrino DM, as well as Z5M, a two scalar singlets with a Z5
symmetry model.
3.2
Compilation of CalcHEP and micrOMEGAs routines.
CalcHEP and micrOMEGAs are compiled by gmake. Go to the micrOMEGAs directory and
launch
gmake
If gmake is not available, then make should work like gmake. In principle micrOMEGAs
defines automatically the names of C and Fortran compilers and the flags for compilation. If you meet a problem, open the file which contains the compiler specifications,
CalcHEP_src/FlagsForSh, improve it, and launch [g]make again. The file is written in
sh script format and looks like
5
# C compiler
CC="gcc"
# Flags for C compiler
CFLAGS="-g -fsigned-char"
# Disposition of header files for X11
HX11=
# Disposition of lX11
LX11="-lX11"
# Fortran compiler
FC="gfortran"
FFLAGS="-fno-automatic"
........
After a successful definition of compilers and their flags, micrOMEGAs rewrites the file
FlagsForSh into FlagsForMake and substitutes its contents in all Makefiles of the package.
[g]make clean deletes all generated files, but asks permission to delete FlagsForSh.
[g]make flags only generates FlagsForSh. It allows to check and change flags
before compilation of codes.
3.3
Module structure of main programs.
Each model included in micrOMEGAs is accompanied with sample files for C and Fortran
programs which call micrOMEGAs routines, the main.c, main.F files. These files consist of
several modules enclosed between the instructions
#ifdef XXXXX
....................
#endif
Each of these blocks contains some code for a specific problem
#define
#define
#define
#define
#define
#define
#define
#define
#define
#define
MASSES_INFO
CONSTRAINTS
HIGGSBOUNDS
LILITH
SModelS
OMEGA
FREEZEIN
INDIRECT_DETECTION
LoopGAMMA
RESET_FORMFACTORS
#define CDM_NUCLEON
#define CDM_NUCLEUS
#define NEUTRINO
#define DECAYS
#define CROSS_SECTIONS
#define CLEAN
//Displays information about mass spectrum
//Displays B_>sgamma, Bs->mumu, etc
//calls HiggsBounds/HiggsSignal to constrain the Higgs sector
//calls LiLith to constrain the Higgs sector
//calls SModelS to constrain the new physics sector
//Calculates the relic density
//Calculates the relic density in the freeze-in mechanism
//Signals of DM annihilation in galactic halo
//Gamma-Ray lines - available only in some models
//Redefinition of Form Factors and other
//parameters
//Calculates amplitudes and cross-sections
//for DM-nucleon collisions
//Calculates number of events for 1kg*day
//and recoil energy distribution for various nuclei
//Calculates flux of solar neutrinos and
//the corresponding muon flux
//Calculates decay widths and branching ratios
//Calculates cross sections
//Removes intermediate files.
6
There is a flag
#define SHOWPLOTS
//which switches on graphic facilities of \micro.
Note that HiggsBounds and HiggsSinals are no longer included in the micrOMEGAs’s
distribution, they are uploaded from our website at first compilation when the option is
activated.
All these modules are completely independent. The user can comment or uncomment
any set of define instructions to suit his/her need.
3.4
Compilation of codes for specific models.
After the compilation of micrOMEGAs one has to compile the executable to compute DM
related observables in a specific model. To do this, go to the model directory, say MSSM,
and launch
[g]make
It should generate the executable main using the main.c source file. In general
gmake main=filename.ext
generates the executable filename based on the source file filename.ext. For ext we support 3 options: ’c’ , ’F’, ’cpp’ which correspond to C, FORTRAN and C++ sources. [g]make
called in the model directory automatically launches [g]make in subdirectories lib and
work to compile
lib/aLib.a - library of auxiliary model functions, e.g. constraints,
work/work_aux.a - library of model particles, free and dependent parameters.
3.5
Command line parameters of main programs.
The default versions of main.c/F programs need some arguments which have to be specified in command lines. If launched without arguments main explains which parameter
are needed. As a rule main needs the name of a file containing the numerical values of
the free parameters of the model. The structure of a file record should be
Name
Value # comment ( optional)
For instance, an Inert Doublet model (IDM) input file contains
Mh
MHC
MH3
MHX
la2
laL
125
200
200
63.2
0.01
0.01
#
#
#
#
#
#
mass of SM Higgs
mass of charged Higgs ~H+
mass of odd Higgs ~H3
mass of ~X particle
\lambda_2 coupling
0.5*(\lambda_3+\lambda_4+\lambda_5)
In other cases, different inputs can be required. For example, in the MSSM with input
parameters defined at the GUT scale, the parameters have to be provided in a command
line. Launching ./main will return
7
This program needs 4 parameters:
m0
common scalar mass at GUT scale
mhf common gaugino mass at GUT scale
a0
trilinear soft breaking parameter at GUT scale
tb
tan(beta)
Auxiliary parameters are:
sgn +/-1, sign of Higgsino mass term (default 1)
Mtp top quark pole mass
MbMb Mb(Mb) scale independent b-quark mass
alfSMZ strong coupling at MZ
Example: ./main 120 500 -350 10 1 173.1
4
Global Parameters and constants
The list of the global parameters and their default values are given in Tables 1, 2.
The numerical value for any of these parameters can be simply reset anywhere in the
code. The numerical values of the scalar quark form factors can also be reset by the
calcScalarQuarkFF routine presented below. Some physical values evaluated by micrOMEGAs
also are presented as global variables, see Table 3.
Table 1: Global input parameters of micrOMEGAs
Name
deltaY
K dif
L dif
Delta dif
Tau dif
Vc dif
Fermi a
Fermi b
Fermi c
Rsun
Rdisk
rhoDM
Vearth
Vrot
Vesc
default value units
0
0.0112
kpc2 /Myr
4
kpc
0.7
1016
s
0
km/s
0.52
fm
-0.6
fm
1.23
fm
8.5
kpc
20
kpc
0.3
GeV/cm3
225.2
km/s
220
km/s
600
km/s
comments
Difference between DM/anti-DM abundances
The normalized diffusion coefficient
Vertical size of the Galaxy diffusive halo
Slope of the diffusion coefficient
Electron energy loss time
Convective Galactic wind
nuclei surface thickness
parameters to set the nuclei radius with
RA = cA1/3 + b
Distance from the Sun to the center of the Galaxy
Radius of the galactic diffusion disk
Dark Matter density at Rsun
Galaxy velocity of the Earth
Galaxy rotation velocity at Rsun
Escape velocity at Rsun
All physical constants used in relic density calculations are defined in the file
include/micromegas_aux.h, they are listed in Table 4.
8
Table 2: Global parameters of micrOMEGAs: nucleon quark form factors
Proton
Name
ScalarFFPd
ScalarFFPu
ScalarFFPs
pVectorFFPd
pVectorFFPu
pVectorFFPs
SigmaFFPd
SigmaFFPu
SigmaFFPs
Neutron
value
0.0191
0.0153
0.0447
-0.427
0.842
-0.085
-0.23
0.84
-0.046
Name
ScalarFFNd
ScalarFFNu
ScalarFFNs
pVectorFFNd
pVectorFFNu
pVectorFFNs
SigmaFFNd
SigmaFFNu
SigmaFFNs
value
0.0273
0.011
0.0447
0.842
-0.427
-0.085
0.84
-0.23
-0.046
comments
Scalar form factor
Axial-vector form factor
Tensor form factor
Table 3: Evaluated global variables
Name
CDM1
CDM2
Mcdm1
Mcdm2
Mcdm
dmAsymm
fracCDM2
Tstart, Tend
units
character
character
GeV
GeV
GeV
GeV
comments
name of first DM particle
name of second DM particle
Mass of the first Dark Matter particle
Mass of the second DM particles
min(Mcdm1,Mcdm2) if both exist
Asymmetry between relic density of DM - DM
fraction of CDM2 in relic density.
Temperature interval
for solving the differential equation
9
Evaluated by
sortOddParticles
sortOddParticles
sortOddParticles
sortOddParticles
sortOddParticles
darkOmega[FO]
darkOmega2
darkOmega[2]
Name
MPlank
EntropyNow
RhoCrit100
Value
1.22091 × 1019
2.8912 × 109
10.537
Units
GeV
m−3
GeVm−3
Description
Planck mass
Present day entropy, s0
ρc /h2 or ρ for H = 100km/s/Mpc
Table 4: Some useful constants included in micrOMEGAs.
5
Setting of model parameters, spectrum calculation,
parameter display.
The independent parameters that characterize a given model are listed in the file
work/models/vars1.mdl. Three functions can be used to set the value of these parameters:
• assignVal(name,val)
• assignValW(name,val)
assign value val to parameter name. The function assignVal returns a non-zero value if
it cannot recognize a parameter name while assignValW writes an error message.
• readVar(fileName)
reads parameters from a file. The file should contain two columns with the following
format (see also Section 3.5)
name
value
readVar returns zero when the file has been read successfully, a negative value when the
file cannot be opened for reading and a positive value corresponding to the line where a
wrong file record was found.
Note that in Fortran, numerical constants should be specified as Real*8, for example
call assignValW(’SW’, 0.473D0)
A common mistake is to use Real*4.
The constrained parameters of the model are stored in work/models/func1.mdl.
Some of these parameters are treated as public parameters. The public parameters include
by default all particle masses and all parameters whose calculation requires external functions (except simple mathematical functions like sin, cos, ... ). The parameters needed
for the calculation of any public parameters in work/models/func1.mdl are also treated
as public. It is possible to enlarge the list of public parameters. There are two ways to do
this. One can type * before a parameter name to make it public or one can add a special
record in work/models/func1.mdl
%Local! |
Then all parameters listed above this record become public.
The calculation of the particle spectrum and of all public model constraints is done
with:
• sortOddParticles(txt)
which also sorts the odd particles with increasing masses. This routine fills the text
parameters CDM1 and CDM2 with the names of the lightest odd particle starting with one
10
and two tildes respectively and assigns the value of the mass of the lightest odd particle
in each sector to the global parameters Mcdm1 and Mcdm2. For models with only one
DM candidate, micrOMEGAs will set CDM2=NULL and Mcdm2=0 (in Fortran the string is
filled by space symbols). This routine returns a non zero error code for a wrong set of
parameters, for example parameters for which some constraint cannot be calculated. The
name of the corresponding constraint is written in txt. This routine has to be called after
a reassignment of any input parameter.
• qNumbers(pName, &spin2,&charge3,&cdim)
returns the quantum numbers for the particle pName. Here spin2 is twice the spin of
the particle; charge3 is three times the electric charge; cdim is the dimension of the
representation of SU (3)c , it can be 1, 3, −3, 6, −6 or 8. The parameters spin2, charge3,
cdim are variables of type int. The value returned is the PDG code. If pName does not
correspond to any particle of the model then qNumbers returns zero.
• pdg2name(nPDG)
returns the name of the particle which PDG code is nPDG. If this particle does not
exist in the model the return value is NULL. In the FORTRAN version this function is
Subroutine pdg2name(pName,nPDG) and the character variable pName consists of white
spaces if the particle does not exist in the model.
• antiParticle(pName)
returns the name of the anti-particle for the particle pName.
• pMass(pName)
returns the numerical value of the particle mass.
• nextOdd(n, &pMass)
returns the name and mass of the nth odd particle assuming that particles are sorted according to increasing masses. For n = 0 the output specifies the name and the mass of the
CDM candidate. In the FORTRAN version this function is Subroutine nextOdd(pName,n,pMass)
• findVal(name,&val)
finds the value of variable name and assigns it to parameter val. It returns a non-zero
value if it cannot recognize a parameter name.
• findValW(name)
returns the value of variable name and writes an error message if it cannot recognize a
parameter name.
The variables accessible by these two commands are all free parameters and the constrained parameters of the model (in file model/func1.mdl) treated as public.
The following routines are used to display the value of the independent and the constrained public parameters:
• printVar(FD)
prints the numerical values of all independent and public constrained parameters into FD
• printMasses(FD, sort)
prints the masses of ’odd’ particles (those whose names started with ~). If sort 6= 0 the
masses are sorted so the mass of the CDM is given first.
• printHiggs(FD, sort)
prints the masses and widths of ’even’ colorless scalars.
11
6
Relic density calculation.
6.1
Switches and auxilary routines
•VWdecay,VZdecay
Switches to turn on/off processes with off-shell gauge bosons in the final state for DM annihilation and particle decays. If VW/VZdecay=1, the 3-body final states will be computed
for annihilation processes only while if VW/VZdecay=2 they will be included in coannihilation processes as well. By default the switches are set to (VW/VZdecay=1).1 Note
that micrOMEGAs calculates the width of each particle only once and stores the result in
Decay Table. A second call to the function pWidth (whether an explicit call or within the
computation of a cross section) will return the same result even if the user has changed
the VW/VZdecay switch. We recommend to call
•cleanDecayTable()
after changing the switches to force micrOMEGAs to recalculate the widths taking into
account the new value of VW/VZdecay. In Fortran, the subroutine
• setVVdecay(VWdecay,VZdecay) changes the switches and calls cleanDecayTable().
The sortOddParticles command which must be used to recompute the particle spectrum after changing the model parameters also clears the decay table.
If the particle widths were stored in a SLHA file (Susy Les Houches Accord [21])
downloaded by micrOMEGAs, then the SLHA value will be used, the widths then do not
depend on the VW/VZdecay switches. To avoid downloading particle widths, one can use
slhaRead(fileName,mode=4) to read the content of the SLHA file, see the description in
Section 14.5.
The temperature dependence of the effective number of degrees of freedom can be set
with
• loadHeffGeff(char*fname)
allows to modify the temperature dependence of the effective number of degrees of freedom
by loading the file fname which contains a table of hef f (T ), gef f (T ) . A positive return
value corresponds to the number of lines in the table. A negative return value indicates
the line which creates a problem (e.g. wrong format), the routine returns zero when the
file fname cannot be opened. The default file is std_thg.tab and is downloaded automatically if loadHeffGeff is not called in the user’s main program [22]. Five other files are
provided in the sources/data directory: HP_A_thg.tab, HP_B_thg.tab, HP_B2_thg.tab,
HP_B3_thg.tab, and HP_C_thg.tab. They correspond to sets A, B, B2, B3, C in [23]. The
user can substitute his/her own table as well, if so, the file must contain three columns
containing the numerical values for T , hef f , gef f , the data file can also contain comments
for lines starting with #.
These functions are accessed via
• gEff(T)
returns the effective number of degrees of freedom for the energy density of radiation at
a bath temperature T, only SM particles are included.
• hEff(T)
returns the effective number of degrees of freedom for the entropy density of radiation at
a bath temperature T.
1
Including the 3-body final states can significantly increase the execution time for the relic density
computation.
12
• hEffLnDiff(T)
d log(hef f (T ))
.
returns the derivative of hef f with respect to the bath temperature, d log(T
)
•Hubble(T)
returns the Hubble expansion rate at a bath temperature T. This applies for the radiationdominated era and is valid for T & 100eV.
•improveCrossSection( p1,p2,p3,p4, Pcm, &cs)
allows to substitute a new cross-section for a given process. Here p1,p2 are the names
of particles in the initial state and p3,p4 those in the final state. Pcm is the center of
mass momentum and cs is the cross-section in [pb]. This function is useful if for example
the user wants to include her/his one-loop improved cross-section calculation in the relic
density computation.This function has to be written by the user. The corresponding code
can be placed in the directory lib. micrOMEGAs will call this routine substituting &cs
by the new calculated cross-section. The dummy version of this routine contained in
micrOMEGAs does not change the default cross section.
6.2
Calculation of relic density for one-component Dark Matter
models.
All routines to calculate the relic density in version 3 are available in further versions.
For these routines, the difference between the two dark sectors is ignored. These routines
are intended for models with either a Z2 or Z3 discrete symmetry.
• vSigma(T,Beps,fast)
calculates the thermally averaged cross section for DM annihilation times velocity at a
temperature T [GeV],
T
σv (T ) =
4
8π n(T )2
X
Z
√ X
s
p2α̃β̃ (s)gα̃ gβ̃
ds sK1
T
√
α̃,β̃
σα̃β̃→xy (s) +
x≥y
n(T ) =
1X
2
xγ̃
T X
mα̃
gα̃ m2α̃ K2 (
),
2
2π α̃
T
σα̃β̃→xγ̃ (s)
(1)
(2)
Here α̃, β̃, γ̃ is used for Odd particles and x,y for Even particles. σα̃β̃→x[γ̃/y] is the cross
section for the corresponding process averaged over the spins of incoming particles and
summed over the spins of outgoing particles. σv should represent the rate of disappearance of Odd particles, therefore when a final state particle has a non-zero decay branching
ratio to odd particles, the annihilation cross section for this process is multiplied by the
corresponding branching ratio of the decays into SM particles. For the same reason cross
sections for semi - annihilation processes contribute to vSigma with a factor 21 . K1 , K2 are
modified Bessel functions of the second kind, and mα̃ and gα̃ stand for the mass and the
number of degrees of freedom of particle α̃. Note, that if α̃ 6= β̃ then each σα̃β̃ term will
be presented twice. The value for σv is expressed in [pb · c]. The parameter Beps defines
the criteria for including coannihilation channels as for darkOmega described below. The
f ast = 1/0/ − 1 option switches between the fast/accurate/very accurate calculation.
13
The global array vSigmaTCh contains the contribution of different channels to vSigma.
vSigmaTCh[i].weight specifies the relative weight of the ith channel,
vSigmaTCh[i].prtcl[j] (j=0, 4) defines the particles names for the ith channel.
The last record in vSigmaTCh array has zero weight and NULL particle names. In the Fortran version, the function vSigmaTCh(i,weight,pdg,process) serves the same purpose.
This function returns 0 if i exceeds the number of annihilation channels and 1 otherwise,
i ≥ 1. The variable real*8 weight gives the relative contribution of each annihilation
channel and integer pdg(5) contains the codes of incoming and outgoing particles in the
annihilation process. character*40 process contains a textual description of annihilation
processes.
• vSigmaCC(T,cc,mode)
calculates the thermally averaged cross section × velocity for 2 → 2, 2 → 3, and 2 → 4
processes. T is the temperature in [GeV], cc is the address of the code for each process.
This address can be obtained by the function newProcess presented in Section 10. The
returned value is given in [c·pb].
If mode 6= 0, vSigmaCC calculates the contribution of a given process to the total
annihilation cross section, see Eq.1. The incoming particles should belong to the odd
sector. For 2 → 2 processes the result after summation over all subprocesses should
be identical to the one obtained via vSigma above. For this mode, vSigmaCC includes
combinatoric factors: 2 if α̃ 6= β̃, an additional factor 2 if the incoming state is not
self-conjugated, and a factor 21 for semi-annihilation.
If mode = 0, vSigmaCC is defined by the integral
√
Z
√
1
s
α̃β̃→X
ds sK1 ( )p2cm (s)σ α̃β̃→X (pcm (s))
< vσ
>T =
mβ̃
mα̃
2 2
T
2T mα̃ m K2 ( T )K2 ( T )
β̃
where pcm is the center of mass momentum of incoming particles. Note that
lim vSigmaCC(T, cc, 0) = lim σ(pcm )vrel (pcm )
pcm →0
T →0
where vrel (pcm ) is the relative velocity of incoming particles. The result of vSigmaCC can
be different from that of vSigma described above when there is an important contribution
from NLSP’s to the total number density of DM particles.
• darkOmega(&Xf,fast,Beps,&err)
calculates the dark matter relic density Ωh2 . This routine solves the differential evolution
equation using the Runge-Kutta method. Xf = M cdm/Tf characterizes the freeze-out
temperature which is defined
p by the condition Y (Tf ) = 2.5Yeq (Tf ). For asymmetric
DM this condition reads 2 Y + (Tf )Y − (Tf ) = 2.5Yeq (Tf ). The value of Xf is given for
information and is also used as an input for the routine that gives the relative contribution
of each channel to Ωh2 , see printChannels below. The f ast = 1 flag forces the fast
calculation (for more details see Ref. [3]). This is the recommended option and gives
an accuracy around 1%. The parameter Beps defines the criteria for including a given
coannihilation channel in the computation of the thermally averaged cross-section, [3].
The recommended value is Beps = 10−4 − 10−6 whereas if Beps = 1 only annihilation of
the lightest odd particle is computed. Non-zero error code means that the temperature
where thermal equilibrium between the DM and SM sectors is too large M cdm/T < 2 or
T > 105 GeV.
14
darkOmega solves the differential equation for the abundance Y (T ) in the temperature
interval [Tend,Tstart] defined by the conditions Y (Tstart ) ≈ 1.1Yeq (Tstart ), Y (Tend ) ≈
10Yeq (Tend ). For temperatures below Tend, the contribution for Yeq is neglected and the
differential equation is integrated explicitely. The solution in the interval [Tend,Tstart]
interval is tabulated and can be displayed via the function YF(T). The equilibrium abundance can be accessed with the function Yeq(T).
• darkOmegaFO(&Xf, fast, Beps)
calculates the dark matter relic density Ωh2 using the freeze-out approximation.
• printChannels(Xf,cut,Beps,prcnt,FD)
writes into FD the contributions of different channels to (Ωh2 )−1 . Here Xf is an input
parameter which should be evaluated first in darkOmega[FO]. Only the channels whose
relative contribution is larger than cut will be displayed. Beps plays the same role as in
the darkOmega[FO] routine. If prcnt 6= 0 the contributions are given in percent. Note
that for this specific purpose we use the freeze-out approximation.
• oneChannel(Xf,Beps,p1,p2,p3,p4)
calculates the relative contribution of the channel p1,p2 → p3,p4 to (Ωh2 )−1 . p1,...,p4
are particle names. To sum over several channels one can write "*" instead of a particle
name, e.g "*" in place of p1.
• omegaCh is an array that contains the relative contribution and particle names for each
annihilation channel. In the Fortran version one uses instead the function
omegaCh(i,weight,pdg,process). These array and function are similar to vSigmaTCh
described above. The array omegaCh if filled after calling either darkOmegaFO or printChannels.
There is an option to calculate the relic density in models with DM -DM asymmetry.
In this case we assume that the number difference DM -DM is conserved in all reactions.
Thus a small difference in initial abundances can lead to a large DM asymmetry after
freeze-out as is the case for the baryon asymmetry.
•deltaY
describes the difference between the DM and anti-DM abundances for the models where
the number of DM particles minus the number of anti-DM is conserved in decays and
collisions. In such models deltaY is a constant during the thermal evolution of the
Universe, see Ref. [7].
•dmAsymm
is defined by the equation
1 ± dmAsymm
Ω± = Ω
2
and evaluated by micrOMEGAs while calculating the relic density with an initial asymmetry
deltaY, see [7]. This parameter can also be reset after the relic density computation and
will then be taken into account for direct and indirect detection rates.
•darkOmegaExt(&Xf, vs_a, vs_sa)
calculates the dark matter relic density Ωh2 for annihilation cross sections provided by
external functions. Here vs_a is the annihilation cross section in [c·pb] as a function of
the temperature in [GeV] units while vs_sa is the semi-annihilation cross section. vs_a
is required for all models, while vs_sa is relevant only for models where semi-annihilation
occurs. The user can substitute NULL for vs_sa when semi-annihilation is not possible.
darkOmegaExt can also be used if 2 → 2 processes do not contribute to DM annihila15
tion. In this case the appropriate annihilation or semi-annihilation cross sections can be
calculated by vSigmaCC and the tabulated results stored in vs_a and vs_sa. Note that
if the user substitute some function which is not in tabular form, darkOmegaExt can be
slow as it has not been optimized.
darkOmegaExt solves the Runge-Kutta equation in the interval [Tstart, Tend] where
Tstart is defined automatically while Tend has a fixed value 10−3 GeV. darkOmegaExt is
sensitive to effect of DM asymmetry.
6.3
Calculation of relic density for two-component Dark Matter
models.
•darkOmega2(fast, Beps)
Calculates Ωh2 for either one- or two-components DM models. In the former case it
should give the same result as darkOmega. The parameters fast and Beps have the same
meaning as for the darkOmega routine. The returned value corresponds to the sum of the
contribution of the two DM components to Ωh2 . darkOmega2 also calculates the global
parameter fracCDM2 which represents the mass fraction of CDM2 in the total relic density
Ω = Ω1 + Ω2
Ω2
fracCDM2 =
Ω
(3)
(4)
This parameter is then used in routines which calculate the total signal from both DM
candidates in direct and indirect detection experiments, nucleusRecoil, calcSpectrum,
and neutrinoFlux. The user can change the global fracCDM2 parameter before the
calculation of these observables to take into account the fact that the value of the DM
fraction in the Milky Way could be different than in the early Universe.
The routines that were described in section 6.2 are not available for two-component
DM models. In particular the individual channel contribution to the relic density cannot
be computed and DM asymmetry is ignored. After calling darkOmega2 the user can check
the cross sections for each class of reactions (but not for individual processes) which were
tabulated during the calculation of the relic density. The functions
•vsabcd F(T)
computes the sum of the cross sections for each class of reactions (a, b, c, d = 0, 1, 2)
tabulated during the calculation of the relic density. Here T is the temperature in [GeV]
and the return value is vσ in [c·pb]. These functions are defined in the interval [Tstart
, Tend] where Tstart is a global parameter defined by darkOmega2, Tend=10−3 GeV.
Specifically the functions available are
vs1100F
vs1220F
vs1110F
vs1222F
vs1120F
vs2200F
vs1112F
vs2210F
vs1122F
vs2220F
vs1210F
vs2211F
vs1211F
vs2221F
The temperature dependence of the abundances can also be called by the user, the
functions are named Y1F(T) and Y2F(T) and are defined only in the interval T ∈ [Tend, Tstart].
The equilibrium abundances are accessible via the Yeq1(T), Yeq2(T) functions and the
deviation from equilibrium by the functions dY1F(T)= Y1F(T)-Y1eq(T) and
dY2F(T)=Y2F(T)-Y2eq(T).
16
6.4
Calculation of relic density for freeze-in.
Several routines are provided in micrOMEGAs to compute the DM abundance in freeze-in
scenarios. These can be found in the file sources/freezein.c. The first line of this file
contains the statement
//#define NOSTATISTICS
This statement can be uncommented for micrOMEGAs to compute the relic density assuming a Maxwell-Boltzman distribution. This option is faster.
The auxiliary functions that are needed for the computation of the factors from statistical
quantum mechanics are
• Stat2(P/T, xY , x1 , x2 , η1 , η2 ),
returns the S function defined in Eq. (5), that takes into account particle statistical
distributions for the decay of a mediator Y → a, b of fixed momentum P .
1
S (P/T, xY , xa , xb , ηa , ηb ) =
2
Z1
dcθ
−1
(eEa (cθ )/T
eEY /T
.
− ηa )(eEb (cθ )/T − ηb )
(5)
where xi ≡ mi /T , Ea , Eb are the energies of the outgoing particles and ηi ≡ ±eµi /T and
µi is the chemical potential.
• K1to2(x1 , x2 , x3 , η1 , η2 , η3 ),
returns the K̃1 function defined in Eq. (6), that takes into account particle statistical
distributions.
Z Y
3 3
1
1
d pi
K̃1 (x1 , x2 , x3 , η1 , η2 , η3 ) ≡
eE1 /T δ 4 (P1 − P2 − P3 )
Ei /T − η
(4π)2 pCM T
E
e
i
i
i=1
(6)
The code does not check whether or not a particle is in thermal equilibrium with the SM
thermal bath and that it is the responsability of the user to specify which particles belong
to the bath, B, or are out of equilibrium, F. This can be done through the function
• toFeebleList(particle_name)
which assigns the particle particle_name to the list of feebly interacting ones (i.e. those
which belong to F). Feebly interacting particles can be odd or even. This function can
be called several times to include more than one particle. All odd or even particles that
are not in this list are assumed to be in thermal equilibrium with the SM and belong to
B. The treatment of the particles that belong to F for the computation of Ωh2 within
the freeze-in routines is described below. Calling toFeebleList(NULL) will reassign all
particles to B.
The actual computation of the freeze-in DM abundance can be performed with the
help of three functions:
• darkOmegaFiDecay(TR, Name, KE, plot)
calculates the DM abundance from the decay of the particle Name into all odd FP. Here
we assume that all odd FP will decay into the lightest one which is the DM. TR is the
reheating temperature and KE is a switch to specify whether the decaying particle is in
17
kinetic equilibrium (KE=1) or not (KE=0) with the SM. The equations used for the three
different cases are described in [8]. Numerically, the latter two methods give very similar
results, however the function with KE=1 is faster. The switch plot=1 displays on the
screen Y (T ) for the decaying particle and dY /d log(T ) for DM.
• darkOmegaFi22(TR, Process, vegas, plot, &err)
calculates the DM abundance taking into account only DM production via the 2 → 2
process defined by the parameter Process. For example "b,B ->~x1,~x1" for the production of DM (here ~x1) via bb̄ scattering. This routine allows the user to extract the
contribution of individual processes. TR is the reheating temperature. When the switch
vegas=1, the collision term is integrated directly as described in [8]. The execution
time for this option is quite long, it is intended mostly for precision checks. The switch
plot=1 displays on the screen dY /d log(T ) for DM. Note that the temperature profile
for DM production obtained by darkOmegaFiDecay and darkOmegaFi22 can be different.
For example, in the case of a s-channel resonance, the temperature for DM production
corresponds to the one of the mediator decay for darkOmegaFiDecay and the temperature
at which the mediator is created for darkOmegaFi22. err is the returned error code, it
has the following meaning
1: the requested processes does not exist
2: 2 → 2 process is expected
3: can not calculate local parameters // some constrain parameters can not be calculated.
4: the reheating temperature is too small, (TR < 1keV )
5: one of the incoming particles belong to F.
6: None of the outgoing particles are odd and feeble.
7: There is an on-shell particle in t/u channel
8: A numerical integral cannot be performed precisely.
When substituting NULL for the error code, the error message is displayed on the screen
and in this case the message does not appear in the address of the variable used for passing
the error code.
• darkOmegaFi(TR, &err)
calculates the DM abundance after summing over all 2 → 2 processes involving particles
in the bath B in the initial state and at least one particle in F in the final state. The routine checks the decay modes of all bath particles and if one of them has no decay modes
into two other bath particles, the 2 → 2 processes involving this particle are removed
from the summation and instead the contribution to the DM abundance computed from
the routine darkOmegaFiDecay is included in the sum. This is done to avoid appearance
of poles in the corresponding 2 → 2 cross-section. We recommend for such models to
compute individual 2 → 2 processes with darkOmegaFi22 described above. As before, we
assume that all odd FIMPs will decay into the lightest one which is the DM. TR has the
same meaning as above. err is the returned error code, err=1 if feeble particles have not
been defined.
• printChannelsFi(cut,prcnt,filename)
writes into the file filename the contribution of different channels to Ωh2 . The cut param18
eter specifies the lowest relative contribution to be printed. If prcnt 6= 0, the contributions
are given in percent. The routine darkOmegaFi fills the array omegaFiCh which contains
the contribution of different channels ( 2 → 2 or 1 → 2) to Ωh2 . omegaFiCh[i].weight
specifies the relative weight of the ith channel, omegaFiCh[i].prtcl[j] (j=0, 4) defines the
particles names for the ith channel. The last record in the array omegaFiCh has zero
weight and NULL particle names.
Note that if no particle has been declared as being feebly interacting, the freeze-out routines darkOmega, darkOmegaFO, and darkOmega2 [24] will work exactly like in previous
versions of micrOMEGAs. A non-empty list of FIMPs, however, will affect these routines
since micrOMEGAs will exclude all odd particles in this list from the computation of the
relic density via freeze-out. For example, if the DM is a bath particle, excluding FPs can
impact the freeze-out computation of Ωh2 when they are nearly degenerate in mass with
the DM, hence could potentially contribute to coannihilation processes. Note also that
if, for example, the lightest odd particle (the DM) belongs to F and the user computes
the freeze-out abundance of the lightest odd bath particle, the resulting value of ΩLBP h2
will be automatically rescaled by a factor MLF P /MLBP , where MLBP (MLF P ) is the mass
of the lightest odd particle in B (F). Conversely, if the DM belongs to B and the user
computes the freeze-in abundance of the lightest feeble odd particle, the corresponding
result for ΩLF P h2 will be rescaled by MLBP /MLF P . In other words, the answer obtained
in both of these cases corresponds to the predicted density of the dark matter particles
and not the heavier ones. Besides, micrOMEGAs does not check whether the decay rate of
an odd particle to the feeble particles is much smaller than H(TF O ) which would justify
the fact that they should not be included in the freeze-out computation.
7
7.1
Direct detection.
Amplitudes for elastic scattering
• nucleonAmplitudes(CDM,pAsi,pAsd,nAsi,nAsd)
calculates the amplitudes for CDM-nucleon elastic scattering at zero momentum. pAsi(nAsi)
are spin independent amplitudes for protons(neutrons) whereas pAsd(nAsd) are the corresponding spin dependent amplitudes. Each of these four parameters is an array of
dimension 2. The zeroth (first) element of these arrays gives the χ-nucleon amplitudes
whereas the second element gives χ-nucleon amplitudes. Amplitudes (in GeV−2 ) are
normalized such that the total cross section for either χ or χ cross sections is
σtot =
4Mχ2 MN2
(|ASI |2 + 3|ASD |2 )
π(Mχ + MN )2
(7)
nucleonAmplitudes depends implicitly on form factors which describe the quark contents in the nucleon. These form factors are global parameters (see Table 1 for default
values)
T ypeFFPq T ypeFFNq
where T ype is either ”Scalar”, ”pVector”, or ”Sigma”, FFP and FFN denote proton and
neutron and q specifies the quark, d, u or s. Heavy quark coefficients are calculated
automatically.
19
micrOMEGAs automaticaly takes into account loop contributions from box diagrams
as calculated in [25] (DM spin 1/2 case) and [26] (DM spin 0 and 1 cases).
• calcScalarQuarkFF(mu /md ,ms /md ,σπN ,σs )
computes the scalar coefficients for the quark content in the nucleon from the quark mass
ratios mu /md , ms /md as well as from σπN and σs . The default values given in Table 2
are obtained for σs = 42MeV, σπN = 34MeV, mu /md = 0.56, ms /md = 20.2 [27]. The
function calcScalarQuarkFF(0.553,18.9,55.,243.5) will reproduce the default values of the
scalar quark form factors used in micrOMEGAs2.4 and earlier versions.
7.2
Scattering on nuclei
• nucleusRecoil(f,A,Z,J,Sxx,dNdE)
This is the main routine of the direct detection module. The input parameters are:
⋄ f - the DM velocity distribution normalized such that
Z ∞
vf (v)dv = 1
0
The units are km/s for v and s2 /km2 for f(v).
⋄ A - atomic number of nucleus;
⋄ Z - number of protons in the nucleus, predefined values for a wide set of isotopes
are called with Z {N ame};
⋄ J - nucleus spin, predefined values for a wide set of isotopes are called with
J {N ame}{atomic number}.
⋄ Sxx - is a routine which calculates nucleus form factors for spin-dependent interactions (S00,S01,S11), it depends on the momentum transfer in f m−1 . The available
form factors are
SxxF19
SxxK39
SxxI127
SxxXe131
SxxNa23
SxxGe73
SxxI127A
SxxXe131A
SxxNa23A SxxAl27
SxxSi29
SxxGe73A SxxNb92
SxxTe125
SxxXe129 SxxXe129A
SxxXe131B
SxxSi29A
SxxTe125A
The last character is used to distinguish different implementations of the form factor
for the same isotope, see details in [5].
The form factors for the spin independent (SI) cross section are defined by a Fermi distribution and depend on the global parameters Fermi_a, Fermi_b, Fermi_c.
The returned value gives the number of events per day and per kilogram of detector
material. The result depends implicitly on the global parameter rhoDM, the density of
DM near the Earth. The distribution over recoil energy is stored in the array dNdE which
by default has N step = 200 elements. The value in the ith element corresponds to
dN dE[i] =
dN
|E=i∗keV ∗step
dE
20
in units of (1/keV/kg/day). By default step is set to 1.
dNdERecoil(E,dNdE) interpolates the dNdE table.
For a complex WIMP, nucleusRecoil averages over χ and χ. For example for
a call to this routine will be:
73
Ge,
nucleusRecoil(Maxwell,73,Z_Ge,J_Ge73,SxxGe73,dNdE);
• setRecoilEnergyGrid(step,Nstep)
changes the values of step and Nstep for the computation of dNdE.
• Maxwell(v)
returns
Z
cnorm
(~v − VEarth )2
3
f (v) =
d ~v exp −
δ(v − |~v |)
v
∆v 2
|~v | 10 Emax .
• addSpectrum(Spect,toAdd)
sums the spectra toAdd and Spect and writes the result in Spect. For example, this
routine can be useful for summing spectra with different maximal energy.
• spectrMult(Spec, func)
allows to multiply the spectrum Spec by any energy dependent function func
• spectrInt(Emin,Emax,Spec)
integrates a spectrum/flux, Spec from Emin to Emax.
• spectrInfo(Emin,Spec,&Etot)
provides information on the spectra. The returned value and Etot corresponds respectively to
Ntot =
E
Zmax
SpectdN dE(E, Spec)dE = spectrInt(Emin , Emax , Spec)
Emin
Etot =
E
Zmax
E SpectdN dE(E, Spec)dE
Emin
where the first element of the table Spec contains the value of Emax .
8.2
Annihilation spectra
• calcSpectrum(key,Sg,Se,Sp,Sne,Snm,Snl,&err)
calculates the spectra of DM annihilation at rest and returns σv in cm3 /s . The calculated spectra for γ, e+ , p̄, νe , νµ , ντ are stored in arrays of dimension NZ as described
above: Sg, Se, Sp, Sne, Snm, Snl. To remove the calculation of a given spectra, substitute NULL for the corresponding argument. key is a switch to include the polarisation
of the W, Z bosons (key=1) or photon radiation (key=2). Note that final state photon
radiation (FSR) is always included. When key=2 the 3-body process χχ′ → XX + γ
is computed for those subprocesses which either contain a light particle in the t-channel
(of mass less than 1.2 Mcdm) or an outgoing W when Mcdm>500GeV. The FSR is then
subtracted to avoid double counting. Only the electron/positron spectrum is modified
with this switch. When key=4 the contibutions for each channel to the total annihilation
22
rate are written on the screen. More than one option can be switched on simultaneously
by adding the corresponding values for key. For example both the W polarization and
photon radiation effects are included if key=3. A problem in the spectrum calculation
will produce a non zero error code, err 6= 0. calcSpectrum interpolates and sums spectra
obtained by Pythia. The spectra tables are provided only for Mcdm> 2GeV. The results
for a dark matter mass below 2 GeV will therefore be wrong, for example an antiproton
spectrum with kinematically forbidden energies will be produced. A warning is issued for
Mcdm< 2GeV.
• vSigmaCh
is an array that contains the relative contribution and particle names for each annihilation channel. It is similar to vSigmaTCh described in Section 6.2. Note that the list of
particles contains five elements to allow to include gamma radiation. For 2->2 processes
vSigmaCh[n].prtcl[4]=NULL. The array vSigmaCh is filled by calcSpectrum. In the Fortran version one uses instead the function
vSigmaCh(i,weight,pdg,process)
which is similar to the Fortran vSigmaTCh described in Section 6.2.
8.3
Distribution of Dark Matter in Galaxy.
The indirect DM detection signals depend on the DM density in our Galaxy. The DM
density is given as the product of the local density at the Sun with the halo profile function
ρ(r) = ρ⊙ Fhalo (r)
(8)
In micrOMEGAs ρ⊙ is a global parameter rhoDM and the Zhao profile [28]
Fhalo (r) =
R⊙
r
γ
α
rcα + R⊙
rcα + rα
β−γ
α
(9)
with α = 1, β = 3, γ = 1, rc = 20[kpc] is used by default. R⊙ , the distance from the
Sun to the galactic center, is also a global parameter, Rsun. The parameters of the Zhao
profile can be reset by
• setProfileZhao(α,β,γ,rc)
The function to set another density profile is
• setHaloProfile(Fhalo (r))
where Fhalo (r) is any function which depends on the distance from the galactic center,
r, defined in [kpc] units. For instance, setHaloProfile(hProfileEinasto) sets Einasto
profile
α
2
r
Fhalo (r) = exp −
−1
α
R⊙
where by default α = 0.17, but can be changed by
• setProfileEinasto(α)
The command setHaloProfile(hProfileZhao) sets back the Zhao profile. Note that
both setProfileZhao and setProfileEinasto call setHaloProfile to define the corresponding profile.
Dark matter annihilation in the Galaxy depends on the average of the square of the
DM density, < ρ2 >. This quantity can be significantly larger than < ρ >2 when clumps
23
of DM are present [29]. In micrOMEGAs, we use a simple model where fcl is a constant
that characterizes the fraction of the total density due to clumps and where all clumps
occupy the same volume Vcl and have a constant density ρcl . Assuming clumps do not
overlap, we get
< ρ2 >= ρ2 + fcl ρcl ρ.
(10)
This simple description allows to demonstrate the main effect of clumps: far from the
Galactic center the rate of DM annihilation falls as ρ(r) rather than as ρ(r)2 . The parameters ρcl and fcl have zero default values. The routine to change these values is
• setClumpConst(fcl ,ρcl )
To be more general, one could assume that ρcl and fcl depend on the distance from the
galactic center. The effect of clumping is then described by the equation
f
< ρ2 > (r) = ρ(r)(ρ(r) + ρef
clump (r)),
(11)
and the function
f
• setRhoClumps(ρef
clump )
allows to implement a more sophisticated clump structure. To return to the default
treatment of clumps call setRhoClumps(rhoClumpsConst) or setClumpConst.
8.4
Photon signal
The photon flux does not depend on the diffusion model parameters but on the angle
φ between the line of sight and the center of the galaxy as well as on the annihilation
spectrum into photons
• gammaFluxTab(fi,dfi,sigmav,Sg,Sobs)
multiplies the annihilation photon spectrum with the integral over the line of sight and
over the opening angle to give the photon flux. fi is the angle between the line of sight
and the center of the galaxy, dfi is half the cone angle which characterizes the detector
resolution (the solid angle is 2π(1 − cos(df i)) , sigmav is the annihilation cross section, Sg
is the DM annihilation spectra. Sobs is the spectra observed in 1/(GeV cm^2 s ) units.
The function gammaFluxTab can be used for the neutrino spectra as well.
• gammaFluxTabGC(l,b,dl,db,sigmav,Sg,Sobs)
is similar to gammaFluxTab but uses standard galactic coordinates. Here l is the galactic
longitude (measured along the galactic equator from the galactic center, and b is the latitude (the angle above the galactic plane). Both l and b are given in radians. The relation
between the angle fi used above and the galactic coordinates is fi= cos−1 (cos(l) cos(b)).
gammaFluxTabGC integrates the flux over a rectangle [(l, b) − (l + dl, b + db)].
• gammaFlux(fi,dfi,dSigmavdE)
computes the photon flux for a given energy E and a differential cross section for photon
production, dSigmavdE. For example, one can substitute dSigmavdE=σvSpectdNdE(E,SpA)
where σv and SpA are obtained by calcSpectrum. This function can also be used to compute the flux from a monochromatic gamma-ray line by substituting the cross section at
fixed energy (in cm3 /s) instead of dSigmavdE, for example the cross sections obtained
with the loopGamma function in the MSSM, NMSSM, CPVMSSM models (vcsAA and
vcsAZ). In this case the flux of photons can be calculated with
gammaFlux(fi,dfi,2*vcsAA+vcsAZ).
• gammaFluxGC(l,b,dl,db,vcs)
is the analog of gammaFlux when using standard galactic coordinates.
24
8.5
Propagation of charged particles.
The observed spectrum of charged particles strongly depends on their propagation in the
Galactic Halo. The propagation depends on the global parameters
K_dif, Delta_dif, L_dif, Rsun, Rdisk
as well as
Tau_dif (positrons), Vc_dif (antiprotons)
• posiFluxTab(Emin,sigmav, Se, Sobs)
computes the positron flux at the Earth. Here sigmav and Se are values obtained by
calcSpectrum. Sobs is the positron spectrum after propagation. Emin is the energy cut
to be defined by the user. Note that a low value for Emin increases the computation time.
The format is the same as for the initial spectrum. The function SpectrdNdE(E,Sobs)
described above can also be used for the interpolation, in this case the flux is returned in
(GeV s cm2 sr)−1 .
• pbarFlux(E,dSigmavdE)
computes the antiproton flux for a given energy E and a differential cross section for
antiproton production, dSigmavdE. For example, one can substitute
dSigmavdE=σvSpectdNdE(E,SpP)
where σv and SpP are obtained by calcSpectrum. This function does not depend on
the details of the particle physics model and allows to analyse the dependence on the
parameters of the propagation model.
• pbarFluxTab(Emin,sigmav, Sp, Sobs)
computes the antiproton flux, this function works like posiFluxTab,
• solarModulation(Phi, mass, stellarTab, earthTab)
takes into account modification of the interstellar positron/antiproton flux caused by the
electro-magnetic fields in the solar system. Here Phi is the effective Fisk potential in
MeV, mass is the particle mass, stellarTab describes the interstellar flux, earthTab is
the calculated particle flux in the Earth orbit.
Note that for solarModulation and for all *FluxTab routines one can use the same
array for the spectrum before and after propagation.
9
Neutrino signal from the Sun and the Earth
This module does not work yet in case of 2DM
After being captured, DM particles concentrate in the center of the Sun/Earth and
then annihilate into Standard Model particles. These SM particles further decay producing neutrinos that can be observed at the Earth. The neutrino spectra originating
from different annihilation channels into SM particles and taking into account oscillations
and Sun medium effects were computed both in WimpSim [30] and in PPPC4DMν [31].
We use the set of tables provided by these two groups as well as those from DMν [32]
which were included in previous versions of micrOMEGAs. The new global parameter
25
WIMPSIM allows to choose the neutrino spectra. The default value WIMPSIM=0 2 corresponds to the PPPC4DMν spectra while WIMPSIM=1 corresponds to the WimpSim spectra
and WIMPSIM=-1 to the DMν spectra.
• neutrinoFlux(f,forSun,nu, nu_bar)
calculates the muon neutrino/anti-neutrino fluxes near Rthe surface of the Earth. Here
∞
f is the DM velocity distribution normalized such that 0 vf (v)dv = 1. The units are
km/s for v and s2 /km2 for f(v). For example, one can use the same Maxwell function
introduced for direct detection. This routine implicitly depends on the WIMPSIM switch.
If forSun==0 then the flux of neutrinos from the Earth is calculated, otherwise this
function computes the flux of neutrinos from the Sun. The calculated fluxes are stored
in nu and nu bar arrays of dimension NZ=250. The neutrino fluxes are expressed in
[1/Year/km2 ].
• muonUpward(nu,Nu,muon)
calculates the muon flux which results from interactions of neutrinos with rocks below the
detector. Here nu and Nu are input arrays containing the neutrino/anti-neutrino fluxes
calculated by neutrinoFlux. muon is an array which stores the resulting sum of µ+ , µ−
fluxes. SpectdNdE(E,muon) gives the differential muon flux in [1/Year/km2 /GeV] units.
The muon flux weakly depends on the propagation medium, e.g. rock or ice. The energy
lost during propagation is described by the equation [33]
dE
= −(α + βE)ρ
dx
(12)
For propagation in ice (the switch forRocks=0), micrOMEGAs substitutes ρ = 1.0 g/cm3 ,
α = 0.00262 GeVcm2 /g, β = 3.5 × 10−6 cm2 /g [34], while for propagation in rocks, ρ =
2.6 g/cm3 , α = 0.002 GeVcm2 /g, β = 3.0 × 10−6 cm2 /g [33]. The result depends on the
ratio α/β .
• muonContained(nu,Nu,rho, muon) calculates the flux of muons produced in a given
detector volume.This function has the same parameters as muonUpward except that the
outgoing array gives the differential muon flux resulting from neutrinos converted to muons
in a km3 volume given in [1/Year/km3 /GeV] units. rho is the density of the detector in
g/cm3 .
• atmNuFlux(nu,cs,E)
returns the atmospheric muon neutrinos (nu > 0) and anti-neutrinos spectrum (nu < 0)
in [1/Year/km^2] units for a given cosine of the zenith angle, cs. This function is based
on [35].
Two functions allow to estimate the background from atmospheric neutrinos creating
muons after interaction with rocks below the detector or with water inside the detector.
• ATMmuonUpward(cosFi,E) calculates the sum of muon and anti-muon fluxes resulting
from the interaction of atmospheric neutrinos with rocks in units of [1/Year/km2 /GeV/Sr].
cosFi is the energy between the direction of observation and the direction to the center
of Earth. E is the muon energy in GeV. The result depends on the forRock switch.
• ATMmuonContained(cosFi, E, rho) calculates the muon flux caused by atmospheric
neutrinos produced in a given (detector) volume. The returned value for the flux is given
2
Since PPPC4DMν does not provide neutrino specrtra produced at the center of the Earth, in this
case and for WIMPSIM=0 micrOMEGAs uses the DMν spectra.
26
in 1/Year/km3 /GeV/Sr. rho is the density of the detector in g/cm3 units. cosFi and E
are the same as above.
9.1
Comparison with IceCube results
These functions are described in [36] and allow to compare the predictions for the neutrino
flux from DM captured in the Sun with results of IceCube22.
• IC22nuAr(E)
effective area in [km2 ] as a function of the neutrino energy, Aνµ (E)
• IC22nuBarAr(E)
effective area in [km2 ] as a function of the anti-neutrino energy, Aν̄µ (E)).
• IC22BGdCos(cs)
dN
angular distribution of the number of background events as a function of cos φ, d cosbgϕ .
• IC22sigma(E)
neutrino angular resolution in radians as a function of energy.
• exLevIC22( nu_flux, nuB_flux,&B)
calculates the exclusion confidence level for number of signal events generated by given νµ
and ν̄µ fluxes, [36]. The fluxes are assumed to be in [GeV km2 Year]−1 . This function uses
the IC22BGdCos(cs) and IC22sigma(E) angular distribution for background and signal
as well as the event files distributed by IceCube22 with φ < φcut = 8◦ . The returned
parameter B is a Bayesian factor representing the ratio of likelihood functions for the
model with given fluxes and the model with null signal. See details in [36].
• fluxFactorIC22(exLev, nu,nuBar)
For given neutrino, nu, and anti-neutrino fluxes, nuBar, this function returns the factor
that should be applied to the fluxes (neutralino-proton cross sections) to obtain a given
exclusion level exLev in exLevIC22. This is used to obtain limits on the SD cross section
for a given annihilation channel.
10
Cross sections and decays.
The calculation of particle widths, decay channels and branching fractions can be done
by the functions
• pWidth(particleName,&address)
returns directly the particle width. If the 1->2 decay channels are kinematically accessible then only these channels are included in the width when VWdecay,VZdecay= 0. If
not, pWidth compiles all open 1->3 channels and use these for computing the width. If
all 1->3 channels are kinematically forbidden, micrOMEGAs compiles 1->4 channels. If
VWdecay(VZdecay)6= 0, then micrOMEGAs also computes the processes with virtual W (Z)
and adds these to the 1->2 decay channels. Note that 1->3 decay channels with a virtual
W will be computed even if the mass of the decaying particle exceeds the threshold for
1->2 decays by several GeV’s. This is done to ensure a proper matching of 1->2 and
1->3 processes. For particles other than gauge bosons, an improved routine with 3-body
processes and a matching between the 1->2 and 1->3 calculations is kept for the future.
The returned parameter address gives an address where information about the decay
channels is stored. In C, the address should be of type txtList. For models which read a
27
SLHA parameter file, the values of the widths and branchings are taken from the SLHA
file unless the user chooses not to read this data, see (Section 14.5) for details.
• printTxtList(address,FD)
lists the decays and their branching fractions and writes them in a file. address is the
address returned by pWidth.
• findBr(address,pattern)
finds the branching fraction for a specific decay channel specified in pattern, a string
containing the particle names in the CalcHEP notation. The names are separated by
commas or spaces and can be specified in any order.
• slhaDecayPrint(pname,delVirt,FD)
uses pWidth described above to calculate the width and branching ratios of particle pname
and writes down the result in SLHA format. The return value is the PDG particles code.
In case of problem, for instance wrong particle names, this function returns zero. This
function first computes 1 → 2 decays. If such decays are kinematically forbidden then
1 → 3 decay channels are computed. Decays via virtual W/Z bosons will be listed via
their decay products when delV irt 6= 0.
• newProcess(procName)
compiles the codes for any 2 → 2 or 1 → 2 reaction. The result of the compilation
is stored in the shared library in the directory work/so-generated/. The name of the
library is generated automatically.
The newProcess routine returns the address of the compiled code for further usage.
If the process can not be compiled, then a NULL address is returned. 3
Note that it is also possible to compute processes with polarized massless beams, for
example for a polarized electron beam use e% to designate the initial electron.
• procInfo1(address,&ntot,&nin,&nout)
provides information about the total number of subprocesses (ntot) stored in the library
specified by address as well as the number of incoming (nin) and outgoing (nout) particles for these subprocesses. Typically, for collisions (decays), nin=2(1) and nout=2,3.
NULL can be substitute if this information is not needed.
• procInfo2(address,nsub,N,M)
fills an array of particle names N and an array of particle masses M for the subprocess nsub
(1 ≤ nsub ≤ ntot) . These arrays have size nin + nout and the elements are listed in the
same order as in CalcHEP starting with the initial state, see the example in MSSM/main.c.
• cs22(address,nsub,P,c1,c2,&err)
calculates the cross-section for a given 2 → 2 process, nsub, with center of mass momentum P (GeV). All model parameters except the strong coupling GG can be specified with
the functions findVal[W]/assignVal[W] described in Section 5. The strong coupling GG
is defined via the scale parameter GGscale. The differential cross-section is integrated
within the range c1 < cos θ < c2. θ is the angle between p~1 and p~3 in the center-of-mass
frame. Here p~1 (~p3 ) denote respectively the momentum of the first initial(final) particle.
err contains a non zero error code if nsub exceeds the maximum value for the number of
subprocesses (given by the argument ntot in the routine procInfo1). To set the polarization of the initial massless beam, define Helicity[i] where i = 0, 1 for the 1th and
2nd particles respectively. The helicity is defined as the projection of the particle spin on
3
The Fortran version of newProcess returns integer*8.
28
the direction of motion. It ranges from [-1,1] for spin 1 particles and from [-0.5,0.5] for
spin 1/2 particles. By definition a left handed particle has a positive helicity.
• hCollider(Pcm,pp,nf, Qren,Qfac, pName1,pName2,Tmin,wrt) calculates the cross
section for particle production at hadron colliders. Here Pcm is the beam energy in the
center-of-mass frame. pp is 1(−1) for pp(pp̄) collisions, nf ≤ 5 defines the number of quark
flavors taken into account. The parameters Qren and Qfac define the renormalisation and
factorization scales respectively. pName1 and pName2 are the names of outgoing particles.
If Tmin ≤ 0 then hCollider calculates the total cross section for the 2-body final state
process. Otherwise it calculates the cross section for
proton, [a]proton -> pName1, pName2, jet
where Tmin > 0 defines the cut on the jet transversee momentum. The jet content is
defined by the parameter nf. If Qfact ≤ 0, then running ŝ is used for the factorization
scale. If Qren ≤ 0, ŝ is used for the renormalization scale for a 2 → 2 process and pT of
the jet is used for the renormalization scale for a process with a jet in the final state. The
last argument in the hCollider routine allows to switch on/off (wrt=1/0) the printing of
the contribution of individual channels to the total cross section. The value returned is
the total cross section in [pb].
One of the arguments pName1,pName2 can be NULL. Then the cross section for 2 → 1
or 2 → 1 + jet process will be calculated. In Fortran, one should pass a blank string
instead of NULL.
By default hCollider uses the cteq6l structure function built-in the micrOMEGAs code.
One can set any other parton distribution included in either micrOMEGAs or LHAPDF.
The list of structure functions in micrOMEGAs can be obtained with the command
• PDFList
and one of these can be activated by
• setPDF(name)
To work with other PDF’s available in LHAPDF one should first define the environment
variable LHAPDFPATH which specifies the path to the LHAPDF library. Then micrOMEGAs
Makefile links it automatically. The list of available LHAPDF distributions can be obtained with the command
• LHAPDFList
and one of these can be activated by
• setLHAPDF(nset,name)
where nset specifies the subset number. Note, that if a wrong input is provided, setLHAPDF
terminates the execution.
The parton densities are defined by the function
• parton_distr(pdg,x,q)
where pdg is the PDG code of a particle. Note that parton_distr defines the parton
number density and contains a factor 1/x with respect to the definition used in LHAPDF.
11
Tools for model independent analysis
A model independent calculation of the DM observables is also available. After specifying
the DM mass, the cross sections for DM spin dependent and spin independent scattering on proton and neutron, the DM annihilation cross section times velocity at rest and
29
the relative contribution of each annihilation channel to the total DM annihilation cross
section, one can compute the direct detection rate on various nuclei, the fluxes for photons, neutrinos and antimatter resulting from DM annihilation in the galaxy and the
neutrino/muon fluxes in neutrino telescopes.
All the routines presented here depend implicitly on the global parameter Mcdm except
for basicSpectra and basicNuSpectra . These routines do not take into account the
multi-component structure of DM and, in particular, possible differences between DM
and anti-DM. To use these for multi-component DM the user has to perform a summation
over the different DM components.
• nucleusRecoilAux(f,A,Z,J,Sxx,csIp,csIn,csDp,csDn,dNdE)
This function is similar to nucleusRecoil. The additional input parameters include
csIp(csIn) the SI cross sections for WIMP scattering on protons(neutrons) and csDp(csDn)
the SD cross sections on protons(neutrons). A negative value for one of these cross sections
is interpreted as a destructive interference between the proton and neutron amplitudes.
Note that the rate of recoil depends implicitly on the WIMP mass, the global parameter
Mcdm. The numerical value for the global parameter has to be set before calling this function.
• nucleusRecoil0Aux(f,A,Z,J,Sp,Sn,csIp,csIn,csDp,csDn,dNdE) is the corresponding modification of nucleusRecoil0.
For indirect detection, we also provide a tool for model independent studies
• basicSpectra(Mass,pdgN,outN,Spectr)
computes the spectra of outgoing particles and writes the result in an array of dimension
250, Spectr, pdgN is the PDG code of the particles produced in the annihilation of a pair
of WIMPs. To get the spectra generated by transverse and longitudinal W’s substitute
pdgN= 24 +′ T ′ and 24 +′ L′ correspondingly. In the same manner pdgN= 23 +′ T ′ and
23+′ L′ provides the spectra produced by a polarized Z boson. outN specifies the outgoing
particle,
outN = {0, 1, 2, 3, 4, 5} for {γ, e+ , p− , νe , νµ , ντ }
The Mass parameter defines the mass of the DM particle. Note that the propagation
routines for e+ , p− , γ can be used after this routine as usual. Note that the result
of basicSpectra are not valid for Mcdm < 2GeV as explained in the description of
calcSpectrum.
To get indirect detection signals one can substitute the obtained spectra in the
[photon/posi/pbar]FluxTab routines. As long as one keeps the default setting CDM1=CDM2=NULL
these routines will use the Mcdm parameter to calculate the number density of DM particles.
• captureAux(f,forSun,Mass,csIp,csIn,csDp,csDn)
calculates the number of DM particles captured per second assuming the cross sections
for spin-independent and spin-dependent interactions with protons and neutrons csIp,
csIn, csDp, csDn are given as input parameters (in [pb]). A negative value for one of
the cross sections is interpreted as a destructive interference between the proton and neutron amplitudes. The first two parameters have the same meaning as in the neutrinoFlux
routine Section 9. The result depends implicitly on the global parameters rhoDM and Mcdm
in Table 1.
• basicNuSpectra(forSun,Mass,pdg, pol, nu_tab, nuB_tab)
calculates the νµ and ν̄µ spectra corresponding to the pair annihilation of DM in the
center of the Sun/Earth into a particle-antiparticle pair with PDG code pdg. Mass is the
30
DM mass. Note that this routine depends implicitly on the global parameter WIMPSIM
(1,0,-1) which selects the neutrino spectra computed by WimpSim [30], PPPC4DMν [31]
and DMν [32]. The parameter pol selects the spectra for polarized particles available in
PPPC4DMν. pol=-1(1) corresponds to longitudinal (transverse) polarisation of vector bosons or to left-handed (right-handed) polarisation of fermions, pol=0 is used for
unpolarized spectra. When polarized spectra are not available, the unpolarized ones are
generated irrespective of the value of pol. The parameter outN is 1 for muon neutrino and
-1 for anti-neutrino. The resulting spectrum is stored in the arrays nu_tab and nuB_tab
with NZ=250 elements.
The files main.c/F in the directory mdlIndep contain an example of the calculation of
the direct detection, indirect detection and neutrino telescope signals using the routines
described in this section. The numerical input data in this sample file corresponds to
’MSSM/mssmh.dat’.
12
12.1
Constraints from colliders
The Higgs sector
To obtain the limits on the Higgs sector for models with one or several Higgs bosons,
the predictions for the signal strengths of the 125 GeV Higgs can be compared to the
latest results of the LHC, for this an interface to the public code HiggsSignals [37] or
Lilith [38] is provided. Moreover the exclusion limits obtained from Higgs searches in
different channels at LEP, Tevatron and the LHC can be applied to other Higgses in the
model using HiggsBounds [39].
12.1.1
Lilith
Lilith [38] is a Python library that is used to construct a global likelihood function
L from the latest ATLAS and CMS results on the 125 GeV Higgs. 4 The Lilith inputs are the set of reduced couplings of the 125 GeV state, i.e., couplings normalized
by the SM ones, as well as the branching ratios of Higgs decays
P to invisible, BRinv , or
to undetected non-SM final states, BRundetected = 1 − BRinv − BR(H → SM SM). By
default, the automatic generation of the input file assumes that only DM contributes to
the invisible width. Note that the reduced couplings of the 125 GeV Higgs are defined
for all the models provided with micrOMEGAs with the exception of the Zγ coupling. The
latter is computed within Lilith assuming that only SM particles run in the loop. The
file include/Lilith.inc (or Lilith.inc_f) contains the instructions to launch Lilith
using a system call. The input file Lilith_in.xml for Lilith can be created by two
commands
LilithMDL("Lilith_in.xml")
LilithMO("Lilith_in.xml")
4
Lilith can test Higgs bosons with masses within the [123, 128] GeV interval, a warning will be issued
if no such state can be found. In the case where two or more states have masses within this interval,
their signal strengths will be summed incoherently and an effective Higgs state will be tested against the
LHC measurements.
31
both also return the number of neutral Higgs particles. LilithMDL requires that the
reduced couplings be defined in lib/lilith.c of the model. These files are defined for all
models provided with micrOMEGAs and should be written by the user for new models. On
the other hand LilithMO generates automatically the input file required by Lilith. The
functionality of LilithMO is described in Section 12.1.3. As input parameters, Lilith
also requires the number of free parameters, n_par, and a reference likelihood point,
m2logL_reference. Those are defined in the main.c file of each model and set to 0 by
default.
The SLHA output file, Lilith_out.slha, consists of six entries which are respectively
the log-likelihood evaluated at a parameter space point P, −2 log L(P), the number of
experimental degrees of freedom, nexp , the reference likelihood point, the number of degrees of freedom, ndf , the p-value, p and the database version. For a detailed description
of the various output and their interpretation see [38]. For example, one could flag or
exclude points with too low p-value. In this context, a point with a p-value smaller than
0.3173, 0.0455, 0.0027 could be excluded at more than the 1σ, 2σ, 3σ levels, respectively.
12.1.2
HiggsBounds and HiggsSignals
Constraints on the properties of the 125 GeV Higgs boson can also be obtained with
HiggsSignals. Moreover exclusion limits provided by the experimental LHC and Tevatron collaborations on additional Higgs bosons are obtained through an interface to
HiggsBounds [40]. Theses codes are no longer distributed with micrOMEGAs but are downloaded when required. The file include/hBandS.inc contains the instructions to call both
HiggsBounds and HiggsSignals, in particular the options for running HiggsSignals are
set in this file by fixing the Dataset, Method and PDF parameters. Moreover the interface
uses the effective coupling option for specifying the input see [39] for more details. Two
functions can be used to generate the input SLHA file, HB.in
hbBlocksMDL("HB.in",&NchHiggs)
hbBlocksMO("HB.in",&NchHiggs)
Both return the number of neutral Higgses and NchHiggs gives the number of charged
Higgs particles. The function hbBlocksMDL is located in lib/hbBlocks.c for each model
and contains the appropriate definition of the reduced couplings, it needs to be rewritten
for new models implemented by the user. The function hbBlocksMDL generates automatically the required input file for any model and is thus very convenient to use for
new models. The content of this function is described in Section 12.1.3. Note that the
theoretical uncertainty on the mass of the Higgs boson should be specified in the SLHA
BLOCK DMASS, see the example given in main.[c/F] to set the uncertainty at 2 GeV.
The complete outputs of HiggsBounds and HiggsSignals are stored in the files HB.out
and HS.out respectively and can be accessed and read by the user using the slhaval function [41]. The screen output of micrOMEGAs contains the following information
HB(version number): result
obsratio
channel
where result = 0, 1, −1 denotes respectively whether a parameter point is excluded at
95% CL, not excluded, or invalid; obsratio gives the ratio of the theoretical expectation
32
relative to the observed value for the most constraining channel specified in channel.
The HiggsSignals output displayed on the screen is simply
HS(version number): Nobservables chi^2
pval
where Nobservables gives the number of observables used in the fit, chi^2 the associated
χ2 and pval the p-value. The interpretation of these values is left to the user, see [37]
for a detailed description.
12.1.3
Automatic generation of interface files
The functions LiLithMO and hbBlocksMO contain two routines that allow to extract the
couplings contained in the model file, lgrng1.mdl. The first routine returns a description
of a given vertex. The format used is
lVert* vv=getLagrVertex(name1,name2,name3,name4);
where name1,..,name4 are the names of the particles included in the vertex; for vertices
with three particles, name4 should be replaced by NULL. The return parameter vv is the
memory address of a structure which contains information about the vertex:
• vv->GGpower
- power of strong coupling included in vertex
• vv->nTerms
- number of different Lorentz structures in vertex
• vv->SymbVert[i] - text form of Lorentz structures i ∈ [0, nT erms]
The second routine allows to obtain the numerical coefficients corresponding to each
Lorentz structure. The command is
getNumCoeff(vv,coeff)
with coeff[i] the numerical coefficient for SymbVert[i]. Note that the strong coupling
is factored out of the coefficients.
All the QCD-neutral scalars belonging to the even sector (not designated with ~) are
considered as Higgs particles. For each of these, micrOMEGAs calculates the couplings to
SM fermions and massive bosons and writes down into the interface file the ratio of these
couplings to the corresponding SM Higgs coupling. Note that the couplings of the Higgs
to SM fermions can significantly depend on the QCD scale; micrOMEGAs assumes that the
quark masses entering the vertices are obtained at the same scale in both the SM and the
new model, thus the scale dependence in the reduced couplings to fermions disappears.
The loop-induced couplings of the Higgs to gluons and photons are calculated by the
LiLithMO/hbBlocksMO routines whether or not the Hgg and Hγγ vertices are already
implemented in the Lagrangian. This includes NLO-QCD corrections and is performed
as described in [7, 42].
The various branching ratios and widths of the Higgs required by HiggsBounds,HiggsSignals
or Lilith are also written automatically in the interface file. When these values are provided in an SLHA file, they will be used. Otherwise LiLithMO/hbBlocksMO check the existence of H → gg and H → γγ in the table of decays generated from the model Lagrangian.
33
If found, the branching ratios and total widths are written in the interface file without
comparing with the internal calculations. If not found, then LiLithMO/hbBlocksMO add
these channels and recompute the total widths and all branching ratios.
12.2
Searches for New particles
12.2.1
SModelS
LHC limits on new (odd) particles can be obtained using SModelS [43, 44], a code which
tests Beyond the Standard Model (BSM) predictions against Simplified Model Spectra (SMS) results from searches for R-parity conserving SUSY by ATLAS and CMS.
SModelS v1.1.1 decomposes any BSM model featuring a Z2 symmetry into its SMS components using a generic procedure where each SMS is defined by the vertex structure and
the SM final state particles; BSM particles are described only by their masses, production cross sections and branching ratios. The underlying assumption is that differences
in the event kinematics (e.g. from different production mechanisms or from the spin of
the BSM particle) do not significantly affect the signal selection efficiencies. Within this
assumption, SModelS can be used for any BSM model with a Z2 symmetry as long as all
heavier odd particles decay promptly to the dark matter candidate. Note that due to the
Z2 symmetry only pair production is considered, and missing transverse energy (MET)
is always implied in the final state description. This code is no longer distributed with
micrOMEGAs but is downloaded when required.
SModelS needs three input files:
• an SLHA-type input file, containing the mass spectrum, decay tables5 and production cross sections for the parameter point under investigation;
• particles.py defining the particle content of the model, specifically which particles
are even (“R-even”) and which ones are odd (“R-odd”) under the Z2 symmetry;
• a file for setting the run parameters, parameters.ini.
The first two are located in the same directory as main.c and are automatically written
by micrOMEGAs by calling the function
smodels(Pcm, nf, csMinFb, fileName, wrt)
where Pcm is the proton beam energy in GeV and nf is the number of parton flavors used
to compute the production cross sections of the Z2 -odd particles. (Note that u, d, ū,
d¯ and gluons are always included while s, c, and b quarks are included for nf = 3, 4, 5
respectively.) csMinFb defines the minimum production cross section in pb for Z2 -odd
particles; processes with lower cross sections are not added to the SLHA file passed to
SModelS, here denoted by fileName. Finally, wrt is a steering flag for the screen output;
if wrt 6= 0 the computed cross sections will be also written on the screen. The call of
smodels() is by default done in the file micromegas_5.0.X/include/SMODELS.inc. For
this functionality, the user first needs to specify which series of LHC results should be
considered. This is done by setting in main.c the switch
LHCrun = LHC8 and/or LHC13
5
Note that all decay products in the decay table need to be on-shell.
34
for Run 1 (8 TeV ), Run 2 (13 TeV) or both. This determines whether 8 TeV or 13 TeV
cross sections or both will be appended to the SLHA file, which in turn determines
which simplified model results will be considered from the SModelS database. Parameters
needed for running SModelS can be set in the file parameters.ini, for example the
parameters sigmacut which is the cutoff cross section for topologies to be considered
in the decomposition. For more information on the different options and parameters
available, see the SModelSv1.1 manual [45].
An SLHA-type output is written to smodels.in.smodelsslha, (or an alternative
name selected by the user). This output consists of the blocks,
• SModelS Settings lists the SModelS code and database versions as well as input
parameters for the decomposition.
• SModelS Exclusion contains as the first line the status information if a point is
excluded (1), not excluded (0), or not tested ( −1). The latter can occur in scenarios
with long-lived charged particles or in scenarios where no matching SMS results are
found.
If a point is excluded (status 1), this is followed by a list of all results with R > 1,
sorted by their R values, R is defined as the ratio of the predicted theory cross
section and the corresponding experimental upper limit. For each of these results,
the SMS topology identifier (entry 0) (so-called Tx-name, see [46] for an explanation
of the terminology), the R value (entry 1), a measure of condition violation (entry
2), and the analysis identifier (entry 3) are listed.
If the point is not excluded (status 0), the result with the highest R value is given
instead to show whether a point is close to the exclusion limit or not.
• SModelS Missing Topos lists up to ten missing topologies sorted
by their weights
√
(= σ × BR). Each entry consists of the line number, the s in TeV, the weight
and a description of the topology in the SModelS bracket notation. Note that this
information is useful mainly for points that are not excluded.
Additional blocks SModelS Outside Grid, SModelS Long Cascade,
SModelS Asymmetric Branches provide information about the coverage by Simplified
Models. In order to exploit decay channels involving a SM-like Higgs for which the
experimental collaborations assume SM branching ratios for the h with the mass fixed to
mh = 125 GeV, micrOMEGAs checks whether neutral scalar particles with a mass in the
range 123–128 GeV have branching ratios to W W, ZZ, τ τ, bb̄ within 15% of those of a SM
Higgs of the same mass. The corresponding particle will be identified as a SM Higgs by
an entry of type
25 : "higgs",
-25 : "higgs"
in the rEven dictionary in the file particles.py. Note that the name higgs is reserved
for a SM-like Higgs and should not be assigned generically. If no particle of that name is
identified in particles.py, then the corresponding SMS results requiring a Higgs in the
final state are not used by SModelS to constrain the parameter point.
35
12.2.2
Other limits
Limits from searches for a new massive Abelian gauge boson at the LHC, from LEP on
an invisible Z as well as limits on light neutralinos from LEP are provided through the
functions:
• Zinvisible()
returns 1 and prints a WARNING if the invisible width of the Z boson of the Standard
Model is larger than 0.5 MeV ( [47]) and returns 0 if this constraint is fulfilled. This
routine can be used in any model with one DM where the Z boson is uniquely defined by
its PDG=23 and whether the neutral LSP is its own antiparticle or not.
• Zprimelimits()
returns 0 if the scenario considered is not excluded by Z ′ constraints, 1 if the point is
excluded and 2 if both subroutines dealing with Z ′ constraints cannot test the given
scenario. The routine can be used for any Z ′ uniquely defined by the PDG code 32.
Currently two types of searches defined in different subroutines √
of Zprimelimits() are
′
implemented. The latest Z search in the dilepton final state at s = 13 TeV from ATLAS [48] is considered in the first subroutine Zprime_dilepton. If the scenario considered
is allowed or not tested by Zprime_dilepton, a second subroutine √
called Zprime_dijet
analyses√the point using constraints from LHC dijet searches at s = 8 TeV [49–51]
and at s = 13 TeV [52, 53]. This subroutine uses the recasting performed in [54] for a
combination of ATLAS and CMS searches. Zprimelimits() returns 1 if MZ ′ < 0.5 TeV
and 2 for points for which the narrow-width approximation is not valid, i.e. Γ/MZ ′ > 0.3.
• LspNlsp_LEP()
checks the ocmpatibility with the upper limit [55] on the cross section for the production of neutralinos σ(e+ e− → χ̃01 χ̃0i ), i 6= 1, when the heavier neutralino
P decays into
quark pairs and the LSP, χ̃0i → χ̃01 q q̄. The function returns σ × BR = i σ(e+ e− →
χ̃01 χ̃0i ) × BR(χ̃0i → χ̃01 q q̄) in pb as well as a flag greater than one if σ × BR > 0.1(0.5) pb
if mNLSP > (<)100 GeV [55]. This function can also be applied for non-SUSY models
which feature the same signature, in this case the function will compute the cross section
for production of the LSP and any other neutral particle from the odd sector which can
decay into the LSP and a Z boson.
• monoJet(pName1,pName2)
√
computes the cross section for p, p → pName1,pName2+jet at s = 8 TeV where pName1,
pName2 are the names of neutral outgoing particles and jet includes light quarks (u,d,s)
and gluons. The function returns the resulting CL obtained with the CLs technique [56,57]
for each signal region (SR) of the 8 TeV mono-jet CMS analysis [58] and chooses the most
constraining one.
13
Additional routines for specific models
The models included in micrOMEGAs contain some specific routines which we describe
here for the sake of completeness. The current distribution includes the following models: MSSM, NMSSM, CPVMSSM, UMSSM, IDM (inert doublet model), LHM(little Higgs model),
Z3IDM (Inert doublet model with Z3 symmetry), Z4IDSM (Inert doublet and singlet model
36
with Z5 symmetry) as well as three freeze-in models for which there are currently no
additional routines (SingletDM, ZpPortal, LLL_scalar).
Some of these models contain a special routine for reading the input parameters:
• readVarMSSM, readVarNMSSM, readVarCPVMSSM, readVarlHiggs, readVarRHNM.
These routines are similar to the general readVar routine described in Section 5 but they
write a warning when a parameter is not found in the input file and display the default
values for these parameters.
The supersymmetric models contain several additional routines to calculate the spectrum and compute various constraints on the parameter space of the models. Some
functions are common to the MSSM,NMSSM,CPVMSSM,UMSSM models:
• o1Contents(FD)
prints the neutralino LSP components as the B̃, W̃ , h˜1 , h˜2 fractions. For the NMSSM the
fifth component is the singlino fraction S̃ and for the UMSSM the sixth component is the
bino’ fraction B̃ ′ . The sum of the squares of the LSP components should add up to 1.
13.1
MSSM
The MSSM has a long list of low scale independent model parameters, those are specified
in the SLHA file [21, 59]. They are directly implemented as parameters of the model.
For EWSB scenarios the input parameters are the soft parameters, the names of these
parameters are given in the MSSM/mssm[1/2].par files. The user can assign new values
to these parameters by means of assignVal or readVarMSSM.
• spectEwsbMSSM()
calculates the masses of Higgs and supersymmetric particles in the MSSM including oneloop corrections starting from weak scale input parameters.
In these functions spect stands for one of the spectrum calculators suspect, isajet,
spheno, or softSusy. The default spectrum calculator package is SuSpect. To work with
another package one has to specify the appropriate path in MSSM/lib/Makefile. For this,
the environment variables ISAJET, SPHENO or SOFTSUSY must be redefined accordingly.
Note that we also provide a special interface for ISAJET to read a SLHA file. This means
that the user must first compile the executable isajet_slha which sets up the SLHA
interface in ISAJET. Specific instructions are provided in the README file.
For other MSSM scenarios, the parameters at the electroweak symmetry breaking scale
are derived from an input at high scale. The same codes suspect, isajet, spheno, or
softSusy are used for this.The corresponding routines are:
• spectSUGRA(tb,MG1,MG2,MG3,Al,At,Ab,signMu,MHu,MHd,Ml1,Ml3,Mr1,Mr3,Mq1,Mq3,
Mu1,Mu3,Md1,Md3)
assumes that all input parameters except tb and signMu are defined at the GUT scale.
The SUGRA/CMSSM scenario is a special case of this general routine.
• spectSUGRAnuh(tb,MG1,MG2,MG3,Al,At,Ab,Ml1,Ml3,Mr1,Mr3,Mq1,Mq3,
Mu1,Mu3,Md1,Md3,mu,MA)
realizes a SUGRA scenario with non universal Higgs parameters. Here the Mhu, MHd
parameters in the Higgs potential are replaced with the mu parameter defined at the
EWSB scale and MA, the pole mass of the CP-odd Higgs. The signMu parameter is
omitted because mu is defined explicitly.
•spectAMSB(am0,m32,tb,sng).
37
does the same as above within the AMSB model.
We have an option to directly read a SLHA input file, this uses the function
• lesHinput(file_name)
which returns a non-zero number in case of problem.
The routines for computing constraints are (see details in [3]).
• deltarho()
calculates the ∆ρ parameter in the MSSM. It contains for example the stop/sbottom
contributions, as well as the two-loop QCD corrections due to gluon exchange and the
correction due to gluino exchange in the heavy gluino limit.
• bsgnlo(&SMbsg)
returns the value of the branching ratio for b → sγ, see Appendix A. We have included
some new contributions beyond the leading order that are especially important for high
tan β. SMbsg gives the SM contribution.
• bsmumu()
returns the value of the branching ratio Bs → µ+ µ− in the MSSM. It includes the loop
contributions due to chargino, sneutrino, stop and Higgs exchange. The ∆mb effect relevant for high tan β is taken into account.
• btaunu()
computes the ratio between the MSSM and SM branching fractions for B̄ + → τ + ντ .
• gmuon()
returns the value of the supersymmetric contribution to the anomalous magnetic moment
of the muon.
• Rl23()
computes the ratio of the MSSM to SM value for Rl23 in K + → µν due to a charged higgs
contribution, see Eq.70 in [7].
• dtaunu(&dmunu)
computes the branching ratio for Ds+ → τ + ντ . dmunu gives the branching ratio for Ds+ →
µ+ ν µ
• masslimits()
returns a positive value and prints a WARNING when the choice of parameters conflicts
with a direct accelerator limits on sparticle masses from LEP. The constraint on the light
Higgs mass from the LHC is included.
We have added a routine for an interface with superIso [60]. This code is not included in micrOMEGAs so one has first to define the global environment variable superIso
to specify the path to the package.
• callSuperIsoSLHA()
launches superIso and downloads the SLHA file which contains the results. The return
value is zero when the program was executed successfully. Results for specific observables can be obtained by the command slhaValFormat described in section (14.5). Both
superIso and callSuperIsoSLHA use a file interface to exchange data. The delFiles
flag specifies whether to save or delete the intermediate files.
• loopGamma(&vcs_gz,&vcs_gg)
calculates σv for loop induced processes of neutralino annihilation into γZ and into γγ.
3
The result is given in cms . In case of a problem the function returns a non-zero value.
38
13.2
The NMSSM
As in the MSSM there are specific routines to compute the parameters of the model as
specified in SLHA. The spectrum calculator is NMSPEC [9] in the NMSSMTools_5.0.1 package [61].
• nmssmEWSB()
calculates the masses of Higgs and supersymmetric particles in the NMSSM starting from
the weak scale input parameters [62]. These can be downloaded by the readVarNMSSM
routine.
• nmssmSUGRA(m0,mhf,a0,tb,sgn,Lambda,aLambda,aKappa)
calculates the parameters of the NMSSM starting from the input parameters of the
CNMSSM.
The routines for computing constraints are taken from NMSSMTools (see details
in [4]).
• bsgnlo(&M,&P), bsmumu(&M,&P), btaunu(&M,&P), gmuon(&M,&P)
are the same as in the MSSM case. Here the output parameters M and P give information
on the lower/upper experimental limits [63]
• deltaMd(&M,&P),deltaMs(&M,&P)
0
compute the supersymmetric contribution to the Bd,s
− B 0 d,s mass differences, ∆Md and
∆Ms .
• NMHwarn(FD)
is similar to masslimits except that it also checks the constraints on the Higgs masses,
returns the number of warnings and writes down warnings in the file FD.
• loopGamma(&vcs_gz,&vcs_gg)
calculates σv for loop induced processes of neutralino annihilation into γZ and into γγ.
3
The result is given in cms . In case of a problem the function returns a non-zero value.
13.3
The CPVMSSM
The independent parameters of the model include, in addition to some standard model
parameters, only the weak scale soft SUSY parameters. The independent parameters are
listed in CPVMSSM/work/models/vars1.mdl. Masses, mixing matrices and parameters
of the effective Higgs potential are read directly from CPsuperH [11, 64], together with
the masses and the mixing matrices of the neutralinos, charginos and third generation
sfermions. Masses of the first two generations of sfermions are evaluated (at tree-level)
within micrOMEGAs in terms of the independent parameters of the model.
The routines for computing constraints are taken from CPsuperH, [65]
• bsgnlo(), bsmumu(), btaunu(), gmuon()
are the same as in the MSSM case.
• deltaMd(),deltaMs()
are the same as in the NMSSM case.
• Bdll()
computes the supersymmetric contribution to the branching fractions for Bd → τ + τ − in
the CPVMSSM.
39
• ABsg()
computes the supersymmetric contribution to the asymmetry for B → Xs γ.
• EDMel(),EDMmu(),EDMTl()
return the value of the electric dipole moment of the electron, de , the muon,dµ , and of
Thallium, dT l in units of ecm.
13.4
The UMSSM
The independent parameters of the UMSSM are the standard model parameters and
weak scale soft SUSY parameters listed in UMSSM/work/models/vars1.mdl. All masses,
mixing matrices and parameters of the different sectors of the model are computed by
micrOMEGAs [13, 14].
Some routines for computing constraints were taken from the MSSM and were adapted
to the UMSSM. For example
• masslimits() which is essentially the same as in the MSSM except that the constraint
on the light Higgs mass from the LHC was removed as other routines in the UMSSM
include this constrain, or
• deltarho()
calculates the ∆ρ parameter in the UMSSM where in addition to MSSM contributions
a pure UMSSM tree-level contribution from the extended Abelian gauge boson sector is
included. Two other routines in C are included, Zinvisible() and Zprimelimits() that
were defined above.
The remaining routines for computing B-physics, K-physics and Higgs observables as
well as the anomalous magnetic moment of the muon were taken from NMSSMTools_5.0.1
and adapted to the UMSSM [13, 66]. To call these routines use umssmtools(PDG_LSP)
where PDG_LSP is the PDG code of the LSP. The result is contained in four files (UMSSM_inp.dat,
UMSSM_spectr.dat, UMSSM_decay.dat and SM_decay.dat) and the WARNING messages
from these routines can be displayed with slhaWarnings(stdout).
As in the NMSSM the following routines are available :
• bsg(&M,&P), bsmumu(&M,&P), btaunu(&M,&P), gmuon(&M,&P),deltamd(&M,&P),
deltams(&M,&P).
as well as the routines
• bdg(&M,&P), bdmumu(&M,&P), bxislllow(&M,&P), bxisllhigh(&M,&P),
bxisnunu(&M,&P), bpkpnunu(&M,&P), bksnunu(&M,&P),
rdtaul(&M,&P), rdstaul(&M,&P)
to compute respectively b → dγ, Bd → µ+ µ− , the b → sl+ l− transition in the low
([1, 6] GeV2 ) and high (≥ 14.4 GeV2 ) m2l+ l− ranges, B → Xs ν ν̄, B + → K + ν ν̄, B → K ∗ ν ν̄,
+ →Dτ + ν )
+ →D ∗ τ + ν )
τ
τ
RD ≡ BR(B
and RD∗ ≡ BR(B
.
BR(B + →Dl+ νl )
BR(B + →D ∗ l+ νl )
+
+
The K-physics observables K → π ν ν̄, KL → π 0 ν ν̄, ∆MK and ǫK can be computed
respectively with
• kppipnunu(&M,&P), klpi0nunu(&M,&P), deltamk(&M,&P), epsk(&M,&P).
See the README of the UMSSM for further details.
40
14
Tools for new model implementation.
It is possible to implement a new particle physics model in micrOMEGAs. For this the
model must be specified in the CalcHEP format. micrOMEGAs then relies on CalcHEP to
generate the libraries for all matrix elements entering DM calculations. Below we describe
the main steps and tools for implementing a new model.
14.1
Main steps
• The command ./newProject MODEL
MODEL launched from the root micrOMEGAs directory creates the directory MODEL.
This directory and the subdirectories contain all files needed to run micrOMEGAs
with the exception of the new model files.
• The new model files in the CalcHEP format should then be included in the subdirectory MODEL/work/models. The files needed are vars1.mdl, func1.mdl, prtcls1.mdl,
lgrng1.mdl, extlib1.mdl. For more details on the format and content of model
files see [1].
• For odd particles and for the Higgs sector it is recommended to use the widths
that are (automatically) calculated internally by CalcHEP/micrOMEGAs. For this
one has to add the ’ !’ symbol before the definition of the particle’s width in the file
prtcls1.mdl, for example
Full name
Higgs 1
| P | aP|PDG
|h1 |h1 |25
|2*spin| mass |width |color|
|0
|Mh1
|!wh1 |1
|
• Some models contain external functions, if this is the case they have to be compiled
and stored in the MODEL/lib/aLib.a library. These functions should be written in
C and both functions and their arguments have to be of type double. The library
aLib.a can also contain some functions which are called directly from the main
program. The MODEL/Makefile automatically launches make in the lib directory
and compiles the external functions provided the prototypes of these external functions are specified in MODEL/lib/pmodel.h. The user can of course rewrite his own
—lib/Makefile if need be.
If the new aLib.a library needs some other libraries, their names should be added
to the SSS variable defined in MODEL/Makefile.
The MODEL directory contains both C and FORTRAN samples of main routines. In
these sample main programs it is assumed that input parameters are provided in a separate
file. In this case the program can be launched with the command:
./main data1.par
Note that for the direct detection module all quarks must be massive. However the cross
sections do not depend significantly on the exact numerical values for the masses of light
quarks.
41
14.2
Automatic width calculation
Automatic width calculation can be implemented by inserting the ’ !’ symbol before the
name of the particle width in the CalcHEP particle table (file prtcls1.mdl). In this case
the width parameter should not be defined as a free or constrained parameter. Actually
the pWidth function described in section 10 is used for width calculation in this case. We
recommend to use the automatic width calculation for all particles from the ’odd’ sector
and for Higgs particles. For models which use SLHA parameter transfer (Section 14.5),
the automatic width option will use the widths contained in the SLHA file unless the user
chooses the option to ignore this data in the SLHA file, see section 14.5.
14.3
Using LanHEP for model file generation.
For models with a large number of parameters and various types of fields/particles such
as the MSSM, it is more convenient to use an automatic tool to implement the model.
LanHEP is a tool for Feynman rules generation. A few minor modifications to the default
format of LanHEP have to be taken into account to get the model files into the micrOMEGAs
format.
• The lhep command has to be launched with the -ca flag
lhep
-ca source_file
• The default format for the file prtcls1.mdl which specifies the particle content has
to be modified to include a column containing the PDG code of particles. For this,
first add the following command in the LanHEP source code, before specifying the
particles
prtcformat fullname:
’Full
Name ’, name:’ P ’, aname:’ aP’, pdg:’ number
spin2,mass,width, color, aux, texname: ’> LaTeX(A) <’,
atexname:’> LateX(A+) <’ .
’,
Then for each particle define the PDG code. For instance:
vector ’W+’/’W-’: (’W boson’, pdg 24, mass MW, width wW).
• LanHEP does not generate the file extlib1.mdl. micrOMEGAs works without this
file but it is required for a CalcHEP interactive session. The role of this file is to
provide the linker with the paths to all user’s libraries needed at compilation. For
example for the lib/aLib.a library define
$CALCHEP/../MODEL/lib/aLib.a
For examples see the extlib1.mdl files in the directory of the models provided.
42
14.4
QCD functions
Here we describe some QCD functions which can be useful for the implementation of a
new model.
• initQCD(alfsMZ,McMc,MbMb,Mtp)
This function initializes the parameters needed for the functions listed below. It has to
be called before any of these functions. The input parameters are the QCD coupling at
the Z scale, αs (MZ ), the quark masses, mc (mc ), mb (mb ) and mt (pole).
• alphaQCD(Q)
calculates the running αs at the scale Q in the M S scheme. The calculation is done using
the NNLO formula in [67]. Thresholds for the b-quark and t-quark are included in nf at
the scales mb (mb ) and mt (mt ) respectively.
• MtRun(Q), MbRun(Q), McRun(Q)
calculates top, bottom and charm quarks running masses evaluated at NNLO.
• MtEff(Q), MbEff(Q), McEff(Q),
calculates effective top, bottom and charm quark masses using [67]
2
2
2
Mef
f (Q) = M (Q) 1 + 5.67a + (35.94 − 1.36nf )a
+ (164.14 − nf (25.77 − 0.259nf ))a3
(13)
where a = αs (Q)/π, M (Q) and αs (Q) are the quark masses and running strong coupling
in the M S-scheme. In micrOMEGAs, we use the effective quark masses calculated at the
scale Q = 2Mcdm. In some special cases one needs a precise treatment of the light quarks
masses. The function
• MqRun(M2GeV, Q)
returns the running quark mass defined at a scale of 2 GeV. The corresponding effective
mass needed for the Higgs decay width is given by
• Mqeff(M2GeV, Q)
14.5
SLHA reader
Very often the calculation of the particle spectra for specific models is done by some
external program which writes down the particle masses, mixing angles and other model
parameters in a file with the so-called SLHA format [21, 59]. The micrOMEGAs program
contains routines for reading files in the SLHA format. Such routines can be very useful
for the implementation of new models.
In general a SLHA file contains several pieces of information which are called blocks. A
block is characterized by its name and, sometimes, by its energy scale. Each block contains
the values of several physical parameters characterized by a key. The key consists in a
sequence of integer numbers. For example:
BLOCK MASS
# PDG Code
25
37
# Mass spectrum
mass
1.15137179E+02
1.48428409E+03
particle
# lightest neutral scalar
# charged Higgs
BLOCK NMIX # Neutralino Mixing Matrix
1 1
9.98499129E-01
# Zn11
1 2
-1.54392008E-02
# Zn12
43
BLOCK Au Q= 4.42653237E+02 # The trilinear couplings
1 1
-8.22783075E+02
# A_u(Q) DRbar
2 2
-8.22783075E+02
# A_c(Q) DRbar
• slhaRead(filename,mode)
downloads all or part of the data from the file filename. mode is an integer which
determines which part of the data should be read form the file, mode= 1*m1+2*m2+4*m4
where
m1 = 0/1 m2 = 0/1 m4 = 0/1 -
overwrites all/keeps old data
reads DECAY /does not read
DECAY
reads BLOCK/does not read
BLOCK
For example mode=2 (m1=0,m2=1) is an instruction to overwrite all previous data and read
only the information stored in the BLOCK sections of filename. In the same manner
mode=3 is an instruction to add information from DECAY to the data obtained previously.
slhaRead returns the values:
0
-1
-2
-3
n>0
-
successful reading
can not open the file
error in spectrum calculator
no data
wrong file format at line n
• slhaValExists(BlockName, keylength, key1, key2,...)
checks the existence of specific data in a given block. BlockName can be substituted
with any case spelling. The keylength parameter defines the length of the key set
{key1,key2,...}. For example slhaValExists("Nmix",2,1,2) will return 1 if the neutralino mass mixing element Zn12 is given in the file and 0 otherwise.
• slhaVal(BlockName,Q, keylength, key1, key2,......)
is the main routine which allows to extract the numerical values of parameters. BlockName
and keylength are defined above. The parameter Q defines the scale dependence. This
parameter is relevant only for the blocks that contain scale dependent parameters, it will
be ignored for other blocks, for example those that give the particle pole masses. In
general a SLHA file can contain several blocks with the same name but different scales
(the scale is specified after the name of the block). slhaVal uses the following algorithm
to read the scale dependent parameters. If Q is less(greater) than all the scales used in
the different blocks for a given parameter slhaVal returns the value corresponding to
the minimum(maximum) scale contained in the file. Otherwise slhaVal reads the values
corresponding to the two scales Q1 and Q2 just below and above Q and performs a linear
interpolation with respect to log(Q) to evaluate the returned values.
Recently it was proposed to use an extension of the SLHA interface to transfer Flavour
Physics data [68]. Unfortunately the structure of the new blocks is such that they cannot
be read with the slhaVal routine. We have added two new routines for reading such data
• slhaValFormat(BlockName, Q, format)
where the format string allows to specify data which one would like to extract from the
given block BlockName. For instance, to get the b → sγ branching ratio from the block
44
Block FOBS # Flavour observables
# ParentPDG type value
q
NDA ID1 ID2 ID3 ... comment
5
1
2.95061156e-04 0
2
3 22
# BR(b->s gamma)
521
4
8.35442304e-02 0
2 313 22
# Delta0(B->K* gamma)
531
1
3.24270419e-09 0
2 13 -13
# BR(B_s->mu+ mu-)
...
one has to use the command slhaValFormat("FOBS", 0., "5 1 %E 0 2 3 22"). In
this command the format string is specified in C-style. The same routine can be used to
read HiggsBound SLHA output.
A block can also contain a textual information. For example, in HIGGSBOUNDS a block
contains the following records,
Block HiggsBoundsResults
#CHANNELTYPE 1: channel with the highest statistical sensitivity
1
1
328
# channel id number
1
2
1
# HBresult
1
3 0.72692779334500290
# obsratio
1
4
1
# ncombined
1
5 ||(p p)->h+..., h=1 where h is SM-like (CMS-PAS-HIG-12-008)|| # text description of channel
In particular, the last record contains the name of the channel which gives the strongest
constraint on the Higgs. To extract the name of this channel one can use the new function
• slhaSTRFormat("HiggsBoundsResults","1 5 || %[^|]||",channel);
which will write the channel name in the text parameter channel.
• slhaWarnings(fileName)
writes into the file the warnings or error message stored in the SPINFO block and returns
the number of warnings. If FD=NULL the warnings are not written in a file.
• slhaWrite(fileName)
writes down the information stored by readSLHA into the file. This function can be used
for testing purposes.
SLHA also describes the format of the information about particle decay widths. Even
though micrOMEGAs also performs the width calculation, one might choose to read the
information from the SLHA file.
• slhaDecayExists(pNum)
checks whether information about the decay of particle pNum exists in the SLHA file.
pNum is the particle PDG code. This function returns the number of decay channels. In
particular zero means that the SLHA file contains information only about the total width,
not on branching ratios while -1 means that even the total width is not given.
• slhaWidth(pNum)
returns the value of particle width.
• slhaBranch(pNum,N, nCh)
returns the branching ratio of particle pNum into the N-th decay channel. Here
0
Source Exif Data:
File Type : PDF
File Type Extension : pdf
MIME Type : application/pdf
PDF Version : 1.4
Linearized : No
Page Count : 56
Page Mode : UseOutlines
XMP Toolkit : XMP toolkit 2.9.1-13, framework 1.6
About : uuid:8e2a1d83-9058-11f3-0000-f5848ea360a7
Producer : dvips + GPL Ghostscript 9.07
Keywords : ()
Modify Date : 2018:05:15 15:32:18+03:00
Create Date : 2018:05:15 15:32:18+03:00
Creator Tool : LaTeX with hyperref package
Document ID : uuid:8e2a1d83-9058-11f3-0000-f5848ea360a7
Format : application/pdf
Title : ()
Creator : ()
Description : ()
Subject :
Author :
EXIF Metadata provided by EXIF.tools