DarkSUSY 6 : An Advanced Tool to Compute Dark Matter Properties Numerically

The nature of dark matter remains one of the key science questions. Weakly Interacting Massive Particles (WIMPs) are among the best motivated particle physics candidates, allowing to explain the measured dark matter density by employing standard big-bang thermodynamics. We introduce here a radically new version of the widely used DarkSUSY package, which allows to compute the properties of such dark matter particles numerically. With DarkSUSY 6 one can accurately predict a large variety of astrophysical signals from dark matter, such as direct detection rates in low-background counting experiments and indirect detection signals through antiprotons, antideuterons, gamma rays and positrons from the Galactic halo, or high-energy neutrinos from the center of the Earth or of the Sun. For thermally produced dark matter like WIMPs, high-precision tools are provided for the computation of the relic density in the Universe today, as well as for the size of the smallest dark matter protohalos. Furthermore, the code allows to calculate dark matter self-interaction rates, which may affect the distribution of dark matter at small cosmological scales. Compared to earlier versions, DarkSUSY 6 introduces many significant physics improvements and extensions. The most fundamental new feature of this release, however, is that the code has been completely re-organized and brought into a highly modular and flexible shape. Switching between different pre-implemented dark matter candidates has thus become straight-forward, just as adding new - WIMP or non-WIMP - particle models or replacing any given functionality in a fully user-specified way. In this article, we describe the physics behind the computer package, along with the main structure and philosophy of this major revision of DarkSUSY. A detailed manual is provided together with the public release at www.darksusy.org.

Interacting Massive Particles (WIMPs) are among the best motivated particle physics candidates, allowing to explain the measured dark matter density by employing standard big-bang thermodynamics. Examples include the lightest supersymmetric particle, though many alternative particles have been suggested as a solution to the dark matter puzzle. We introduce here a radically new version of the widely used DarkSUSY package, which allows to compute the properties of such dark matter particles numerically. With DarkSUSY 6 one can accurately predict a large variety of astrophysical signals from dark matter, such as direct detection rates in low-background counting experiments and indirect detection signals through antiprotons, antideuterons, gamma rays and positrons from the Galactic halo, or high-energy neutrinos from the center of the Earth or of the Sun. For thermally produced dark matter like WIMPs, high-precision tools are provided for the computation of the relic density in the Universe today, as well as for the size of the smallest dark matter protohalos. Furthermore, the code allows to calculate dark matter self-interaction rates, which may affect the distribution of dark matter at small cosmological scales. Compared to earlier versions, DarkSUSY 6 introduces many significant physics improvements and extensions. The most fundamental new feature of this release, however, is that the code has been completely re-organized and brought into a highly modular and flexible shape. Switching between different pre-implemented dark matter candidates has thus become straight-forward, just as adding new -WIMP or non-WIMP -particle models or replacing any given functionality in a fully user-specified way. In this article, we describe the physics behind the computer package, along with the main structure and philosophy of this major revision of DarkSUSY. A detailed manual is provided together with the public release at www.darksusy.org.

Introduction
The problem of finding the identity of dark matter (DM), one of the most interesting problems in our picture of the universe, still lacks a solution despite an impressive improvement of detection capabilities. Indications have, however, grown even stronger that an unknown form of matter contributes to the gravitationally attractive matter density of the universe, as shown, for example, by analyses of the recent high-precision data from the PLANCK satellite [1]. Also, the discovery [2] of the "bullet cluster", two colliding galaxy clusters where the mass distribution, as extracted from gravitational lensing, and the baryonic mass, as detected by its X-ray emission, are clearly separated, which makes it very difficult to modify gravity to circumvent the problem. We are thus led to believe, as already Zwicky conjectured in the 1930's [3] (and where Lundmark reached a similar conclusion even a few years earlier [4]), that most of the mass of the observable Universe is invisible, or 'dark'. A broad consensus has been reached that the most likely explanation for this excess mass is a so far unknown elementary particle [5]. Analyses combining high-redshift supernova luminosity distances, microwave background fluctuations and baryon acoustic oscillations in the galaxy distribution [1,6] give tight constraints on the present mass density of matter in the Universe. This is usually expressed in the ratio Ω M = ρ M /ρ crit , normalized to the critical density ρ crit = 3H 2 0 /(8πG N ) = h 2 × 1.9 · 10 −29 g cm −3 . The value obtained for the most recent Planck data ( Table 4 in [1]) for cold DM for the (unknown) particle, here called χ, is Ω χ h 2 = 0.1188 ± 0.0010, which is around 5 times higher than the value obtained for ordinary baryonic matter, Ω B h 2 = 0.02225 ± 0.00016. Here h = 0.6727 ± 0.0066 is the derived [1] present value of the Hubble constant in units of 100 km s −1 Mpc −1 . In addition, the Planck data is consistent with a flat universe (Ω tot = 1) and a value for the dark energy component, e.g. a cosmological constant Λ, of Ω Λ = 0.6844 ± 0.0091. Also from the point of view of structure formation, non-baryonic DM seems to be necessary, and the main part of it should consist of particles that were non-relativistic at the time when structure formed (cold dark matter, CDM). This excludes the light standard model neutrinos, although there is perhaps a window remaining for sterile neutrinos in the keV range which would act like "warm" DM (for a review, see [7]). The Planck collaboration [1], limits the contribution of standard model neutrinos to Ω ν h 2 ∼ < 0.005. A well-motivated particle physics candidate which has the required properties has long been the lightest supersymmetric particle, assumed to be a neutralino [8][9][10]. (For thorough reviews of supersymmetric DM, see [11][12][13].) Although supersymmetry has long been generally accepted as a very promising enlargement of the Standard Model of particle physics (for instance it would solve the so-called hierarchy problem of understanding why the electroweak scale is protected against Planck-scale corrections), the details of a viable model are largely unknown. In fact, simplified templates with only a few free parameters, like the so-called constrained minimal supersymmetric model (CMSSM) or mSUGRA, are already very strongly constrained by early LHC data [14,15]. All these results seem to push the mass scale of the lightest supersymmetric particle, the lightest neutralino, up to the TeV scale. For a detailed and comprehensive global analysis of simple supersymmetric scenarios, we refer to the recent results [16,17] of the GAMBIT [18] collaboration.
Of course, there is no compelling reason why the actual model, if nature is supersymmetric at all, should be of this simplest kind. However, the minimal supersymmetric standard model (MSSM) has served and still serves as a useful template with which to test current ideas about detection, both in particle physics accelerators and in DM experiments, and contains many features which are expected to be universal for any weakly interacting massive particle (WIMP) model. Such WIMP DM candidates appear almost inevitably in many nonsupersymmetric extensions to the standard model (e.g. universal extra dimensions [19] or little Higgs models [20]) and have the additional advantage that the observed DM abundance can be understood from a simple thermal production mechanism in the early universe. Although the tendency at present is that supersymmetric (SUSY) models are getting slightly out of fashion, we keep such models in DarkSUSY, not the least to allow comparison of this well-studied candidate with other DM candidates that the user of DarkSUSY may want to supply as input to the program. It should be noted that there is also a large number of DM models that cannot be described within the WIMP paradigm [21][22][23], which will become more and more important alternatives if the classical WIMP searches continue to report the absences of any (uncontroversial) signals. We stress that many routines in DarkSUSY can be used, without any modification, even for DM candidates that are not WIMPs.
The main new feature of DarkSUSY 6, as compared to earlier DarkSUSY [24] versions, is a radically new modular and flexible structure of the code. This allows to handle a wide range of DM particle models, included in the release or supplied by the user, and quite in general makes it much easier to include any user-designed changes or additions. Besides this, we have in DarkSUSY 6 further developed the analytical and numerical tools for computing the relic density of WIMP-like particles from the early universe, and the ensuing direct and indirect detection rates today. Among the completely newly added features are the possibility to compute the cutoff in the power spectrum of matter density perturbations, related to the kinetic decoupling of WIMP-like particles in the early universe, the possibility to compute DM self-interaction rates as well as a detailed treatment of the impact of radiative corrections on the spectra of DM-induced cosmic-ray spectra relevant for indirect detection. Overall, we provide a program package with a rather high level of sophistication, aiming to benefit the scientific community working with problems related to DM. This paper describes the basic structure and the underlying physical and astrophysical formulas contained in DarkSUSY, as well as examples of its use. Even if it has the same name as previous versions, DarkSUSY 6 constitutes a major revision in that it introduces a conceptual change to the very structure of the package. Let us stress that these changes in the structure of the code are physicsdriven, greatly enlarging the range of potential physics applications. In particular, the often hard-coded focus on supersymmetric neutralino DM, as implemented in previous versions, has now been abandoned. In this sense, one may argue that also the meaning has changed, and DarkSUSY should now be interpreted as 1

Dark SUsy Samt Ytterligare modeller
To demonstrate this new interpretation, we explicitly include in this release modules for a generic WIMP and a generic decaying DM particle, as well as for the scalar singlet model [25] (besides, obviously, supersymmetry). Further particle modules will be supplied with future code releases, but can at any time also easily be added by the user.
This article is organized as follows. We start, in Section 2, by summarizing the main new physics features since the last main release of DarkSUSY (describing version 4.0 of the code [24]). This is closely intertwined with the general structure and guiding principles of the code, which we describe in Section 3. In the coming sections we then focus in more detail on the various physics aspects behind the implementation, with special emphasis on improvements or additions since the previous publication [24]. We first describe the thermal production and decoupling of WIMPs (Section 4), the direct detection of DM particles (Section 5) and models for the DM distribution in the Milky Way or other halos (Section 6). As discussed in Section 7, the structure of such DM halos is expected to be affected by strong DM self-interactions, thus providing potential DM observables. We then move on to discuss the particle yield from DM annihilation or decay (Section 8) that are needed for the indirect detection of DM signals using gamma rays and neutrinos from the Galactic halo (Section 9), charged cosmic rays (Section 10) or neutrinos from the interior of the Sun or Earth (Section 11). We summarize, and provide our conclusions in Section 12. In four Appendices, we describe the specific particle physics models that are currently implemented, namely the generic WIMP case (Appendix A), a generic decaying DM scenario (Appendix B), the MSSM (Appendix C), the scalar singlet model (Appendix D) and a module for self-interacting DM confined to a dark sector (Appendix E). Finally, in Appendix F, we provide a quick-start guide and, for users already familiar with DarkSUSY 5 or earlier, list the main technical differences of this release compared to earlier code versions.
The version of the package described in this paper is DarkSUSY 6.1. To download the latest version of DarkSUSY, and for a more technical manual, please visit the official DarkSUSY website at www.darksusy.org. This webpage also contains an updated list of contributed particle modules and other extensions of general interest that have been externally developed by users of the code.

Physics highlights
Let us start with a quick reading guide, mostly for users already familiar with earlier code versions, where we list some of the most important new or considerably improved physics features since the last DarkSUSY publication [24], which described version 4.0 of the code (see also Appendix F.4 for a list of the most important technical changes).
• As described in more detail in Section 3, DarkSUSY can now directly be used with particle models different from supersymmetry. An example of a fully implemented, UV-complete new model is Scalar Singlet DM (Appendix D). Another example is selfinteracting DM in a secluded dark sector (Appendix E).
• The relic density routines have been made more general (allowing for different particle models, including situations where the dark sector and photon temperatures differ), while retaining their ability to treat co-annihilations, thresholds and resonances with high numerical accuracy; partial parallelisation of the code has led to a speed-up that allows to compute the relic density even for 'critical' models (Section 4.1).
• A newly added feature is the calculation of the kinetic decoupling decoupling of WIMPs, along with the size of the smallest protohalos (Section 4.2).
• The routines for direct detection have been significantly revised, and generalised to include both spin-dependent and spin-independent scattering, in the corresponding limit, as well as descriptions using nonrelativistic effective operators (Section 5).
• Another completely newly added feature is the possibility to calculate DM self-interaction rates (Section 7).
• Extensive runs of Pythia [26] have been tabulated to provide the yields of stable particles relevant for DM indirect detection, improving both the statistics of previous tables and extending them to new channels (e.g. light quarks) and final states (anti-deuterons) (Section 8).
• The routines for the numerical propagation of cosmic rays have been completely rewritten; they are now both much more flexible and stable than in previous versions (Section 10). This allows in particular to efficiently scan over propagation parameters, as well as to treat cuspy DM density profiles.
• New capture rate calculations for WIMP capture in the Sun. We now include much more detailed models of the solar composition from Ref. [27]. We include up to 289 different isotopes for spin-independent capture and 112 elements for spin-dependent capture. We also perform the form factor integration numerically, which allows us to use better form factors than the usual Gaussian ones.
• Astrophysics and particle physics have now been completely disentangled, allowing to link the master library to any particle physics module. This release, in particular, ships with complete modules for a generic WIMP (App. A), a generic decaying DM scenario (App. B), the MSSM (App. C), the scalar singlet model (App. D) and a range of simplified DM models with velocity-dependent self-interactions (App. E).
• For the mssm module, the internal bremsstrahlung of both U (1), SU (2) and SU (3) gauge bosons has been fully implemented, resulting in highly accurate predictions of all cosmic-ray spectra relevant for indirect DM detection (Section C.3.2).
• For the relic density calculations in the mssm module, a new effective framework to take into account the effect of QCD corrections to neutralino-neutralino annihilation has been implemented; an interface to DM@NLO [28] for the full NLO result is in preparation (Section C.3.3).
• The mssm module has an updated SLHA reader, which is more versatile and flexible than in previous versions of the code. For example, DarkSUSY can now handle more supersymmetric models.
• The core library contains a new set of functionalities to compute limits [29,30] on DM models from the formation of ultracompact minihalos [31][32][33]. 2 Apart from these additions included in the released version 6 of the code, there currently exist several publicly available packages to be used with DarkSUSY. The full list will be constantly updated at www.darksusy.org, but let us point out in particular dedicated code to compute the Sommerfeld effect in the MSSM [34] as well as detailed likelihoods for neutrino signals from the IceCube experiment [35,36]. After this overview of the basic structure and guiding principles of the code, we will next turn to the physics that is implemented in DarkSUSY 6. In the remaining Sections of the main text of this article, we will describe in some detail those parts of DarkSUSY that are independent of any specific particle-physics framework beyond the standard model (BSM), as implemented in the main library ds_core, and refer to the Appendix for any functionality contained in the particle modules. Throughout, most of our focus will be on new aspects compared to earlier code versions. Before we start, and as a quick reading guide to users familiar with earlier code versions, let us briefly list some of the most important new or considerably improved physics features since the last DarkSUSY publication [24], which described version 4.0 of the code (see also Appendix F.4 for a list of the most important technical changes).

Guiding principles
DarkSUSY 6 has a new structure compared to earlier versions of the code. The most striking difference is that we have split the particle physics model dependent parts from the rest of the code. This means in practice that we have separated DarkSUSY into one set of routines, ds_core, which contains no reference to any specific particle model, as well as distinct sets of routines for each implemented model of particle physics. For supersymmetry, for example, all routines that require model-dependent information now reside in the mssm module. The advantage with this setup is that ds_core and the particle physics modules may be put in separate FORTRAN libraries, which implies that DarkSUSY can now have several particle physics modules side by side. The user then simply decides at the linking stage, i.e. when making the main program, which particle physics module to include.
For this to work, the main library needs each particle physics module to supply so-called interface functions (or subroutines), with pre-defined signatures and functionalities. Note that  Figure 1. Conceptual illustration of how to use DarkSUSY 6. The main program links to both the main library, ds_core, and to one of the available particle physics modules. User-replaceable functions are optional and may be linked to directly from the main program, or indirectly by including them in the various libraries. See examples/dsmain_wimp.F for an example of a main program that demonstrates typical usage of DarkSUSY for different particle physics modules. a particle physics module does not have to provide all of these predefined functions: which of them are required is ultimately determined only when the user links the main program to these libraries. Assume for example that the main program wants to calculate the gammaray flux from DM (a functionality provided by ds_core). This is only possible if the particle module provides an interface function for the DM contribution to the local cosmic ray source function; if it does not, the main program will not compile and a warning is issued that points to the missing interface function. If, on the other hand, interface functions required by direct detection routines would be missing in this example, this would not create any problems at either runtime or the compile stage.
In addition, we have added the concept of replaceable functions, which allows users to replace essentially any function in DarkSUSY with a user-supplied version. DarkSUSY ships with dedicated tools to assist you setting up both replaceable functions and new particle physics modules (see also Appendix F). Fig. 1 illustrates these concepts by showing how a typical program would use DarkSUSY; below, we describe each of them in more detail.

The DarkSUSY core library
As introduced above, the library ds_core is in some sense the heart of the new DarkSUSY 6 , offering all the functionality that a user typically would be interested in without having explicitly to refer to specific characteristics of a given particle physics model (after initialization of such a model). The main library thus contains routines for, e.g., cosmic ray propagation, solar models, capture rates for the Sun/Earth, a Boltzmann solver for the relic density calculation, yield tables from annihilation/decay etc. None of the routines in the main library contains any information about the particle physics module. Instead any information needed is obtained by calling an interface function that resides in the particle physics module which is linked to (see below). The source code for all functions and subroutines in the main library can be found in the src/ directory of the DarkSUSY installation folder, with subdirectory names indicating subject areas as summarized in Table 1.

Particle physics modules
The particle physics modules are collected in src_models/ and contain the parts of the code that depend on the respective particle physics model. Examples are cross section calculations, yield calculations etc. The routines in the particle physics module have access to all routines in ds_core, whereas the reverse is in general not true (with the exception of a very limited set of interface functions that each particle module provides).
DarkSUSY 5 and earlier was primarily used for supersymmetric, and specifically neutralino DM, and those parts of the code now reside in the MSSM module mssm. However, many people used DarkSUSY even before for e.g. a generic WIMP setup, which was doable for parts of the code, but not all of it. We now provide a generic WIMP module generic_wimp that can be used for these kinds of calculations in a much more general way. In a similar spirit, we also provide a module generic_decayingDM for phenomenological studies of decaying DM scenarios. As an example for an actual particle physics model other than supersymmetry, DarkSUSY 6 furthermore includes a module silveira_zee which implements the DM candidate originally proposed by Silveira and Zee [40] and which now is often referred to as Scalar Singlet DM [25,41,42]. We furthermore provide a set of simplified renormalizable models with light mediators that are fully confined to a dark sector, in the vdSIDM module; the characteristic feature of these models are Sommerfeld-enhanced DM annihilation rates, strong self-interactions and late kinetic decoupling. We also include an empty model, empty, which is of course not doing any real calculations, but contains (empty versions of) all interface functions that the core library is aware of -which is very useful for debugging and testing purposes. Designing new particle modules is a straight-forward exercise, see Appendix F.2, and it is generally a good idea to start with the most similar module that is already available.
The source code for the particle modules can be found in the directory src_models/ of the DarkSUSY installation folder, each of the subdirectories (e.g. src_models/mssm/, src_models/generic_wimp/) typically reflecting a (sub)subdirectory structure analogous to what is shown in Table 1 for the core library. A given model, defined by its model parameters, is then in general initialized with a call to a routine like dsgivemodel_decayingDM for decaying DM, or dsgivemodel_25 for a pMSSM model with 25 parameters (see the Appendix for more examples), followed by a call to dsmodelsetup. Many dark matter models, furthermore, constitute only relatively simple extensions to the standard model, inheriting most of its structure. For convenience, we therefore also provide various auxiliary routines, in src_models/common/sm, that each particle module automatically has access to, and which return basic standard model quantities like, e.g., the masses of standard model particles and their running (thus, additional BSM effects have to be implemented in the respective particle module).

Interface functions
Interface functions (or subroutines) are the limited set of functions that are provided by the particle physics module and which routines in ds_core might need to call; if they are not provided, a main program that makes use of those routines in ds_core will not compile. For this inter-library communication to work seamlessly, the behaviour of each interface function must be uniquely defined by its usage in the core library, and is explicitly specified in the function header. Furthermore, interface functions must have a well defined signature; the number and types of all their arguments, in other words is fixed.
Examples of interface functions include dsddsigma that returns the equivalent DM nucleus cross section (relevant for direct detection), dscrsource that returns the source term for DM-induced cosmic rays (relevant for indirect detection), and dsanxw that returns the invariant annihilation rate (in the case of WIMPs). All interface functions contain the keyword 'interface' in the function or subroutine header. A complete list of all interface functions known to the core library, about a dozen for this version 6 of the code, can be found by looking in the empty module.

Replaceable functions
The concept of replaceable functions introduces the possibility to replace a DarkSUSY routine with one of your own. The way it works is that the user-provided function will be linked when you make your main program instead of the DarkSUSY default one. If, for example, you want to replace the yields from a typical final state of DM annihilation or decay (likebb) with a new function of your own -e.g. because you are interested in comparing the tabulated Pythia [26] yields with those provided by PPPC [43] -you create your own version of the routine dsanyield_sim and let DarkSUSY use this one instead. Note that both routines in ds_core and in any of the particle physics modules can be replaced in this way, including interface functions. 3 To help you with this setup, we provide tools that can create (or delete) a replaceable function from any DarkSUSY function, and update the makefiles to use this user-supplied function instead. We also provide a simple way of managing large 'libraries' of user-supplied functions, via a list imported by the makefiles, where the user can determine on the fly which user-replaceable functions should be included and which ones should not.

Halo models
Several routines in the core library need information about which DM halo should be adopted for the calculations. With this DarkSUSY version, we introduce a new and flexible scheme that avoids pre-defined hardcoded functions to describe the DM density profiles, and allows to consistently define different DM targets at the same time. For convenience, we still provide several pre-defined options for such halo parameterizations, and the user can choose between, e.g., the Einasto [44] and the Navarrow-Frenk-White profile [45], or read in any tabulated axi-symmetric (or spherically symmetric) profile. On a technical level, the halo models are handled by the dsdmsdriver routine which contains a database of which halo profiles the user has set up, and consistently passes this information to all parts in the code where it is needed. For more details, see Section 6.

Thermal production and decoupling
The evolution of the DM phase-space density f (p) is described by the Boltzmann equation, which in an expanding Friedmann-Robertson-Walker universe reads [46,47] (4.1) Here, H =ȧ/a is the Hubble parameter, and all interactions between DM and standard model particles are contained in the collision term C[f ]. For WIMPs, for example, the righthand side is sufficiently large at high temperatures that these DM candidates have been in full equilibrium with the thermal bath in the very early universe. In general, one can distinguish two types of interactions: i) processes that change the number of the involved DM particles, or other new particles close in mass, and ii) proper scattering events that are particle number-conserving. DarkSUSY provides advanced routines to handle both, and to compute the corresponding observables with high precision; the DM relic density for the first class of processes (Section 4.1) and the mass of the smallest protohalos for the second (Section 4.2).

Chemical freeze-out and relic density
To compute the DM relic density, one is typically not only interested in the interactions of the DM particle itself, but also in the possibility that other new particles may be thermally produced as well. To take into account the effect of such coannihilations, one needs to consider an equation like (4.1) for each particle that is involved, and the collision terms will in general also involve interactions between those new particles. In most scenarios the stability of DM is guaranteed by a discrete parity, like R-parity in supersymmetry, implying that all other particles will eventually decay into the DM particle. In DarkSUSY we therefore currently only keep track of the total number density of particles n ≡ N i=1 n i , where n 1 denotes the number density of the lightest (i.e. DM) particle. The evolution equation for this quantity is obtained by summing over all individual Boltzmann equations of the form (4.1), and then integrating over the DM momenta [48]: where n eq denotes the total number density in thermal equilibrium with the heat bath and ... denotes a velocity average. Note that this assumes kinetic decoupling (see below) to happen well after chemical decoupling. Motivated by the WIMP case, DarkSUSY provides this velocity average under the assumption that the DM particles follow a Maxwell-Boltzmann distribution. In that case, the thermal average of the effective cross section weighted by the relative velocity of the incoming particles is given by where K 1 (K 2 ) is the modified Bessel function of the second kind of order 1 (2), Tγ is the heat bath temperature (which can, but does not have to coincide with the photon temperature T ), and s is one of the Mandelstam variables. The effective invariant rate introduced above is defined as where the g i denote internal degrees of freedom and the relative momenta are given by with p eff ≡ p 11 . The relation between the invariant rate W ij and the more familiar cross section σ ij for the annihilation of two particles i and j is 4 where v ij is the invariant relative velocity, i.e., the magnitude of the velocity of one particle in the rest frame of the other, (4.8) 4 Further useful relations between Wij and σij include Wij = 4pij √ sσij = 4EiEjσijv Møl,ij , where v Møl,ij = [|vi − vj| 2 − |vi × vj| 2 ] 1/2 is the Møller velocity. For the special case of self-annihilation or particle-antiparticle annihilation of particles of mass m, this becomes where v = s(s − 4m 2 )/(s − 2m 2 ) = 2vcm/(1 + v 2 cm ) and vcm = (1 − 4m 2 s ) 1/2 .
Note that the sum in Eq. (4.3) only involves processes that change the total number n of non-SM particles. For more details on the standard way of calculating the relic density, as implemented in DarkSUSY, see Refs. [48,49]. For complex particle models involving many coannihilations, like the MSSM, the calculation of the invariant rate W eff can be rather time-consuming. Given that it is independent of temperature, DarkSUSY therefore tabulates it once and for all for every set of model parameters (with a particularly dense sampling of potentially critical regions like thresholds and resonances). The thermal average is then performed by means of an adaptive gaussian integration, using a spline routine for interpolation between the tabulated points. In order to integrate the density evolution equation (4.2), we follow Ref. [48] and rephrase it as an evolution equation for the abundance, Y = n/s (with s being the entropy density); instead of time, we use the dimensionless quantity x = m χ /T as our integration variable. The numerical integration is still subtle, because the differential equation is stiff, and we developed a dedicated implicit trapezoidal method with adaptive stepsize to guarantee a high precision result. For more details on the numerical implementation of these routines, see Ref. [24]. As a default for the energy and entropy degrees of freedom in the early universe, finally, we implemented the results by Drees et al [50].
Despite those similarities at the detailed implementation level, the relic density routines are now considerable more general than in previous versions of the code. This allows to link them to a given particle physics model in a simple and straightforward way. The effective invariant rate W eff in Eq. (4.4), in particular, must be supplied by the particle module in the form of the interface function dsanwx. In order to calculate the thermal average in Eq. (4.3), the main library also needs to know the masses and internal degrees of freedom of the coannihilating particles and, in order to obtain an as precise result as possible, the location of possible resonances or thresholds in W eff (p eff ). This information is provided by the only other interface function in the relic density part, dsrdparticles. Another new feature since DarkSUSY 6.1 is that the relic density routines now can handle heat bath temperatures different from the photon temperature, parameterized in terms of the ratio This ratio is returned by the replaceable function dsrdxi which simply returns a constant factor of unity in its default version -as supplied by the core library -but can take any user-supplied form more appropriate for the particle model in question. 5 The DM relic density (of both components, in case the DM particle is not its own antiparticle) in terms of Ωh 2 is provided by a call to the function dsrdomega, which manages the calls to the DarkSUSY relic density routines. In their full version, these routines return the relic density with a precision beyond per-cent accuracy -assuming that W eff is provided with corresponding precision -but there are also options for considerably faster calculations (for example by taking into account only coannihilations up to mass differences of f co = 1.4). This 5 For example, if DM is confined to a dark sector that has been in thermal contact with the visible sector until some decoupling temperature T dc , entropy is afterwards typically conserved separately in the two sectors. This implies that the temperature ratio changes with the ratios of the degrees of freedom (see, e.g., Ref. [51,52]): where g SM,DS * are the entropy degrees of freedom in the visible and dark sector, respectively. In Appendix E, we will present a module already implemented in DarkSUSY where we make use of this relation.
can be very useful in situations where an exceptional error of a few per cent is an acceptable trade-off in view of significantly improved performance with respect to computation time (e.g. in large global scans). For comparison, the DM density today is now observationally determined with a precision of about 1% [1].

Kinetic decoupling and the smallest protohalos
Even after the DM particles have chemically decoupled from the thermal bath, they are generally still in local equilibrium through frequent scattering processes (typically with SM particles). As long as this is the case, the DM temperature T χ ≡ gχ 3mχnχ is the same as the temperature of the scattering partners. After kinetic decoupling, the DM 'temperature' will simply decrease due to the expansion of the universe and, for non-relativistic DM, scale like T χ ∝ a −2 . It is thus natural to define the moment of decoupling as the transition between these two regimes [53]. Allowing for the scattering partners again to have a different temperature (Tγ) than the photons (T ), this implies The evolution of T χ (T ) is obtained by solving the second moment of the full Boltzmann equation 4.1. Introducing x ≡ m χ /T and y ≡ m χ T χ s −2/3 , where s is the total entropy density, this is to leading order in p 2 /m 2 χ given by [54] Here, g * S is the number of effective entropy degrees of freedom, y eq is the value of y in thermal equilibrium and γ denotes the momentum transfer rate, where the sum runs over all DM scattering partners (counting, e.g., electrons and positrons separately), k ≡ |k| is their momentum and ω their energy. The scattering amplitude squared in this expression, |M| 2 , is understood to be summed over all internal (spin or color) degrees of freedom, including initial ones. If it is not Taylor expandable around t = 0, one has to make the replacement [57, 58] in the above expression.
As long as DM is in local thermal equilibrium, the frequent scattering processes will wash out any small-scale perturbations in the DM fluid as soon as they enter the horizon. For temperatures close to T kd , first a remaining viscous coupling and then the free-streaming of the DM particles generates an exponential cutoff in the power spectrum of density fluctuations [59], which translates to the smallest mass a protohalo can have once structure formation enters the non-linear regime. Acoustic oscillations of the DM fluid also generate an effective exponential cutoff in the power spectrum [60,61]. Which of the two mechanisms dominates, i.e. leads to a larger cutoff mass M cut , depends on the combination of decoupling temperature T kd and DM mass m χ [53].
In DarkSUSY, the kinetic decoupling temperature is obtained with a call to the routine dskdtkd, which solves the Boltzmann equation and determines the asymptotic solution of T χ /T to return T kd as defined in Eq. (4.11). A call to dskdmcut then takes T kd and the DM mass as input, and returns M cut (as given in Ref. [53]). These functions reside in the core library. There are three interface functions that a particle physics module must provide if the main program uses this functionality: dskdm2 returns the full scattering matrix element squared, evaluated at t = 0 or averaged over t, while dskdm2simp returns only the leading contribution for small ω (expressed as a simple power-law in ω, in which case there exists an analytic solution for T kd [47]). Lastly, the particle physics module must provide a routine dskdparticles which, in analogy to the routine dsrdparticles for the case of chemical decoupling, sets the location of resonances in |M| 2 . This allows the corresponding routine in the core library to perform the energy integration in Eq. (4.13) to a much higher precision. Non-standard heat bath temperatures Tγ = T are automatically handled by supplying a user-defined function dsrdxi, in the same way as described in the previous Section.

Direct detection
The rate for nuclear recoil events in direct detection experiments (often expressed in dru, or differential rate unit, i.e. counts/kg/day/keV) can be written as The sum here runs over nuclear species in the detector, c T being the detector mass fraction in nuclear species i. m T is the nuclear (target) mass, and µ χT = m χ m T /(m χ + m T ) is the reduced DM-nucleus mass. Furthermore, ρ 0 χ is the local DM density, v the DM velocity relative to the detector, v = |v|, and f (v, t) is the (3D) DM velocity distribution. In order to impart a recoil energy E R to the nucleus, the DM particle needs a minimal speed of v min = M T E R /2µ 2 χT . Finally, dσ χT /dE R describes the unpolarized differential scattering cross section for DM scattering off a target nucleus of mass m T , differential with respect to the nucleus recoil energy E R in the initial rest frame of the target nucleus. This is of the form Here v is the WIMP-nucleus relative velocity, and S is a Lorentz-invariant scattering factor (the overline indicates a sum over final polarizations and an average over initial polarizations). If the WIMP and the nucleus are treated as elementary particles, the scattering factor is given by S = |M/(4k · p)| 2 , where k and p are the initial four-momenta of the WIMP and nucleus and M is the invariant amplitude in the normalization of the Review of Particle Physics [62]. For a nucleus of finite size, with nonrelativistic single-nucleon interactions governed by WIMP-nucleon contact operators O i , one can separate the scattering factor S into the contributions of single operators and their interference. One can write The sum in Eq. (5.3) is over the operator indices i and j. A "square" term S ii is the contribution from a single operator O i , and a "cross" term S ij with i = j is half of the contribution from the interference of operators O i and O j (the other half is the term S ji = S ij ). The sum in Eq. (5.4) is over proton and neutrons (N = p, n), the G N i 's are effective WIMPnucleon coupling constants generalizing the Fermi coupling constant G F , given in the same units ( 3 c 3 GeV −2 ), and the P N N ij 's are dimensionless factors containing the nuclear structure functions for specific nucleon currents (see below for details).
In place of the differential cross section dσ/dE R , DarkSUSY provides an equivalent cross sectionσ χT , and partial cross sectionsσ ij , defined by the relatioñ Here µ χT = m χ m T /(m χ +m T ) is the WIMP-nucleus reduced mass, and the prefactor is equal to the maximum recoil energy for elastic scattering, Notice thatσ χT is in general not equal to the total WIMP-nucleus cross section σ χT . The latter is obtained by integration of dσ χT /dE R over the recoil energy E R and the dependence of the nuclear structure functions on E R may be important. In the limit of vanishing E R and v, however, the two quantities coincide. Moreover, the definition ofσ χT applies to both elastic and inelastic scattering, and is merely a convention.
In the usual spin-independent/spin-dependent treatment, the equivalent cross sectioñ σ χT contains only two terms, in this case without interference terms, the spin-independent term σ SI and the spin-dependent term σ SD , In the convention of [63], these two terms arise from the operators O 1 and O 4 , thus we will use the indices 1 and 4 forσ SI =σ 11 andσ SD =σ 44 . The coupling constants G p 1 , G n 1 , G p 4 , G n 4 are hence equal to the previous DarkSUSY constants G p s , G n s , G p a , G n a , respectively. They are related to the coupling constants f p , f n , a p and a n in [64] by the following equations.
Here s χ is the spin of the WIMP, J is the spin of the nucleus, A and Z are the mass number and atomic number of the nucleus, q = √ 2m T E R is the magnitude of the momentum transfer, F (q) is the spin-independent nuclear form factor (e.g., the Helm form factor) normalized to F (0) = 1, and the S pp (q), S pn (q), and S nn (q) are the nucleus spin structure functions as defined in [65]. The factor 4s χ (s χ + 1)/3 has been introduced in P N N 44 because the original spin structure functions S N N (q) were computed for s χ = 1 2 . In the nonrelativistic effective operator approach of [63], there are several operators O i with corresponding coupling constants G p i and G n i . As physical quantities, the constants G N i are equal to the constants c N i in [66], (20)). So numerically in Eq. (77) of [63] and are given in terms of the WIMP and nuclear response functions of [66] by Alternative formulation (where the factor q 2 /m 2 N in front of the Rs in the equation above is included in the R s): in Eqs. (73) and (74) of [63], which are related to the functions W N N α in Eq. (41) of [66] by As examples of this formalism, for the spin-independent operator O 1 one has and for the spin-dependent operator O 4 one has Comparison of these expressions with Eqs. (5.10)-(5.15) shows that both the usual SI/SD approach and the nonrelativistic effective operator approach of [63] are included in the same formalism of Eq. (5.4) and that they differ only by the choice of nuclear structure functions appearing in P N N ij and the number of nonzero WIMP-nucleon coupling constants G N i . Specific choices of nuclear structure functions can be selected by calling dsddset('sf', label), where the character variable label indicates the set of structure functions. For the default option ('best'), e.g., the code automatically picks the best currently available structure function (depending on the nucleus). The function returning the value ofσ χT is to be provided by an interface function dsddsigma(v,Er,A,Z,sigij, ierr) residing in the particle physics module, where on input v=v, Er=E R , A=A, Z=Z and on output the 27×27 array sigij contains the (partial) equivalent cross sectionsσ ij in cm 2 and the integer ierr contains a possible error code. The order of the entries in sigij corresponds to that of the independent nonrelativistic operators O i ; for the first 11 entries, we use the same operators and convention as in Ref. [63], while for the last 16 entries we add the additional operators discussed in Ref. [67]. In particular, sigij(1,1) is the usual spin-independent cross section and sigij(4,4) is the usual spin-dependent cross section. In addition, the direct detection module in DarkSUSY provides utility functions that can be used in the computation of the cross section. For example, the subroutine dsddgg2sigma(v, er,A,Z,gg,sigij,ierr) computes the (partial) equivalent scattering cross sectionsσ ij for nucleus (A, Z) at relative velocity v and recoil energy E R starting from values of the G N i constants in gg, with nuclear structure functions set by the previous call to dsddset. The actual nuclear recoil event rate as given in Eq. (5.1), finally, is computed by the function dsdddrde. The latter two functions are independent of the specific particle physics implementation and hence are contained in the core library.

Halo models
The distribution of the DM density ρ χ (r) in the Milky Way and external objects, like nearby dwarf spheroidal galaxies, is a crucial information for calculating rates both in highly local direct detection experiments (as described in the previous Section) and for indirect DM searches that probe the DM distribution inside a much larger volume (as described in the following Sections). Observationally, however, the distribution of DM on scales relevant for DM searches is often only poorly constrained (with the exception of classical dwarf spheroidal galaxies -though even here uncertainties in halo properties relevant for DM searches may be somewhat larger [68] than what is typically quoted [69]). The situation is somewhat improved when instead referring to the results of large N -body simulations of gravitational clustering, which consistently find that DM halos on average are well described by Einasto [44] or Navarro-Frenk-White profiles [45], with the halo mass being essentially the only free parameter (after taking into account that the halo concentration strongly correlates with the halo mass [70]). On the other hand, there is a considerable halo-to-halo scatter associated to these findings, so that it remains challenging to make concrete predictions for individual objects -in particular if they are located in cosmologically somewhat 'special' environments like in the case of the Milky Way and its embedding in the Local Group. Even worse, baryonic physics can have a large impact on the DM profiles, especially on their inner parts most relevant for indirect detection, and even though hydrodynamic simulations taking into account such effect have made tremendous progress in recent years [71][72][73][74], there is still a significant uncertainty related to the modelling of the underlying processes. In light of this situation, there is a considerable degree of freedom concerning halo models and the DM density profiles, and a code computing observables related to DM should be able to fully explore this freedom. One can also note that even in the Milky Way, the local halo density has some uncertainty. If we assume a given halo profile, the uncertainty goes down as we get constraints from the galaxy as a whole, but the fluxes of cosmic ray particles locally depend on the density not only here, but also further away which cause a dependence on the halo profile and not just the local density.
In DarkSUSY 6, the implementation of dark matter halo models and related quantities in the library ds_core follows a new and highly flexible scheme, avoiding pre-defined hardcoded functions. For convenience a few pre-defined options are provided, however these can be either complemented by other profiles eventually needed, or the entire sample configuration can be simply replaced by linking to a user-defined setup -in both cases without editing routines provided in the release version of the code. A further improvement compared to previous versions of DarkSUSY is that different dark matter density profiles, possibly referring to different dark matter detection targets, can be defined at the same time: E.g., one can easily switch back and forth from a computation of the local antiproton flux induced by dark matter annihilations or decays in the Milky Way halo to the computation of the gammaray flux from an external halo or a Milky Way satellite within the same particle physics scenario. A special label is reserved to indicate which profile corresponds to the Milky Way, such that rates that are Milky Way specific, such as the local contribution to antimatter fluxes, can be computed for this (and only this) profile. Finally the present implementation simplifies the task of keeping track of consistent definitions for related quantities, such as, e.g., a proper connection between the dark matter profile and the source function for a given dark matter yield (see Section 8), or calling within an axisymmetric coordinate system a spherically symmetric function (and preventing the opposite).
At the heart of the implementation, there is the subroutine dsdmsdriver which contains the complete set of instructions necessary to act as an interface to all quantities related to the DM density profiles. 7 For example, it checks the scaling of the various DM source functions (Section 8) with the DM density ρ χ -namely ρ for decaying dark matter and ρ 2 for dark matter pair annihilations -and passes this information to the routines for propagations of charged cosmic rays in the Galaxy (Section 10); the same goes for line of sight integration routines (needed, e.g. for the computation of gamma-ray fluxes, see Section 9). dsdmsdriver also contains information on whether the dark matter density profile is spherically symmetric or axisymmetric, and consistently handles requests of functions assuming one or the other. Furthermore, it can be used for initialization calls, for instance to set parameters for a given parametric density profile, and test calls, for instance to print which dark matter profile is currently active within a set of available profiles.
While a specific dsdmsdriver routine should match the user needs in the problem at hand, a sample version is provided with the DarkSUSY release, illustrating the flexibility of the setup. In particular, the simple version included in DarkSUSY 6 assumes that the DM density profile is spherically symmetric and does not include dark matter substructure; it allows to choose as dark matter profile one among three parametric profiles, namely the Einasto [44], the Navarro-Frenk-White [45] and the Burkert [75] profiles, or a profile interpolated from a table of values of the dark matter density at a given radius. The present version also implements a halo profile database for later use, where the code allows to associate a given input label to every full set of entries specifying a halo profile; such a profile can be reloaded at any time when needed (but profiles can also be declared as 'temporary', in which case only the latest set temporary profile is available at any given time). For example, all indirect detection flux routines have this label among their input parameters, so that in case of several dark matter detection targets or several profiles for the same targets it becomes unambiguous which profile is being considered. For objects in this database, one can furthermore save or read tabulated quantities from the disc, such as the Green's function needed for the computation of the local positron flux. In Appendix F.3, we provide more details on how to use this driver routine in praxis, illustrated by example programs shipped with the DarkSUSY release, as well as how to expand it such as to include different setups and profile parameterizations.
Let us finally mention that the local DM velocity profile that enters in the direct detection rate, Eq. (5.1), is in principle not independent of the chosen DM density profile. For a spherically symmetric and isotropic system, e.g., the two profiles are related by the Eddington equation [76,77]. A fully self-consistent implementation of phase-space distributions will be available with a later DarkSUSY version; until then, the user can freely choose a DM velocity distribution among those provided in src/dmd_vel -but should keep in mind this consistency requirement when comparing direct detection rates to, e.g., the gamma-ray flux from the galactic center (which requires choosing a density profile). Concretely, it is the function dshmuDF that returns the 3D distribution function f (v) needed by the direct detection routines. It allows to switch between various pre-implemented functional forms, including tabulated velocity profiles, but can of course also be replaced by an arbitrary function supplied by the user (c.f. Section 3.4).

Dark matter self-interactions
The general expectations for DM profiles described in the previous section apply to cold, collisionless DM. While remarkably successful at large cosmological scales, however, the CDM paradigm is less well tested at smaller scales. In fact, DM self-interactions are phenomenologically still allowed to be as strong as the interaction between nucleons, which is much stronger than current limits on DM interacting with SM particles. Partially also triggered by the absence of undisputed DM signals in traditional searches, the possibility that DM could be self-interacting [78] and thereby leave imprints on cosmological observables related to structure formation has thus seen greatly increased interest in recent years, both from an astrophysical and a model-building perspective. In fact, it has been pointed out [55,[79][80][81][82][83][84][85] that self-interacting DM (SIDM) could address the most pressing potential small-scale problems of ΛCDM cosmology [86], referred to as 'core-vs.-cusp" [87,88], 'too-big-to-fail' [89,90], 'diversity' [91,92] and (for late kinetic decoupling) 'missing satellites' problems [93,94]. While those specific discrepancies with the ΛCDM paradigm may well turn out to be either due to poorly modelled baryonic effects or observational uncertainties, it is clear that SIDM can leave observable imprints which, in case of an unambiguous detection, would significantly narrow down the range of possible DM particle models.
In the literature, one typically uses the momentum transfer cross section, as the relevant quantity to measure the impact of DM self-interactions on the halo structure, where σ is the standard cross section for DM-DM scattering (which equals σ T for isotropic scattering). This has the advantage of regulating large forward-scattering amplitudes, which should not affect the DM distribution. 8 In particular, is very roughly the typical scale of DM self-interactions of cosmological relevance: cross sections in this ballpark may leave observable imprints and possibly address the afore-mentioned small-scale structure problems, while much smaller cross sections have no impact on the structure of DM halos; much larger values of σ T /m χ are ruled out, on the other hand, in particular from observations of galaxy clusters. For a detailed review that summarizes both various constraints on DM self-interactions and models that have been discussed in the literature, see Ref. [97]. In DarkSUSY, σ T /m χ is computed by an interface function dssigtm provided by the particle module. In general, this function can depend on the relative velocity of the scattering DM particles -which from a model-building point of view has the advantage that one can evade the strong cluster constraints but still have sizeable interaction rates at the scale of dwarf galaxies (where the supposed discrepancies with ΛCDM expectations have been reported). In such a situation, one needs to average over the velocity distribution of DM particles in the halo; this is done by dssigtmav provided by the core library. The core library furthermore provides several auxiliary routines, to be used by any particle module, for the commonly encountered specific transfer cross sections in the presence of (attractive or repulsive) Yukawa potentials. These result from the exchange of single mediator particles of mass m med with coupling α χ ≡ g 2 /4π. In particular, dssisigmatborn returns the resulting σ T in the Born limit, α χ m χ m med , where the cross section can be calculated perturbatively [98], and dssisigmatclassical returns σ T in the classical limit, m χ v m med , for which we use the parameterizations from Ref. [99]. For the intermediate (resonant) regime, dssisigmatres implements the analytical expressions from Ref. [95] that result when approximating the Yukawa potential by a Hulthén potential in the limit where scattering dominantly proceeds via an s-wave (which is only guaranteed for m χ v m med ). We stress that these specific auxiliary routines are only useful for models where the self-interaction is mediated by a single scalar or vector particle. For other cases (see e.g. Ref. [100] for a classification according to effective operators) the corresponding self-interaction cross sections are currently not preimplemented and thus have to be provided along with the corresponding particle module.

Particle yields
The annihilation or decay of DM particles in places with large astrophysical DM densities (like in our Milky Way halo) produces SM particles which either propagate on straight lines like gamma rays or neutrinos (Section 9) or, in the case of charged particles, are deflected by stochastically distributed inhomogeneities in the Galactic magnetic fields (Section 10). In either case, the DM signal ultimately only depends on the local injection rate of some stable (cosmic ray) particle f , per volume and energy, Here, ρ χ is the DM density (of the respective component, in case of multi-component DM) and the ensemble average ... is taken over the DM velocities; in principle, it depends on the spatial location x, but in many applications of interest this can be neglected. The particle source terms S n (E f ) encode the full information about the DM particle physics model. For a typical WIMP DM candidate, e.g., only the annihilation part (n = 2) contributes, where σ i is the annihilation cross section of two DM particles into final state i and dN i /dE f is the resulting number of stable particles of type f per such annihilation and unit energy. N χ is a symmetry factor that depends on the nature of DM; if DM is (not) its own antiparticle we have N χ = 2 (N χ = 4). For decaying DM, on the other hand, we have where Γ i denotes the partial decay widths. Let us stress, however, that Eq. (8.1) is much more general in that it encapsulates also DM that is both annihilating and decaying, multicomponent SM, as well as DM models with an internal Z n symmetry [101][102][103][104].
In DarkSUSY 6, the particle source term S n (E f ) therefore takes a central rôle as essentially the only interface between the particle physics modules and the indirect detection routines provided by the core library (apart from the routines providing neutrino yields from the Sun or Earth, see Section 11, which rely on a slightly different setup). Concretely, there are two interface functions that provide this quantity, dscrsource and dscrsource_line. The former returns S n (E f ) for any continuous energy distribution of the particle f , while the latter returns location, strength and width of (almost) monochromatic contributions -which result if (at least) one of the two particles in a two-body channel i is the stable particle f itself.
In general, it is up to the particle module how to provide S n (E f ). For convenience, however, the core library contains a function dsanyield_sim to provide dN i /dE f for any two-particle SM model final state. Those are based on Pythia 6.426 [26] runs. The Pythia simulations are run for 30 different masses from 3 GeV to 20 TeV and for annihilation into dd, uū, ss, cc, bb, tt, gg, W + W − , Z 0 Z 0 , µ + µ − and τ + τ − . 9 For each of these annihilation channels we store the resulting spectra of γ's, positrons, antiprotons, anti-deuterons, ν e /ν e , ν µ /ν µ , ν t au/ν τ . For the muon neutrinos we also simulate neutrino-nucleon interactions at neutrino detector on Earth and calculate the muon yields at the neutrino-nucleon vertex (i.e. a yield/volume) and at a plane at the detector (i.e. a yield/area). These latter neutrino-nucleon interactions are simulated with nusigma [105]. All these yields are available in dsanyield_sim which interpolates the supplied tables in both energy and mass.
For the anti-deuterons, we have performed a simulation with Pythia 6.426 where we from the event record connect anti-neutrons and anti-protons if they are within a coalescence momentum p 0 . These yields can then in DarkSUSY be retrieved for different coalescence momenta p 0 from 0 to 300 MeV. For the same Pythia setup, we have also performed simulations of the yields of anti-deuterons on the Z 0 resonance that have been compared with ALEPH measurements [106]. The best fit coalescence momentum and its errors are available in DarkSUSY (and for the default simulation the best fit value is 202 MeV). The parameter dbp0bar determines which p 0 value to use. If zero, the best-fit p 0 value is used, if non-zero it determines how much above or below the best fit p 0 to use (in units of standard deviations from the fit, σ).
It is important to realize though that these yield functions in DarkSUSY are a prime example of a replaceable function. If the user wants to supply yields calculated in a specific scenario or adding a specific new process, or just using a different simulation/analytic calculation, it is fairly easy to replace either the yield function dsanyield in the chosen particle physics module or the simulation yields in dsanyield_sim in src/.
The setup described above is for WIMP annihilations in e.g. the galactic halo, a dwarf galaxy or any other similar environment where the background density is low. For annihilations in the Sun and the Earth, we have similar routines that will be described in Section 11 below.
9 Gamma rays and Neutrinos from the halo Gamma rays produced in the galactic halo propagate on straight lines and hence point directly back to their sources. Along with the fact that they are expected to be copiously produced in many DM models, and can carry distinct spectral features that would allow to relatively easily distinguish signals from astrophysical backgrounds, this is the reason they are sometimes referred to as the golden channel of indirect DM searches [107]. Unless produced in celestial bodies like the Sun or the Earth (see Section 11), the propagation of neutrinos follows the same simple pattern. Indeed, while not as common as for gamma rays, and not as easily distinguished because of the poorer energy resolution of neutrino telescopes, neutrino spectra may also carry characteristic features related to their DM origin [108][109][110][111].
For a telescope pointing in the direction ψ, the expected DM-induced differential flux in gamma rays or neutrinos -i.e. the expected number of particles per unit area, time and energy -from a sky-region ∆ψ is thus given by a line-of-sight integral where the local injection rate dQ/dE was earlier introduced in Eq. (8.1). For decaying DM, the above line-of-sight integral always factorizes into the particle source term S 1 given in Eq. (8.3) and a term that only depends on the DM distribution, For annihilating DM, the corresponding factorization strictly speaking only holds if the annihilation rate is independent of velocity: While notable exceptions exist (in particular for resonances [112], p-wave annihilation [113] and Sommerfeld-enhanced annihilation [114]), this is a commonly encountered situation and hence of general interest.
DarkSUSY therefore provides the functions dscrgaflux_dec and dscrgaflux_v0ann that take J dec (or J ann ) as input and return the fluxes given in the above two equations for gamma rays (as well as corresponding routines dscrnuflux_dec and dscrnuflux_v0ann for neutrinos). Here, the subscript _v0ann refers to the fact that, for the purpose of those routines, S 2 is evaluated in the limit of vanishing relative velocity of the annihilating DM pair. While this is the only situation of practical interest in many DM models, future DarkSUSY versions will offer support for a full velocity dependence of this quantity. In analogy with dscrsource_line, the core library furthermore provides routines dscrgaflux_line_dec and dscrgaflux_line_v0ann (and correspondingly for neutrinos) that return strength, width and location of monochromatic ('line') signals.
The values for J dec and J ann can be supplied by the user, e.g. as taken from the literature, or be calculated with the routines dsjfactor and dsdfactor, respectively. Those latter routines take as input a halo label as described in Section 6, and only take into account the contribution from the smooth halo profile (thus neglecting a possible enhancement of J ann due to the existence of DM substructure, which would have to be added by hand). If such a halo label is available, one can also more conveniently call dsgafluxsph instead. This function directly returns the gamma-ray flux, fully automatically calculating the required line-of-sight integrals and adding decaying and annihilating DM components, depending on which particle model is initialized.

Cosmic ray propagation and antimatter signals
There is a clean asymmetry between particles and antiparticles in the standard cosmic ray picture: The bulk of cosmic rays -protons, nuclei and electrons -are mainly "primary" species, i.e. particles accelerated in astrophysical sources and then copiously injected in the interstellar medium; "secondary" components, including antimatter, are instead produced in the interaction of primaries with the interstellar medium during the propagation. It follows that there is a pronounced deficit of antimatter compared to matter in the locally measured cosmic ray flux (about 1 antiproton in 10 4 protons). When considering instead a source term due to DM annihilations or decays, a significant particle-antiparticle asymmetry is in general not expected, and antiprotons, positrons and antideuterons turn out to be competitive indirect DM probes.
Charged particles propagate diffusively through the regular and turbulent components of Galactic magnetic fields. This makes it more involved for local measurements to track spectral and morphological imprints of DM sources than, e.g., for the gamma-ray and neutrino channels (though searches for spectral features in CR positron fluxes still lead to very competitive limits [115]). In fact the transport of cosmic rays in the Galaxy is still a debated subject: Most often one refers to the quasi-linear theory picture (with magnetic inhomogeneities as a perturbation compared to regular field lines) in which propagation can be described in terms of a (set of) equation(s) linear in the density of a given species, containing terms describing diffusion in real space, diffusion in momentum space (the so-called reacceleration), convective effects due to Galactic winds, energy or fragmentation losses and primary and secondary sources (see, e.g., Ref. [116] for a review). Dedicated codes have been developed to solve numerically this transport equation, including GALPROP [117], DRAGON [118] and PICARD [119].
Here we follow instead a semi-analytical approach, analogous to that developed for the USINE code [120]. In particular, we model the propagation of antiprotons and antideuterons by considering the steady state equation [121] We solve it for situations where i) the diffusion coefficient D can have an arbitrary dependence on the particle rigidity but can at most take two different values in the Galactic disc and in the diffusive halo, ii) the convective velocity u has a given fixed modulus and is oriented outwards and perpendicular to the disc, iii) the loss term due to inelastic collisions has an interaction time τ N which is energy dependent but spatially constant in the disc and going to infinity in the halo (corresponding to a constant target gas density in the disc and no gas in the halo), iv) the DM source Q is spatially axisymmetric and has a generic energy dendence. Under these approximations and assuming, as is usually done, that the propagation volume is a cylinder centred at the disc and that particles can freely escape at the boundaries of the diffusion region, Eq. (10.1) can be solved analytically by expanding N in a Fourier-Bessel series; the computation of the flux involves, at each energy, a sum over the zeros of a Bessel function of first kind and order zero, and a volume integral of the spatially dependent part in the axisymmetric source term Q, as introduced in Section 8, times a weight function depending on the given zero in the series (see [121] for further details). We are neglecting in this equation energy changing operators, such as the reacceleration term and the adiabatic flow term connected to the convection operator; while these are know to have a relevant impact at low energy, their effect can be approximately mimicked by a proper reshaping of the scaling of the diffusion coefficient at low rigidities. Since the path lengths for antiprotons and antideuterons are rather large, of the order of a few kpc, taking average values for parameters in the transport equation rather than the more realistic modelling that can be implemented in full numerical solutions, has no large impact in case of extended and rather smooth sources such as for DM. Eq. (10.1) neglects reacceleration effects, which may in general be relevant at rigidities below a few GV; however even this does not have a large impact in case of the species at hand, see, e.g., [122] for a comparison of results with numerical and semi-analytical solutions for cosmic-ray antiprotons. The power of our semi-analytic approach is that one can store the solution of the transport equation for a given mono-energetic source -provided by the functions dspbtdaxi and dsdbtdaxi for, respectively, antiprotons and antideuterons -and then apply these as weights to any particle source term S n (E f ) as introduced above. For antiprotons and antideuterons, this latter step is done in the functions that compute the local galactic differential fluxes from DM annihilation and decay, dspbdphidtaxi and dsdbdphidtaxi, respectively.
The structure we implemented gives a particularly clear advantage when the code is used to scan over many particle physics DM models, but only over a limited number of propagation parameters and DM density profiles. For such an application, it is useful to tabulate the 'confinement times' (returned by dspbtdaxi and dsdbtdaxi) over a predefined range of energy; this is done in the functions dspbtdaxitab and dsdbtdaxitab when calling the flux routines with an appropriate option. Such tabulations are time-consuming at first run but can be saved and re-loaded for later use; here the proper table association is ensured by a propagation parameter label setting system in analogy to the one implemented for the halo profile database.
The case for positrons is treated analogously, except that energy losses and spatial diffusion are the most important effects for propagating cosmic ray leptons. The transport equation we solve semi-analytically therefore has the form [123] ∂N ∂t = 0 = ∇ · (D ∇N ) + ∂ ∂p dp dt N + Q , where the functional form of D and Q can be chosen as for antiprotons and antideuterons, and the energy loss rate dp/dt can have a generic momentum dependence but needs again to be spatially constant. Assuming the same topology for the propagation volume and free escape conditions at the vertical boundaries (the radial boundary, being much farther away, is actually irrelevant when computing local densities), the solution of the propagation equation is given in terms of a Green's function in energy (the function dsepgreenaxi in the code) to be convoluted over the source energy spectrum at emission for a given particle DM candidate. This last step is performed by the function dsepdndpaxi returning the local positron number density, while dsepdphidpaxi converts it to a flux and is the function which should be called from the main file. The method to implement this solution is a slight generalization of the one described [123] and generalizes the one introduced in [124] for a spherically symmetric configuration to an axisymmetric system. The computation of the Green's function involves a volume integral over the spatially dependent part of the DM source function Q and implements the so-called method of image charges. It can again be CPU expensive for singular halo profiles, but its tabulation is always needed since the Green function appears in a convolution. The main limiting factor with respect to full numerical solutions is that one is forced to assume an (spatially) average value for dp/dt. However this has not a severe impact on our results for the local DM-induced positron flux since, especially for energies above 10 GeV, the bulk of the DM contribution to the local flux stems from a rather close-by emission volume; one thus just has to make sure to normalize dp/dt to the mean value for local energy losses, as opposed to the mean value in the Galaxy, which are mainly due to the synchrotron and inverse Compton processes.
While the transport equations (10.1) and (10.2) are essentially the same as considered in previous releases of the code, their implementation in the present release is completely new and appears to be numerically more stable. In particular cases with very singular DM profiles still give numerically accurate results and converge faster. Note however that the case of very singular DM profiles is also the one in which any propagation model is subject to a significant uncertainty related to the underlying physics, since propagation in the Galactic center region is difficult to model and probably rather different from what can be tested in the local neighbourhood by measurements of primary and secondary cosmic rays. Finally, the implementation in DarkSUSY 6 is more flexible regarding parameter choices, such as for the rigidity scaling of the diffusion coefficient and energy scaling of energy losses, in a framework which is now fully consistent for antiprotons, antideuterons and positrons.

Neutrinos from the Sun and Earth
We can write the capture rate in the Sun/Earth as (see e.g. [125]) where Ω v,i (w) is the capture probability on element i, f (u) is the velocity distribution of the WIMPs, u is their velocity (far away from the Sun/Earth), w = √ u 2 + v 2 is the velocity at the interaction radius, v is the escape velocity at that point and the sum is over all the elements in the Sun/Earth. The capture probability Ω v,i (w) depends on the scattering cross sections, including form factors of the different nuclei, the kinematics of the scattering process and the element abundance (as a function of position) in the Sun/Earth. In DarkSUSY we have different ways to calculate the capture rate and different solar models available. For the Sun, the most general calculation available contains a full numerical integration over the solar radius, the velocity distribution and the momentum exchange and sums over all the elements as given by the solar model by Serenelli et al [27] (up to 289 isotopes for spin-independent scattering and 112 for spin-dependent scattering). In Fig. 2 we show the capture rate in the Sun for spin-independent (left) and spin-dependent (right) scattering and the most important elements that contribute to the total capture rate. The example program to calculate the curves in Fig. 2 is available in examples/aux. We include a simple setup where the user can easily choose how many elements to include. To speed up the calculation when scanning over particle physics models, we write the total capture rate as so that we can tabulate the C i :s as a function of mass once and for all for a given halo and solar composition setup. The default is to use these tabuled capture rates, and if not available for the chosen halo model, they are recreated. DarkSUSY also contains simpler (and faster) routines where e.g. the momentum transfer is integrated analytically if the form factor is Gaussian. The main routine to calculate the capture rate in the Sun including the full form factor integration is dssenu_capsunnumff, whereas the tabulated version is dssenu_capsuntabff (which will call dssenu_capsunnumff to set up the tables if they do not exist). DarkSUSY obtains the couplings G p/n i above from a call to the interface routine dsddgpgn. In principle, these routines could use the more general routine dsddsigma as the direct detection routines do, but this would make tabulation very difficult, hence we currently require the particle physics module to provide dsddgpgn for the best capture rate calculation.
For the Earth, we use essentially the same set of routines, but currently assume that the form factor is Gaussian. This is motivated by the escape velocity being much lower in the Earth and hence the momentum transfer is also smaller, meaning that the dependence on the form factors is smaller as well. Also for the Earth, DarkSUSY tabulates the capture rates to speed up the calculation when the halo model does not change. The main routine to calculate the capture rate in the Earth is dssenu_capearthnum.
To get the number of WIMPs in the Sun/Earth, we have to solve the evolution equation The annihilation term C A N 2 depends on the potential and temperature in the Sun/Earth and is proportional to the annihilation cross section. The evaporation term C E N can be neglected for masses above a few GeV. Eq. (11.3) can then be solved and the WIMP annihilation rate Γ A is then The annihilation rate today is then obtained for t = t 4.5 · 10 9 years. The annihilation rates DarkSUSY return are the ones given from the solution above. It is instructive though to consider what happens when t /τ 1. We then have equilibrium between capture and annihilation, i.e. dN dt = 0, and we then have In e.g. MSSM models we typically have equilibrium in the Sun, but not in the Earth. Once we have the annihilation rate, we can calculate the flux of neutrinos from WIMP annihilations in the Sun/Earth. This is simulated with WimpSim [126,127] which uses Pythia [26] to simulate the yields from WIMP annihilations in the Sun/Earth, and then takes care of interactions and oscillations on the way out of the Sun/Earth and propagates the neutrinos to the detector. Once at the detector, interactions are again simulated and we obtains both neutrino-induced lepton volumetric fluxes (particles created per volume element) at the neutrino-nucleon vertex and lepton fluxes (particles passing per area element) at the detector. We also obtain hadronic showers both from charged and neutral current interactions. The neutrino-nucleon interactions are simulated with nusigma [105]. For neutrino oscillations we use the normal ordering best fit points of [128,129]. In the same way as for halo yields, these yields (i.e. how many particles we get at the detector per annihilation in the Sun/Earth) are tabulated for a range of masses and annihilation channels. These simulation results can be accessed via the function dsseyield_sim. The particle physics module is expected to provide a function dsseyield that returns the yield for the current particle physics model. This routine can of course use the simulation results from dsseyield_sim summing over different branching fractions and include cascade decays. For the MSSM model this is for example done including yields from Higgs bosons decaying in flight.
The routine dssenu_rates performs all the steps of the calculation including the calculation of the capture rates, solution of the time evolution equation (11.3) and inclusion of the yields from dsseyield. The end result is then a differential or integrated rate of events. The default in DarkSUSY is to use the full numerically integrated capture rates in dssenu_rates, but tabulated to speed up the calculation. For a more sophisticated likelihood analysis for IceCube data, the external package nulike [35,36] can be used in conjunction with DarkSUSY.

Conclusions
In this article, we have presented a fully revised new version 6 of the numerical package DarkSUSY to compute dark matter (DM) properties and observables. Compared to earlier versions of the code, the main new feature is a manifestly modular and flexible structure. In particular, DarkSUSY is no longer restricted to supersymmetry and neutralino DM, but allows to handle a large variety of DM candidates from particle physics. Each such particle model is contained in a separate DarkSUSY library, or module, and communicates with the (particle-physics independent) core library via a small set of interface functions.
With DarkSUSY 6, one can compute a large variety of very accurate predictions for DM observables, using state-of-the-art techniques. This includes the relic density of thermally produced DM and the associated cutoff scale in the spectrum of matter density perturbations, predictions for DM self-interaction rates, as well as comprehensive direct detection routines that can handle not only the traditional spin-dependent and spin-independent cross sections, but arbitrary effective non-relativistic operators describing DM scattering on nuclei. A large focus of DarkSUSY is furthermore on indirect DM detection, where the code can in principle simultaneously handle a large number of halo objects without re-computing time-consuming quantities, in particular cosmic ray fluxes in positrons, antiprotons or anti-deuterium after propagation. For the same general halo setup, efficient line-of-sight integration routines are set up to ensure accurate predictions for the DM-induced flux in gamma rays and cosmic neutrinos. For some particle models, in particular the MSSM, the spectra from DM annihilation fully include radiative corrections beyond the simplifying, but often adopted, 'model-independent' approach. Sophisticated DM capture rate routines, furthermore, take into account the individual scatterings of DM on the different elements in the Sun and Earth, resulting in highly reliable predictions for the flux of neutrinos from the center of these objects.
In this article we have briefly described the main ingredients and structure of DarkSUSY. In the Appendix, we have provided details on the currently implemented particle moduleswhose number will increase with upcoming releases -and provided various examples of what can be calculated with DarkSUSY. For further details, we encourage the reader to download DarkSUSY, take a look at the manual, and start using the code. We believe that this comprehensive package can be of great use to the physics community, complementary to codes that are similar in scope (like micrOMEGAs [130]). We also note that precision observables computed with DarkSUSY are particularly useful to feed into accurate experimental likelihoods, as provided e.g. by DarkBit [131], to reconstruct or constrain individual signals, or to be used in global scans like in the GAMBIT [18] framework.

A.1 Model parameters
The module generic_wimp provides the simplest example of a particle physics library that the DarkSUSY core library can link to. Rather than being based on an actual particle physics model, it mostly serves to provide an illustration of how the functionalities of DarkSUSY can be used in phenomenological studies of 'vanilla' WIMP DM, when only providing the absolute minimum of input parameters. A generic WIMP model in DarkSUSY is set up by a call to dsgivemodel_generic_wimp, and hence fully defined by the input parameters of that routine: the mass m χ of the DM particle and a flag stating whether the DM particle is its own anti-particle or not; a constant annihilation rate σv, 10 along with the dominant annihilation channel into SM particles; the spin-independent scattering cross section σ SI of DM with nucleons.
While the purpose of this example module is to provide the most simple realization of a generic WIMP, we note that extensions in various directions are straightforward to implement (see Appendix F.2 for technical instructions). How to provide a list of final states, rather than only one dominant annihilation channel, is for example illustrated in the generic_decayingDM module described in Appendix B.

A.2 DM annihilation
The interface function dssigmav0tot simply returns the input model parameter for σvbut only if the stated dominant annihilation channel is kinematically accessible, otherwise it returns zero. In a similar fashion, the module provides the interface function dsanwx which calculates W eff according to Eq. (4.6), if kinematically accessible, from the input value for σv. The module also provides the interface functions dscrsource, dscrsource_line and dsseyield, which link to the tabulated yields provided by the core library. As explained in Section 8, these interface functions are sufficient to allow full access to all indirect detection routines of DarkSUSY.
As an example of a simple result obtained with this particle module, we show in Fig. 3 the annihilation rate σv that is required to obtain a DM relic abundance in accordance with current observations, where we use the results from Planck [1], Ω χ h 2 = 0.1193 ± 0.0028 (2σ). We show separately the cases where DM annihilates only toνν, µ + µ − , W + W − andtt. As introduced above, the generic_wimp particle module does not allow off-shell final state particles. For DM masses below the kinematic threshold, the cross section thus sharply drops to zero, and the thermal average in Eq. (4.3) only picks up contributions from the tail of the thermal distribution. This leads to a rapid increase in the resulting relic density, explaining the almost vertical lines in the figure (for better visibility, we zoom in on the threshold region for W + W − final states). Apart from this effect, the non-constant nature of the 'thermal' annihilation cross section is entirely due to changes in the effective number of relativistic degrees of freedom after decoupling. As a default, DarkSUSY adopts tabulated values from Drees et al [50], but it is straight-forward for the user to choose different options or provide own tables. The program to produce the curves in Fig. 3

is available in examples/aux.
As an illustration of the concept of replaceable functions, we also used the identical program in examples/aux, but at compile time replaced dsanwx with a version of W eff that Figure 3. 'Thermal' annihilation rate σv of a generic WIMP, as a function of its mass m χ , assuming that it only couples to the indicated standard model particles. The width of the band corresponds to the accuracy, at the 2σ level, with which the relic density is determined by Planck [1]. The dashed line shows the results for a model with virtual final state particles (see text for details). For comparison, we also show the results by Drees et al [50].
takes into account the effect of allowing for virtual particles in the final states. The result for W + W − final states is indicated in Fig. 3 as dashed lines. Concretely, we assumed that for m χ m W the constant annihilation cross section of generic_wimp is recovered, (σv) χχ→W + W − = (σv) gen.WIMP , but that close to threshold we still have the phase-space suppression of (σv) χχ→W + W − ∝ λ 1 2 (s, m 2 W , m 2 W ) = s 1 − 4m 2 W /s that is unavoidable in any realistic model. Around the threshold of unstable final states, general expressions for the corresponding s-wave annihilation cross section can be found in Ref. [132]. Assuming W to decay to massless particles, with Γ W ∝ m W and the above scaling of the two-body cross section, in particular, we have where µ W ≡ m 2 W /s and γ ≡ Γ W /m W . For m χ > m W , this ratio is smaller than unity due to the phase-space suppression and the required annihilation rate to obtain the observed relic density is thus larger than in the generic_wimp case. Even for m χ < m W , on the other hand, the ratio does not vanish completely, such that the 'thermal' cross section becomes smaller than in the generic_wimp case for m χ 72 GeV.

A.3 DM scattering
The interface function required for direct DM detection, dsddsigma, returns the usual 27×27 array sigij containing the equivalent cross sectionsσ ij . Here, sigij(1,1) contains the parameter value of σ SI used to initialize the model and the other elements ofσ ij equal zero (note in particular that the spin-dependent nucleon coupling is assumed to vanish in this simple model.). The module also provides a routine dsddgpgn, which translates the input parameter σ SI to the corresponding couplings to nucleons, g a/s , assuming a four-point fermion interaction.
For the kinetic decoupling routines, the interface function dskdm2simp provides the squared amplitude for DM-SM particle scattering. Here, the implementation assumes that DM couples only to one type of SM particles, as in the case of annihilation, and that the same constant amplitude describes both annihilation and scattering.

B.1 Model parameters
Similar in spirit to the generic WIMP case, the module generic_decayingDM provides an illustration of how to use the DarkSUSY core library with a minimal phenomenological framework where DM is not stable but assumed to decay on cosmological time scales. In this case, there is no generic way of thermally producing DM, and we choose to be agnostic about the genesis of DM in the early universe.
Phenomenological studies for such a minimal decaying DM candidate are thus restricted to indirect searches, see below, and fully characterized by the decay rate into standard model particles. More concretely, a generic decaying DM candidate in this module is initialized by a call to dsgivemodel_decayingDM, which takes as input the DM mass m χ , the total decay rate Γ, as well as a list of decay channels (specified by branching ratios and PDG [133] codes of the final state standard model particles).

B.2 Indirect detection routines
The only interface routines provided by the module generic_decayingDM are the source functions needed by the various indirect detection routines, namely dscrsource and dscrsource_line. The former loops over the list of decay channels provided in the model setup and then adds the tabulated yields dN f /dE from the corresponding SM final states f by a call to dsanyield_sim_ls, which is provided as an auxiliary routine in the core library. The quantity returned by dscrsource is then simply S 1 given in Eq. (8.3), i.e.
where E is the energy of the messenger requested by the corresponding indirect detection routine (e.g. photons for gamma rays, or antiprotons for the charged cosmic ray propagation routines). If a final state particle in a two-particle decay channel f is stable, this results in a monochromatic contribution to the corresponding cosmic-ray yield. These cases are caught by a call to dscrsource_line, which returns ΓBR f /m χ (times 2 if there are two identical particles in the final state), along with the energy and width of the 'line' signal (where the width equals the total decay rate of the other final-state particle). Note that for decaying DM no interface function is provided for neutrino yields from the interior of the Sun or the Earth, because the capture rate routines require a non-vanishing DM-SM scattering rate.

C.1 Model framework
Linking the core library to the particle module mssm allows DarkSUSY 6 to access the same particle-physics specific functionality as earlier versions of the code. In particular, the conventions for the superpotential and soft supersymmetry-breaking potential are the same as implemented in [24] (following [134] and similar to [135,136]). The full set of input parameters to be provided at the weak scale thus consists of the pseudoscalar mass (m A ), the ratio of Higgs vacuum expectation values (tan β), the Higgsino (µ) and gaugino (M 1 , M 2 , M 3 ) mass parameters, trilinear couplings (A Eaa , A U aa , A Daa , with a = 1, 2, 3) as well as soft sfermion Eaa , with a = 1, 2, 3). 11 Internally, those values are stored in mssm common blocks. The user may either provide them directly or by setting up pre-defined phenomenological MSSM models with a reduced number of parameters through a call to a routine like dsgive_model or dsgive_model25 (followed by a call to dsmodelsetup). The former sets up the simplest of those models, defined by the input parameters µ, M 2 , m A , tan β, a common scalar mass m 0 , and trilinear parameters A t and A b ; M 1 and M 3 are then calculated by assuming the GUT condition, and the remaining MSSM parameters are given by sets up a pMSSM model with 25 free parameters (see the header of that file for details). Alternatively, all those values can be set by reading in an SLHA file, or providing GUT scale parameters in the case of cMSSM models (via an interface to the ISASUGRA code, as included in ISAJET [137,138]).
Compared to previous versions of the code, DarkSUSY 6 has a new interface to read and write SUSY Les Houches Accord (SLHA) files [139,140]. To use the general structure of DarkSUSY to its fullest, SLHA2 files are preferred, but DarkSUSY will also work with SLHA1 files, but then of course within the limitations of those (no inter-generation sfermion mixing e.g.). DarkSUSY can both read and write SLHA files from/to other codes and uses the SLHA I/O library of Ref. [141] as supplied with FeynHiggs [142][143][144][145][146] to perform this.
All particle and sparticle masses are stored in a common block array mass(). For neutralino masses, we include the leading loop corrections [147][148][149] but neglect the relatively small corrections for charginos [147] (in both cases, masses cannot be negative in our convention). 12 Likewise, all mixing matrices and decay widths are available as common block arrays. The latter are currently only computed for the Higgs particles (via an interface to the FeynHiggs [142][143][144][145][146] package), while the other sparticles have fictitious widths of 0.5% of the sparticle mass (for the sole purpose of regularizing annihilation amplitudes close to poles). Again, the conventions for masses and mixings follow exactly those of Ref. [24], to which we refer for further details.

C.2 Experimental constraints on supersymmetry
There are various experimental constraints on supersymmetric models that are not directly connected to the DM observables that can be computed with DarkSUSY. Furthermore, Dark-SUSY generally focusses on theoretical predictions for such observables, given a DM model realization, rather than on the implementation of experimental likelihoods and the possibility to derive statistically well-defined limits from those. For the latter, we instead refer to packages like DarkBit [131] (or ColliderBit [150] for accelerator-based constraints), which are particularly useful in conjunction with global scans to determine the viable parameter space of a given particle physics model, for example with advanced tools like GAMBIT [18], Still, it is useful to have at least a rough idea of whether a particle model is already clearly excluded before performing advanced calculations to determine DM-related observables. For this purpose, the mssm module allows to perform simple checks on i) the theoretical viability of a given combination of model parameters, ii) a very approximate implementation of current accelerator bounds, as well as to calculate two traditional observables that have turned out to be particularly useful in constraining supersymmetry, namely iii) the anomalous decay b → sγ and iv) the anomalous magnetic moment of the muon (g − 2) µ .
The theoretical viability of a given model is immediately indicated after the initializing call to the subroutine dsmodelsetup(unphys,warning), where a non-zero flag unphys on return indicates a problem of type i), e.g. no electroweak symmetry breaking or the appearance of tachyonic particles (while a non-zero flag warning in the case of the mssm module indicates the breakdown of approximations used in radiative corrections in the Higgs sector). The subroutine dsacbnd allows to check whether any (and which) of the above-mentioned observables ii) -iv) most likely violates current bounds. Compared to previous DarkSUSY versions, we use in particular updated limits from HiggsBounds [151] on the mass of the MSSM Higgs bosons, as well as approximate bounds on squark and gluino masses from LHC 8 TeV data [152]. For b → sγ, we keep our genuine routines for this rare decay (see Ref. [24] for a more detailed description) but now use as a default the result from SuperIso [153]; we compare this to the current limit of B(B → X s γ) = (3.27 ± 0.14) × 10 −4 as adopted in FlavBit [154], based on data from BarBar and Belle [155][156][157]. SuperIso also computes the rate for the rare leptonic decay B 0 s → µ + µ − , which we compare to the LHCb measurement of B(B 0 s → µ + µ − ) = (3.0 ± 0.6 +0.3 −0.2 ) × 10 −9 [158]. Finally, a µ ≡ (g − 2) µ /2 is calculated by dsgm2muon, based on [159]; in dsacbnd, this is compared to the observed value of a µ, obs = (11659208.9 ± 6.3) × 10 −10 [160] after subtracting the SM expectation as specified in PrecisionBit [161].

C.3 Annihilation rates C.3.1 The invariant rate at tree level
In order to allow an efficient numerical integration, all tree-level diagrams contributing to the invariant rate W eff introduced in Eq. (4.4) have been classified according to their topology (s-, t-or u-channel) and to the spin of the involved particles. For initial states with two fermions, all corresponding helicity amplitudes are calculated analytically and then summed to give where θ is the angle between particles k and i. The strength of this method lies in obtaining relatively compact analytic expressions that are summed at the amplitude level and then squared numerically. For initial states containing bosons, the squared amplitude is directly calculated in the standard way. Finally, a numerical integration over cos θ is performed by means of an adaptive method [37]. An important new feature of this DarkSUSY version is that we have parallelized this integration, leading to a spead-up by a factor of up to about six in particular for models with many coannihilations. We have also generalized the annihilation rate routines to that they now allow for a fully general flavour changing structure (i.e. we always include the full 6 × 6 sfermion mass matrices and do not treat them generation by generation.) In DarkSUSY, all coannihilations between neutralinos, charginos and sfermions as calculated in [162] are included. For advanced users, there exists the possibility to manually decide exactly which coannihilation processes to include.

C.3.2 Internal bremsstrahlung
The emission of an additional boson in the final state can greatly enhance neutralino annihilation rates for small DM velocities [163,164]. For photons, this is fully implemented in the mssm module, using analytic expressions for |M| 2 for all processesχ 1 0χ 1 0 → X + X − γ in the limit of vanishing relative velocity of the annihilating neutralinos, where X = q, , H ± , W ± [165]. While this may lead to striking spectral features in the spectrum of gamma rays [166] or positrons [167] relevant for indirect detection, the effect on relic density calculations is very small. In DarkSUSY these contributions are thus not added to the interface function dsanwx, but only to dscrsource. More specifically, only the primary contributions from final state γ and e + are added, because the secondary contributions from the hadronization or decay of the other final states X ± are subdominant. A gauge-invariant procedure subtracts the soft/collinear (final state radiation) part that is already included in the tabulated yields from two-body final states [165]. For performance reasons, these processes are per default only calculated for the most relevant models -e.g. in the case of sfermions almost degenerate in mass with the lightest neutralino -but this behaviour can be steered with a call to a dedicated subroutine dsibset.
For the internal bremsstrahlung of gluons,χ 1 0χ 1 0 →qqg, DarkSUSY uses the above mentioned full analytic expressions for photon IB and rescales the results for quark final states by Q 2 α em → (4/3)α s . In this case, the much larger coupling strength implies that even the integrated annihilation rate can be significantly affected; both gluon IB and QCD loop corrections (which contribute at the same order in α s , see below) are therefore added to dsanwx. While the energy distribution of the final state gluons can be sharply peaked just as in the case of photons, the resulting spectrum of gamma-rays or antiprotons (after gluon and quark fragmentation) is relatively feature-less -though noticeably different to that resulting from two-body quark final states [168]. Rather than simulating these spectra for each MSSM model point with Pythia, DarkSUSY interpolates between pre-tabulated versions of the most extreme 3-body spectra by following the method developed in Ref. [168]; for typical MSSM models, this procedure approximates the true gamma-ray and antiproton spectra to an accuracy of a few percent. These contributions are per default included in the particle source term provided by dscrsource.
The implementation of electroweak internal bremsstrahlung is considerably more involved, not the least due to the shear number of contributing diagrams. For all processes with a fermion pair and either a Higgs or an electroweak gauge boson in the final state, the helicity amplitudes have been determined fully analytically and implemented in DarkSUSY [111,132], extending Eq. (C.1) to three-body final states. Another complication is that for these processes -unlike in the case of the U (1) and SU (3) corrections discussed above -intermediate particles can go on-shell, so a sophisticated subtraction scheme has been implemented both at the cross section and at the yield level to avoid double-counting when including con-  Figure 4. Cosmic ray spectra for an MSSM benchmark model with a degenerate sfermion spectrum, featuring a lightest neutralino with mass m χ ∼ 3.4 TeV and the correct thermal relic abundance (model D2 in [111]). The left panel shows leptonic cosmic ray particles e ± , ν µ , and ν τ (green, blue and cyan line respectively), and the right panel the case of photons (red) and antiprotons (orange). Solid lines indicate the total spectrum, while dashed lines represent the spectrum expected when only taking into account 2-body final states. Dotted lines show the effect of adding photon bremsstrahlung to the 2-body result, so shaded areas highlight the difference due to the inclusion of electroweak corrections.
tributions fromtt or bosonic two-body final states [111]. In DarkSUSY 6 electroweak internal bremsstrahlung is now fully implemented for all cosmic ray yields; while disabled per default due to performance reasons (it takes O(10 s) to compute a full cosmic ray spectrum), these contributions are automatically added to dscrsource after a corresponding call to dsib2set.
We stress that the resulting continuum spectra can differ significantly from two-body final states, or the spectra that result from including 'model-independent' electroweak corrections as implemented in [43,169], and even may lead to striking spectral signatures in positron or neutrino yields. For any detailed phenomenological study we thus strongly recommend to switch on these corrections, despite the additional computation time that this requires. For illustration, we show in Fig. 4 the example of cosmic-ray spectra computed with DarkSUSY 6 , showing separately the impact of electromagnetic and electroweak IB (for a pMSSM-7 benchmark model taken from [111]).

C.3.3 Loop and higher-order corrections
In DarkSUSY, the full analytic one-loop expressions for the processesχ 1 0χ 1 0 → γγ andχ 1 0χ 1 0 → gg [170,171], as well asχ 1 0χ 1 0 → γZ [172], are implemented in the limit of vanishing relative velocity of the annihilating neutralino pair, both as contributions to the particle source function dscrsource (and dscrsource_line) and the invariant annihilation rate dsanwx. The annihilation to a gluon pair, e.g., can be the dominant process in setting the relic density for light squarks [168]. For implementation details, we refer to Ref. [24]. Let us mention that these one-loop expressions formally violate unitarity for diagrams with electroweak gaugeboson exchange because they do not take into account the nonperturbative effects related to multiple electroweak gauge boson exchange ('Sommerfeld effect') as described in [173,174]. In practice, one can still trust the result for sub-TeV neutralinos -but should keep in mind that the implemented cross sections become unphysical in the limit m W /m χ → 0.
For neutralino annihilation to quark final states, we implement the prescription of Ref. [168] to capture the leading effects of QCD corrections in the limit of vanishing relative velocity of the incoming DM particles. This essentially amounts to modelling the incoming neutralino pair as a decaying pseudo-scalar, for which we compute QCD corrections by resumming leading logarithms as in Refs. [175,176], and to add separately the potentially large correction from gluon IB discussed above.

C.4 Neutralino scattering rates C.4.1 Direct detection
Neutralino-nucleus elastic scattering is provided as in the previous version of DarkSUSY apart from small corrections reflecting the fact that the neutralino scattering matrix does not have a pole at mb = m χ + m b [177]. We note that a very similar reasoning also applies to all other quarks. 13 The default option in the mssm module is therefore not to include poles in direct detection amplitudes, but otherwise adopt the prescription by Drees and Nojiri [178]. This can be changed with a call to dsddset_mssm, allowing both to switch to tree-level amplitudes -with or with our without poles -and to adopt the original prescription by Drees and Nojiri (which includes poles except for b and t quark scattering) [178].
Presently, only the traditional spin-independent and spin-dependent cross sections are implemented, i.e. the interface function dsddsigma returns the (partial) equivalent cross sectionsσ ij off a target with all elements but sigij(1,1) and sigij(4,4) set to zero. While it is typically preferable to directly call dsddsigma from a main program, the module also provides a function dsddgpgn that returns the individual couplings to nucleons as given in Eqs. (5.8, 5.9).

C.4.2 Kinetic decoupling
The interface function dskdm2 provides the full amplitude for neutralinos scattering with fermions in the limit of vanishing momentum transfer, as determined in Ref. [53]; dskdm2 implements analytical expressions for the same quantity in the simplified limit where the scattering partners are relativistic, and their energies are well below any kinematically accessible threshold. The subroutine dskdparticles which is usually called directly from dskdtkd, finally, determines all relevant resonances for the scattering processes (stemming from s-channel sleptons).

D The Silveira-Zee module (scalar singlet)
As an example of implementing a nonsupersymmetric model in DarkSUSY 6.0 we present the scalar singlet model of Silveira and Zee [40]. Dark matter particles in this model go under many names, among them: scalar phantoms (which is the original name in [40]) and singlet Higgs dark matter. To avoid confusion in the naming of DarkSUSY modules, we have adopted the unambiguous name silveira_zee for the particle-physics module containing it. 13 For heavy quarks the argument is indeed identical: given that the probability to find a heavy quark in the nucleon is negligible, it has to be co-produced with the squark -which leads to a pole in the unphysical region, at mχ = mq + mq. The same reasoning applies to light quarks that are not valence quarks. For valence quarks, on the other hand, the parton distribution function at the momentum transfer required to produce a squark in the s-channel, Q 2 ∼ m 2 χ , is negligible, and so is the probability of forming a resonance.

D.1 Model parameters
The Silveira-Zee model [40] adds a gauge-singlet real scalar field S to the standard model onto which it imposes a Z 2 symmetry S → −S. Its Lagrangian is where H is the Standard Model Higgs doublet. After electroweak symmetry breaking, the S boson acquires a tree-level mass GeV is the Higgs vacuum expectation value. The module uses the S mass m S and the S-Higgs coupling constant λ as model parameters. These parameters are set with a call to dsgivemodel_silveira_zee, followed as usual by a call to dsmodelsetup to initialize a given model.

D.2 DM annihilation
We have included the following annihilation channels in the annihilation of a pair of S bosons: SS → + − , qq, γγ, W + W − , ZZ, Zγ, gg, HH. All these annihilation channels but SS → HH are mediated exclusively by Higgs exchange in the s-channel. We have computed the invariant annihilation rate and obtained  [25] once it is recognized that their quantity σv rel = W/s (notice that their σv rel = 2σv/(1 + √ 1 + v 2 ) = σv). For Γ H→XY ( √ s) we follow the procedure in [25], i.e. we use tabulated values for √ s < 300 GeV and analytic expressions at higher √ s, but instead of the tables in [179] we have produced our own tables using HDecay 6.51, and we use the analytic expressions in [25,180].
The annihilation channel SS → HH is a sum of Higgs-mediated s-channel, S-mediated t-channel, S-mediated u-channel, and contact SSHH diagrams. We have computed the invariant annihilation rate and find When comparing Eq. (D.5) with the analogous expression in [25], there appears to be a typo in their Eq. (A4): the sign in front of their log term is incorrect. Besides dsanwx, which returns the invariant rate W , the module also provides the annihilation cross section times relative velocity (as it appears e.g. in indirect detection calculations), in terms of interface functions dssigmav and dssigmavpartial. This allows in the usual way access to all indirect detection routines of the core library, via interface functions dscrsource, dscrsource_line and dsseyield.

D.3 DM scattering
For the scattering cross section of S bosons off nuclei, we follow [25] and use the following nonzero S-nucleon coupling constants from Higgs-boson exchange with quarks and gluons, where we take f N = 0.30 as in the Erratum of [25]. These G's are computed in a function dsddgpgn. Based on this, the interface function dsddsigma returns the (partial) equivalent cross sectionsσ ij off a target, as explained in Section 5. Presently, only the traditional spin-independent and spin-dependent cross sections are implemented, i.e. sigij(1,1) and sigij (4,4) are the only non-vanishing components of sigij as returned by dsddsigma in this module. For the kinetic decoupling routines, we need the full scattering amplitude with fermions, averaged over the transferred momentum as stated in Eq. (4.14). This was recently calculated in [56], and we provide it in terms of the interface function dskdm2.
Here, m f is the mass of the fermionic scattering partner and the colour factor is N f = 3 for quarks and N f = 1 for leptons. We note that this is an example where the t = 0 prescription for the scattering amplitude (see the description in Section 4.2) would fail rather badly, given that only a t-channel process contributes. In particular, we find while |M| 2 t=0 = 2N f λ 2 (m f /m H ) 4 in the same limit, which features a different scaling with energy (k ω for relativistic scattering partners). Given that there are no s-channel resonances in the scattering amplitude with fermions, the interface function dskdparticles is trivial for the silveira_zee module.
As an example of a result obtained with DarkSUSY using this specific particle module, we plot in Fig. 5 the value of the coupling λ that is required to obtain the correct relic density, as a function of the DM mass m S . For these values of λ, we also show the mass of the smallest protohalos in this model, given by the cutoff M cut in the power spectrum λ Figure 5. Value of λ that results in Ωh 2 = 0.112, along with the corresponding mass scale of the smallest protohalos, as a function of the Scalar Singlet mass m S . The shaded band reflects the uncertainty in M cut due to quark scattering: the lower curve follows [58] and assumes that all quarks contribute as free particles for T > T QCD ≡ 154 MeV, with no DM-quark scattering for T < T QCD , while the upper curve follows the conservative treatment of [53] and assumes that only light quarks contribute at T > 4T QCD .
of matter density perturbations. We note that close to the resonance, m S ∼ m H /2, the assumption of kinetic equilibrium during chemical freeze-out (which the standard treatment [48,49] discussed in Section 4.1 relies on) is not satisfied; as a result, the actual value of λ resulting in the correct relic density differs by a factor of up to about 2 from what is shown in this figure [56]. Extended relic density routines to capture this effect will be available with a future DarkSUSY release.

E A dark sector module with velocity-dependent self-interactions
With the vdSIDM module, we provide an example for a class of models confined to a dark sector with no or only little interaction with the SM. These simplified models contain a CDM particle, a relativistic dark radiation (DR) particle that is thermally distributed, and one additional -typically light -particle that mediates a renormalizable interaction between them. Such models have been extensively discussed in the literature as they affect structure formation in a way that leads to interesting cosmological observables for DM particles that are complementary to the traditional direct, indirect or collider searches for DM. In fact, it has been argued that such simple setups may be sufficient to simultaneously mitigate (combinations of) the main known observational discrepancies with the ΛCDM concordance model [52,55,79,80,[181][182][183].

E.1 Model parameters
Similar in spirit to the simplified model analysis presented in Ref. [54], the module provides a structure that allows any of the three particles just mentioned to be either a scalar, a fermion or a vector, with arbitrary masses and couplings. In the present version 6.1 of the code, the module has two concrete example models fully implemented, where DM is a massive Dirac fermion ψ χ , DR is a massless fermion ψγ, and the mediator is either a vector V or a scalar φ. The interaction parts of the respective Lagrangians are thus given by and DM mass, mediator mass and couplings are set by calls to dsgivemodel_vdSIDM_vector or dsgivemodel_vdSIDM_scalar, followed as usual by a call to dsmodelsetup. In these models, the mediator thus only decays into invisible particles. Additional decay channels into SM particles are in principle straightforward to add, in a similar fashion as done in the generic_decayingDM module, and will be included in a later version of the code. We caution, however, that such decay channels are strongly constrained from direct detection experiments (for the case of scalar mediators [184]) or cosmology (for the case of vector mediators [185]).
A dark radiation component with temperature equal to that of the photons would not be cosmologically viable, as the new DS particles would result in an unacceptably large additional contribution to the radiation density. This is conventionally stated in terms of where ρ 1ν is the energy density contributed by one massless neutrino species, and T ν is the neutrino temperature (which differs from the photon temperature after e ± annihilation). In the vdSIDM module, we provide a function dsrdxi which assumes that the dark sector has been in thermal contact with the SM heat bath at very early times, but then decoupled at a temperature when all SM particles were still relativistic. This function thus automatically replaces the trivial version of dsrdxi provided by the core library. In the expression in footnote 5, we thus have g SM * (T dc ) = 106.75 and g DS * (T dc ) = 10 (8) for a vector (scalar) mediator. Aggressively assuming that the mediators are still relativistic during BBN, this results in for the effective number of relativistic degrees of freedom, which in both cases respects current limits at 2-3σ [186,187]. At late times, when only the DR particles are present, we find ∆N eff = 0.38 (0.28), largely compatible with the current CMB bound of ∆N eff < 0.35 (95% C.L.) from the Planck collaboration [1]. For convenience, the module provides a function dsrddeltaneff to compute ∆N eff as a function of (photon) temperature.

E.2 DM annihilation
At tree-level, DM annihilates via the t-and u-channel to a pair of mediators, which in the early universe are kept in thermal equilibrium with the DR particles. We have calculated and implemented the full expressions for the invariant rate, and only quote here the result in the limit of small CMS energies √ s → 2m χ : where α χ ≡ g χ /(4π) and p CM is the initial CMS momentum of (each of) the DM particles. Just as in Eq. (D.9), the invariant rate can be translated to the annihilation cross section times relative velocity (which is returned by dssigmav). We note that the resulting expressions, in the above limit, differ from the corresponding lowest-order results stated in Eq. (28) of Ref. [95]. Furthermore, DM can annihilate to a pair of DR particles, by a mediator exchange in the s-channel. For this we find the following contributions to the invariant rate: where D V,φ are defined in analogy to Eq. (D.4) as the inverse of the denominator of the respective propagator. For the total widths that enter in D V,φ we add to the decay rate of the mediator into dark radiation the partial width of the mediator decaying into DM particles (if this is kinematically allowed). If the mediator is much lighter than the DM particle, the annihilation rates for small relative velocities are strongly enhanced by the Sommerfeld effect [188,189]. The interface function dsanwx therefore adds to the full tree-level expressions an enhancement factor dsansommerfeld to multiply the leading-order expressions in the v → 0 limit. This factor depends on whether the process is s-wave dominated (as in the vector mediator case) or p-wave dominated (as in the scalar mediator case), and implements in the present version the analytic expressions from Ref. [95,190], which result from approximating the Yukawa potential with a Hulthén potential.
Since the mediator in this simplest realization is assumed to exclusively decay into invisible DR particles, all source functions for indirect detection routines in the vdSIDM module currently return zero.

E.3 DM scattering
Since the mediators do not couple to SM particles, the DM scattering cross section off nuclei is zero in this module. There is, however, significant scattering of the DM particles with the DR particles, through t-channel exchange of mediator particles, which can lead to kinetic decoupling much later than in the case of ordinary WIMPs. For the momentum-averaged scattering amplitude returned by the interface function dskdm2, we implement the expressions from Ref. [54]: where ω is the energy of the DR particle. We also add the amplitudes for DM scattering directly off heat bath vector mediators, |M| 2 t = 64g 4 χ /3, and scalar mediators, |M| 2 t = 16g 4 χ /3 [54], though these are typically subdominant for the parameter ranges of interest. The fact that the mediators generate a Yukawa potential does not only lead to the already mentioned Sommerfeld effect enhancing the annihilation of DM particles, but also implies strong and velocity-dependent DM self-interactions if the mediators are light (which is the reason for the 'vd' in the module name). The interface function dssisigtm returns the resulting momentum transfer cross section per unit mass, σ T /m χ , as introduced in Section 7, by determining whether scattering takes place in the classical, Born or resonant regime (following Ref. [95]). For a vector mediator, it averages between the corresponding scattering rates for a repulsive and an attractive potential to properly take into account both components of the Dirac DM particles χ andχ (for asymmetric DM, only the repulsive part would contribute). For a scalar mediator, an attractive potential is implemented.
In Fig. 6 we provide a simple example of a DarkSUSY result obtained using the vdSIDM module. For each point in the DM mass vs. mediator mass plane, the coupling constant g χ is determined by requiring that thermal production leads to a relic density of (Ω χ + Ωχ)h 2 = 0.112, thus matching the observed DM density (fully including the Sommerfeld enhancement). The blue band shows σ T /m χ , averaged over a Maxwellian distribution with a most probable velocity of v 0 = 30 km/s. This corresponds to typical velocities in dwarf spheroidal galaxies, where it has been argued that a value of σ T /m χ ∼ 1 cm 2 /g can mitigate both the cuspcore problem, the too-big-to-fail problem and the diversity problem (see section 7). One can clearly see the resonant regime of the self-interaction cross section -the main reason for the different structure being that in the vector mediator case we average between attractive and repulsive potentials (the latter of which does not result in resonances). The green band, finally, indicates the cutoff mass M cut resulting from late kinetic decoupling for a range that may be relevant in mitigating the missing satellite problem (while too large values are excluded by Lyman-α data, see Refs. [54,182] for a recent discussion). For the purpose of this figure, we set gγ = g χ , noting that M cut ∝ g 3/2 γ m .3 γ for a given DM mass m χ m φ,V [54]. As explained above, the temperature of the dark sector is not a constant in this model; during DM freezeout, we have ξ ≈ 0.54 (0.57) for vector (scalar) mediators, but ξ ≈ 0.47 (0.44) during late kinetic decoupling.
We note that the relic density calculation that enters in fixing the coupling constant g χ assumes that the combinations of the model parameters (g χ , m χ , m V,φ ) do not result in a Sommerfeld enhancement very close to a resonance. In that case, DM would re-enter chemical equilibrium after kinetic decoupling, with a subsequent second freeze-out period [191][192][193][194]. The necessary coupled set of Boltzmann equations to accurately describe such a scenario [194] will be available with a future release of DarkSUSY. While requiring some fine-tuning in the model parameters, such a scenario would provide additional phenomenological interest as it may not only address the common ΛCDM small-scale issues, but additionally reconcile the observed tension in σ 8 and H 0 measurements between high-and low-redshift cosmological observables [183,195].

F.1 Getting started
To get started, first download DarkSUSY from www.darksusy.org and unpack the tar file. To compile, run the following in the folder where you unpacked it ./configure make You will then have compiled the main DarkSUSY library, as well as all supplied particle physics modules. To test whether the installation was successful, type cd examples/test ./dstest The program will take up to about a minute to run and reports if there are any problems. 14 Even if you now have DarkSUSY running, it is more fun to start doing some calculations on your own, so the next steps might be to • Look at the code examples/dsmain_wimp.F which contains an example main program to calculate various observables for WIMP DM candidates. It can be the starting point if you want to write your own programs (simply copy both the program and the makefile in examples to a user-defined directory, and run make again in that new directory) • Use another particle physics module, e.g. the generic WIMP or the scalar singlet (Silveira-Zee) model. To do this, just change your makefile to select the particle physics module you wish to use. Again, a good starting point to test this is the example program examples/dsmain_wimp.F: In this case, you can run the same program for different particle modules simply by changing the first line in examples/dsmain_wimp.driver. Alternatively, you can call make -B DS_MODULE=<your_module_choice> to override the entry in the driver file.
• If you want to modify some existing DarkSUSY function or subroutine, please don't do it! Instead, add your own routine as a replaceable function, by running the script scr/make_replaceable.pl on the routine you wish to have a user-replaceable version of (just run the script without arguments for more details). You will then find a dummy version of the routine in the corresponding user_replaceables folder, where you can edit it to your liking. 15 Happy running! F.2 How to add a new particle physics module To create a new particle physics module, the easiest way is to start from an existing one as a template and create a new one from that one. To help you in this process we provide a script scr/make_module.pl that takes two arguments, the module you want to start from and the new one you want to create (for further instructions, just call the script without arguments). It will then copy the module to a new one, change its name throughout the module and make sure that it is compiled by the makefiles and included properly when requested by the main programs. If you specify the option -i only interface functions will be copied (which creates a cleaner starting point, but also will most likely not compile without modifications). When creating a new module this way, the best is to copy from a module that is as similar as possible to your new model. If you want a clean setup, you can always copy from the empty module. A general advice is to view the modules we provide as a starting point as inspiration for your new modules. Even though a particle physics module does not need to include all interface functions (which ones are needed only depends on the observables you try to calculate in your main programs), it needs to provide an initialization routine dsinit_module.f. This routine should set a global variable moduletag to the name of the module so that routines that need to check if the correct module is loaded can do so. When using the script scr/make_module.pl this routine is always created and moduletag set as it should.

F.3 How to choose pre-defined halo profiles and add new ones
As described in more detail in Section 6, the halo profiles are handled by the dsdmsdriver function in src/dmd_mod. This driver function is like a container, or an interface, for all presently active halo profiles. It ensures that all parts of the code use the halo parameterization in a consistent and optimized way (e.g. routines using spherical vs. axial symmetry for improved performance, storage of computationally intensive quantities, etc.). Given the complexity of this function, we provide various example main program in examples/aux to explicitly demonstrate possible use cases. The example DMhalo_predef.f illustrates how to use the default version of dsdmsdriver (provided with the DarkSUSY 6 release) to load additional profiles into the currently active halo database by using its pre-defined halo parameterizations, while the example DMhalo_table.f shows how to load a profile from a table. In DMhalo_new.f, we demonstrate instead how to correctly extend dsdmsdriver when adding a new profile parameterization in order to consistently make it available to all DarkSUSY routines that rely on the DM density (in this concrete example, we add the spherical Zhao profile [196], aka αβγ profile). Given the complexity of the previous example, finally, the program DMhalo_bypass.f demonstrates a work-around of completely bypassing the default dsdmsdriver setup when switching to a user-provided new DM density profile; in this approach, the advanced DarkSUSY system of automatic tabulation of quantities related to DM rates cannot be easily exploited.

F.4 Main technical changes compared to previous DarkSUSY releases
For those familiar with DarkSUSY 4 and 5, we list here some of the most important technical changes introduced in DarkSUSY 6. Most physics differences/ improvements are described earlier in this document, with a summary of highlights provided in Section 2, and for the general changes and the overall structure see section 3.
• The way DM halos are initialized and used has fundamentally changed. See Sections 6 and Appendix F.3.
• The separation between the particle physics model independent ds_core and the particle physics modules. This has the effect that e.g. the routines that calculate yields of particles from WIMP annihilation, the former dshayield and dswayield, are now split into two parts. The most important routines are the interface functions that reside in the particle physics module (dscrsource[_line] and dsseyield). Those can call auxiliary routines such as dsanyield_sim in ds_core to get results from the simulated yield tables. Cascade decays (e.g. decaying Higgs bosons) are handled in the particle physics module as this is particle physics dependent.
• The yield functions have changed name from dshayield to dsan_yield for annihilation in the halo and from dswayield to dsse_yield for annihilation in the Sun/Earth.
• In general, many routines have changed names to clarify what they do and which belong together.
• The standard model resides in the particle physics library. The reason for this is that one most often wants to implement a particle physics model and the standard model at the same time, and separating it out to ds_core would not be very practical. At the same time, the standard model is obviously not a BSM theory in itself, with a viable DM candidate. Particle definitions and basic properties have therefore, for convenience, simply been collected in the directory src_models/common/ and can be included by any (BSM) particle physics module if the user so wishes.
• The particle codes (k-variables) are now treated as internal codes. They can be used by the particle physics module if the module so wishes, but the interface functions and routines in ds_core instead use PDG [133] codes when referring to particles.
• A more stream-lined make system where most makefiles and the configuration script can be updated via a simple script. This makes it easier if the user wants to add particle physics modules or other contributions to the code.
• A new test routine dstest that tests the code and compares with expected results for a set of given observables.
• dsmain.F is now particle physics module independent in the sense that it can be linked to different particle physics modules. This assumes that the module has implemented the observables one is asking for -which is why dsmain.F presently comes in two versions (dsmain_wimp.F and dsmain_decay.F).

F.5 Interfaces to other codes
For some of its calculations, DarkSUSY uses external codes. We here briefly mention which codes we use and how they are interfaced: • FeynHiggs [142] is used for Higgs boson mass calculations in the MSSM model. We interface with FeynHiggs by calling their setup routines directly.
• HiggsBounds [151] and HiggsSignals [197] are used for LEP and LHC limits and constraints on the observed Higgs boson. We interface with both these codes by piggybacking on the FeynHiggs interface following the example programs in HiggsBounds and HiggsSignals.
• SuperIso [153] is used for rare decays and is interfaced via our SLHA interface (i.e. by writing an SLHA file to disk and read this in from SuperIso).
• IsaSugra is used for RGE running in CMSSM models and is interfaced by manually extracting all the low-energy parameters from the [137,138] common blocks.