\documentclass[11pt]{cernrep}
\usepackage{graphicx,epsfig,amsmath}
\bibliographystyle{lesHouches}
\begin{document}

\title{DLHA: Dark Matter Les Houches Agreement}

\author{{\rm Coordinator}: 
Csaba Bal\'azs$^{1, 2, 3}$,
David G. Cerde\~no$^{4, 5}$,
Rebecca Leane$^{1, 3}$,
Mitsuru Kakizaki$^{6}$,
Sabine Kraml$^{7}$,
Christopher Savage$^{8}$,
Pat Scott$^{9}$,
Sezen Sekmen$^{10}$,
\\{\bf AMIDAS}:
Chung-Lin Shan$^{11}$,
\\{\bf DarkSUSY}: 
Lars Bergstr\"om$^{8}$,
Torsten Bringmann$^{12}$,
Gintaras Duda$^{13}$,
Joakim Edsj\"o$^{8}$,
Paolo Gondolo$^{8, 14}$,
\\{\bf DMFIT and DMMW}:
Stefano Profumo$^{15}$,
Tesla Jeltema$^{15}$,
\\{\bf DRAGON}: 
Carmelo Evoli$^{16}$,
Daniele Gaggero$^{17}$,
Dario Grasso$^{17}$,
Luca Maccione$^{18}$,
\\{\bf FeynRules}: 
Claude Duhr$^{19}$,
Benjamin Fuks$^{20}$,
\\{\bf GALPROP}:  
Igor Moskalenko$^{21}$,
Elena Orlando$^{21}$,
Troy Porter$^{21}$,
Andrew Strong$^{22}$,
\\{\bf ISAJET}:
Howie Baer$^{23}$,
\\{\bf micrOmegas}:  
Genevieve B\'elanger$^{24}$,
Fawzi Boudjema$^{24}$,
Alexander Pukhov$^{25}$,
\\{\bf PPPC4DMID}: 
Marco Cirelli$^{26, 27}$,
Mario Kadastik$^{28}$,
Paolo Panci$^{29}$,
Martti Raidal$^{28}$,
\\{\bf SLHA}:  
Sven Kreiss$^{30}$,
Frank S. Thomas$^{31}$,
\\{\bf Semi-analytical propagation models}:
Timur Delahaye$^{4}$,
Fiorenza Donato$^{32}$,
Julien Lavalle$^{33}$,
David Maurin$^{34}$,
Pierre Salati$^{24}$,
Richard Taillet$^{29}$,
\\{\bf SuperISO}:
Alexandre Arbey$^{26, 35, 36}$,
Farvah Mahmoudi$^{26, 37}$
}
 
 \institute{
  $^{1}$ARC Centre of Excellence for Particle Physics at the Tera-scale, Monash University, Victoria 3800 Australia
\\$^{2}$Monash Centre for Astrophysics, Monash University, Victoria 3800 Australia
\\$^{3}$School of Physics, Monash University, Melbourne, Victoria 3800 Australia
\\$^{4}$Instituto de F\'isica Te\'orica UAM/CSIC, Universidad Aut\'onoma de Madrid, 28049 Madrid, Spain
\\$^{5}$Departamento de F\'isica Te\'orica, Universidad Aut\'onoma de Madrid, 28049 Madrid Spain
\\$^{6}$Department of Physics, University of Toyama, 3190 Gofuku, Toyama 930-8555, Japan
\\$^{7}$Laboratoire de Physique Subatomique et de Cosmologie, UJF Genoble 1, CNRS/IN2P3, INPG, F-38026 Grenoble, France
\\$^{8}$The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden
\\$^{9}$Department of Physics, McGill University, 3600 rue University, Montreal QC H3A 2T8 Canada
\\$^{10}$Department of Physics, Florida State University, Tallahassee, FL 32306 USA
\\$^{11}$Institute of Physics, Academia Sinica, Taipei 11529, Taiwan, R.O.C
\\$^{12}$II. Institut f\"ur Theoretische Physik, Hamburg University, Luruper Chaussee 149, 22761 Hamburg Germany
\\$^{13}$Department of Physics, Creighton University, 2500 California Plaza, Omaha, NE 68178 USA
\\$^{14}$Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112 USA
\\$^{15}$Department of Physics and Santa Cruz Institute for Particle Physics, University of California, 1156 High St., Santa Cruz, CA 95064, USA
\\$^{16}$II. Institut f\"ur Theoretische Physik, Hamburg University, Luruper Chaussee 149, 22761 Hamburg, Germany
\\$^{17}$INFN sezione di Pisa, Largo Bruno Pontecorvo 3, Pisa Italy
\\$^{18}$Arnold Sommerfeld Center, Ludwig-Maximilians-Universitaet and Max-Planck-Institut fuer Physik, Theresienstr. 37, 80333 Munich Germany
\\$^{19}$Institut f\"ur theoretische Physik, ETH Z\"urich, Wolfgang-Paulistr. 27, CH-8093, Switzerland
\\$^{20}$Institut Pluridisciplinaire Hubert Curien, D\'epartement Recherches Subatomiques, Universit\'e de Strasbourg, CNRS-IN2P3, 23 Rue du Loess, F-67037 Strasbourg, France
\\$^{21}$Hansen Experimental Physics Laboratory and Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Stanford CA 94309 USA
\\$^{22}$Max Planck Institut f\"ur extraterrestrische Physik, Postfach 1312, 85741 Garching, Germany
\\$^{23}$Department of Physics and Astronomy, University of Oklahoma, Norman, OK 73019 USA
\\$^{24}$LAPTH, Universit\'e de Savoie, CNRS, BP110, F-74941 Annecy-le-Vieux Cedex, France
\\$^{25}$Skobeltsyn Institute of Nuclear Physics, Moscow State University, Moscow 119992 Russia
\\$^{26}$CERN Theory Division, CH-1211 Gen\`eve, Switzerland
\\$^{27}$Institut de Physique Th\'eorique, CNRS, URA 2306 \& CEA/Saclay, F-91191 Gif-sur-Yvette, France
\\$^{28}$National Institute of Chemical Physics and Biophysics, Ravala 10, Tallinn 10143, Estonia
\\$^{29}$CP3-origins and the Danish Institute for Advanced Study DIAS, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark
\\$^{30}$New York University, Department of Physics, 4 Washington Place, New York NY 10003 USA
\\$^{31}$Institut f\"ur Theoretische Physik und Astrophysik, Universit\"at W\"urzburg, Am Hubland, 97074 W\"urzburg, Germany
\\$^{32}$Dipartimento di Fisica Teorica, Universit\`a di Torino \& INFN, 10122 Torino Italy
\\$^{33}$Laboratoire Univers \& Particules de Montpellier [LUPM] CNRS-IN2P3, Universit\'e Montpellier II [UMR-5299] Place Eug\`ene Bataillon - Bat.13 CC.72 F-34095 Montpellier Cedex 05, France
\\$^{34}$Laboratoire de Physique Subatomique et de Cosmologie, Universit\'e Joseph Fourier Grenoble 1, CNRS/IN2P3, Institut Polytechnique de Grenoble, 53 avenue des Martyrs, Grenoble, 38026, France
\\$^{35}$Centre de Recherche Astrophysique de Lyon, Observatoire de Lyon, Saint-Genis Laval Cedex, F-69561, France; CNRS, UMR 5574; Ecole Normale Sup\'erieure de Lyon, Lyon, France
\\$^{36}$Universit\'e de Lyon, France; Universit\'e Lyon 1, CNRS/IN2P3, UMR5822 IPNL, F-69622~Villeurbanne Cedex, France
\\$^{37}$Clermont Universit\'e, Universit\'e Blaise Pascal, CNRS/IN2P3, LPC, BP 10448, 63000 Clermont-Ferrand, France
}

\maketitle


\begin{abstract}
This work presents a set of conventions and numerical structures that aim to provide a universal interface between computer programs calculating dark matter related observables.  It specifies input and output parameters for the calculation of observables such as abundance, direct and various indirect detection rates.  These parameters range from cosmological to astrophysical to nuclear observables.  The present conventions lay the foundations for defining a future Les Houches Dark Matter Accord. 
\end{abstract}


\section{INTRODUCTION}

% DM related calculations

Over the last decade a burst of activity has surrounded various dark matter related problems.  As a byproduct of this activity, a number of robust numerical tools have been created by particle, astro-particle and astrophysicists.  These tools predict values for observables that, when contrasted with observation, can shed light on the properties of dark matter.  The dark matter related computer codes address a wide range of physical problems, of which the most important are:
\begin{enumerate}
\item deriving Feynman rules from Lagrangians,
\item calculating scattering amplitudes from Feynman rules,
\item building cross sections or decay rates from amplitudes,
\item computing dark matter abundance from cross sections or decay rates,
\item obtaining dark matter-nucleon scattering rates from cross sections,
\item evaluating cosmic ray yields based on annihilation cross sections or decay rates.
\end{enumerate}

% DM related calculators

A certainly incomplete but representative list of such computer codes is:
\begin{center}
\begin{tabular}[h]{lll}
 \cite{Pukhov:2004ca}           CalcHEP        (2, 3)          ~~~ & 
 \cite{Gondolo:2005we}          DarkSUSY       (3, 4, 5, 6)    ~~~ &
 \cite{Jeltema:2008hf,DMMWRef}  DMFIT and DMMW (6)             ~~~ \\  
 \cite{DRAGONRef}               DRAGON         (6)             ~~~ & 
 \cite{Christensen:2008py}      FeynRules      (1)             ~~~ & 
 \cite{Strong:2009xj}           GALPROP        (6)             ~~~ \\
 \cite{Paige:2003mg}            ISATools       (4, 5)          ~~~ & 
 \cite{Semenov:2008jy}          LanHEP         (1)             ~~~ & 
 \cite{Belanger:2010pz}         micrOmegas     (2, 3, 4, 5, 6) ~~~ \\
 \cite{PPPC4DMID}               PPPC4DMID      (6)             ~~~ &  
 \cite{Arbey:2011zz}            SuperIso Relic (4)             
\end{tabular}
\end{center}
where the numbers in parentheses after the names indicate what these codes do in the context of the task-list above.  These two lists show that 
% \renewcommand{\labelitemi}{$-$}
% \begin{itemize}
% \item 
today there exists no code that covers all the dark matter related calculations (although some come very close), and
% \item 
there are significant overlaps between the capabilities of these codes.
% \end{itemize}
% \renewcommand{\labelitemi}{$\bullet$}
%
What these lists do not reveal is how general and sophisticated the various programs are.  Not surprisingly, each of these tools has its strengths and weaknesses.  Necessarily, all these codes contain hard-wired assumptions that potentially limit their capabilities. 

% LHAs

When particle physicists found themselves in a similar situation regarding collider related calculations, they created a series of interfaces that allowed their various tools to communicate with each other \cite{Boos:2001cv, Whalley:2005nh, Skands:2003cj, Allanach:2008qq, Mahmoudi:2010iz}.  These interfaces allow the users of the codes to easily and selectively exploit the features that each of the tools offers.  The success of the Supersymmetric Les Houches Accord (SLHA) and the other accords lies in the fact that they significantly increase ease and flexibility for collider related calculations \cite{Skands:2003cj, Allanach:2008qq}. 

% DLHA interface

Following the successful path of previous Les Houches accords, the Dark Matter Les Houches agreement (DLHA) proposes a format for storing and exchanging information relevant to dark matter calculations.  In the spirit of the Les Houches accords, DLHA aims to interface various calculators to provide increased flexibility to users of these tools.  Presently these calculators work semi-independently from each other and their consecutive use typically requires tedious interfacing.  However careful this interfacing is, it may jeopardise the integrity of the tools and, in turn, the results of the calculation.  

% DLHA transparency

A related problem DLHA targets is the transparency of the dark matter calculations.  Present codes input and output limited amounts of information while they may contain large amounts of implicit, hard-wired assumptions that affect their results.  The more of this information is accessible, the more control a user has over the calculation.  Additionally, making more assumptions explicit gives the user a chance to change them and this could lead to a more diverse and productive use of these codes.  Since a DLHA file can be part of the input or the output of a numerical code calculating dark matter related observables, it aids in making implicit assumptions explicit and, potentially, it allows for changing some of those assumptions.  

% Why do we need DLHA?

Here we give some examples of how such an agreement enables the user to easily interface various codes with different capabilities and thereby calculate quantities that none of these codes could calculate alone.
\begin{itemize}
\item Assume that code $A$ has the capability to calculate decay rates of a dark matter candidate in an exotic particle model, but code $A$ can only handle standard cosmology when calculating relic abundance of the dark matter.  Assume that code $B$ does not implement the exotic particle model of interest but can calculate relic abundance for non-standard cosmologies.  With the appropriate DLHA interface the decay rate can be communicated from code $A$ to $B$ and relic abundance of an exotic candidate can be calculated in a non-standard cosmology.
\item Assume that code $C$ has the capability to calculate the annihilation cross section of a dark matter candidate with next-to-leading order (NLO) corrections but it has no indirect detection routines built into it.  Meanwhile code $D$ can only do one thing, but it does it excellently: calculate indirect detection rates.  An obvious task for DLHA is to pass the annihilation cross section from $C$ to $D$, thereby making it possible to calculate a cosmic ray yield from annihilating dark matter with NLO corrections.
\end{itemize}
These examples are just a limited sample of the possibilities that DLHA offers for a resourceful user.


\section{CALCULATION OF DARK MATTER RELATED OBSERVABLES}

This section summarizes the relevant details of dark matter related calculations so that we can fix our notation and define the accord in the next section.


\subsection{Relic abundance}

Numerous microscopic models have been proposed to describe the identity of dark matter particles.  Depending on the details of these models, and conditions in the early Universe, dark matter may have been produced in various different ways.  For simplicity, here we only address the thermal production mechanism of weakly interacting massive particles.


\subsubsection{Thermal relic abundance with standard cosmological assumptions}

In this section we recapture some details of the calculation of thermally produced relic abundance of dark matter following standard cosmological assumptions.  According to Big Bang cosmology the Universe cooled from temperatures substantially higher than the mass of a typical dark matter particle.  At those temperatures dark matter particles were in chemical equilibrium with their environment and consequently their (co-moving) number density was unchanged.  At temperatures comparable to the typical weakly interacting dark matter particle, at $T = {\cal O}$(100 GeV/$k_B$), the Universe was radiation dominated.  (Here $k_B$ is the Boltzmann constant.)  As the Universe cooled and expanded further, dark matter particles came in thermal equilibrium with the cooling radiation and consequently their number density decreased exponentially.  Later scattering between dark matter particles became rare until the scattering rate fell below the expansion rate of the Universe for the dark matter particles to stay in equilibrium.  At that stage dark matter particles froze-out, that is their (co-moving) number density ceased to change. 

The change of the average number density of dark matter particles, $n_\chi$, is described by Boltzmann's equation \cite{Goldberg:1983nd,Krauss:1983ik,Srednicki:1988ce,Gondolo:1990dk,Edsjo:1997bg}. This change occurs due to (co-)annihilation to Standard Model particles and to the expansion of the early Universe. Boltzmann's equation reads:
\begin{equation}
\frac{dn_\chi}{dt} = -\langle\sigma v\rangle\left(n_{\chi}^2-n_{\chi,eq} ^2\right) - 3 H n_{\chi} .
\label{DLHA_boltzbest}
\end{equation}
On the right hand side the thermally averaged product of the annihilation cross section and relative velocity of the two annihilating dark matter particles is defined by
\begin{equation}
\langle\sigma v\rangle = \frac{\int{\sigma v\ dn_{1,eq}dn_{2,eq}}}{\int{dn_{1,eq} dn_{2,eq}}}
= \frac{\int{\sigma v \ f(E_1)f(E_2)d^3p_1 d^3p_2}}{\int{f(E_1)f(E_2)d^3p_1 d^3p_2}} .
\label{DLHA_sigmav}
\end{equation}
Here $p_i$ are the three-momenta and $E_i = \sqrt{m_i^2 + p_i^2}$ are the energies of the colliding particles with $i=1,2$.%
\footnote{Throughout this section we use natural units with $c = \hbar = 1$, unless stated otherwise.}
%\footnote{Throughout this section we use natural units with $G_N = c = \hbar = k_B = 1$, unless stated otherwise.}
The total (co-)annihilation cross section $\sigma$ depends on the microscopic properties of the dark matter particles and can be calculated assuming a particular dark matter model.  The equilibrium number density is given by thermodynamics
\begin{equation}
n_{\chi,eq} = \frac{g_\chi}{(2\pi)^3}\int{f(E)d^3E} ,
\end{equation}
where $g_\chi$ is the number of internal degrees of freedom of dark matter particles.  Their energy distribution is the familiar
\begin{equation}
f(E)=\frac{1}{\exp(E/k_B T)\pm1} .
\label{DLHA_eq:f(E)}
\end{equation}
In Eq. (\ref{DLHA_boltzbest}) the Hubble parameter describes the expansion rate, and in a nearly flat Friedmann-Lemaitre-Robertson-Walker (FLRW) Universe it is given by
\begin{equation}
H^2 = \frac{8\pi G_N}{3} \rho .
\label{DLHA_hubs}
\end{equation}
As the Universe was radiation dominated, the total energy density was the standard function of the temperature of the radiation $T$:
\begin{equation}
\rho = \rho_{rad}(T) = g_{eff}(T)\frac{\pi^2}{30}T^4 ,
\label{DLHA_energydens}
\end{equation}
where $g_{eff}(T)$ is the number of internal degrees of freedom of the (effectively) massless particles in the thermal bath.\footnote{Here we assume that the various species contributing to the total energy density are all in thermal equilibrium with each other at a common temperature $T$.} 

The decrease in number density due to expansion of the Universe can be made implicit by rewriting Eq. (\ref{DLHA_boltzbest}) in terms of the co-moving number density 
\begin{equation}
Y = \frac{n_{\chi}}{s} .
\label{DLHA_y}
\end{equation}
For the radiation dominated Universe the total entropy density is
\begin{equation}
s = s_{rad}(T) = h_{eff}(T)\frac{2\pi^2}{45}T^3 .
\label{DLHA_entropy}
\end{equation}
The number of internal degrees of freedom of the particles contributing to the entropy is $h_{eff}(T)$.\footnote{For a relativistic particle with one internal degree of freedom, such as a spin zero boson, $g_{eff} = h_{eff} = 1$.}
Recasting Boltzmann's equation for $Y$ as a function of $x=m_{\chi}/T$, where $m_\chi$ is the mass of the dark matter particle, we arrive at
\begin{equation}
\frac{dY}{dx} = \frac{1}{3H}\frac{ds}{dx}\langle\sigma v\rangle(Y^2-Y^2_{eq}) .
\label{DLHA_boltznew}
\end{equation}
By utilising Eqs. (\ref{DLHA_hubs}), (\ref{DLHA_energydens}) and (\ref{DLHA_entropy}) this may also be written as
\begin{equation}
 \frac{dY}{dx} = 
 -\sqrt{\frac{\pi g_\star}{45 G_N}}\frac{m_\chi}{x^2}\langle\sigma v\rangle(Y^2-Y^2_{eq}) .
\label{DLHA_boltznew2}
\end{equation}
The parameter $g_\star$ incorporates the degrees of freedom arising from the energy and entropy densities
\begin{equation}
 g_\star ^{1/2} = \frac{h_{eff}}{g_{eff}^{1/2}}\left(1+\frac{1}{3}\frac{T}{h_{eff}}\frac{dh_{eff}}{dT}\right) .
\label{DLHA_eq:gstar}
\end{equation}
A value for $Y_{eq}$ of a cold gas of dark matter is obtained for (non-)relativistic dark matter particles by taking the (low) high energy limit of the Maxwell-Boltzmann distribution.  For non-relativistic dark matter particles
\begin{equation}
Y_{eq}=\frac{45g_{eff}}{4\pi^4} \frac{x^2 K_2(x)}{h_{eff}(m_\chi/x)} ,
\label{DLHA_yeq}
\end{equation}
where $K_2(x)$ is the modified Bessel function of second order.

In various calculations it can be useful to introduce the freeze-out temperature $T_f$ based on the condition 
\begin{equation}
Y-Y_{eq} = \alpha Y_{eq} ,
\label{DLHA_freezecond}
\end{equation}
where $\alpha > 0$ is a number of order 1. Incorporating this condition into Eq. (\ref{DLHA_boltznew2}) a statement for freeze-out is obtained
\begin{equation}
 \sqrt{\frac{\pi g_\star}{45G_N}}\frac{m_\chi}{x^2}\langle\sigma v\rangle \alpha(\alpha+2)Y_{eq} = -\frac{d(\ln Y_{eq})}{dx} .
\label{DLHA_boltznew3}
\end{equation}
Upon substitution of Eq. (\ref{DLHA_yeq}) into Eq. (\ref{DLHA_boltznew3}) this allows for a numerical solution for $x=x_f$.

Equipped with the co-moving number density after freeze-out, $Y_0$, the relic energy density of dark matter may be calculated.  This is typically given in units of the critical energy density
\begin{equation}
 \Omega_\chi=\frac{\rho_\chi}{\rho_c}=\frac{s_0Y_0m_\chi}{\rho_c} ,
\label{DLHA_reldens}
\end{equation}
where $\rho_c$ is the energy density of a flat FLRW Universe.
Here $s_0$ is the current entropy density of the Universe.  Combining Eqs. (\ref{DLHA_hubs}) and (\ref{DLHA_reldens}), the relic density may also be expressed as
\begin{equation}
 \Omega_\chi h^2 = 
 \frac{8 \pi G_N}{3} s_0 Y_0 m_\chi \times 10^{-4} {\rm (s~Mpc/km)}^2 .
\label{DLHA_omegachi}
\end{equation}
The present day normalized Hubble expansion rate $h$ is determined via 
\begin{equation}
 H_0 = 100~h~{\rm km/(s~Mpc)} .
\label{DLHA_eq:def_h}
\end{equation}


\subsubsection{Thermal relic abundance with non-standard cosmological assumptions}

In this section, we consider the generic scenario described in Refs.~\cite{Arbey:2008kv,Arbey:2009gt} and implemented in {\tt SuperIso Relic} \cite{Arbey:2009gu,Arbey:2011zz} and {\tt AlterBBN} \cite{Arbey:2011nf}. In this scenario the total energy density and entropy density of the Universe are modified: 
\begin{equation}
\rho=\rho_{rad}+\rho_{D} ,  ~~~ s=s_{rad}+s_{D} .
\label{DLHA_update1}
\end{equation}
The temperature dependence of the additional components $\rho_{D}$ and $s_{D}$ can be parametrized as
\begin{equation}
\alpha_{D}(T)=\kappa_\alpha \alpha_{rad}(T_{BBN})\left(\frac{T}{T_{BBN}}\right)^{n_\alpha} ,
\label{DLHA_alphadark}
\end{equation}
where $\alpha=\rho$ or $s$, 
\begin{equation}
\kappa_{\alpha}=\frac{\alpha_{dark}(T_{BBN})}{\alpha_{rad}(T_{BBN})} ,
\end{equation}
and $T_{BBN}$ is the Big Bang Nucleosynthesis temperature.

The time evolution of the total entropy density is given by
\begin{equation}
 \frac{ds}{dt}=-3H s+\Sigma_{D},
 \label{DLHA_eq:dsdt}
\end{equation}
where $\Sigma_{D}$ describes the entropy production of the additional component. Instead of parametrizing the entropy density $s_D$, it is sometimes better (for example for reheating models) to write $\Sigma_{D}$ in the form of Eq. (\ref{DLHA_alphadark}) with $\alpha=\Sigma$, and it is then related to the entropy density by
% \begin{equation}
% \Sigma_{rad}(T_{BBN})=T_{BBN}^2 s_{rad}(T_{BBN})\sqrt{\frac{4\pi^3 G}{5}g_{eff}(T_{BBN})}.
% \end{equation}
% It can be brought to the form
\begin{equation}
 \Sigma_{D} = T^2\sqrt{\frac{4\pi^3G}{5}\left(1+\frac{\rho_{D}}{\rho_{rad}}\right)}\left(s_{D} g_{eff}^{1/2}-\frac{1}{3}\frac{h_{eff}}{g_\star^{1/2}}T\frac{ds_{D}}{dT}\right) .
\label{DLHA_en1}
\end{equation}
A possible generalization of the above parameterizations consists of relaxing Eq. (\ref{DLHA_alphadark}) and letting $\rho_D$, $s_D$, $\Sigma_D$ and/or $H$ be general functions of the temperature.

Introducing the co-moving number density $Y = n_{\chi}/s$, Eq.(\ref{DLHA_boltzbest}) becomes
\begin{equation}
 \frac{dY}{dx} = 
  -\sqrt{\frac{\pi g_\star}{45}}\frac{m_\chi}{x^2}\left( \frac{1+s_D/s_{rad}}{\sqrt{1+\rho_D/\rho_{rad}}} \right)
  \left( \langle\sigma v\rangle (Y^2-Y^2_{eq}) +
  \frac{Y \Sigma_D}{s_{rad}^2(1+s_D/s_{rad})^2} \right) ,
\label{DLHA_eq:dYdxD}
\end{equation}
where 
\begin{equation}
 Y_{eq}=\frac{45 g_{eff}}{4\pi^4} \frac{1}{1+s_D/s_{rad}} \frac{x^2 K_2(x)}{h_{eff}(m_\chi/x)} .
\label{DLHA_eq:YeqD}
\end{equation}
The relic abundance can then be calculated based on Eq. (\ref{DLHA_omegachi}).

% \subsubsection{Non-thermal relic abundance}

% We postpone this issue for later DLHA updates.


\subsection{Direct detection}

If dark matter and Standard Model particles interact with strength comparable to that of the electroweak interactions, and if dark matter is massive and fast enough, then it may be detected by observing dark matter scattering on atoms.  This is the aim of the direct dark matter detection experiments.  To assess detection rates we calculate elastic scattering cross sections of dark matter on nuclei.  These cross sections are derived from a Lagrangian describing the effective interaction of dark matter with nuclei.  Nuclear form factors provide a link between the partons, the nucleons and the nucleus of the target species.  
%
For simplicity, in this first agreement, we only consider elastic dark matter-nucleus scattering. Inelastic scattering may be included in a future edition of this agreement.

The dark matter-nucleus recoil rate per unit detector mass, unit time, and unit recoil energy $E$ is written in the form
\begin{equation}
\frac{dR}{dE} = \frac{\sigma_A(E)}{2m_{\chi} \mu_{A}^2} \, \rho_{\chi} \, \eta(E,t),
\end{equation}
where
\begin{equation}
\mu_{A} = \frac{m_{\chi}m_{A}}{m_{\chi}+m_{A}}
\end{equation}
is the dark matter-nucleus reduced mass with $m_A$ the nucleus mass, $\sigma_A(E)$ is the dark matter-nucleus scattering cross section,  $\rho_{\chi}$ is the local dark matter density, and
\begin{equation}
\eta(E,t) = \int_{|{\bf v}|>v_{\rm min}(E)} \frac{f({\bf v},t)}{|{\bf v}|} \, d^3 v
\label{DLHA_eq:etaEt}
\end{equation}
is an average inverse dark matter speed, also called the velocity integral. Here $f({\bf v},t)$ is the dark matter velocity distribution in the reference frame of the detector, which is expected to be time dependent, and 
\begin{equation}
v_{\rm min}(E) = \sqrt{ \frac{m_A E}{2\mu_A^2} } 
\end{equation}
is the minimum dark matter speed that can impart a recoil energy $E$ to the nucleus. Notice that the dark matter-nucleus differential cross section in recoil energy $E$ for a dark matter particle of initial velocity $v$ with respect to the nucleus has the expression
\begin{equation}
\frac{d\sigma_A}{dE} = \frac{\sigma_A(E)}{E_{\rm max}(v)} \, \Theta\!\left( E_{\rm max}(v) - E\right) ,
\end{equation}
where
\begin{equation}
E_{\rm max}(v) = \frac{2 \mu_A^2 v^2}{m_A}
\end{equation}
is the maximum recoil energy that a particle of velocity $v$ can impart to the nucleus.

Only the non-relativistic limit is relevant for dark matter-nucleus scattering. In this limit, only two kinds of nucleon currents survive: the spin-independent current $\psi^\dagger \psi$ and the spin-dependent current $\psi^\dagger {\bf S} \psi$. Here $\psi$ is the nucleon wave function (a non-relativistic two-component Pauli spinor) and ${\bf S}=\frac{1}{2} \boldsymbol\sigma$ is the nucleon spin operator. 
Therefore one splits the differential dark matter-nucleus scattering cross section $\sigma_{A}(E)$  into its spin-independent (SI) and spin-dependent (SD) contributions,
\begin{equation}
\sigma_{A}(E) = \sigma_{A}^{SI}(E) + \sigma_{A}^{SD}(E).
\label{DLHA_sigman}
\end{equation}
Correspondingly, one often separates the spin-independent and spin-dependent contributions to the  recoil rate $dR/dE$ as
\begin{equation}
\frac{dR}{dE} = \left({\frac{dR}{dE}}\right)^{SI} + \left({\frac{dR}{dE}}\right)^{SD}.
\end{equation}

The spin-independent part $\sigma_{A}^{SI}(E)$ is written as
\begin{equation}
\sigma_{A}^{SI}(E) = \frac{4\mu_{A}^{2}}{\pi}\, {\Big|Z \, f_{\rm p} \, F^{(Z,A)}_{\rm p}(E) + (A - Z) \, f_{\rm n} \, F^{(Z,A)}_{\rm n}(E) \Big|}^{2},
\label{DLHA_sigmanSI}
\end{equation}
where 
$Z$ is the number of protons in the nucleus (atomic number), $A$ is the mass number of the nucleus, $F^{(Z,A)}_{\rm p}(E)$ and $F^{(Z,A)}_{\rm n}(E)$ are proton and nucleon number density form factors for the nucleus $(Z,A)$, normalized to $F^{(Z,A)}_{\rm p}(0)=F^{(Z,A)}_{\rm n}(0)=1$, and finally $2f_{\rm p}$ and $2f_{\rm n}$ are the dimensionless four-particle vertices for the SI dark matter-proton and dark matter-neutron interactions, respectively. 

One also introduces the pointlike dark matter-proton and dark matter-neutron cross sections, which by convention are used when reporting or interpreting experimental results,
\begin{equation}
\sigma_{\rm p}^{SI} = \frac{4\mu_{\rm p}^{2}}{\pi}\, {|f_{\rm p}|}^{2},  
~~~~~~
% \label{DLHA_eq:sigmaSIp} \\
\sigma_{\rm n}^{SI} = \frac{4\mu_{\rm n}^{2}}{\pi}\, {|f_{\rm n}|}^{2}. \label{DLHA_eq:sigmaSIn}
\end{equation}
Here $\mu_{\rm p}$ and $\mu_{\rm n}$ are the reduced dark matter-proton and dark matter-neutron masses.

\begin{table}
\caption{Four-particle effective vertices for dark matter-proton elastic scattering. Dark matter-neutron vertices are obtained by changing ${\rm p} \to {\rm n}$. $s_{\chi}$ is the spin of the dark matter particle.}
\begin{center}
\setlength{\unitlength}{0.8in}
\begin{tabular}{cccc}
\hbox to 1in{~}& \hbox to 1.5in{\hfil $s_{\chi}=0$ \hfil} & \hbox to 1.5in{\hfil $s_{\chi}=\frac{1}{2}$ \hfil}  &\hbox to 1.5in{\hfil $s_{\chi}=1$ \hfil} \\ &&& \\
\hline
\begin{picture}(1.3,1)
\put(0,.45){spin-independent}
\end{picture}
&
\begin{picture}(1,1)
\put(.2,.2){\line(1,1){0.6}}
\put(.2,.8){\line(1,-1){0.6}}
\put(.09,.1){${\rm p}$}
\put(.825,.1){${\rm p}$}
\put(.09,.82){${\chi}$}
\put(.825,.82){${\chi}$}
\put(.62,.45){$$}
\end{picture}
&
\begin{picture}(1,1)
\put(.2,.2){\line(1,1){0.6}}
\put(.2,.8){\line(1,-1){0.6}}
\put(.09,.1){${\rm p}$}
\put(.825,.1){${\rm p}$}
\put(.09,.82){${\chi}$}
\put(.825,.82){${\chi}$}
\put(.62,.45){$2f_{\rm p}$}
\end{picture}
&
\begin{picture}(1,1)
\put(.2,.2){\line(1,1){0.6}}
\put(.2,.8){\line(1,-1){0.6}}
\put(.09,.1){${\rm p}$}
\put(.825,.1){${\rm p}$}
\put(.09,.82){${\chi}$}
\put(.825,.82){${\chi}$}
\put(.62,.45){$$}
\end{picture}
\\
\begin{picture}(1.3,1)
\put(0,.45){spin-dependent}
\end{picture}
&
\begin{picture}(1.3,1)
\put(0.45,.45){N/A}
\end{picture}
&
\begin{picture}(1,1)
\put(.2,.2){\line(1,1){0.6}}
\put(.2,.8){\line(1,-1){0.6}}
\put(.09,.1){${\rm p}$}
\put(.825,.1){${\rm p}$}
\put(.09,.82){${\chi}$}
\put(.825,.82){${\chi}$}
\put(.62,.45){$2\sqrt{2}G_Fa_{\rm p} \boldsymbol\sigma_{\chi} \!\cdot\!\boldsymbol\sigma_{\rm p}$}
\end{picture}
&
\begin{picture}(1,1)
\put(.2,.2){\line(1,1){0.6}}
\put(.2,.8){\line(1,-1){0.6}}
\put(.09,.1){${\rm p}$}
\put(.825,.1){${\rm p}$}
\put(.09,.82){${\chi}$}
\put(.825,.82){${\chi}$}
\put(.62,.45){$$}
\end{picture}
\\
\hline
\end{tabular}
\end{center}
\end{table}

One often assumes, and we will do so in this first agreement, that 
\begin{equation}
 F^{(Z,A)}_{\rm p}(E) = F^{(Z,A)}_{\rm n}(E) \equiv F_A(E).
\label{DLHA_eq:F_A}
\end{equation}
In this case, one sometimes introduces the pointlike dark matter-nucleus cross section
\begin{equation}
\sigma_{A,0}^{SI} = \frac{4\mu_{A}^{2}}{\pi}\, {\Big|Z \, f_{\rm p} + (A - Z) \, f_{\rm n} \Big|}^{2},
\end{equation} 
which is $\sigma_{A}^{SI}$ with $F^{(Z,A)}_{\rm p}(E)=F^{(Z,A)}_{\rm n}(E)=1$.

The spin-dependent part $\sigma_{A}^{SD}(E)$ is written as
\begin{equation}
\sigma_{(Z,A)}^{SD}(E) = \frac{32\mu_{A}^{2}G_{F}^{2}}{(2J_{A}+1)}\,\big[a_{\rm p}^{2}\,S_{{\rm pp}}^{(Z,A)}(E) + a_{\rm n}^{2}\,S_{{\rm nn}}^{(Z,A)}(E) + a_{\rm p}\,a_{\rm n}\,S_{{\rm pn}}^{(Z,A)}(E)\big].
\label{DLHA_sigmanSD}
\end{equation}
Here 
%$\hbar$ is the reduced Planck constant, 
$G_F=1.16637\times10^{-5}$ $(\hbar c)^3/{\rm GeV}^{2}$ is the Fermi coupling constant, $J_{A}$ is the nucleus total angular momentum in units of $\hbar$, and finally $2\sqrt{2}\,G_{F}a_{\rm p}$ and $2\sqrt{2}\,G_{F}a_{\rm n}$ are the effective four-particle vertices for the SD interaction of DM particles with protons and neutrons.
 
In Eq.~(\ref{DLHA_sigmanSD}), the dimensionless functions $S_{{\rm pp}}^{(Z,A)}(E)$, $S_{{\rm nn}}^{(Z,A)}(E)$, and $S_{{\rm pn}}^{(Z,A)}(E)$ are the proton-proton, neutron-neutron, and proton-neutron nuclear spin structure functions. They can be written in terms of isoscalar and isovector spin structure functions  $S_{00}^{(Z,A)}$, $S_{11}^{(Z,A)}$, $S_{01}^{(Z,A)}$ as 
\begin{eqnarray}
S_{{\rm pp}}^{(Z,A)} & = & S_{00}^{(Z,A)}+S_{11}^{(Z,A)}+S_{01}^{(Z,A)}\,,
\label{DLHA_Ss1}
\\
 S_{{\rm nn}}^{(Z,A)} & = & S_{00}^{(Z,A)}+S_{11}^{A}-S_{01}^{(Z,A)}\,,
\label{DLHA_Ss2}
\\
S_{{\rm pn}}^{(Z,A)} & = & 2(S_{00}^{(Z,A)}-S_{11}^{(Z,A)})\,.
\label{DLHA_Ss3}
\end{eqnarray}
One similarly introduces
\begin{equation}
a_{0} = a_{\rm p} + a_{\rm n}, ~~~~~~
a_{1} = a_{\rm p} - a_{\rm n} .
\end{equation}

When the nuclear spin is approximated by the spin of the odd nucleon only, one finds
\begin{equation}
S_{{\rm pp}}^{(Z,A)} = \frac{\lambda_{A}^{2}\,J_{A}(J_{A}+1)(2J_{A}+1)}{\pi}\,,\,\,S_{{\rm nn}}^{nA} = 0\,,\,\,S_{{\rm pn}}^{(Z,A)} = 0\,,
\label{DLHA_Spp}
\end{equation}
for a proton-odd nucleus, and
\begin{equation}
S_{{\rm pp}}^{(Z,A)} = 0\,\,,\,S_{{\rm nn}}^{A} = \frac{\lambda_{A}^{2}\,J_{A}(J_{A}+1)(2J_{A}+1)}{\pi}\,,\,\,S_{{\rm pn}}^{(Z,A)} = 0\,,
\end{equation}
for a neutron-odd nucleus. Here $\lambda_{A}$ is conventionally defined through the relation $\langle\psi_A|{\bf S}_A|\psi_A\rangle = \lambda_{A}\, \langle \psi_A|{\bf J}_A|\psi_A \rangle$, where $|\psi_A \rangle$ is the nuclear state, ${\bf S}_A$ is the nucleus spin vector, and ${\bf J}_A$ is its total angular momentum vector.

One sometimes introduces dark matter-nucleon pointlike cross sections. For a single proton or a single neutron, $\lambda_{\rm p}=\lambda_{\rm n}=1$, $J_{\rm p}=J_{\rm n}=\frac{1}{2}$, and
\begin{equation}
\sigma_{\rm p}^{SD} = \frac{36\mu_{\rm p}^{2}G_{F}^{2}}{\pi^2}\,|a_{\rm p}|^{2}, 
~~~~~~
% \label{DLHA_eq:sigmaSDp} \\
\sigma_{\rm n}^{SD} = \frac{36\mu_{\rm n}^{2}G_{F}^{2}}{\pi^2}\,|a_{\rm n}|^{2}. 
\label{DLHA_eq:sigmaSDn}
\end{equation}

We may list the expressions of the effective four-particle dark matter-nucleon couplings $f_{\rm p}$, $f_{\rm n}$, $a_{\rm p}$ and $a_{\rm n}$ in terms of elementary couplings to quarks and gluons for various kinds of dark matter particles, via an effective four-particle lagrangian with spin $s$
\begin{equation}
 {\cal L}_{eff}(s) = 
 \frac{1}{2}\left(\sum_{i=e,o} f_{N,i}(s) {\cal L}_{N,i}^{SI}(s) +
 \sum_{i=e,o} a_{N,i}(s)     {\cal L}_{N,i}^{SD}(s)\right)   .
\label{DLHA_eq:L_eff}
\end{equation}
The operators ${\cal L}_e^{SI}(s)$ and ${\cal L}_o^{SI}(s)$ describe spin independent even and odd interactions, while ${\cal L}_e^{SD}(s)$ and ${\cal L}_o^{SD}(s)$ are their spin dependent counterparts.\footnote{Even and odd couplings are introduced to trace symmetries under particle-antiparticle interchange.  Majorana fermions, for example, have even couplings and non-Majorana fermions have odd couplings.}
Table \ref{DLHA_tableop} shows the explicit forms of these operators for various values of $s$.  
%
The dark matter-nucleon couplings can be obtained as
\begin{equation}
 f_{\rm N} = \sum_{i=e,o} f_{N,i}(s), ~~~ 
 a_{\rm N} = \sum_{i=e,o} a_{N,i}(s),
\end{equation}

\begin{table}
\begin{center}
\caption{Even and odd operators ${\cal L}_{e,o}^{SI,SD}(s)$ for dark matter interactions with standard quarks $q$.  A scalar field only interacts in a spin independent manner \cite{micromegas22}.}
\vspace{3mm}
\begin{tabular}{|l|c|c|c|}
\hline

 & s & Even operators ($i=e$) & Odd operators ($i=o$) \\ \hline

 Spin independent & 0 & $2M_\chi \chi\chi^\dagger\bar{q}q$ & $i(\partial_\mu\chi\chi^\dagger-\chi\partial_\mu\chi^\dagger)\bar{q}\gamma^\mu q$ \\

 ($SI$) & 1/2 & $\bar{\chi}\chi \bar{q}q$ & $\bar{\chi}\gamma_\mu\chi\bar{q}\gamma^\mu q$ \\
 
 & 1 & $2M_\chi \chi_\mu\chi^\mu\bar{q}q$ & $i(\chi^{\dagger\alpha}\partial_\mu\chi,\alpha-\chi^\alpha\partial_\mu\chi^\dagger_\alpha)\bar{q}\gamma_\mu q$ \\ \hline

 Spin dependent & 1/2 & $\bar{\chi}\gamma_5\gamma_\mu\chi\bar{q}\gamma_5\gamma^\mu q$ & $-\frac{1}{2}\bar{\chi}\sigma_{\mu\nu}\chi\bar{q}\sigma^{\mu\nu}q$ \\
 ($SD$) & 1 & $\sqrt{6}(\partial_\mu\chi^{\dagger\beta}\chi_\gamma-\chi^\dagger_\beta\partial_\mu\chi_\gamma)\epsilon^{\alpha\beta\gamma\mu}\bar{q}\gamma_5\gamma_\mu q$ & $\frac{\sqrt{3}}{2}i(\chi_\mu\chi^{\dagger}_\nu-\chi^\dagger_\mu\chi_\nu)\bar{q}\sigma^{\mu\nu}q$ \\
 
\hline
\end{tabular}\label{DLHA_tableop}
\end{center}
\end{table}

Dark matter-nucleus scattering amplitudes can be calculated based on parton level amplitudes after relating nuclear couplings to dark-matter quark couplings via nucleon form factors
\begin{eqnarray}
 2f_{N,e}=\sum_{q}\frac{m_p}{m_q} {f_q^N s_{q,e}} , ~~~
 2f_{N,o}=\sum_{q}\frac{m_p}{m_q}{f_{V_q}^N s_{q,o}} , \label{DLHA_nuc1form} \\
 2\sqrt{2}G_Fa_{N,e}=\sum_{q}{\Delta_q^N a_{q,e}} , ~~~
 2\sqrt{2}G_Fa_{N,o}=\sum_{q}{\delta_q^N a_{q,o}}  .
\label{DLHA_nuc2form}
\end{eqnarray}
These form factors capture the distribution of quarks within the nucleons.  The light quark flavor scalar form factors are related to the pion-nucleon sigma term and the nucleon and quark masses as
\begin{equation}
 f_u^{p,n} = \frac{m_u}{m_d} \alpha^{+1,-1} f_d^{p,n} , ~~~
 f_d^{p,n} = \frac{2 \sigma_{\pi N}}{(1 + \frac{m_u}{m_d}) m_{p,n}} \frac{\alpha^{0,1}}{1 + \alpha} , ~~~
 f_s^{p,n} = \frac{\sigma_{\pi N} y}{(1 + \frac{m_u}{m_d}) m_{p,n}} \frac{m_s}{m_d} .
\label{DLHA_eq:f_q^N}
\end{equation}
Here $\alpha=B_u/B_d$, $\sigma_{\pi N}=(m_u+m_d)(B_u+B_d)/2$, $y=2B_s/(B_u+B_d)$ and $B_q=\langle N|{\bar q}q|N \rangle$.
The heavy flavor scalar form factors are typically calculated as 
\begin{equation}
 f_Q^N = \frac{2}{27} \left( 1-\sum_{q=u,d,s} f_q^N \right) .
\label{DLHA_eq:f_Q^N}
\end{equation}

\subsection{Indirect detection}

Perhaps the most challenging way to discover dark matter is to detect its annihilation or decay products in the astrophysical environment.  This is called indirect detection and involves the description of an initial (source) distribution of standard particles originating from dense dark matter concentrations.  After the source properties are fixed, the propagation of the secondary particles has to be followed through.  Here we review only the simplest cases: electron, positron, antiproton or photon production via dark matter annihilation and their subsequent propagation to us through our Galaxy using a simplified but effective treatment. We also comment on the modifications which would need to be introduced for a more detailed treatment.

Charged cosmic ray propagation through the Galaxy can be usually described by the diffusion-convection model \cite{Ginzburg:1990sk}. This model assumes homogeneous propagation of charged particles within a certain diffusive region (similar to one of the simplest models of propagation called the leaky box model), but it also takes into account cooling (energy loss) effects. The diffusive region is usually assumed to have the shape of a solid flat cylinder of half-height $L$ and radius $R$ that sandwiches the Galactic plane: inside it, charged cosmic rays are trapped by magnetic fields; outside, they are free to stream away. The (cylindrical) coordinates of the solar system correspond to $\vec r_\odot = (8.33 \ {\rm kpc}, 0  \ {\rm kpc})$ \cite{Gillessen:2008qv}. The phase-space density $\psi_a (\vec r,p,t)$ of a particular cosmic ray species $a$ at an instant $t$, at a Galactic position $\vec r$ and with momentum $p$ can be calculated by solving the cosmic ray transport equation, which has the general form \cite{Strong:2007nh}
\begin{eqnarray}
 \frac{\partial \psi_a (\vec r,p,t)}{\partial t} 
 & = &
 Q_a(\vec r, p, t) 
 + \nabla \cdot ( D_{xx}\nabla\psi_a - \vec V\psi_a ) \nonumber \\
 & + & \frac{\partial}{\partial p}\, \left( p^2 D_{pp} \frac{\partial}{\partial p}\, \frac{1}{p^2}\, \psi_a \right)  
 - \frac{\partial}{\partial p} \left( \dot{p} \psi_a
 - \frac{p}{3} \, (\nabla \cdot \vec V )\psi_a \right)
 - \frac{1}{\tau_f}\psi_a - \frac{1}{\tau_r}\psi_a .
\label{DLHA_eq:transport}
\end{eqnarray}
Usually one assumes that steady state conditions hold, as they do if the typical time scales of the dark matter galactic collapse and of the variation of propagation conditions are much longer than the time scale of propagation itself (which is of the order of 1 Myr at 100 GeV energies). In this case, the l.h.s. can be equated to zero and the dependence on time is dropped for all quantities.

We concentrate on a version of the propagation equation which is sufficient to describe in first approximation the electron or positron and proton or antiproton flux through the Galaxy: 
\begin{equation}
 0 = 
 Q_a(\vec r,E) +
K(E)\ \nabla^2\psi_a +
 \frac{\partial}{\partial E}\Big(b(E)\ \psi_a - \ K_{\rm EE}(E)\ \psi_a\Big) - 
 \frac{\partial}{\partial z}\left({\rm sign}(z)\, V_C \ \psi_a\right),
\label{DLHA_evol_energy}
\end{equation}
where now $E$ is the energy of the secondary particle species $a$. Boundary conditions are imposed such that the CR density vanishes on the outer surface of the cylinder, outside of which the particles are supposed to freely propagate and escape, consistently with the physical picture described above. At $r=0$, one imposes a symmetric condition $\partial \psi_a / \partial r (r=0) = 0$. In momentum space one imposes null boundary conditions. 
The terms containing the spatial diffusion coefficient $K(E)$, the energy loss rate $b(E)$ and the diffusive reacceleration coefficient $K_{\rm EE}(E)$ describe respectively the transport of cosmic ray species through turbulent magnetic fields, their cooling due to different phenomena (such as Inverse Compton scattering (ICS), synchrotron radiation, Coulomb scattering or bremsstrahlung) and their reacceleration due to hits on moving magnetized scattering targets in the Galaxy. The term with the convective velocity $V_C$ describes the characteristics of the Galactic winds emanating vertically from the stars in the disk.  A source term resulting from dark matter annihilation can be written as
\begin{equation}
 Q_a(\vec r,E)=\frac{1}{2}\frac{dN_a}{dE}\langle\sigma_a v\rangle_0 \left(\frac{\rho_g(\vec r)}{m_\chi}\right)^2 .
\label{DLHA_source}
\end{equation}
Here $\langle\sigma_a v \rangle_0$ is the value of the thermally averaged annihilation cross section into the relevant species, and $\rho_g(r)$ is the dark matter energy density in the Galaxy.
The energy distribution of the secondary particle $a$ is $dN_a/dE$, normalized per annihilation.  This formula applies to self-conjugated annihilating dark matter. 
%
% If dark matter is not composed of self-conjugated particles, and $n$ indicates the number density of particles and ${\overline n}$ the number density of antiparticles, the factor $n^2/2$ in the formula above has to be replaced by $n {\overline n}$.
In the case of non-self-conjugated dark matter, or of multicomponent dark matter, the quantities in Eq.~(\ref{DLHA_source}) should be replaced as follows, where an index $i$ denotes a charge state and/or particle species (indeed any particle property, collectively called "component") and $f_i=n_i/n$ is the number fraction of the $i$-th component:
\begin{eqnarray}
m_\chi&  \to&   \sum_i f_i m_i \qquad\qquad\qquad\qquad\qquad \text{(mean mass)},\\
\langle \sigma_a v \rangle&  \to&  \sum_{ij} f_i f_j \langle \sigma_{a,ij} v_{ij} \rangle \qquad\qquad\qquad\,\,\,\text{(mean cross section times relative velocity)},\\
dN_a/dE&  \to&  \frac{\sum_{ij} f_i f_j \sigma_{a,ij} v_{ij} \, (dN_{a,ij}/dE)}{\sum_{ij} f_i f_j \sigma_{a,ij} v_{ij} }, \quad \text{(annihilation spectrum per annihilation).}
\end{eqnarray}

The spatial diffusion coefficient $K(E)$ is generally taken to have the form
\begin{eqnarray}
K(E) =  K_0 \ v^\eta \left(\frac{R}{\rm GeV}\right)^\delta ,
\label{DLHA_eq:KE}
\end{eqnarray}
where $v$ is the speed (in units of $c$) and $R=p/eZ$ is the magnetic rigidity of the cosmic ray particles.  Here $Z$ is the effective nuclear charge of the particle and $e$ is the absolute value of its electric charge (of course the quantities are different from 1 only in the case in which particles other than electrons, positrons, protons or antiprotons are considered). The parameter $\eta$ controls the behavior of diffusion at low energy: recently, in departure from the traditional choice $\eta = 1$, other values (possibly negative) have been advocated. In more detailed treatments, the diffusion coefficient can be considered as space dependent ($K(\vec r, E)$) and a possible dependence on the particle direction of motion, leading to anisotropic diffusion, can be introduced. 

The energy loss rate can be parametrized as
\begin{eqnarray}
b(E) = b_0 \ E^2 .
\label{DLHA_eq:bE}
\end{eqnarray}
This form holds as long as one neglects the fact that energy losses are position dependent in the Galactic halo (e.g. synchrotron radiation depends on the intensity of the magnetic field, which varies in the Galaxy, and ICS depends on the density of the background light distribution, which also varies in the Galaxy), and as long as one assumes that all energy loss phenomena are proportional to $E^2$ (which is true only if one neglects Coulomb losses and bremsstrahlung and one considers ICS only in the Thomson scattering regime, i.e. at relatively low electron or positron energy). In a more detailed treatment, the energy loss rate can also be considered as space dependent, $b(\vec r, E)$, and with a more general dependence on the energy. Coulomb losses ($dE/dt \sim$ const) and bremsstrahlung losses ($dE/dt \sim bE$) can also be taken into account.  Some codes compute these using detailed formulae as a function of position and energy, based on the gas, interstellar radiation and magnetic field distributions \cite{Strong:2009xj}.  The ultimate form of the energy loss rate may be a complex function which, as will be described below, can be given under \verb|FUNCTION EnerLoss|. 

Finally, the diffusive reacceleration coefficient $K_{\rm EE}(E)$ is usually parameterized as 
\begin{eqnarray}
K_{\rm EE}(E) = \frac{2}{9} v_{\rm A}^2 \frac{v^4 \, E^2}{K(E)},
\label{DLHA_eq:KEE}
\end{eqnarray}
where $v_{\rm A}$ is the Alfv\'en speed.

A propagator, or Green's function $G$, is used to evolve the flux which originates from the source $Q$ at $\vec r_S$ with energy $E_S$ through the diffusive halo, to reach the Earth at point $\vec r$ with energy $E$.  This allows the general solution for Eq. (\ref{DLHA_evol_energy}) to be written as
\begin{equation}
 \psi_a(\vec r,E)=\int_E^{m_\chi}dE_S \int{d^3r_S \, G(\vec r,E;\vec r_S,E_S)\, Q_a(\vec r_S,E_S)} .
\label{DLHA_energysoln}
\end{equation}
The differential flux is related to the solution in Eq. (\ref{DLHA_energysoln}) via
\begin{equation}
 \frac{d\Phi_a}{dE}=\frac{v(E)}{4\pi}\psi_a(\vec r,E) .
\label{DLHA_eq:dPhi_a/dD}
\end{equation}
For proton or antiproton propagation in the Galactic halo, additional terms in Eq. (\ref{DLHA_evol_energy}) should be introduced to account for spallations on the gas in the disk. 

\medskip

We next turn to photons. In this case, propagation is much simpler since no scattering processes take place. This leads to a straight line propagation without any energy loss, which, in the formalism above, corresponds to a trivial propagator. Thus, for photons $a = \gamma$, by utilising Eqs. (\ref{DLHA_source}) and (\ref{DLHA_energysoln}), Eq. (\ref{DLHA_eq:dPhi_a/dD}) can be approximated with
\begin{equation}
 \frac{d\Phi_\gamma}{dE}(\phi) =  \langle\sigma_\gamma v\rangle_0 \frac{dN_\gamma}{dE}\frac{1}{8\pi \, m_\chi^2}\int_0^\infty{\rho_g^2(r(s,\phi))\, ds} .
\label{DLHA_flux}
\end{equation}
Here the square of the Galactic dark matter density profile $\rho_g^2$ is integrated over the line of sight, parameterized by the coordinate $s$.  The angle $\phi$ is the aperture between the direction of the line of sight and the axis connecting the Earth to the Galactic Center. Explicitly, the coordinate $r$, centered on the Galactic Center, reads 
\begin{equation}
 r(s,\phi) = (r_\odot^2+s^2-2\,r_\odot\,s\cos\phi)^{1/2} . 
\label{DLHA_rsphi}
\end{equation}
As for Eq.~(\ref{DLHA_source}), Eq.~(\ref{DLHA_flux}) applies to self-conjugated annihilating dark matter.  
%
If dark matter is not composed of self-conjugated particles, and $n$ indicates the number density of particles and ${\overline n}$ the number density of antiparticles, the factor $n^2/2$ in the formula above has to be replaced by $n {\overline n}$.

%Piece added by Lars Bergstrom:
This formula will depend very sensitively on the direction, especially near the galactic center. However, a given gamma-ray instrument will only measure an averaged value, smeared over the angular resolution. Therefore, a more useful quantity to compute is that average value in the direction given by $\phi$ \cite{Bringmann:2008kj},
\begin{equation}
\frac{d\Phi_ \gamma}{dE}(\phi;\Delta\Omega)= \int d\Omega' \frac{d\Phi_ \gamma}{dE}(\varphi',\theta')
R_{\Delta\Omega}(\theta'), 
\end{equation}
where
\begin{equation}
R_{\Delta\Omega}(\theta')= \frac{1}{2\pi\theta_r^2}\exp\left(\frac{-\theta'^2}{2\theta_r^2}\right), 
\end{equation}
describes the angular resolution $\theta_r$ with $d\Omega'=d\varphi' d\cos \theta'$. Here $\phi'$ and $\theta'$ are polar coordinates centered on the direction $\phi$, and $\Delta\Omega = \pi\theta_r^2$ (for small $\theta_r$). For the different, slightly varying lines of sight entering the integral, $\cos\phi$ in Eq.~(\ref{DLHA_rsphi}) is replaced by
$\cos\psi=\cos\phi\cos\theta'-\cos\varphi'\sin\phi\sin\theta'$.


\subsubsection{Dark matter substructures}

If dark matter particles are packed inside dense clumps, their annihilations are enhanced, and so are their indirect signatures at the Earth.  The boost factor by which the yield of a smooth dark matter halo has to be multiplied depends in a simple, but not obvious, way on the spatial distribution of the clumps and on their inner structure.  A population ${\cal P}$ of substructures $i$ generates at the Earth the cosmic ray density 
\begin{equation}
 \psi_a^{sub}({\vec r},E)=
 \left( {\cal S} \equiv \frac{1}{2} \langle\sigma_a v\rangle_0 \left(\frac{\rho_{\odot}}{m_\chi}\right)^2 \right) \,
 {\displaystyle \sum_{i \in {\cal P}}} \; \tilde{G}_{i} \; \xi_{i} \;\; ,
\end{equation}
where the effective propagator $\tilde{G}_{i}$, defined as
\begin{equation}
 \tilde{G}_{i}=\int_E^{m_\chi}dE_S \, G({\vec r},E;{\vec r}_i,E_S) \, \frac{dN_a}{dE} ,
\end{equation}
takes into account the propagation from clump $i$ located
at ${\vec r}_{i}$, and the injection spectrum at the source.
The annihilation volume $\xi_{i}$ would be the volume of
clump $i$ should its density be equal to the Milky Way
dark matter density $\rho_{\odot}$ at the Sun. It is
defined as the integral over the volume of the $i$-th clump
\begin{equation}
\xi_{i}=\int_{\rm clump \, i} d^3r_S \left( \frac{\rho_{\rm DM}(\vec r_S)}{\rho_{\odot}} \right)^2 .
\end{equation}
%
Because we have no idea of the actual population of
dark matter substructures inside which we are embedded,
a statistical analysis needs to be performed on the
ensemble of all possible realizations of galactic
clump distributions. A population of ${\cal N}_{\rm H}$
substructures inside the Milky Way dark matter halo
yields, on average, the cosmic ray density
\cite{Lavalle:2006vb, Brun:2007tn, Lavalle:1900wn}
\begin{equation}
 \langle \psi_a^{sub}({\vec r},E) \rangle=
 {\cal N}_{\rm H} \, {\cal S} \,
%\frac{1}{2} \langle\sigma_a v\rangle_0 \left(\frac{\rho_{\odot}}{m_\chi}\right)^2 \, {\cal N}_{\rm H} \,
 \int\!d^3r_S \int\!d\xi \; {\cal D}({\vec r}_S,\xi) \, \tilde{G}({\vec r},E;{\vec r}_S) \, \xi .
\end{equation}
The propability to find a dark matter clump at location
${\vec r}_S$ with annihilation volume $\xi$ is denoted
by ${\cal D}({\vec r}_S,\xi)$. The variance associated
to the average substructure signal $\langle \psi_a^{sub}({\vec r},E) \rangle$
can be expressed as
\begin{equation}
 \sigma_{\psi}^{2}({\vec r},E)=
 {\cal N}_{\rm H} \, {\cal S}^{2} \,
 \int\!d^3r_S \int\!d\xi \; {\cal D}({\vec r}_S,\xi) \, \tilde{G}^{2}({\vec r},E;{\vec r}_S) \, \xi^{2}
 \, - \,
 {\displaystyle \frac{\langle \psi_a^{sub}({\vec r},E) \rangle^{2}}{{\cal N}_{\rm H}}} .
\end{equation}


\section{DESCRIPTION OF THE DLHA BLOCKS AND FUNCTIONS}

Information in a DLHA file is organized into blocks.  The general properties of the DLHA blocks follow those of the Supersymmetry Les Houches Accord (SLHA) \cite{Skands:2003cj}.  Similarly to SLHA, the entries within the blocks are identified by the first numerical value(s) within the block.  This feature allows for a flexible order of entries within a block.  Most block entries are optional and can be omitted at writing.  A missing entry typically signals the lack of a calculation within the program that wrote the block.

A DLHA file may contain the following blocks or statements (listed here alphabetically):
\begin{quote}
\verb|BLOCK ABUNDANCE|       \\
\verb|BLOCK ANNIHILATION|    \\
\verb|BLOCK ASTROPROPAG|     \\
\verb|BLOCK COSMOLOGY|       \\
\verb|BLOCK DETECTOR_NUCLEI| \\
\verb|BLOCK DMCLUMPS|        \\
\verb|BLOCK DMSPADIST|       \\
\verb|BLOCK DMVELDIST|       \\
\verb|BLOCK DOFREEDOM|       \\
\verb|BLOCK EFFCOUPLING|     \\
\verb|BLOCK FORMFACTS|       \\
\verb|BLOCK INDIRDETSPECTRUM|\\
\verb|BLOCK MASS|            \\
\verb|BLOCK NDMCROSSSECT|    \\
\verb|BLOCK QNUMBERS|        \\
\verb|BLOCK STRUCTFUN| 
\end{quote}
Blocks \verb|ASTROPROPAG|, \verb|COSMOLOGY|, \verb|DETECTOR_NUCLEI|, \verb|DMCLUMPS|, \verb|DMSPADIST|, \verb|DMVELDIST|, \verb|DOFREEDOM|, \verb|FORMFACTS| and \verb|STRUCTFUN| depend on the cosmological, astrophysical, standard particle and nuclear physics assumptions but are kept independent from the microscopic properties of dark matter.  
% Blocks \verb|DMSPADIST|, \verb|DMVELDIST| depend on the particular galaxy under consideration and thus carry an index \verb|j| identifying that galaxy.
The rest of the blocks depend on the microscopic physics describing dark matter.  A given dark matter candidate is identified in the block \verb|MASS| by the PDG number and mass of the particle.  If the PDG code of a particle does not exist then an arbitrary code identifying the candidate can be supplied.  Further microscopic properties of the dark matter candidate are given in block \verb|QNUMBERS|. 

For dark matter models containing multiple dark matter candidates a separate DLHA file has to be created for each candidate.  If multiple dark matter candidates contribute simultaneously to the present abundance and direct or indirect detection signals, multiple sets of blocks may appear for each candidate in separate files.  For example, a DLHA file may contain 
\begin{quote}
\verb|BLOCK MASS| \\
\verb|#  PDG code   mass                particle name | \\
\verb|   1000022    1.29098165E+00   #  1st neutralino| \\
\verb|...|             \\
\verb|BLOCK QNUMBERS|  \\
\verb|...|             \\
\verb|BLOCK ABUNDANCE| \\
\verb|...| 
\end{quote}
while another DLHA file might contain
\begin{quote}
\verb|BLOCK MASS| \\
\verb|#  PDG code   mass                particle name | \\
\verb|   1000039    2.35019093E+00   #  gravitino     | \\
\verb|BLOCK QNUMBERS|  \\
\verb|...|             \\
\verb|BLOCK ABUNDANCE| \\
\verb|...| 
\end{quote}
Here the ellipses denote entries irrelevant to this discussion.  A user of DLHA is responsible for knowing that the two files contain information in the context of multiple dark matter candidates. 

For decaying dark matter particles a standard SLHA decay file can be used to read and write the total decay width of the dark matter particle and its branching ratios into various final states.


\subsection{The FUNCTION object}

Departing from the tradition of previous Les Houches accords, DLHA introduces a new structure that specifies a function.  A function definition is facilitated by DLHA using the following construction
\begin{quote}
\verb|FUNCTION <name> type=<type> args=<number of arguments>| \\
\verb| ...                                                  | \\
\verb|END_FUNCTION| 
\end{quote}
The \verb|FUNCTION| heading denotes the beginning and \verb|END_FUNCTION| the end of the structure.  Each function is identified by a name, which follows the \verb|FUNCTION| heading.  In the various block descriptions below we fix the names of the possible functions and specify their content.  The content of a function can be given in several different ways in DLHA.  The \verb|<type>| variable differentiates between these methods:
\begin{description}
\item[] \verb|type = P| for a predefined function,
\item[] \verb|type = C| for a C function,
\item[] \verb|type = F| for a Fortran function,
\item[] \verb|type = T| for tabular information.
\end{description}
The number of independent variables of the function is given by the numerical value of the last argument \verb|<number of arguments>| of the \verb|FUNCTION| structure.  

While a \verb|FUNCTION| depends on one or more independent variables, it may also carry information about related parameters inside the body of the function.  These parameters can be listed as follows:
\begin{quote}
\verb|FUNCTION <name> type=<type> args=<number of arguments>| \\
\verb| PARAMETERS                                           | \\
\verb|  <parameter name 1>=<value 1>                        | \\
\verb|  <parameter name 2>=<value 2>                        | \\
\verb|  <parameter name 3>=<value 3>                        | \\
\verb|  ...                                                 | \\
\verb| END_PARAMETERS                                       | \\
\verb| <function body>                                      | \\
\verb|END_FUNCTION                                          | 
\end{quote}
Parameter names are fixed by DLHA in the block descriptions similarly to names of functions.  If a parameter value is specified both outside and inside of a function, the value given inside the FUNCTION construction overrides the one appearing outside. 

Predefined functions, typically the most commonly used functions for a given quantity, are specified by DLHA under the description of the various functions.  A predefined function which is listed in the DLHA write-up can be referred to by the following construction following the \verb|FUNCTION| heading:
\begin{quote}
\verb|FUNCTION <name> type=<type> args=<number of arguments>| \\
\verb| DLHA <name> <identifier>                             | \\
\verb| ...                                                  | \\
\verb|END_FUNCTION| 
\end{quote}
Within the block descriptions DLHA fixes the function choices corresponding to various numerical values of identifiers.

As an example, the following structure specifies the use of the Einasto profile for the dark matter galactic halo profile $\rho_g(r)$:
\begin{quote}
\verb|FUNCTION rho_g type=P args=2| \\
\verb| DLHA rho_g 5               | \\
\verb|END_FUNCTION                |
\end{quote}
The Einasto profile depends on two arguments.  One of them may also be fixed as:
\begin{quote}
\verb|FUNCTION rho_g type=P args=1| \\
\verb| DLHA rho_g 5               | \\
\verb| PARAMETERS                 | \\
\verb|  alpha=1                   | \\
\verb| END_PARAMETERS             | \\
\verb|END_FUNCTION                |
\end{quote}

In the present agreement, we only consider functions tabulated on a rectangular (but not necessary equidistant) grid.  In this case the list of the independent variables and the function value is given as an \verb|n|+1 column table.  For tabular functions the line after the function name gives the names of the independent variables and the dimensions of the rectangular grid on which the function is specified.  A schematic example of a function given in tabular format is the following:
\begin{quote}
\verb|FUNCTION rho_g type=T args=2                |\\
\verb| r: 2  alpha: 3                             |\\
\verb|#  r             alpha         rho_g        |\\
\verb|   0.000000E+00  1.000000E+00  2.345678E+00 |\\
\verb|   1.000000E-05  1.000000E+00  1.234567E+00 |\\
\verb|   0.000000E-05  2.000000E+00  4.567890E+00 |\\
\verb|   1.000000E-05  2.000000E+00  3.456789E+00 |\\
\verb|   0.000000E-05  3.000000E+00  6.789012E+00 |\\
\verb|   1.000000E-05  3.000000E+00  5.678901E+00 |\\
\verb|END_FUNCTION                                |
\end{quote}

A \verb|FUNCTION| given in the form of a C language function could appear as:
\begin{quote}
\verb|FUNCTION rho_g type=C args=2              |\\
\verb| #include<math.h>                         |\\
\verb| double Einasto(double r, double alpha)   |\\
\verb| { return exp(-2*(pow(r,alpha)-1)/alpha);}|\\
\verb|END_FUNCTION                              |
\end{quote}

The \verb|FUNCTION| construction also allows passing function names that are included in a pre-compiled library:
\begin{quote}
\verb|FUNCTION <name> type=<type> args=<number of arguments>|\\
\verb| libName=<name of compiled library>                   |\\
\verb| funcName=<name of function in library>               |\\
\verb|END_FUNCTION                              |
\end{quote}
An example of this is shown in Section \ref{DLHA_subsec:BDN}.

% In the rest of this writeup we describe the content of each of the DLHA blocks in detail.


\subsection{BLOCK COSMOLOGY}

The \verb|COSMOLOGY| block specifies the values of the cosmological parameters that enter into the calculation of dark matter related observables such as the abundance.  The cosmology block may contain numerical entries such as those below.  In these entries, the temperature is given in units of GeV$/k$, where $k$ is the Boltzmann constant.

\begin{description}
\item[1] The value of the current entropy density of the Universe $s_0$ in units of $({\rm GeV}/\hbar c)^3/k$, appearing in Eq. (\ref{DLHA_omegachi}).
\end{description}

The cosmology block also accommodates various entries that are defined as functions as described in the previous subsection.  The list of these functions is given below.
 
Non-standard energy density  $\rho_{D}(T)$, as described in Eq. (\ref{DLHA_update1}), as a function of the temperature and in units of $({\rm GeV}/\hbar c)^4$.
\begin{quote}
\verb|FUNCTION rho_D type=<type> args=1|\\
\verb| ...                             |\\
\verb|END_FUNCTION                     |
\end{quote}

Non-standard entropy density  $s_{D}(T)$, as described in Eq. (\ref{DLHA_update1}), as a function of the temperature and  in units of $({\rm GeV}/\hbar c)^3$.
\begin{quote}
\verb|FUNCTION s_D type=<type> args=1|\\
\verb| ...                           |\\
\verb|END_FUNCTION                   |
\end{quote}

Non-standard normalized entropy production rate $\Sigma_{D}(T)/\sqrt{G}$, as described in Eq. (\ref{DLHA_update1}), as a function of the temperature and in units of $({\rm GeV}/\hbar c)^5$.
\begin{quote}
\verb|FUNCTION Sigma_D type=<type> args=1|\\
\verb| ...                               |\\
\verb|END_FUNCTION                       |
\end{quote}

The value of the normalized Hubble expansion rate $H(T)/\sqrt{G}$, as per Eq. (\ref{DLHA_eq:def_h}), as a function of the temperature.  The normalized Hubble expansion rate is in units of $({\rm GeV}/\hbar c)^2$.
\begin{quote}
\verb|FUNCTION Hubble type=<type> args=1|\\
\verb| ...                              |\\
\verb|END_FUNCTION                      |
\end{quote}

Note that the functions \verb|s_D| and \verb|Sigma_D| are not to be given simultaneously, as they refer to two different parameterizations of the entropy content of the Universe. Also, the functions \verb|rho_D| and \verb|Hubble| do not need to be given simultaneously.

A specific cosmology block, describing a standard cosmological scenario may appear as this:
\begin{quote}
\verb|BLOCK COSMOLOGY| \\
\verb|# identifier(s)  parameter value   comment         | \\
\verb|  1              2.18421585E-03  # s_0 [GeV^3]     | \\
\verb|# -------------------------------------------------| \\
\verb|FUNCTION rho_D type=C args=1     # rho_D(T) [GeV^4]| \\
\verb| ...                                               | \\
\verb|END_FUNCTION                                       | \\
\verb|FUNCTION s_D type=F args=1       # s_D(T) [GeV^3]  | \\
\verb| ...                                               | \\
\verb|END_FUNCTION                                       | 
\end{quote}


\subsection{BLOCK DOFREEDOM}

This block contains the various degrees of freedom entering into the abundance calculation.

The effective degrees of freedom $g_{eff}(T)$ in Eq. (\ref{DLHA_energydens}) as a function of the temperature are given by   
\begin{quote}
\verb|FUNCTION g_eff type=<type> args=1         |\\
\verb| ...                                      |\\
\verb|END_FUNCTION                              |
\end{quote}
The units of temperature are GeV$/k$.

As an alternative to $g_{eff}(T)$, as a function of the temperature, $g_\star$ can also be given.  The latter is defined in Eq. (\ref{DLHA_eq:gstar}).  
\begin{quote}
\verb|FUNCTION gstar type=<type> args=1         |\\
\verb| ...                                      |\\
\verb|END_FUNCTION                              |
\end{quote}

The effective degrees of freedom $h_{eff}(T)$ in Eq. (\ref{DLHA_entropy}) as a function of the temperature (measured in GeV$/k$) are given by
\begin{quote}
\verb|FUNCTION h_eff type=<type> args=1         |\\
\verb| ...                                      |\\
\verb|END_FUNCTION                              |
\end{quote}

This block containing 57 entries of $g_{eff}$ and 48 entries for $h_{eff}$ in tabular format might appear as follows.  
\begin{quote}
\verb|BLOCK DOFREEDOM| \\
\verb|FUNCTION g_eff type=T args=1     | \\
\verb| T: 57                           | \\
\verb|# T [GeV/k]         g_{eff}(T)   | \\
\verb|  1.00000000E-05  2.00000000E-00 | \\
\verb|  5.00000000E-05  2.00000000E-00 | \\
\verb|  1.00000000E-04  4.00000000E-00 | \\
\verb|...                              | \\
\verb|END_FUNCTION                     | \\
\verb|FUNCTION h_eff type=T args=1     | \\
\verb| T: 48                           | \\
\verb|# T [GeV/k]         h_{eff}(T)   | \\
\verb|  2.00000000E-05  1.00000000E-00 | \\
\verb|  4.00000000E-05  1.00000000E-00 | \\
\verb|  6.00000000E-05  2.00000000E-00 | \\
\verb|...                              | \\
\verb|END_FUNCTION                     | 
\end{quote}


\subsection{BLOCK DMSPADIST}

Block \verb|DMSPADIST| contains information regarding the dark matter energy density distribution in various astrophysical objects, such as galaxies or nebulae.  These enter the calculation of direct and indirect detection rates.  
%
Before discussing the block attributes, we list the most commonly used halo profiles and give their corresponding mathematical expressions.  
%
% \subsubsection{Dark matter galactic halo distributions}
%
Each profile has a set of parameters which are dependent on the particular halo of interest.  Unless otherwise stated $\rho_0$ sets the normalization of the profile and $r_s$ is the scale radius.

Hernquist-Zhao density profile:
\begin{equation}
\rho_g(r)=\frac{\rho_0}{(r/r_s)^\gamma\left(1+(r/r_s)^{(\alpha-\gamma)/\beta}\right)^\beta} .
\label{DLHA_eq:generalnfw}
\end{equation}

The Navarro-Frenk-White (NFW) distribution is a special case of this with $\alpha = 3$, $\beta = 2$, and $\gamma = 1$ \cite{Navarro:1995iw}.
% \begin{equation}
% \rho_g(r)=\frac{\rho_0}{(r/r_s)\left(1+r/r_s\right)^2} .
% \label{DLHA_eq:NFW}
% \end{equation}
The Kravtsov et al. profile can also be obtained from Eq. (\ref{DLHA_eq:generalnfw}) by setting $\alpha = 3$ and $\beta = 2$ \cite{kra}. 
% \begin{equation}
% \rho_g(r)=\frac{\rho_0}{(r/r_s)^\gamma\left(1+(r/r_s)^{(3-\gamma)/2}\right)^2} .
% \label{DLHA_eq:Kravtsov}
% \end{equation}

Modified isothermal profile \cite{isothermal}:
\begin{equation}
\rho_g(r)=\frac{\rho_0}{1+(r/r_s)^2} .
\label{DLHA_eq:isotherm}
\end{equation}

Einasto profile \cite{eina}:
\begin{equation}
 \rho_g(r) = 
 \rho_0 \exp\left( -\frac{2}{\alpha} \left( \frac{r^\alpha}{r_s^\alpha} -1 \right) \right) .
\label{DLHA_eq:Einasto}
\end{equation}

Moore et al. profile \cite{mooreetal}:
\begin{equation}
\rho_g(r)=\frac{\rho_0}{(r/r_s)^{3/2}\left(1+r/r_s\right)^{3/2}} .
\label{DLHA_eq:Moore}
\end{equation}

Burkert profile \cite{burkert}:
\begin{equation}
\rho_g(r)=\frac{\rho_0 r_s^3}{(r+r_s)(r^2+r_s^2)} .
\label{DLHA_eq:Burkert}
\end{equation}


\subsubsection{The structure of BLOCK DMSPADIST}

Since dark matter energy density distributions can be specified for various astrophysical objects, such as external galaxies and halo substructures, more than one \verb|BLOCK DMSPADIST| may appear in a DLHA file.  In this case multiple \verb|BLOCK DMSPADIST| blocks carry an index differentiating them from each other.  
For example
\begin{quote}
\verb|BLOCK DMSPADIST_MilkyWay| \\
\verb|... | \\
\verb|    | \\
\verb|BLOCK DMSPADIST_NGC6388| \\
\verb|... | \\
\verb|    | \\
\verb|BLOCK DMSPADIST_VIRGOHI21| \\
\verb|... | 
\end{quote}

The shape of the dark matter density distribution $\rho_g(r)$, featured in Eqs. (\ref{DLHA_source}) and (\ref{DLHA_flux}), is defined by the function
\begin{quote}
\verb|FUNCTION rho_g type=<type> args=<number of arguments>|\\
\verb| ...                                                 |\\
\verb|END_FUNCTION                                         |
\end{quote}
The galactic halo profile is assumed to be given in units of GeV/$c^2$/cm$^3$, and its radial argument in units of kpc.

Identifiers for numerical parameters featured in this function for predefined density profiles are the following: 
\begin{description}
\item[1] the normalization of the profile $\rho_0$ in units of GeV/$c^2$/cm$^3$,
\item[2] the scale radius $r_s$ in units of kpc,
\item[3] $\alpha$ for the Hernquist-Zhao and Einasto profiles,
\item[4] $\beta$ for the Hernquist-Zhao  profile,
\item[5] $\gamma$ for the Hernquist-Zhao  profile,
\item[6] the value of the dark matter mass density near Earth $\rho_\oplus$ in units of GeV/$c^2$/cm$^3$, 
\item[7] our Galactocentric distance, $R_0$ in units of kpc. 
\end{description}
These variables can also be used within predefined functions as \verb|PARAMETERS| with the following names
\begin{description}
\item[] \verb|rho_0    |: $\rho_0$      ,
\item[] \verb|r_s      |: $r_s$         ,
\item[] \verb|alpha    |: $\alpha$      ,
\item[] \verb|beta     |: $\beta$       ,
\item[] \verb|gamma    |: $\gamma$      ,
\item[] \verb|rho_oplus|: $\rho_\oplus$ ,
\item[] \verb|R_0      |: $R_0$         . 
\end{description}

Predefined values of the most commonly used halo distributions are
\begin{description}
\item[] \verb|DLHA rho_g 1|: Hernquist-Zhao profile, as given in Eq. (\ref{DLHA_eq:generalnfw}),
\item[] \verb|DLHA rho_g 2|: NFW density profile, defined by Eq. (\ref{DLHA_eq:generalnfw}),
\item[] \verb|DLHA rho_g 3|: Kravtsov et al. profile, defined by Eq. (\ref{DLHA_eq:generalnfw}),
\item[] \verb|DLHA rho_g 4|: Modified isothermal profile, as given in Eq. (\ref{DLHA_eq:isotherm}),
\item[] \verb|DLHA rho_g 5|: Einasto profile, as given in Eq. (\ref{DLHA_eq:Einasto}),
\item[] \verb|DLHA rho_g 6|: Moore et al. profile, as given in Eq. (\ref{DLHA_eq:Moore}),
\item[] \verb|DLHA rho_g 7|: Burkert profile, as given in Eq. (\ref{DLHA_eq:Burkert}).
\end{description}

For example, in a DLHA file a predefined NFW profile can be referred to as
\begin{quote}
\verb|# identifier(s)  parameter value    comment               | \\
\verb|  1              0.385954823E+00  # rho_0 [GeV/cm^3]      | \\
\verb|  2              2.000000000E+01  # scale radius r_s [kpc]| \\
\verb|  3              3.000000000E+00  # parameter $\alpha$    | \\
\verb|  4              2.000000000E+00  # parameter $\beta$     | \\
\verb|  5              1.000000000E+00  # parameter $\gamma$    | \\
\verb|# --------------------------------------------------------| \\
\verb|FUNCTION rho_g type=P args=1                              | \\
\verb| DLHA rho_g 2 # NFW profile                               | \\
\verb|END_FUNCTION                                              |
\end{quote}
Alternatively, the following definition can also be used:
\begin{quote}
\verb|FUNCTION rho_g type=P args=1| \\
\verb| DLHA rho_g 2 # NFW profile | \\
\verb| PARAMETERS                 | \\
\verb|  rho_0 = 0.385954823E+00   | \\
\verb|  r_s   = 2.000000000E+01   | \\
\verb|  alpha = 3.000000000E+00   | \\
\verb|  beta  = 2.000000000E+00   | \\
\verb|  gamma = 1.000000000E+00   | \\
\verb| END_PARAMETERS             | \\
\verb|END_FUNCTION                |
\end{quote}

A sample \verb|BLOCK DMSPADIST| for a predefined Kravtsov et al. halo profile and a tabulated velocity distribution may appear as follows:
\begin{quote}
\verb|BLOCK DMSPADIST_MilkyWay| \\
\verb|# identifier(s)  parameter value    comment              | \\
\verb|  1              0.385954823E+00  # rho_0 [GeV/cm^3]     | \\
\verb|  2              2.000000000E+01  # r_s [kpc]            | \\
\verb|  3              3.000000000E+00  # alpha                | \\
\verb|  4              2.000000000E+00  # beta                 | \\
\verb|  6              0.385954821E+00  # rho_Earth [GeV/cm^3] | \\
\verb|# -------------------------------------------------------| \\
\verb|FUNCTION rho_g type=P args=1      # predefined Kravtsov  | \\
\verb| DLHA rho_g 3                                            | \\
\verb|END_FUNCTION                                             | 
\end{quote}


\subsection{BLOCK DMVELDIST}

This block contains information about velocity distributions of dark matter needed in direct and indirect detection calculations.  
% 
% \indent \indent Truncated Maxwellian velocity distribution:
% \begin{equation}
% \widetilde{f}(\mathbf{v}) = \frac{1}{N_{esc}} (\pi v_0^2)^{-3/2} \, e^{-\mathbf{v}^2\/v_0^2} \, \theta(v_{esc}-v) \, . 
% \label{DLHA_eq:TMaxwellian}
% \end{equation}
% \indent Subtracted Maxwellian distribution:
% \begin{equation}
% \widetilde{f}(\mathbf{v}) = \frac{1}{N_{esc}} (\pi v_0^2)^{-3/2} \, \left[ e^{-\mathbf{v}^2\/v_0^2} - e^{-v_{esc}^2/v_0^2}\right] \, \theta(v_{esc}-v) \, . 
% \label{DLHA_eq:SMaxwellian}
% \end{equation}
% where $v_{esc}$ is the escape speed, a finite cutoff in the distribution expected due to high-speed dark matter being able to escape from the galaxy's gravitational potential, and $N_{esc}$ is a normalization factor.  In the large $v_{esc}$ limit, both distributions reduce to the traditional Maxwell distribution.  For an observer moving at a velocity $\mathbf{v}_{obs}$ relative to the rest frame of the dark matter distribution,
% \begin{equation}
% f(\mathbf{v}) = \widetilde{f}(\mathbf{v}_{obs} + \mathbf{v})
% \label{DLHA_eq:shiftedfv}
% \end{equation}
% is the observed velocity distribution.  In the case of direct detection (Eq.~\ref{DLHA_eq:etaEt}), the relative velocity is given by $\mathbf{v}_{obs} = \mathbf{v}_{\odot} + \mathbf{V}_{\oplus}$, where $\mathbf{v}_{\odot}$ is the the velocity of the Sun relative to the rest frame of the dark matter and $\mathbf{V}_{\oplus}$ is the (time-varying) velocity of the Earth relative to the Sun.
%
% \subsubsection{Dark matter velocity distributions}
%
Velocities are given using the Galactic coordinate system, described
below. 
%
Before discussing the block attributes, we specify common velocity distributions.
%
%Commonly-used speed distributions are the following.
%

Truncated Maxwellian distribution:
\begin{equation}
\widetilde{f}(v) = \frac{4\pi v^2}{N_{esc}  (2\pi \sigma_v^2)^{3/2}} \, e^{-v^2/2\sigma_v^2} \, \theta(v_{esc}-v) \, 
\label{DLHA_eq:TMaxwellian}
\end{equation}
\indent\indent where $N_{esc} = \mathop{\rm erf}(z)-2ze^{-z^2}/\sqrt{\pi}$ with $z=v_{esc}/\sqrt{2}\sigma_v$.

\indent Subtracted Maxwellian distribution:
\begin{equation}
\widetilde{f}(v) = \frac{4\pi v^2}{N_{esc} (2\pi \sigma_v^2)^{3/2}} \,  \left[ e^{-v^2/2\sigma_v^2} - e^{-v_{esc}^2/2\sigma_v^2}\right] \, \theta(v_{esc}-v) \, 
\label{DLHA_eq:SMaxwellian}
\end{equation}
\indent\indent where $N_{esc} = \mathop{\rm erf}(z)-2z(1+2z^2/3)e^{-z^2}/\sqrt{\pi}$ with $z=v_{esc}/\sqrt{2}\sigma_v$.

Here $N_{esc}$ is a normalization factor such that $\int_0^\infty \widetilde{f}(v) \,dv=1$, $v_{esc}$ is the escape speed (a finite cutoff in the distribution expected due to high-speed dark matter being able to escape from the object's gravitational potential), and $\sigma_v$ is a velocity dispersion parameter.  In the large $v_{esc}$ limit, the distributions in Eqs. (\ref{DLHA_eq:TMaxwellian}) and (\ref{DLHA_eq:SMaxwellian}) reduce to the same Maxwellian distribution with most probable speed $v_0=\sqrt{2}\, \sigma_v$, 1D velocity dispersion $\sigma_v$, and 3D velocity dispersion $\sqrt{3} \, \sigma_v$.

% \indent Truncated Gaussian distribution:
% \begin{equation}
% f_{\rm Gal}(\mathbf{v}) = \frac{1}{N_{esc}  (2\pi \sigma_v^2)^{3/2}} \, e^{-\mathbf{v}^2/2\sigma_v^2} \, \theta(v_{esc}-v) \, . 
% \label{DLHA_eq:TGaussian}
% \end{equation}
% \indent Subtracted Gaussian  distribution:
% \begin{equation}
% f_{\rm Gal}(\mathbf{v}) = \frac{1}{N_{esc} (2\pi \sigma_v^2)^{3/2}} \,  \left[ e^{-\mathbf{v}^2/2\sigma_v^2} - e^{-v_{esc}^2/2\sigma_v^2}\right] \, \theta(v_{esc}-v) \, . 
% \label{DLHA_eq:SGaussian}
% \end{equation}
% where $v_{esc}$ is the escape speed (a finite cutoff in the distribution expected due to high-speed dark matter being able to escape from the galaxy's gravitational potential), $\sigma_v$ is a velocity dispersion parameter (which becomes the 1D velocity dispersion in the limit of large $v_{esc}$; the 3D velocity dispersion is $\sqrt{3}\sigma_v$), and $N_{esc}$ is a normalization factor (such that $\int f_{\rm Gal}({\bf v}) \,d^3 v=1$). In the large $v_{esc}$ limit, both distributions reduce to the same (traditional) Maxwell distribution.  For an observer moving at a Galactocentric velocity $\mathbf{v}_{\rm obs}$,
% \begin{equation}
% f(\mathbf{v}) = f_{\rm Gal}(\mathbf{v}_{\rm obs} + \mathbf{v})
% \label{DLHA_eq:shiftedfv}
% \end{equation}
% is the velocity distribution in the frame of the observer.  In the case of direct detection (Eq.~\ref{DLHA_eq:etaEt}), the Galactocentric velocity of the observer is given by $\mathbf{v}_{\rm obs} = \mathbf{v}_{\odot} + \mathbf{v}_{\oplus\rm rev}+{\bf v}_{\oplus\rm rot}$, where  $\mathbf{v}_{\odot}$ is the  velocity of the Sun relative to the Galactic rest frame, $\mathbf{v}_{\oplus\rm rev}$ is the (time-varying) velocity of the Earth revolution relative to the Sun, and ${\bf v}_{\oplus\rm rot}$ is the (time-varying) velocity of the Earth rotation at the location of the detector relative to the Earth rest frame. 
% For reference, speed distributions $\widetilde{f}(v)$ in the observer's rest frame are related to $f(\mathbf{v})$ through
% an integral over the velocity directions
% \begin{equation}
% \widetilde{f}(v) = \int f(\mathbf{v}) |\mathbf{v}|^2 d\Omega_v .
% \end{equation}


% \subsubsection{Local dark matter velocity distributions}

For direct detection, one needs the local dark matter velocity distribution multiplied by the local dark matter density.  The latter is given in the \verb|DMSPADIST| block.  Here we address the velocity distribution.  Notice that if the local dark matter density is split into different velocity groups, each group has its own block 
%s \verb|DMDENSITY| and 
\verb|DMVELDIST|.
%, which must be coordinated.

For direct detection (Eq.~\ref{DLHA_eq:etaEt}), one needs the dark matter velocity distribution $f(\mathbf{v},t)$ in the rest frame of the detector. Its time dependence is expected to arise from the motion of the Earth around the Sun and about its axis. Since these motions are known, $f(\mathbf{v},t)$ can be obtained through a Galilean transformation from the Sun's rest frame,
\begin{equation}
f(\mathbf{v},t) = f_{\odot}(\mathbf{v}_{\rm det} + \mathbf{v}).
\label{DLHA_eq:shiftedfv}
\end{equation}
Here  $f_{\odot}(\mathbf{u})$ is the dark matter heliocentric velocity distribution, which is a function of the dark matter velocity $\mathbf{u}$ with respect to the Sun (heliocentric velocity), and $\mathbf{v}_{\rm det}$ is the detector velocity with respect to the Sun given by $ \mathbf{v}_{\rm det} = \mathbf{v}_{\oplus\rm rev}+{\bf v}_{\oplus\rm rot}$, where $\mathbf{v}_{\oplus\rm rev}$ is the (time-varying) velocity of the Earth revolution relative to the Sun, and ${\bf v}_{\oplus\rm rot}$ is the (time-varying) velocity of the Earth rotation at the location of the detector relative to the Earth's rest frame. 

Enough information should be given in \verb|BLOCK DMVELDIST| to
reconstruct the three-dimensional heliocentric velocity distribution 
$f_{\odot}(\mathbf{u})$. This is achieved by giving the velocity
distribution function $\widetilde{f}(\mathbf{v})$ in some specified but
otherwise arbitrary reference frame $\widetilde{S}$, together with the
velocity $\widetilde{\mathbf{v}}$ of $\widetilde{S}$ with respect to the
Sun, i.e.
\begin{equation}
f_{\odot}(\mathbf{u}) = \widetilde{f}(\mathbf{u}-\widetilde{\mathbf{v}}) .
\label{DLHA_eq:shiftedfvS}
\end{equation}
Though not required (except in a case described below), the reference
frame $\widetilde{S}$ will typically be chosen such that the form for
$\widetilde{f}(\mathbf{v})$ becomes more simplified, as occurs in the
rest frame of an isotropic velocity distribution.


\subsubsection{The structure of BLOCK DMVELDIST}

Similarly to the energy density, dark matter velocity distributions can be specified for various astrophysical objects: galaxies, halo substructures, or even different local velocity components.  Thus, more than one \verb|BLOCK DMVELDIST| may appear in a DLHA file.  Multiple \verb|BLOCK DMVELDIST| blocks are differentiated from each other as
\begin{quote}
\verb|BLOCK DMVELDIST_MilkyWay| \\
\verb|... | \\
\verb|    | \\
\verb|BLOCK DMVELDIST_NGC6388| \\
\verb|... | \\
\verb|    | \\
\verb|BLOCK DMVELDIST_VIRGOHI21| \\
\verb|... | 
\end{quote}

% The dark matter galactic velocity distribution $\widetilde{f}(v)$, and parameters therein, are featured in Eqs. (\ref{DLHA_eq:etaEt}) and (\ref{DLHA_eq:shiftedfv}).  The shape of the distribution is given by the function 

The function $\widetilde{f}(\mathbf{v})$ is given by
\begin{quote}
\verb|FUNCTION fv_g type=<type> args=<number of arguments>|\\
\verb| ...                                                |\\
\verb|END_FUNCTION                                        |
\end{quote}
% where an isotropic distribution, $\widetilde{f}(v\!=\!|\mathbf{v}|) = \widetilde{f}(\mathbf{v})$, is assumed for single-argument functions, including the pre-defined distributions. The unit of the velocity distribution is assumed to be s/km.

% The units of the velocity distribution are assumed to be $({\rm km/s})^{-n}$, where $n$ is the number of arguments. The function $\verb|fv_g|$ must be normalized so that $\int \verb|fv_g| \, d^nv=1$. Notice the following peculiarity: if only one argument is given, it is assumed that this argument is positive and represents the speed $v=|{\bf v}|$, i.e.\ the magnitude of the velocity; in this case, the normalization is $\int_0^\infty \verb|fv_g|(v) \, dv = 1$, and $\verb|fv_g|$ is in fact the speed distribution. 

The units of $\widetilde{f}(\mathbf{v})$  are assumed to be $({\rm km/s})^{-n}$, where $n$ is the number of arguments. 
%
Two different conventions apply depending on whether $\widetilde{f}$ is presented as a function of one or more than one variable, i.e., depending on whether $n=1$ or $n>1$. 
%
If $n>1$, $\widetilde{f}(\mathbf{v})$  must be normalized so that $\int \widetilde{f}(\mathbf{v})  \, d^nv=1$. 
%
%If $n=1$, i.e., if $\widetilde{f}(v)$ is a function of only one variable $v$, it is assumed that this variable is positive and represents the speed $v=|{\bf v}|$, i.e., the magnitude of the velocity. In this case, the normalization is $\int_0^\infty \widetilde{f}(v) \, dv = 1$, and $\widetilde{f}(v)$  is the speed distribution. It is also assumed that the frame $\widetilde{S}$ coincides with the rest frame of the distribution (i.e.,  the frame in which the average velocity is zero) and that the velocity distribution is isotropic in the frame $\widetilde{S}$.
%
If $n=1$, i.e., if $\widetilde{f}(v)$ is a function of only one variable $v$, it is assumed that (i)  the frame $\widetilde{S}$ coincides with the frame in which the average velocity is zero, (ii)  in this frame the velocity distribution is isotropic, (iii)  the variable $v=|{\bf v}|$ is positive and represents the speed, i.e., the magnitude of the velocity, in the frame $\widetilde{S}$, and (iv)  $\widetilde{f}(v)$ is the speed distribution in $\widetilde{S}$, normalized so that $\int_0^\infty \widetilde{f}(v) \, dv = 1$.

Predefined values of the most commonly used velocity distributions are
\begin{description}
\item[] \verb|DLHA fv_g 1|: Truncated Maxwellian distribution, as given in Eq. (\ref{DLHA_eq:TMaxwellian}),
\item[] \verb|DLHA fv_g 2|: Subtracted Maxwellian distribution, as given in Eq. (\ref{DLHA_eq:SMaxwellian}).
\end{description}
Names of numerical \verb|PARAMETERS| featured in these predefined functions are the following: 
\begin{description}
\item[] \verb|sigma_v|: the velocity dispersion parameter $\sigma_v$ in units of km/s,
\item[] \verb|v_esc  |: the escape speed $v_{esc}$ in units of km/s.
\end{description}
These parameters can also be given before the \verb|FUNCTION| definition with the following identifiers
\begin{description}
\item[1] $\sigma_v$ in units of km/s,
\item[2] $v_{esc}$ in units of km/s.
\end{description}

% Identifiers for numerical parameters featured in this function for predefined density profiles are the following: 
% \begin{description}
% \item[1] the parameter $v_0$ in units of km/s,
% \item[2] the escape velocity $v_{esc}$ in units of km/s,
% \end{description}
% along with the velocity of the of the Sun $\mathbf{v}_{\odot}$ (in galactic coordinates) relative to the rest frame of the dark matter distribution, in units of km/s
% \begin{description}
% \item[11] the Sun's relative velocity in the direction of the galactic center ($\mathbf{v}_{\odot,x}$),
% \item[12] the Sun's relative velocity in the direction of the local disk rotation ($\mathbf{v}_{\odot,y}$), and
% \item[13] the Sun's relative velocity in the direction of the galaxy's axis of rotation ($\mathbf{v}_{\odot,z}$).
% \end{description}

% The function \verb|fv_g| must be accompanied by the specification of the Sun's velocity with respect to the Galactic rest frame.  This information is encoded by by the following parameters identified as
% \begin{description}
% \item[11] the Sun's Galactocentric velocity $\mathbf{v}_{{\rm DM},x}$ in the direction of the Galactic Center $(l=0,b=0)$  in units of km/s,
% \item[12] the Sun's Galactocentric  velocity $\mathbf{v}_{{\rm DM},y}$ in the direction of the Galactic rotation $(l=90^\circ,b=0)$  in units of km/s, 
% \item[13] the Sun's Galactocentric  velocity $\mathbf{v}_{{\rm DM},z}$ in the direction of the Galaxy's axis of rotation (North Galactic Pole $b=90^\circ$) in units of km/s.
% \end{description}
% Notice that $xyz$ is a right-handed coordinate system.

% Predefined values of the most commonly used speed distributions are
% \begin{description}
% \item[] \verb|DLHA fv_g 1|: Truncated Maxwellian distribution, as given in Eq. (\ref{DLHA_eq:TMaxwellian}),
% \item[] \verb|DLHA fv_g 2|: Subtracted Maxwellian distribution, as given in Eq. (\ref{DLHA_eq:SMaxwellian}).
% \end{description}
% Variable names for use as \verb|PARAMETERS| for the predefined speed distributions are
% \begin{description}
% \item[] \verb|N_esc|: $N_{esc}$ normalization factor                ,
% \item[] \verb|v_esc|: $v_{esc}$ the escape speed in units of km/s ,
% \item[] \verb|v_0  |: $v_0$ in units of km/s                      .
% \end{description}

The  $\widetilde{S}$ heliocentric velocity vector $\widetilde{\mathbf{v}}$ is to be given by its components in the Galactic reference frame. The Galactic reference frame is a right-handed Cartesian coordinate system $xyz$ with $x$ axis in the direction of the Galactic Center $(l=0,b=0)$,  $y$ axis  in the direction of the Galactic rotation $(l=90^\circ,b=0)$, and  $z$ axis in the direction of the Galaxy's axis of rotation (North Galactic Pole $b=90^\circ$). The  velocity $\widetilde{\mathbf{v}}$ is thus given by the following parameters.
\begin{description}
\item[11] the  $x$ component $\widetilde{v}_{x}$ of the $\widetilde{S}$ heliocentric velocity in units of km/s,
\item[12] the  $y$ component $\widetilde{v}_{y}$ of the $\widetilde{S}$ heliocentric velocity in units of km/s,
\item[13] the  $z$ component $\widetilde{v}_{z}$ of the $\widetilde{S}$ heliocentric velocity in units of km/s.
\end{description}

A sample \verb|BLOCK DMVELDIST| for a predefined truncated Maxwellian distribution may appear as follows:
\begin{quote}
\verb|BLOCK DMVELDIST_standard_dark_halo| \\
\verb|# identifier(s)  parameter value    comment        | \\
\verb|  1              1.555634918E+02  # sigma_v [km/s] | \\
\verb|  2              6.500000000E+02  # v_esc [km/s]   | \\
\verb|  11             0.000000000E+00  # vframe_x [km/s]| \\
\verb|  12            -2.200000000E+02  # vframe_y [km/s]| \\
\verb|  13             0.000000000E+00  # vframe_z [km/s]| \\
\verb|# -------------------------------------------------| \\
\verb|FUNCTION fv_g type=P args=1 # predefined v distr.  | \\
\verb| DLHA fv_g 1                # truncated Maxwellian | 
\end{quote}

A sample \verb|BLOCK DMVELDIST| for a tabulated velocity distribution may appear as follows:
\begin{quote}
\verb|BLOCK DMVELDIST_MilkyWay| \\
\verb|  1              1.555634918E+02  # sigma_v [km/s]| \\
\verb|  2              6.500000000E+02  # v_esc [km/s]  | \\
\verb|  11             0.000000000E+00  # vframe_x [km/s]| \\
\verb|  12            -2.200000000E+02  # vframe_y [km/s]| \\
\verb|  13             0.000000000E+00  # vframe_z [km/s]| \\
\verb|FUNCTION fv_g type=T args=1       # fv_g tabulated| \\
\verb| v: 57                                            | \\
\verb|# v [km/s]        f(v) [s/km]                     | \\
\verb|  5.00000000E+01  1.20000000E-01                  | \\
\verb|  7.50000000E+01  2.40000000E-01                  | \\
\verb|  1.00000000E+02  4.90000000E-01                  | \\
\verb|  ...                                             | \\
\verb|END_FUNCTION                                      | 
\end{quote}


\subsection{BLOCK FORMFACTS}

% The structure of the \verb|FORMFACTS| block is as follows.

The nucleon form factors $f^N$ 
%, $f_{V_q}^N$ 
and $\Delta_q^N$ 
%and $\delta_q^N$ 
as introduced in Eqs. (\ref{DLHA_nuc1form}) and (\ref{DLHA_nuc2form}).  
% and defined by Eqs. (\ref{DLHA_eq:f_q^N}) and (\ref{DLHA_eq:f_Q^N}).
The corresponding identifiers are:
\begin{description}
\item[1~~~~~1]      up quark scalar form factor  for the proton $f_u^p$,
\item[1~~~~~2]    down quark scalar form factor  for the proton $f_d^p$,
\item[1~~~~~3] strange quark scalar form factor  for the proton $f_s^p$,
\item[1~~~~~4]   heavy quark scalar form factors for the proton $f_Q^p$,
\item[2~~~~~1]      up quark scalar form factor  for the neutron $f_u^n$,
\item[2~~~~~2]    down quark scalar form factor  for the neutron $f_d^n$,
\item[2~~~~~3] strange quark scalar form factor  for the neutron $f_s^n$,
\item[2~~~~~4]   heavy quark scalar form factors for the neutron $f_Q^n$,
\item[3~~~~~1]      up quark vector form factor for the proton $f_{V_u}^p$,
\item[3~~~~~2]    down quark vector form factor for the proton $f_{V_d}^p$,
\item[4~~~~~1]      up quark vector form factor for the neutron $f_{V_u}^n$,
\item[4~~~~~2]    down quark vector form factor for the neutron $f_{V_d}^n$,
\item[5~~~~~1]      up quark axial-vector form factor for the proton $\Delta_u^p$,
\item[5~~~~~2]    down quark axial-vector form factor for the proton $\Delta_d^p$,
\item[5~~~~~3] strange quark axial-vector form factor for the proton $\Delta_s^p$,
\item[6~~~~~1]      up quark axial-vector form factor for the neutron $\Delta_u^n$,
\item[6~~~~~2]    down quark axial-vector form factor for the neutron $\Delta_d^n$,
\item[6~~~~~3] strange quark axial-vector form factor for the neutron $\Delta_s^n$,
\item[7~~~~~1]      up quark $\sigma_{\mu\nu}$ form factor for the proton $\delta_u^p$,
\item[7~~~~~2]    down quark $\sigma_{\mu\nu}$ form factor for the proton $\delta_d^p$,
\item[7~~~~~3] strange quark $\sigma_{\mu\nu}$ form factor for the proton $\delta_s^p$,
\item[8~~~~~1]      up quark $\sigma_{\mu\nu}$ form factor for the neutron $\delta_u^n$,
\item[8~~~~~2]    down quark $\sigma_{\mu\nu}$ form factor for the neutron $\delta_d^n$,
\item[8~~~~~3] strange quark $\sigma_{\mu\nu}$ form factor for the neutron $\delta_s^n$.
\end{description}

A specific \verb|FORMFACTS| block for a spin-zero dark matter candidate may look like this:
\begin{quote}
\verb|BLOCK FORMFACTS| \\
\verb|# identifier(s)  value              comment| \\
\verb|  1     1        2.30000000E-02  #  f_u^p  | \\
\verb|  1     2        3.30000000E-02  #  f_d^p  | \\
\verb|  1     3        2.60000000E-01  #  f_s^p  | \\
\verb|  2     1        1.80000000E-02  #  f_u^n  | \\
\verb|  2     2        4.20000000E-02  #  f_d^n  | \\
\verb|  2     3        2.60000000E-01  #  f_s^n  | 
\end{quote}


\subsection{BLOCK STRUCTFUN}

The nuclear structure functions $F_A(E)$ and $S_{ij}(E)$ are defined in Eqs. (\ref{DLHA_sigmanSI}) and (\ref{DLHA_eq:F_A}) for SI interactions, and Eqs. (\ref{DLHA_sigmanSD}) and (\ref{DLHA_Ss1}-\ref{DLHA_Ss3}) for SD interactions.  They may be given in the form of functions of the recoil energy $E$, such as
\begin{quote}
\verb|FUNCTION F_A type=<type> args=1| \\
\verb|...| \\
\verb|END_FUNCTION| 
\end{quote}
and 
\begin{quote}
\verb|FUNCTION S_ij type=<type> args=1| \\
\verb|...| \\
\verb|END_FUNCTION| 
\end{quote}
The nuclear form factors should carry no unit, while the transferred energy $E$  is in units of keV. For reference, the transferred energy $E$ is related to the momentum transfer $q$ through $E=q^2/(2m_A)$, where $m_A$ is the nucleus mass.

The following common parametrizations of $F_A$ are used.  The exponential form factor defined as
\begin{equation}
 F_A(q r_n) = e^{-\alpha (q r_n)^2 / 2} ~~~ {\rm with} ~~~ 
 r_n = a_n A^{1/3} + b_n .
\label{DLHA_eq:F_A_exp}
\end{equation} 
The names of the corresponding parameters that appear in the \verb|PARAMETERS| list are
\begin{description}
\item[] \verb|A    |: $A$ the mass number of the nucleus   ,
\item[] \verb|Z    |: $Z$ the atomic number of the nucleus ,
\item[] \verb|alpha|: $\alpha$                             ,
\item[] \verb|a_n  |: $a_n$ in units of fm                 ,
\item[] \verb|b_n  |: $b_n$ in units of fm                 .
\end{description}

The Helm form factor
\begin{equation}
 F_A(q r_n) = 3 \frac{j_1(q r_n)}{q r_n} e^{-(q s)^2/2} ~~~ {\rm with} ~~~
 r_n^2 = c^2 + \frac{7}{3} \pi^2 a^2 - 5 s^2 ~~~ {\rm and} ~~~
 c = c_0 A^{1/3} + c_1 .
\end{equation} 
\label{DLHA_eq:F_A_Helm}
Helm specific names for \verb|PARAMETERS| are
\begin{description}
\item[] \verb|A  |: $A$ the mass number of the nucleus  ,
\item[] \verb|Z  |: $Z$ the atomic number of the nucleus ,
\item[] \verb|c_0|: $c_0$ in units of fm                ,
\item[] \verb|c_1|: $c_1$ in units of fm                ,
\item[] \verb|a  |: $a  $ in units of fm                ,
\item[] \verb|s  |: $s  $ skin thickness in units of fm .
\end{description}

The Fermi distribution
\begin{equation}
 F_A(q) =  \int \rho(r) e^{ i {\mathbf q} \cdot {\mathbf r} } d^3r 
 ~~~ {\rm with} ~~~
 \rho(r) = \rho_0 ( 1 + e^{(r-c)/a} )^{-1} ~~~ {\rm and} ~~~
 c = c_0 A^{1/3} + c_1 .
\label{DLHA_eq:F_A_Fermi}
\end{equation} 
Corresponding names of \verb|PARAMETERS| are
\begin{description}
\item[] \verb|A    |: $A$ the mass number of the nucleus ,
\item[] \verb|Z    |: $Z$ the atomic number of the nucleus ,
\item[] \verb|rho_0|: $\rho_0$                           ,
\item[] \verb|c_0  |: $c_0$ in units of fm               ,
\item[] \verb|c_1  |: $c_1$ in units of fm               ,
\item[] \verb|a    |: $a  $ in units of fm               .
\end{description}

Corresponding to these in DLHA the following predefined $F_A$ form factors can be used:
\begin{description}
\item[] \verb|DLHA F_A 1|: exponential form factor,
\item[] \verb|DLHA F_A 2|: Helm form factor       ,
\item[] \verb|DLHA F_A 3|: Fermi distribution     . 
\end{description}

Examples of specifying these form factors are:
\begin{quote}
\verb|# Exponential form factor:                   |\\
\verb|#   F(q r_n) = e^{-\alpha (q r_n)^2 / 2)}    |\\
\verb|#   r_n = a_n A^{1/3} + b_n                  |\\
\verb|FUNCTION F_A type=P args=1                   |\\
\verb| DLHA F_A 1                # Exponential form factor|\\
\verb| PARAMETERS                                  |\\
\verb|  Z = 11                                     |\\
\verb|  A = 23                                     |\\
\verb|  alpha = 0.20000000E+00   #                 |\\
\verb|  a_n   = 1.15000000E+00   # [fm]            |\\
\verb|  b_n   = 0.39000000E+00   # [fm]            |\\
\verb| END_PARAMETERS                              |\\
\verb|END_FUNCTION                                 |\\
\verb|                                             |\\
\verb|# General spin-independent form factor (Helm):        |\\
\verb|#  F(q r_n) = 3 \frac{j_1(q r_n)}{q r_n} e^{-(qs)^2/2}|\\
\verb|#  r_n^2 = c^2 + \frac{7}{3} \pi^2 a^2 - 5 s^2        |\\
\verb|#  c = c_0 A^{1/3} + c_1                              |\\
\verb|FUNCTION F_A type=P args=1                            |\\
\verb| DLHA F_A 2             # Helm form factor            |\\
\verb| PARAMETERS                                           |\\
\verb|  Z = 11                                     |\\
\verb|  A = 23                                     |\\
\verb|  c_0 =  1.23000000E+00 # [fm]                        |\\
\verb|  c_1 = -0.60000000E+00 # [fm]                        |\\
\verb|  a   =  0.52000000E+00 # [fm]                        |\\
\verb|  s   =  0.90000000E+00 # [fm] skin thickness         |\\
\verb| END_PARAMETERS                                       |\\
\verb|END_FUNCTION                                          |\\
\verb|                                             |\\
\verb|# Fermi distribution:                        |\\
\verb|FUNCTION F_A type=P args=1                   |\\
\verb| DLHA F_A 3             # Fermi distribution |\\
\verb| PARAMETERS                                  |\\
\verb|  Z = 11                                     |\\
\verb|  A = 23                                     |\\
\verb|  c_0 =  1.23000000E+00 # [fm]               |\\
\verb|  c_1 = -0.60000000E+00 # [fm]               |\\
\verb|  a   =  0.52000000E+00 # [fm]               |\\
\verb| END_PARAMETERS                              |\\
\verb|END_FUNCTION                                 |\\
\verb|                                             |\\
\verb|# Sodium (spin-dependent, tabulated)|\\
\verb|FUNCTION S_00 type=T args=1         |\\
\verb| PARAMETERS                         |\\
\verb|  Z = 11                            |\\
\verb|  A = 23                            |\\
\verb| END PARAMETERS                     |\\
\verb| E: 101                             |\\
\verb|#  E             S_00               |\\
\verb|   0.000000E+00  1.000000E+00       |\\
\verb|   1.000000E+00  0.997000E+00       |\\
\verb|   2.000000E+00  0.994000E+00       |\\
\verb|   ...                              |\\
\verb|END FUNCTION                        |
\end{quote}


\subsection{BLOCK DETECTOR\_NUCLEI \label{DLHA_subsec:BDN}}

An alternative way to specify nuclear structure functions $S_{ij}$ is the following.

\begin{quote}
\verb|BLOCK DETECTOR_NUCLEI  # DAMA                                |\\
\verb|# Num  Fraction  A    Z   J    FSD  S00      S01      S11    |\\
\verb|  1    0.153     23   11  1.5  STD  S00Na23  S00Na23  S11Na23|\\
\verb|  2    0.847     127  53  2.5  STD  S00I127  S00I127  S11I127|\\ 
\verb|                                                             |\\
\verb|FUNCTION S00Na23 type=C args=1                                |\\
\verb| libName=libmicromegas.so                                    |\\
\verb| funcName=S00Na23                                            |\\
\verb|END_FUNCTION                                                 |
\end{quote}


\subsection{BLOCK MASS}

Block \verb|MASS| is part of SLHA.  It is used by DLHA to identify the dark matter particle and to specify its mass.  Specific examples of \verb|BLOCK MASS| appear in the introduction of this section.


\subsection{BLOCK QNUMBERS}

This block is part of the Les Houches BSM Generator Accord and defined in Ref. \cite{Alwall:2007mw}.  The block contains information on the spin, self-conjugate nature, the standard $(SU(3)_C$, $SU(2)_W$, $U(1)_Y$) and exotic charges of the particle.  It may be extended to also contain information about quantum numbers corresponding to discrete symmetries such as R-parity, KK-parity, T-parity, Z-parities.  %This extension may be published in the same proceedings as this work. 


\subsection{BLOCK ABUNDANCE}

Abundances of any species of dark matter particles receive their own block, detailing their properties.  The elements of the abundance block are the following. 

\begin{description}
%\item[1] Thermal/non-thermal tag. Possible values are
%\item[~~~~~]0: the DM species is created thermally,
%\item[~~~~~]1: the DM species is created non-thermally,
%\item[~~~~~]2: both thermal and non-thermal mechanisms contribute to the abundance.
\item[1] The freeze-out temperature $T_f$ in units of GeV$/k$ as defined by the condition in Eq. (\ref{DLHA_freezecond}).
\item[2] The chosen value of $\alpha$ for Eq. (\ref{DLHA_freezecond}).
\item[3] The thermal average of the total (co-)annihilation cross section times velocity at freeze-out $\langle\sigma v\rangle(T_f)$ in units of cm$^3$/s as defined in Eq. (\ref{DLHA_sigmav}).
\item[4] The average energy density of the dark matter particle $\Omega_\chi h^2$ due to thermal production, in units of the critical density as given in Eq. (\ref{DLHA_omegachi}).
\item[5] Same as above but for non-thermal production.
\item[6] Percentage contribution to the total cross section by (co-)annihilation channels.  This line will appear multiple times for a cross section with multiple (co-)annihilation channels.  Each line lists the PDG codes of the two final state particles along with the percent at which it contributes to the total (co-)annihilation cross section.
\end{description}

An example of \verb|BLOCK ABUNDANCE| for a thermally produced dark matter candidate with an abundance of $\Omega_\chi h^2 = 0.11018437$ is the following:
\begin{quote}
\verb|BLOCK ABUNDANCE| \\
\verb|# identifier(s)  parameter value    comment                | \\
%\verb|  1              0               #  thermal production     | \\
\verb|  1              5.16392660E+00  #  T_f [GeV/k]            | \\
\verb|  2              1.50000000E+00  #  alpha                  | \\
\verb|  3              3.18452057E-26  #  <sigma v>(T_f) [cm^3/s]| \\
\verb|  4              0.11018437E+00  #  Omega h^2 thermal      | \\
\verb|# annihilation channel contribution to <sigma v>(T_f) [%]  | \\
\verb|# identifier(s) PGD code 1  PGD code 2    %                | \\
\verb|  6             ...         ...           ...              | 
\end{quote}


\subsection{BLOCK EFFCOUPLING}

The entries of this block are the effective dark matter-nucleon couplings as defined in Eqs. (\ref{DLHA_nuc1form}) and (\ref{DLHA_nuc2form}).
% Eq. (\ref{DLHA_eq:L_eff}) and by Table \ref{DLHA_tableop}.
\begin{description}
\item[1] Spin-independent scalar coupling $f_p$ for $s_\chi=1/2$ and the proton,
\item[2] spin-independent scalar coupling $f_n$ for $s_\chi=1/2$ and the neutron,
\item[3] spin-dependent axial-vector coupling $a_p$ for $s_\chi=1/2$ and the proton,
\item[4] spin-dependent axial-vector coupling $a_n$ for $s_\chi=1/2$ and the neutron.
\end{description}
% \begin{description}
% \item[1~~~~~1] Spin-independent even scalar coupling $\lambda_e(0)$,
% \item[1~~~~~2] Spin-independent  odd scalar coupling $\lambda_o(0)$,
% \item[2~~~~~1] Spin-independent even spinor coupling $\lambda_e(1/2)$,
% \item[2~~~~~2] Spin-independent  odd spinor coupling $\lambda_o(1/2)$,
% \item[3~~~~~1] Spin-independent even vector coupling $\lambda_e(1)$,
% \item[3~~~~~2] Spin-independent  odd vector coupling $\lambda_o(1)$,
% %\item[4~~~~~1] Spin-dependent even scalar coupling $\xi_e(0)$,
% %\item[4~~~~~2] Spin-dependent  odd scalar coupling $\xi_o(0)$,
% \item[5~~~~~1] Spin-dependent even spinor coupling $\xi_e(1/2)$,
% \item[5~~~~~2] Spin-dependent  odd spinor coupling $\xi_o(1/2)$,
% \item[6~~~~~1] Spin-dependent even vector coupling $\xi_e(1)$,
% \item[6~~~~~2] Spin-dependent  odd vector coupling $\xi_o(1)$.
% \end{description}

An example for an \verb|EFFCOUPLING| block is the following
\begin{quote}
\verb|BLOCK EFFCOUPLING| \\
\verb|# identifier(s)  parameter value   comment | \\
\verb|  1              0.10000000E+00  # f_p     | \\
\verb|  2              0.10000000E+00  # f_n     | \\
\verb|  3              0.00000000E+00  # a_p     | \\
\verb|  4              0.00000000E+00  # a_n     | 
\end{quote}


\subsection{BLOCK NDMCROSSSECT}

This block contains the spin-independent and spin-dependent nucleon-dark matter elastic scattering cross sections $\sigma^{SI,SD}_N$ found using Eqs. (\ref{DLHA_eq:sigmaSIn}) and (\ref{DLHA_eq:sigmaSDn}).  
\begin{description}
\item[1] The spin-independent  proton-dark matter elastic scattering cross section $\sigma^{SI}_{p}$ in units of cm$^2$, 
\item[2] The spin-independent neutron-dark matter elastic scattering cross section $\sigma^{SI}_{n}$ in units of cm$^2$, 
\item[3] The   spin-dependent  proton-dark matter elastic scattering cross section $\sigma^{SD}_{p}$ in units of cm$^2$, 
\item[4] The   spin-dependent neutron-dark matter elastic scattering cross section $\sigma^{SD}_{n}$ in units of cm$^2$.
\end{description}

An example for a \verb|NDMCROSSSECT| block is the following
\begin{quote}
\verb|BLOCK NDMCROSSSECT                                      | \\
\verb|# identifier(s)  parameter value   comment              | \\
\verb|  1              0.12345678E-41  # \sigma^{SI}_p [cm^2] | \\
\verb|  2              0.23456789E-41  # \sigma^{SI}_n [cm^2] | \\
\verb|  3              0.34567890E-41  # \sigma^{SD}_p [cm^2] | \\
\verb|  4              0.45678901E-41  # \sigma^{SD}_n [cm^2] | 
\end{quote}


\subsection{BLOCK ASTROPROPAG}

This and the following blocks describe aspects of the calculations relevant to indirect detection experiments.  \verb|BLOCK ASTROPROPAG| contains parameters and functions of the propagation model for charged cosmic rays.  The parameters in this block are the following:

\begin{description}
\item[1] normalization of the spatial diffusion coefficient $K_0$, in Eq. (\ref{DLHA_eq:KE}),
\item[2] coefficient $\eta$, controlling the dependence on $v$ of the spatial diffusion coefficient in Eq. (\ref{DLHA_eq:KE}),
\item[3] coefficient $\delta$ parameterizing the steepness of the spatial diffusion coefficient in energy, in Eq. (\ref{DLHA_eq:KE}),
\item[4] normalization of the energy loss coefficient $b_0$, in Eq. (\ref{DLHA_eq:bE}),
\item[5] half-height $L$ of the diffusion box in units of kpc,
\item[6] radius $R$ of the diffusion box in units of kpc,
\item[7] re-acceleration parameter $v_A$ in units of km/s, in Eq. (\ref{DLHA_eq:KEE}),
\item[8] Galactic wind parameter $V_C$ in units of km/s, in Eq. (\ref{DLHA_evol_energy}).
%\item[8] magnetic rigidity of the cosmic ray particles $R_0$ in units of $MV$,
\end{description}

A spatial diffusion coefficient with generic energy or rigidity dependence can be defined by the function
\begin{quote}
\verb|FUNCTION SpatDiff type=<type> args=1| \\
\verb| ...                                | \\
\verb|END_FUNCTION                        | 
\end{quote}
A predefined spatial diffusion coefficient dependence is given by Eq. (\ref{DLHA_eq:KE}) and can be indicated by 
\begin{quote}
\verb|FUNCTION SpatDiff type=P args=1| \\
\verb| DLHA SpatDiff 1               | \\
\verb|END_FUNCTION                   | 
\end{quote}

The energy loss coefficient for charged cosmic rays, $b(E)$, can also be defined by 
\begin{quote}
\verb|FUNCTION EnerLoss type=<type> args=1| \\
\verb| ...                                | \\
\verb|END_FUNCTION                        | 
\end{quote}
Its predefined form, given by Eq. (\ref{DLHA_eq:bE}), is indicated as
\begin{quote}
\verb|FUNCTION EnerLoss type=P args=1| \\
\verb| DLHA EnerLoss 1               | \\
\verb|END_FUNCTION                   | 
\end{quote}


\subsection{BLOCK DMCLUMPS}

This is a recommended block and it is only sketched in this version of DLHA.  It should store information regarding dark matter substructures.  It stores the relevant information on the distribution of dark matter substructures inside the Milky Way halo.  
%Here are a few suggestions borrowed from \cite{Lavalle:1900wn} and extracted from my C code.

The inner structures of clumps are similar to
the dark matter distributions inside galactic halos. Whatever the
favourite density profile, code builders need to compute the
annihilation volume $\xi$ as a function of the substructure mass
$M_{\rm cl}$. For this, they need to compute the virial radius
$R_{\rm vir}$ of a given clump. This radius encompasses an average
density
\begin{equation}
 \bar{\rho}(R_{\rm vir})={\displaystyle \frac{M_{\rm cl}}{(4\pi/3)R_{\rm vir}^{3}}}=
 \Delta_{\rm vir} \, \Omega_{\rm M} \, \rho_{C}^{0} ,
\end{equation}
which is $\Delta_{\rm vir}$ times larger than the average cosmological
matter density $\Omega_{\rm M}\rho_{C}^{0}$. The concentration of the
clump is defined as the ratio
\begin{equation}
c_{\rm vir}=\frac{R_{\rm vir}}{r_{-2}} ,
\end{equation}
where the clump density profile has slope $-2$ at radius $r=r_{-2}$.
The mass-to-concentration relation can be parameterized as
\cite{Lavalle:1900wn}
\begin{equation}
 \ln(c_{\rm vir})={\sum_{i=0}^{4}} \, C_{i} \,
 \left\{ \ln \left({\displaystyle \frac{M_{\rm cl}}{M_{\odot}}}\right) \right\}^{i} .
\end{equation}
The output is the annihilation volume $\xi$ as a function of the clump
mass $M_{\rm cl}$.
%
% \begin{quote}
% \verb|//PARAMETERS C_i FOR THE BULLOCK CONCENTRATION           |\\
% \verb|#define COEF_C0 ( 4.34   )                               |\\
% \verb|#define COEF_C1 (-0.0384 )                               |\\
% \verb|#define COEF_C2 (-3.91e-4)                               |\\
% \verb|#define COEF_C3 (-2.2e-6 )                               |\\
% \verb|#define COEF_C4 (-5.5e-7 )                               |\\
% \verb|                                                         |\\
% \verb|double DELTA_vir = 340.                                  |\\
% \verb|double OMEGA_M   = 0.2646                                |\\
% \verb|double rho_0_C   = 139.87 //[SOLAR MASS PER CUBIC PARSEC]|\\
% \verb|                                                         |\\
% \verb|double c_vir(double M)    //[SOLAR_MASS]                 |\\
% \verb|double x = log(M)                                        |\\
% \verb|double c = COEF_C0 + COEF_C1*x + COEF_C2*x*x +           |\\
% \verb|           COEF_C3*x*x*x + COEF_C4*x*x*x*x               |\\
% \verb|return exp(c)                                            |
% \end{quote}

The substructure probability distribution
${\cal D}({\vec r}_S,\xi)$ and the total number of clumps
${\cal N}_{\rm H}$ inside the Milky Way halo can be inferred
from the space and mass distribution function
${d^{4}{\cal N}_{\rm cl}}/{dM_{\rm cl \,}}{d^{3}r_S}$. As an illustration,
we can assume the space distribution of clumps to be isotropic and independent
of the mass spectrum, so that
\begin{equation}
 {\displaystyle \frac{d^{2}{\cal N}_{\rm cl}}{{4{\pi}r_S^{2}dr_S}\,{dM_{\rm cl}}}}=
 {\displaystyle \frac{d{\cal N}_{\rm cl}}{dM_{\rm cl}}} \times {\cal P}(r_S) .
\end{equation}
The mass distribution is taken in general to be the power law
\begin{equation}
 {\displaystyle \frac{d{\cal N}_{\rm cl}}{dM_{\rm cl}}}=
 {\displaystyle \frac{C}{{M_{\rm cl}}^{\alpha}}} ,
%C \times {M_{\rm cl}}^{-\alpha} ,
\end{equation}
where the normalization constant $C$ is calculated by requiring that
the Milky Way halo contains a certain number of substructures
with mass in the range $M_{1}$ to $M_{2}$. In Ref. \cite{Lavalle:1900wn}
for instance, there are 100 clumps between $10^{8}$ and $10^{10} M_{\odot}$.
The mass spectrum of the clumps extends from $M_{\rm inf}$ to $M_{\rm sup}$
which are inputs of the code.
%
% \begin{quote}
% \verb|//PARAMETERS FOR THE CLUMP MASS SPECTRUM                    |\\
% \verb|#define M_INF   = 1.0e-6 // MINIMUM CLUMP MASS [SOLAR_MASS] |\\
% \verb|#define M_SUP   = 1.0e10 // MAXIMUM CLUMP MASS [SOLAR_MASS] |\\
% \verb|#define ALPHA_M = 1.9    // SLOPE OF THE MASS SPECTRUM      |
% \end{quote}

The last input is the space distribution function
${\cal P}(r_S)$. Substructures do not follow the smooth dark matter
distribution. In particular, they are tidally disrupted near the Galactic
center, hence a deficit with respect to the smooth component. As an
illustration, a possible clump distribution is given by the isothermal
profile with core radius $a$
\begin{equation}
 {\cal P}(r_S)={\displaystyle \frac{\kappa}{a(a^2+r_S^2)}} .
\end{equation}
Assuming that the integral of ${\cal P}(r_S)$ from $r=r_{\rm min}$
to $r=r_{\rm max}$ is normalized to unity leads to the constant
\begin{equation}
 \kappa={\displaystyle \frac{{1}/{4{\pi}}}{(u_2 - u_1) - \left\{\arctan(u_2) - \arctan(u_1)\right\}}} ,
\end{equation}
where $u_1={r_{\rm min}}/a$ and $u_2={r_{\rm max}}/a$. This is just an
example. What is actually needed is the input function ${\cal P}(r_S)$.


\subsection{BLOCK ANNIHILATION}

Input block (from various particle physics model codes) that can be used by other codes to calculate the source spectra from annihilation.  This block contains the total annihilation cross section and partial cross sections into different Standard Model final states.

\begin{description}
\item[1] Annihilation cross section times velocity, $\sigma v$ in the $v\rightarrow 0$ limit and in units of cm$^3$/s, followed by a table of annihilation channels with branching fractions with the following columns:
\item[column 1] branching fraction (BR),
\item[column 2] number of annihilation products (NDA),
\item[column 3] PDG code of annihilation product 1,
\item[column 4] PDG code of annihilation product 2.
\end{description}

An example entry could look like
\begin{quote}
\verb|BLOCK ANNIHILATION| \\
%\verb|# identifier(s)  parameter value  comment                 |\\
\verb|      1          3.000000E-26     # sigma v (v->0) [cm^3/s]|\\
\verb|#     BR         NDA   ID1    ID2   comment                |\\
\verb| 9.50000000E-01   2     24    -24 # X X -> W+ W-           |\\
\verb| 4.00000000E-02   2      5     -5 # X X -> b b-bar         |\\
\verb| 1.00000000E-02   2     25     36 # X X -> H_1^0 H_3^0     |
\end{quote}


% \subsection{BLOCK DECAY}
% 
% This block is utilized for dark matter particles that decay to exchange information about its decay properties.  The block contains the total decay rate and partial decay widths into different standard model final states.
% 
% \begin{description}
% \item[1] Decay life time $\tau$ of the dark matter particle in seconds,
% \item[] Table of decay channels with branching fractions
% \item[] Branching fraction (BR), number of decay products (NDD) and PDG codes of annihilation products
% \end{description}
% 
% An example entry could look like
% \begin{quote}
% \verb|BLOCK DECAY| \\
% \verb|      1          1.00000000E26   # decay rate tau [s] | \\
% \verb|#     BR         NDD     ID1      ID2      channel| \\
% \verb| 9.50000000E-01   2       24      -24    # W+ W-  | \\
% \verb| 5.00000000E-02   2        5       -5    # b b-bar| \\
% \end{quote}


\subsection{BLOCK INDIRDETSPECTRUM $i$ $n$}
The spectra of end products (such as positrons, gamma rays, antiprotons, etc.) in the halo (or vacuum) or in an environment (like the Sun) can be calculated from the information given in the blocks above. These spectra can then be used as input for programs that solve the cosmic ray propagation equations, or calculate the fluxes from particular sources in the sky in case of neutral particles. If for example the block \verb|ANNIHILATION| is given at the same time as any of the spectrum blocks given below, the blocks below should take precedence, i.e. the codes should \emph{not} recalculate the spectra if a spectrum is already given in the DLHA file.

The block definition is
\begin{quote}
\verb|BLOCK INDIRDETSPECTRUM i n| 
\end{quote}
where $i$ is the dark matter particle and $n$ is the spectrum type:
\begin{description}
\item[~] $n=1$: positron spectrum in vacuum (or the halo),  
\item[~] $n=2$: antiproton spectrum in vacuum (or the halo),  
\item[~] $n=3$: antideutron spectrum in vacuum (or the halo),  
\item[~] $n=4$: gamma-ray spectrum in vacuum (or the halo),  
\item[~] $n=5$: electron neutrino spectrum in vacuum (or the halo),  
\item[~] $n=6$: muon neutrino spectrum in vacuum (or the halo),  
\item[~] $n=7$: tau neutrino spectrum in vacuum (or the halo),
\item[~] $n=101$: electron neutrino spectrum at the Earth from annihilations in the Sun,
\item[~] $n=102$: electron anti-neutrino spectrum at the Earth from annihilations in the Sun,
\item[~] $n=103$: muon neutrino spectrum at the Earth from annihilations in the Sun,
\item[~] $n=104$: muon anti-neutrino spectrum at the Earth from annihilations in the Sun,
\item[~] $n=105$: tau neutrino spectrum at the Earth from annihilations in the Sun,
\item[~] $n=106$: tau anti-neutrino spectrum at the Earth from annihilations in the Sun,
%\item[~] $n=109$: $\mu^-$ spectrum at the neutrino-nucleon vertex in the Earth from annihilations in the Sun
%\item[~] $n=110$: $\mu^+$ spectrum at the neutrino-nucleon vertex in the Earth from annihilations in the Sun
\item[~] $n=113$: $\mu^-$ spectrum at the Earth (coming from muon neutrino nucleon interactions) from annihilations in the Sun,
\item[~] $n=114$: $\mu^+$ spectrum at the Earth (coming from muon neutrino nucleon interactions) from annihilations in the Sun,
\item[~] $n=201$: electron neutrino spectrum at the Earth from annihilations in the Earth,
\item[~] $n=202$: electron anti-neutrino spectrum at the Earth from annihilations in the Earth,
\item[~] $n=203$: muon neutrino spectrum at the Earth from annihilations in the Earth,
\item[~] $n=204$: muon anti-neutrino spectrum at the Earth from annihilations in the Earth,
\item[~] $n=205$: tau neutrino spectrum at the Earth from annihilations in the Earth,
\item[~] $n=206$: tau anti-neutrino spectrum at the Earth from annihilations in the Earth,
%\item[~] $n=209$: $\mu^-$ spectrum at the neutrino-nucleon vertex in the Earth from annihilations in the Earth
%\item[~] $n=210$: $\mu^+$ spectrum at the neutrino-nucleon vertex in the Earth from annihilations in the Earth
\item[~] $n=213$: $\mu^-$ spectrum at the Earth (coming from muon neutrino nucleon interactions) from annihilations in the Earth,
\item[~] $n=214$: $\mu^+$ spectrum at the Earth (coming from muon neutrino nucleon interactions) from annihilations in the Earth.
\end{description}

The content of the block is the energy distribution $dN/dE$ as the function of the energy $E$:
\begin{quote}
\verb|E        dN/dE | 
\end{quote}
The spectrum $dN/dE$ is the yield per annihilation or decay at that energy (for $n<100$). For spectra from the Sun/Earth ($n>100$), the units of $dN/dE$ are per annihilation per m$^2$.

A typical spectrum block could then look like
\begin{quote}
\verb|BLOCK INDIRDETSPECTRUM 1 4 # gamma-ray spectrum | \\
\verb|#      E        dN/dE                           | \\
\verb|  1.000000E-02  1.000000E-4                     | \\
\verb|  2.000000E-02  1.200000E-4                     | \\
\vdots
\end{quote}


\subsection{DECAY files}

For decaying dark matter particles a standard SLHA decay file can be used to read and write the total decay width of the dark matter particle and its branching ratios into various final states.  In SLHA DECAY entries are possible for decay channels of various particles, including the dark matter particle. They look similar to the \verb|ANNIHILATION| block above. To be able to calculate the complete spectrum from annihilation and decay of a dark matter particle, we also need the partial decay widths (or branching ratios) for other new physics particles, e.g.\ new Higgs bosons. These DECAY structures should follow the same format as in SLHA2. For example, the decay of a Higgs boson might look like
\begin{quote}
\verb|#      PDG       Width                     | \\
\verb|DECAY  1000039   3.287443E+35              | \\
\verb|#      BR        NDD   ID1   ID2    channel| \\
\verb| 9.50000000E-01    2    24   -24  # W+ W-  | \\
\verb| 5.00000000E-02    2     5    -5  # b b-bar| \\
\verb|# -----------------------------------------| \\
\verb|#      PDG       Width                     | \\
\verb|DECAY  35        3.287443E+00              | \\
\verb|#      BR        NDA   ID1   ID2           | \\
\verb|   9.000000E-01    2     5    -5  # b b-bar| \\
\verb|   6.000000E-02    2    24   -24  # W- W+  | \\
\verb|   4.000000E-02    2    23   -23  # Z0 Z0  | 
\end{quote}


\section{OPEN ISSUES}

A partial list of open issues is addressed in various degrees of detail.


\subsection{Cosmology related open issues}

\vspace{3mm}
\begin{itemize}
\item Other generalizations of the standard cosmological equations would be useful. 
\item The standard inflation scenario should be discussed as well.
\item What about a decaying inflaton scenario? Can be presented as an example. 
\end{itemize}


\subsection{Astrophysics related related open issues}

\vspace{3mm}
\begin{itemize}

\item The notation in BLOCK DMSPADIST is confusing.  $\rho_0$ is not the density at $R_0$.  The local density could be called $\rho_\odot$ and the Galactocentric distance $R_\odot$.  It should be clearly stated that either the density is defined with respect to the change of slope (or the interior mass) and then parameters (1, 2, 3, 4, and 5) are needed, or with respect to the local density, which requires parameters (2, 3, 4, 5, 6, and 7). Giving all the parameters is also possible but one needs to check consistency.
\item Ideally, DLHA should be able to accommodate innately anisotropic distributions.
\item A capability for non-spherical and/or clumped halo distributions are also desirable.
\item Notation: should we relabel $\rho_0$ to $\rho_N$ (such as $\rho$ normalisation without any more specification as N-body codes use different definitions) and $\rho_\oplus$ to $\rho_0$ so that $\rho_0 = \rho(r=R_0)$?
\end{itemize}


\subsection{Direct detection related open issues}

\vspace{3mm}
\begin{itemize}
\item For non-self conjugate dark matter particles there is a need to present nuclear cross sections for antiparticles.  This is not possible within the present setup.
\end{itemize}


\subsection{Indirect detection related open issues}

\vspace{3mm}
\begin{itemize}
\item Photon energy losses need to be improved.
\item Solar modulation needs a unified treatment.
\item The quantities that enter the transport equation (for example, the diffusion coefficient, the re-acceleration term, \ldots) should be standardized, in such a way that free parameters are identified. 
\end{itemize}


\subsection{Other open issues}

\vspace{3mm}
\begin{itemize}
\item Kinetic decoupling of DM particles should be discussed. This sets the small-scale cutoff in the spectrum of density perturbations, viz. the mass of smallest dark matter halos, and can have impact on, e.g., the anisotropy spectrum and the 'boost factor' for indirect searches. 
\item Concerning the calculation of the relic density, the importance of the QCD phase transition should be stressed.  This may impact strongly the calculation if the dark matter candidate is light (10 GeV/c$^2$ or so).  The QCD phase transition temperature should be a parameter to put in BLOCK DOFREEDOM.
\item BLOCK DMPDGCODE: suggestion for a new block more explicitly identifying the DM candidate(s). 
\begin{quote}
\verb| 1   0  1000022 # neutralino|\\
\verb| 1   1  1000012 # sneutrino |\\
\verb| ...                        |\\
\verb| 2   0  9999999 # axion     |
\end{quote}
\item With large tables in large parameter scans, I/O may take a long time.
\item In scenarios/models with $Z_3$, $Z_4$ etc., different interactions and processes appear.
\item What about semi-annihilations?
\item Possible conflict between multi-component dark matter and asymmtry/clumping/non-spherical halo.
\end{itemize}


\section*{ACKNOWLEDGEMENTS}

The authors thank Xerxes Tata for useful comments.
%
C.B. and R.L. were funded in part by the Australian Research Council under Project ID DP0877916 and in part by the ARC Centre of Excellence for Particle Physics at the Terascale.
%
C.B. was partially funded by Project of Knowledge Innovation Program (PKIP) of the Chinese Academy of Sciences, Grant No. KJCX2.YW.W10. 
%
G.B., F.B. and A.P. were supported in part by the GDRI-ACPP of CNRS.
%
T.B. acknowledges support from the German Research Foundation (DFG) through 
Emmy Noether grant BR 3954/1-1.
%
D.G.C. is supported by the Ram\'on y Cajal program of the Spanish MICINN. T.D. was supported by the Spanish MICINN�s  Consolider-Ingenio 2010 Programme under grant CPAN CSD2007-00042. D.G.C. and T.D. also thank the support of the Spanish MICINNs Consolider-Ingenio 2010 Programme under grant MultiDark CSD2009-00064, the Spanish MICINN under grant FPA2009-08958, the Community of Madrid under grant HEPHACOS S2009/ESP-1473, and the European Union under the Marie Curie-ITN program PITN-GA-2009-237920.
%
The work of M.C. in supported in part by the French national research agency ANR under contract ANR 2010 BLANC 041301 and by the EU ITN network UNILHC.
%
J.E. and C.S. are grateful to the Swedish Research Council (VR) for
financial support.
%
P.G. was supported by NSF award PHY-1068111 and acknowledges support by the Oskar Klein Centre during his sabbatical visit.
%
S.K. was supported by the US National Science Foundation grants PHY-0854724 and PHY-0955626.
%
A.P. was supported by the Russian foundation for Basic Research, grant RFBR-10-02-01443-a. 
%
S.P. is partly supported by the US Department of Energy, Contract DE-FG02-04ER41268, and by NSF Grant PHY-0757911. 
%
P.S. is supported by the Lorne Trottier Chair in Astrophysics and an Institute for Particle Physics Theory Fellowship.
%
The development of GALPROP is partially supported by NASA grants NNX09AC15G and NNX10AE78G, and by the Max-Planck Society.

\bibliography{DLHA}

\end{document}
