Effect of polarisation and choice of event generator on spectra from dark matter annihilations

TL;DR

This paper examines how the choice of event generator and particle polarisation affect the predicted energy spectra of secondary particles from dark matter annihilations, highlighting the importance of these factors in indirect detection analyses.

Contribution

It provides a comparative analysis of Pythia8 and Herwig7 for secondary particle yields and investigates the impact of polarisation on neutrino flux predictions from dark matter annihilations.

Findings

Differences between event generators are larger for hadronic yields than leptonic.

Polarisation significantly affects leptonic spectra, especially neutrinos from the Sun.

Propagation effects in the Sun reduce the impact of polarisation differences.

Abstract

If indirect detection searches are to be used to discriminate between dark matter particle models, it is crucial to understand the expected energy spectra of secondary particles such as neutrinos, charged antiparticles and gamma-rays emerging from dark matter annihilations in the local Universe. In this work we study the effect that both the choice of event generator and the polarisation of the final state particles can have on these predictions. For a variety of annihilation channels and dark matter masses, we compare yields obtained with Pythia8 and Herwig7 of all of the aforementioned secondary particle species. We investigate how polarised final states can change these results and do an extensive study of how the polarisation can impact the expected flux of neutrinos from dark matter annihilations in the centre of the Sun. We find that differences between the event generators are larger for yields of hadronic end products such as antiprotons, than for leptonic end products. Concerning polarisation, we conversely find the largest differences in the leptonic spectra. The large differences in the leptonic spectra point to the importance of including polarisation effects in searches for neutrinos from dark matter annihilations in the Sun. However, we find that these differences are ultimately somewhat washed out by propagation effects of the neutrinos in the Sun.

Tables

Table 1. Table 1 : The values of the neutrinos oscillation parameters that we use in our WimpSim simulations. These are the best-fit values for the normal mass hierarchy from the NuFit-3.2 global fit [ 51 , 52 ] .

Parameter	Best-fit value
${\theta }_{12}$	$33.62°$
${\theta }_{23}$	$47.2°$
${\theta }_{13}$	$8.54°$
${\delta }_{\mathrm{CP}}$	$234°$
$\Delta m_{21}^2$	$7.40\text{×}{10}^{-5}{\mathrm{eV}}^2$
$\Delta m_{31}^2$	$2.494\text{×}{10}^{-3}{\mathrm{eV}}^2$

Table 2. Table 2 : The number of events (per annihilation) in IceCube for different masses, annihilation channels and polarisation states. The Pythia6 ( P6 ) results are from a regular WimpSim run and the results for the polarised states are for the runs specified in this section. The IceCube number of events are calculated with the effective areas in ref.  [ 54 ] for a livetime of 532 days (corresponds to three years of data taking). The ratio in the last column is the ratio of the two polarisation results (column 5/column 4).

${𝒎}_{\chi }$ [GeV]	Channel	Events per annihilation			Ratio
		P6 (unpol.)	$W_L^{-}{W}_L^{+}$	$W_T^{-}{W}_T^{+}$	(col 5/col 4)
100	$W^{-}{W}^{+}$	$8.99\text{×}{10}^{-25}$	$8.70\text{×}{10}^{-25}$	$9.58\text{×}{10}^{-25}$	1.13
1000		$3.44\text{×}{10}^{-23}$	$4.04\text{×}{10}^{-23}$	$3.46\text{×}{10}^{-23}$	0.86
		P6 (unpol.)	${\tau }_R^{-}{\tau }_L^{+}$	${\tau }_L^{-}{\tau }_R^{+}$
100	${\tau }^{-}{\tau }^{+}$	$2.13\text{×}{10}^{-24}$	$1.66\text{×}{10}^{-24}$	$2.60\text{×}{10}^{-24}$	1.56
1000		$1.16\text{×}{10}^{-22}$	$1.14\text{×}{10}^{-22}$	$1.24\text{×}{10}^{-22}$	1.08

Equations

$$\frac{d\mathcal{Q}}{d{E}_f}=\frac{1}{N_{\chi }}\frac{{\rho }_{\chi }^2}{m_{\chi }^2}\sum _i⟨{\sigma }_{i}v⟩\frac{d{N}_i}{d{E}_f},$$

$${\tau }_L^{-}{\tau }_R^{+},{\tau }_R^{-}{\tau }_L^{+},t_L{\overline{t}}_R,t_R{\overline{t}}_L,b\overline{b},W_L^{-}{W}_L^{+},W_T^{-}{W}_T^{+}$$

$$E_{\pm }=\frac{m_{\chi }}{2}\left(1\pm \sqrt{1-\frac{m_W^2}{m_{\chi }^2}}\right).$$

Full text

We find that differences between the event generators are larger for yields of hadronic end products such as antiprotons, than for leptonic end products. Concerning polarisation, we conversely find the largest differences in the leptonic spectra. The large differences in the leptonic spectra point to the importance of including polarisation effects in searches for neutrinos from dark matter annihilations in the Sun. However, we find that these differences are ultimately somewhat washed out by propagation effects of the neutrinos in the Sun.

1 Introduction

Despite years of investigation, the identity of dark matter remains an open question. One of the most widely studied ideas is that dark matter consists of a hitherto unknown type of particle called a WIMP (Weakly Interacting Massive Particle), which is well motivated because it allows for the production of the correct abundance of dark matter through a simple thermal mechanism in the early universe [1]. The same annihilation processes that are important in setting the WIMP relic abundance also lead to substantial astronomical fluxes of Standard Model (SM) particles from regions of space where the dark matter density is high. This leads to the programme of so-called indirect detection, where one attempts to discover dark matter particles, such as WIMPs, by observing excesses of these particles, typically charged antiparticles, gamma-rays or neutrinos, compared to the astrophysical background.

To be able to predict the signal expected in indirect detection for a particular dark matter model it is essential to be able to accurately predict the energy spectra of the yield particles (by yield particles we here refer to the particles one looks for in indirect detection experiments) expected directly from annihilations for a given dark matter model. To compare with observations, these production yields are translated to an expected flux at Earth by propagating the produced SM particles from the source to the Earth, including necessary effects such as diffusion by magnetic fields in the case of charged particles, or flavour oscillations in the case of neutrinos. Depending on the astrophysical environment where the dark matter annihilation takes place, other effects can also play an important role, for example accounting for the interactions of neutrinos with the solar material is crucial in the searches for neutrinos from dark matter annihilations in the Sun [2, 3]. The characteristics of the production yields are therefore a vital part of the indirect detection programme, and any uncertainty in the production yields will translate into an uncertainty in experimental constraints on dark matter properties.

To determine the production yields, one often uses a model-independent approach where the specific particle physics model containing the dark matter particle (e.g. the neutralino in a supersymmetric model [4]) is not specified, but rather Monte Carlo simulations are done for various two body final states at a range of centre-of-mass energies. By interpolating between the tabulated yields from these simulations, the energy spectra of the yield particles of interest can be determined for annihilations of a dark matter particle of arbitrary mass to a particular final state (it is usually assumed $E_{\rm c.o.m.}=2m_{\chi}$, where $m_{\chi}$ is the dark matter mass). The total injection rate at the source $\mathcal{Q}$ of a particular yield particle $f$ is then given by the relation

[TABLE]

where $\langle\sigma_{i}v\rangle$ is the thermally averaged annihilation cross section for annihilation to final state $i$, $N_{\chi}=2$ for self-conjugate dark matter and 4 otherwise, and $dN_{i}/dE_{f}$ is the corresponding spectra of yield particle $f$ for the that final state, as determined from the aforementioned interpolation. For parameter space scans, using this approach rather than than running Monte Carlo simulations for every model point leads to overwhelming speed increases. However, there are some models in which this procedure can fail to capture important physics. For example, internal bremsstrahlung, where a gauge boson is emitted from one of the internal propagators of the annihilation process, can change significantly the shape of the obtained photon spectrum [5, 6]. Also, as the final state particles are often assumed to be unpolarised in this approach, it can lead to substantial mismodelling of the yields in models in which the final state particle are produced with a particular polarisation [7, 8, 9, 10, 11, 12].

One typically employs general purpose Monte Carlo event generators [13, 14, 15, 16, 17] to simulate the dark matter annihilations and the resulting energy spectra of yield particles in order to handle the complications of parton showering and hadronisation in the case of coloured final state particles. These event generators include the treatment of parton showering and hadronisation and the subsequent decays into stable particles. Apart from coloured radiation also electromagnetic and sometimes electroweak radiation can be included [18]. There are however differences in how these codes treat the simulation of the involved particle physics. In particular, since hadronisation is a non-perturbative process, one is forced to resort to phenomenological modelling for which the codes follow varying approaches. Also the radiation of coloured particles, the parton showering, can be different from code to code. Therefore, it is important to quantify the differences between the different event generators in order to assess their ultimate impact on indirect dark matter searches.

In this paper we look at two of the most common event generators used for showering and hadronisation, Pythia8 [19, 13] and Herwig7 [15, 14], and look at the differences between them for the energy spectra of the various yield particles of interest in indirect detection experiments. We also compare to results obtained previously with Pythia6 [19]. These spectra are implemented in the current version of the DarkSUSY code [20] and are thus useful as a reference. In addition to comparing the yields between event generators, we consider the effect of spin polarisation on the spectra, looking specifically at how this affects the predictions for event rates in neutrino telescopes, where these effects are expected to have a significant impact.

There are a number of previous studies that consider the yields from dark matter annihilations. In refs. [21, 3] in particular, a large number of annihilation channels, dark matter masses and yield particles are considered, with the results publicly available. In this paper, we perform a complementary analysis while also updating refs. [21, 3] with results from new versions of the event generators and including further details on the effects of final state polarisation. The effect of varying the event generator and its parameters on the gamma-ray flux from dark matter annihilations has been studied in refs. [21, 22, 23]. Of these studies, ours is similar in spirit to refs. [21, 22], but we expand upon these works by considering a broader range of yield particles beyond photons, as well as polarised final states. The approach of ref. [23] is different. There the authors focus mainly on the variance in the photon yield that is derived from the uncertainties of the QCD model parameters in Pythia, while we (and the previously mentioned studies) focus on comparing the results of various event generators with their default tunes. Even with this less comprehensive approach, we still find significant uncertainties in the predicted fluxes of all of the yield particles of interest, hopefully motivating detailed studies along the lines of ref. [23] that consider yield particles other than photons, as well as the effects of varying the tune parameters in Herwig as well as Pythia. Finally, the effect of polarisation in gauge boson final states in dark matter annihilations in the Sun has been studied in ref. [24], where a 5-10% effect is found on the second moment of the neutrino fluxes. For neutralino dark matter, the effect of polarisation on gauge boson final states was also studied in ref. [7] where they emphasise the importance of including the polarisation effects. Here, we study polarisation of gauge bosons and also extend the study to fermion polarisations and re-investigate the effect of polarisation in a more complete framework for the simulation of the neutrino propagation from the centre of the Sun to the Earth.

2 Final state polarisation in dark matter annihilations

In many cases, the final states in the dark matter annihilation processes can be preferentially produced in a specific helicity configuration. This will affect the shapes of the spectra of the yield particles [10, 25], a result of non-isotropic angular distributions of the particles that emerge from the decay of the polarised particle. Consider for example the annihilation of dark matter particles into a polarised pair of final state particles. If the final state particles decay through a two-body decay, the emitted particles are monochromatic in the rest frame of the decaying particle. Their angular decay distributions are however not necessarily isotropic, with some directions preferred over others in the decay, for example when the decay proceeds through a chiral coupling. Therefore, when the decay products are boosted to the rest frame of the annihilation, they will have an energy distribution which is generally non-flat. The flat distribution arises only in case the final state particle decays isotropically (as for example for a decaying scalar) or is unpolarised.

In the case of annihilation into massive gauge bosons, some degree of polarisation is expected from the separation of the longitudinal and transverse degrees of freedom—the longitudinal degrees of freedom are connected to the Goldstone bosons of electroweak symmetry breaking and thus connected to the interactions with the Higgs sector, whereas the transverse degrees of freedom are associated with the pure gauge part [26, 10]. Depending on the strength of the interaction of the dark sector with the Higgs and pure gauge part respectively, the production rate of longitudinal and transverse polarised states can therefore differ. As an example, we note that, in neutralino annihilation into a pair of gauge bosons, the gauge bosons can be produced with purely transverse polarisation [26, 7].

Another case where polarisation can play an important role is when a Dirac dark matter candidate annihilates with its conjugate particle to SM fermions in a spin 1 configuration through a vector mediator. It is possible for the SM fermions to be produced with a net polarisation if this vector mediator $B^{\mu}$ has a chiral coupling to the SM fermions, i.e. if it couples to the fermions via a term of the form $B^{\mu}(g_{L}\bar{f}\gamma_{\mu}P_{L}f+g_{R}\bar{f}\gamma_{\mu}P_{R}f)$. If $g_{R}>g_{L}$, dark matter annihilations to left-handed SM fermions and right-handed anti-fermions will dominate (and vice versa). Such a coupling structure can arise in $Z^{\prime}$ models that are anomaly free [12]. Dirac neutralino dark matter can also lead to a similar excess of one polarisation state over another [11].

In general, the polarisation of the final state particles depends on the nature of the interactions between the dark sector and the SM. It is therefore difficult to make model-independent statements about what polarisation to expect. Here we refrain from committing to a specific particle physics model and see what effect the polarisation can have by considering spectra from pure final state polarisations that should in most cases bracket the spectra in a full particle physics model.

3 Event generators

An event generator uses Monte Carlo methods to simulate the complexity of a particle physics process [17, 27]. In order to obtain fluxes of the stable particles from a particular dark matter annihilation process we have used two of the more popular general purpose event generators, Pythia8 and Herwig7, for hadronisation and showering of hard process events simulated with MadGraph5. For a given final state of the annihilation process, the event generator is responsible for the evolution of the hard process final states into the final list of stable particles from the annihilation. This includes decay of unstable annihilation products, radiation of coloured particles and photons as well as the hadronisation of partons and the subsequent decay of unstable hadrons. We will in the following review the main characteristics of general purpose event generators, restricting ourselves to that which is relevant for the dark matter annihilation fluxes. We then specifically consider the properties of Pythia8 and Herwig7.

3.1 General purpose event generators

Event generators typically divide the modelling of a particle physics process into several parts. The high momentum transfer interaction is called the hard process. If coloured particles are produced in the hard process, they undergo a showering process, successively radiating coloured particles and evolving downwards to lower scale until they reach a cutoff $Q_{0}$, typically around $1\text{\,}\mathrm{GeV}$, which signifies the transition to the non-perturbative regime of QCD. At this point, all partons form colour-neutral hadrons in a process called hadronisation. Unstable hadrons are then decayed, leading to the ultimate list of stable particles that are produced in the interaction.111Note that whether a particle is defined as stable or not depends on the application, in an accelerator all particles that decay outside the detector are denoted stable. In the current situation where we consider dark matter annihilations in an astrophysical environment, we have to make sure that these particles are allowed to decay. Non-coloured particles are decayed to stable states as well, and photon and electroweak radiation can be taken into account in the showering and decay processes.222In general one also has to consider radiation of the incoming partons in the hard process as well as multiple interactions between incoming partons for collisions of composite particles; this is however not necessary in the current case of dark matter annihilations, since the dark matter particle is assumed to be an uncharged elementary particle.

The hard process is represented by the matrix element of the process in question; in the case of annihilation into coloured states it is the partonic amplitude. In this paper we take the hard process to be dark matter annihilating via an $s$-channel mediator to a pair of SM particles. While our results are broadly applicable to a range of dark matter models considered in the literature, it is important to note that there exist models in which final states other than two SM particles dominate, such as secluded dark matter scenarios (e.g. [28, 29]) or models in which radiation off a charged $t$-channel mediator is important (e.g. [5, 30]).

In the case of partonic annihilation channels, or when the annihilation channel particle decays to quarks, the outgoing partons from the hard process radiate coloured particles, losing energy in the process. These parton showers take into account radiation of coloured particles like quarks and gluons, and also photons and electroweak bosons [18] can be included in the showering procedure. This radiation is dominated by collinear and soft emission and therefore showering algorithms have traditionally focused on radiation in this limit. The showering is performed iteratively, starting from the outgoing coloured states from the hard process and evolving down to a scale $\sim$1\text{,}\mathrm{GeV}$$, where hadronisation kicks in.

The choice of evolution variable in the showers differs between generators. Older Pythia versions evolved in terms of the virtuality of the partons [19]. Using the virtuality as evolution variable can be shown to be incorrect when moving away from the collinear and soft limit, due to insufficient modelling of colour coherence effects [31]. To take these effects into account one can instead use an angular ordering, evolving in the angle between the outgoing partons in the shower [32, 33]. Another way of accounting for the colour coherence effects is to consider radiation from colour dipoles instead of partons, so that dipole branchings are considered rather than parton branchings [31, 34, 35].

Hadronisation and hadron decays is the final part of the event generation. This consists of the formation of colour neutral hadrons out of the partons coming from the parton showering. Since this occurs in the non-perturbative regime of QCD one has to resort to phenomenological modelling of the hadronisation process. The two most common models are the Lund string model [36, 37] and the clustering model [38]. The hadronisation part of the event generation includes phenomenological parameters that must be tuned to experimental data.

In the string model, all the partons from the showering are connected by strings that represent the strong force potential between them. These strings are stretched as partons continue to move outwards and splits when the potential energy is sufficient to create a new $q\bar{q}$ pair. The gluons in the shower are represented by kinks on the strings. Further complexities allow for baryon formation. After string evolution stops, unstable hadrons are decayed.

The clustering model starts by enforcing a splitting of all gluons into $q\bar{q}$ pairs. Colour neutral clusters are then formed out of the resulting quarks and antiquarks with the clusters then forming low-mass hadrons, where the specific hadron is determined mainly by phase space and mass. Any unstable hadrons are then decayed.

3.2 Comparing Pythia8 and Herwig7

The history of Pythia goes back to the development of the Lund string model [39]. Today it is one of the most used event generators. The latest major release Pythia8 [13] consisted, apart from physics updates, of a complete rewrite of the old Fortran version Pythia6 [19] into a C++ version.

For the parton showering, Pythia8 uses a modified version of the dipole showering scheme, where a variable closely related to the transverse momentum is used as the evolution variable [40, 41]. The hadronisation is performed with the Lund string model.

Herwig7 [16, 15, 14] continues the development of Herwig++ [14] which in turn was a major C++ rewrite of the the earlier Fortran version. Herwig has traditionally relied on angular-ordered parton showers, this continues to be the case in Herwig7, with additional options for dipole showering also included. Hadronisation is in Herwig7 performed using the clustering model, as in older versions.

By using helicity amplitudes, Herwig7 tracks spin throughout the event generation machinery and therefore generally includes spin correlation effects. This is generally not the case in Pythia8 where the effects are included only for some specific cases. As described in section 2, such effects can significantly change the shape of the energy distributions.

Given the differences between Pythia8 and Herwig7 in the modelling of showers and hadronisation, we expect some differences in the fluxes of the final yield particles, even though the generators are for the most part tuned to the same data sets. This is particularly true in the case of annihilation into quarks or for the fluxes of, for example, anti-protons where QCD effects are expected to be more important.

4 Simulation procedure

To obtain the fluxes of the yield particles of interest we have simulated annihilations of dark matter particles at rest into an $s$-channel uncharged resonance of mass $2m_{\chi}$ that decays into the desired annihilation channel. We then let the annihilation products undergo parton showers, decays, and hadronisation. All unstable particles (including muons, pions, kaons and neutrons except for the case of annihilations inside the Sun) are forced to decay. In each annihilation we count the number of yield particles appearing and record each particle’s kinetic energy. Dividing by the number of annihilations generated we thus obtain $dN/dE_{\rm{kin}}$—the number of yield particles per annihilation as a function of their kinetic energy.

We have considered a mix of quark, lepton, and gauge boson annihilation channels

[TABLE]

with the heavier quark and leptons chosen because they will lead to more complicated showering and decay processes. Polarised $b$ quarks are not considered since these hadronise before decaying and the polarisation is lost in the showering and hadronisation. The missing $Z^{0}Z^{0}$ final state has similar phenomenology to the $W^{+}W^{-}$, with the important difference that there will be no QED bremsstrahlung off the $Z$.333Note that our choice of polarisations for the final state fermions differs from that of the PPPC [21], which instead includes final states in which both the fermion and anti-fermion are in the same polarisation state. The reason behind this is theoretical – while a chiral coupling of the fermions to a vector boson of the form discussed in section 2 can lead to annihilations to fermions in a definite LR or RL helicity state, the need for the Lagrangian to be Hermitian seems to exclude the possibility of preferential annihilations to LL or RR final states. For example two dark matter particles in a $J=0$ configuration could annihilate to two fermions in a RR state if they annihilated via an $s$-channel scalar mediator $\phi$ with couplings of the form $\phi\bar{f}P_{L}f$. However such a term is not Hermitian by itself, so its Hermitian conjugate, $\phi\bar{f}P_{R}f$, must also be included in the Lagrangian, leading to final state fermions that are unpolarised on average. We also do not directly simulate gauge boson final states that are a mix of transverse and longitudinal polarisations. Assuming that $W^{-}_{L}W^{+}_{T}$ is produced at the same rate as $W^{-}_{T}W^{+}_{L}$, the average spectra from an annihilation to these states can be determined by combining our results via the relation $(W^{-}_{L}W^{+}_{T}+W^{-}_{T}W^{+}_{L})/2=(W^{-}_{L}W^{+}_{L}+W^{-}_{T}W^{+}_{T})/2$. For these final states, we have constructed the energy distributions for the yield particles $e^{\pm}$, $p$, $\bar{p}$, $\gamma$ and all $\nu_{\ell}$, $\bar{\nu}_{\ell}$ for $\ell=e,\mu,\tau$. For the annihilation channels, the unpolarised state is used in the cases where there is no subscript, otherwise the subscript indicates the helicities of the particles. The subscript denotes for fermions whether it is right- ($R$) or left-handed ($L$), while for gauge bosons a subscript of $L$ ($T$) indicates longitudinal (transverse) polarisation, with the two transverse polarisation states assumed to occur at equal rates.444To convert our results to the yields from unpolarised final states, the simple relations $f^{-}f^{+}=(f^{-}_{L}f^{+}_{R}+f^{-}_{R}f^{+}_{L})/2$ and $W^{-}W^{+}=(2W^{-}_{T}W^{+}_{T}+W^{-}_{L}W^{+}_{L})/3$ can be used. As described below, we have used MadGraph5 and MadEvent [42] for generation of hard process events and some decays, while Pythia8 or Herwig7 are used to simulate the subsequent showering, hadronisation, and other decays not handled by MadEvent. We use MadGraph5 rather than the internal Pythia matrix elements because of the limited ability of Pythia to handle polarised final state particles. The same hard process events from MadEvent are then used in Herwig for consistency in our comparisons. We have used version 8.240 of Pythia, version 2.6.5 of MadGraph5 and version 7.1.1 of Herwig in our simulations. These will be referred to as Pythia8, MadGraph5 and Herwig7 throughout.

For the MadGraph5/MadEvent simulations, we have set up a simple FeynRules model that includes a scalar resonance (for the gauge boson channels) and a spin-1 resonance (for the fermionic channels), which couple to all fermionic annihilation channel particles through maximally parity violating couplings with left- and right-handed projectors.555A vector or pseudo-vector resonance that couples to the fermions via $\bar{f}\gamma_{\mu}f$ and $\bar{f}\gamma_{\mu}\gamma_{5}f$ respectively would give unpolarised final states on average. We can then control the polarisation of fermion final states by setting these couplings. For the gauge boson final states we achieve a longitudinally polarised final state by coupling the resonance to a term of the form $V_{\mu}V^{\mu}$, for gauge fields $V^{\mu}$, whereas the transverse polarisation is achieved by a coupling to the squared field strength tensor of the gauge field, $F^{\mu\nu}{F}_{\mu\nu}$. A pseudoscalar mediator which couples to $F^{\mu\nu}\widetilde{F}_{\mu\nu}$ would result in the same energy spectrum [8]. We implement the model in MadGraph5 via the UFO [43] interface, allowing us to use it in event generation. From the MadEvent simulations, we produce event files in the Les Houches Event File (LHEF) format [44, 45], which are then read into Pythia8 and Herwig7 to perform the showering and hadronisation.666In Herwig7, the EvtGen code is used for the decay of heavy hadrons containing $b$- or $c$-quarks. In order to improve efficiency in the simulations, we have had to turn off this interface and use the internal Herwig7 decays instead. This affects primarily the $b\overline{b}$ channel, where many of these heavy hadrons are produced. The difference is generally smaller than the difference compared to Pythia8 simulations.

For annihilation channels into unstable particles where we are interested in polarisation effects, and there are no built-in features in the event generator to take into account the polarisation of the final state, we let MadEvent perform the decays. We use the decay chain syntax to decay electroweak gauge bosons and top quarks down to parton or lepton level. The top quarks are decayed as $t\to bW^{+}$ and electroweak gauge bosons are decayed into all possible two-body decays. Since there are functionalities in Pythia8 and Herwig7 to perform polarised $\tau$ decays [46, 47], we do not decay these in MadGraph5 but use Pythia8 or Herwig7. Finally, following the code defaults and in order to make consistent comparisons we do not turn on electroweak radiation in the showers (this is turned off by default in Pythia8 and not available in Herwig7).

5 Results for annihilation fluxes

In this section we look at the fluxes of the various yield particles of interest produced in the annihilations. Broadly speaking, for a given annihilation channel, the flux of a specific yield particle will be harder the earlier it can appear in the decay chain, with the hardest flux produced when the dark matter particles annihilate directly into the yield particle in question. Typically, there is also a soft component of the spectrum from radiation and decays of lower energy particles produced in the hadronisation. In the following, we go into more details on the spectra obtained for the various annihilation channels.

5.1 Differences between Pythia8 and Herwig7

In figures 1, 2, 3 and 4 we show the yields of muon neutrinos, positrons, gamma-rays and antiprotons for the simulated annihilation channels, here for a dark matter mass $m_{\chi}=$1000\text{,}\mathrm{GeV}$$. We show results for simulations using MadGraph5 and Pythia8 (referred to in the following as P8) or MadGraph5 and Herwig7 (H7) and indicate the difference with the coloured band. As reference, we also show results extracted from DarkSUSY in dash-dotted lines, that have been simulated using Pythia6 (P6) with unpolarised final states, and the results of ref. [21] (PPPC) in dotted lines.

In Appendix A, similar plots are shown for a dark matter mass of $m_{\chi}=$100\text{,}\mathrm{GeV}$$. While the shape of the spectra are different at this lower energy, ultimately the ratios between the spectra produced are similar, with the exception of photons from the $W^{+}W^{-}$ final state, as discussed in section 5.1.3.

5.1.1 Neutrino yields

We show in figure 1 the muon neutrino yields. For the neutrinos, the hardest spectra are produced in the gauge boson annihilation channels, since the gauge bosons can decay directly into leptons. In two-body decays such as $W^{-}\to\ell^{-}\overline{\nu}_{\ell}$, the leptons are monochromatic in the rest frame of the decaying gauge boson. When boosted to the annihilation rest frame, the spectrum is spread between energies $E_{-}$ and $E_{+}$, in the annihilation rest frame given by

[TABLE]

The shape of the spectrum between $E_{-}$ and $E_{+}$ depends on the polarisation of the decaying gauge boson, with the longitudinal (transverse) polarisation resulting in a convex (concave) shaped lepton spectrum.

Annihilation into tau leptons and top quarks also produces relatively hard neutrino spectra, although somewhat softer than the gauge boson spectra. The softest neutrino spectra are produced in the $b\overline{b}$ channel, where the neutrinos predominantly come from decays of hadrons.

We note that in general, differences between P8 and H7 are small for neutrinos. Given that the main differences between the event generators lie in the handling of QCD effects, it is expected that there will be less differences in leptonic processes. As the P6 simulations are for unpolarised final states, their yields typically lie between those from polarised final states, as expected. The PPPC results show an excess at high energies relative to our results from both P8 and H7, particularly in the quark channels.

5.1.2 Positron yields

As one can see in figure 2, the situation for positrons is similar as for the muon neutrinos. The gauge boson spectra are very similar at high energies, as expected, since these positrons also come mainly from direct $W$ boson decays. At lower energies there are some differences between P8 and H7, with the P8 simulation giving a larger spectrum. This difference can most likely be traced to the differences in the hadronisation procedure, as the low energy positrons mainly come from decays of particles produced in the hadronisation step.

As expected, the main differences between the event generators are found for the quark channels, primarily for the $b\overline{b}$ channel, where they reach around 20%. Again the PPPC yields show an excess at high energies in the hadronic channels relative to our results.

5.1.3 Gamma-ray yields

For the gamma-ray yields, shown in figure 3, there are large differences for the $W$ boson channels. The difference is related to the bremsstrahlung of photons from the $W$ bosons in the final state. The rate of bremsstrahlung in H7 is significantly lower than in P8.777We have been unable to conclusively trace the source of this difference. We have seen in the simulations that in Herwig7, unlike in Pythia8, no photons at all seem to be emitted from the $W$ bosons. We conclude that QED bremsstrahlung has a significant effect on the high energy part of the gamma-ray spectra. Reducing the energy of the $W^{+}W^{-}$ pair of course reduces the amount of bremsstrahlung radiation, and so at lower dark matter masses, the agreement between the codes is better, as shown in figure 12.

For the other channels, differences are smaller, and mainly concentrated to the tails of the distributions. This is expected, since at high energies the event generators rely to a higher extent on extrapolation away from the tuning regions.

5.1.4 Antiproton yields

As expected, consistent differences between the generators are found in the antiproton yields, shown in figure 4. The ratio between the two are for most energies around a factor of two and the peaks of the distributions are slightly shifted, with P8 generally giving harder spectra than H7, except in the $b\overline{b}$ channel, where the H7 spectrum is instead more energetic.

5.2 Effect of polarisation

We show in figure 5 yields of positrons, muon neutrinos, antiprotons and gamma-rays for different polarisations in the final state, with the difference between the polarisations highlighted by the coloured band, and the average represented by the dotted line of the same colour. An unpolarised final state will follow this line in the centre of the band. We show here results obtained with the P8 simulation, compared with simulations performed with P6 (for unpolarised final states) in the dash-dotted curves. The coloured band in this case represents the effect that the polarisation of the final states can have on the yields.

We see in figure 5 that the differences between different polarisations of the final states are most apparent at higher energies. This is expected, since the polarisation in our simulations affects the decays occurring earliest in the decay chain. For example, we can clearly see the expected concave and convex shape of the neutrino spectrum from $W_{T}^{-}W_{T}^{+}$ and $W_{L}^{-}W_{L}^{+}$ final states respectively, with the spectrum in the $W_{T}^{-}W_{T}^{+}$ final state clearly peaking at the highest energies. We also see, as expected, the unpolarised P6 result following closely the average between the two polarisations apart from the gamma-ray spectrum in the $WW$ channel. The difference between Pythia6 and Pythia8 is here most probably linked to a lack of QED bremsstrahlung in Pythia6.

We can also see differences in the gamma-ray and antiproton spectra. The differences for gamma-rays can likely be attributed to the different energy spectra of the charged particles which radiate the photons. We can also see that, in particular for the antiproton spectra, there are differences between the polarisation average and the P6 result. This likely depends on the fact that the hadronisation and showering has changed between P6 and P8.

6 Impact on dark matter searches with neutrino telescopes

For the neutrino fluxes, polarisation in the final state can have an important impact on dark matter searches. The polarisation of the final state results in differences in the high-energy part of the spectra, as shown in figure 5 above. This difference is further enhanced for the case of neutrino searches since the detectability of neutrinos in a neutrino telescope scales roughly as the flux times the neutrino energy squared, $E_{\nu}^{2}$, where one factor of $E_{\nu}$ comes from the rise of the neutrino interaction cross section with $E_{\nu}$ and the other from the muon range in the detector, which rises approximately linearly with muon energy, and hence is proportional to $E_{\nu}$. In addition, the main background in dark matter searches, the atmospheric neutrino flux, decreases sharply with energy. Therefore, the high end part of the flux becomes highly relevant for dark matter searches with neutrino telescopes. For this reason we focus in this section on the impact of polarisation on the neutrino fluxes specifically.

6.1 Simulation of neutrino event rates with WimpSim

To make predictions of event rates in neutrino telescopes for different polarisations we use the WimpSim package [2, 48, 49]. WimpSim is an event based Fortran code that simulates the full chain from dark matter annihilation to the interaction in a neutrino telescope on Earth. The neutrinos are generated in dark matter annihilations in the centre of the Sun or Earth. The produced neutrinos are then propagated from the production point to a detector on Earth, where interactions in the Sun or Earth and oscillations are handled in a full three-flavour setup. Apart from simulating dark matter annihilations in the centre of the Sun or Earth [2], there are also options to simulate the neutrino flux from cosmic ray interactions in the solar atmosphere [48] as well as the flux of neutrinos, gamma-rays and charged cosmic rays from annihilation of dark matter into long-lived mediators that decay away from the solar core [49].

Traditionally, WimpSim has used Pythia6 for the simulation of the dark matter annihilations. For this study, we primarily want to study other event generators (namely Pythia8 and Herwig7) and have therefore added a new functionality in WimpSim to make it possible to read in event files with neutrino yields at production. This allows us to study the effect the different spectra from our simulations will have on the event rates in neutrino telescopes. We have chosen to do this part of the analysis using only the Pythia8 simulation data, but results using Herwig7 will be similar given that the differences between the generators is generally small for neutrino spectra. For annihilations in the Sun, we have taken care not to include neutrinos coming from the decays of muons, neutrons or charged kaons and pions, as these particles are quickly stopped in the solar interior and do not give rise to high energy neutrinos.

For neutrinos coming from heavy hadrons containing $b$ or $c$-quarks, one should take into account the interactions of the hadrons with the solar material and reduce the neutrino energies correspondingly. Currently we have not included this when reading event files into WimpSim for the production fluxes. This effect will be most important for final states giving many of these heavy hadrons (for example the $b\overline{b}$ channel) and since we focus on $W$ boson and $\tau$ lepton channels for this part of the study, it will not have a large effect.

The annihilation point is sampled from a thermal distribution around the solar core. Assuming a thermal annihilation distribution has been shown to be a very good approximation [50]. The neutrinos are then allowed to interact through charged current (CC) or neutral current (NC) interactions and oscillate (including matter effects) on their way from the annihilation point to the detector. If a $\tau^{\pm}$ is produced in a CC interaction, the $\tau^{\pm}$ is decayed and any neutrinos from the decay are collected and included in the remaining neutrino propagation. For this decay, Pythia6 is used. As we are here interested in polarised final states in the dark matter annihilation, the fact that the Pythia6 $\tau$ decays do not include polarisation effects is not relevant for us.

As soon as $m_{\chi}$ is large enough it is essential to take into account the CC and NC interactions when simulating the neutrino fluxes. In the CC interactions, neutrinos are effectively lost, leading to an attenuation of the flux at high energies (since the cross section rises with the neutrino energy). The main effect of the NC interactions is to shift the neutrino flux to lower energies.

In the case of annihilations in the Sun, neutrinos are propagated to a distance of $1\text{\,}\mathrm{AU}$ from the centre of the Sun. They are then propagated from $1\text{\,}\mathrm{AU}$ to a specified detector location on Earth, taking into account the eccentricity of the Earth’s orbit, and allowed to interact in the detector volume to produce a neutrino signal in the telescope.

For neutrino oscillations and interactions we need values for the neutrino oscillation parameters, and furthermore need a solar model. We have used the best-fit values of the neutrino oscillation parameters from the NuFit-3.2 global fit [51, 52] (summarised in table 1). We consider only the normal mass hierarchy here. For the solar model, we use the model by Serenelli et al. [53] which is the default in DarkSUSY.

7 Neutrino telescope results

In figures 6, 7, 8 and 9 we show neutrino fluxes obtained in the WimpSim simulations for the $W$ and $\tau$ annihilation channels. We show the sum of neutrino and antineutrino fluxes for all neutrino flavours, in the left column of the figures at creation and in the right column after propagation to $1\text{\,}\mathrm{AU}$ from the Sun. We show in figures 6 and  7 the fluxes for a dark matter mass of $100\text{\,}\mathrm{GeV}$, where interactions in the Sun are less important. To highlight the effect of interactions we show in figures 8 and  9 the fluxes for a dark matter mass of $1000\text{\,}\mathrm{GeV}$.

Generally, we can see the effect of oscillations between the annihilation point and the Earth as sinusoidal behaviours in the energy spectra. The oscillations will appear in given ranges of neutrino energies (with one energy range corresponding to each of the neutrino squared mass differences) and are mainly visible above $100\text{\,}\mathrm{GeV}$ in the spectra in the right-hand parts of figs. 8 and 9 for $m_{\chi}=$1000\text{,}\mathrm{GeV}$$. For more details regarding the oscillation effects we refer to ref. [2].

For the fluxes from the $\tau$ lepton channels, we can see that the $\nu_{e}+\overline{\nu}_{e}$ and $\nu_{\mu}+\overline{\nu}_{\mu}$ fluxes are very similar at production, whereas the $\nu_{\tau}+\overline{\nu}_{\tau}$ fluxes are significantly larger, especially at high energy. The explanation from this comes from the fact that the high energy part of the $\nu_{e}+\overline{\nu}_{e}$ and $\nu_{\mu}+\overline{\nu}_{\mu}$ fluxes come from the $W$ bosons in the $\tau$ decays, and the branching ratio into the two flavours is approximately equal whereas the high energy $\nu_{\tau}$ can come directly from the $\tau$ decays, resulting in higher average energies. We can also see that for the polarisation channel where the $\nu_{\tau}+\overline{\nu}_{\tau}$ flux is higher, the $\nu_{e}+\overline{\nu}_{e}$ and $\nu_{\mu}+\overline{\nu}_{\mu}$ flux is lower, and vice versa.

For the $W$ channels at $m_{\chi}=$100\text{,}\mathrm{GeV}$$, we can clearly see the box-shaped spectrum (which here appears tilted, since we show $x^{2}dN/dx$ on the vertical axis), is now either concave, linear or convex depending on the polarisation of the $W$ bosons. In the $\nu_{\tau}+\overline{\nu}_{\tau}$ fluxes there is an additional component in the lower end of the box, coming from $\tau$ decays.

In general, we can see in our results that differences between polarisation states at creation can be to some extent washed out by propagation effects (interactions and oscillations). For example we can see that for the $\tau$ annihilation channels, although there is a rather large difference between the $\nu_{\mu}+\overline{\nu}_{\mu}$ fluxes from the $\tau_{L}^{-}\tau_{R}^{+}$ and $\tau_{R}^{-}\tau_{L}^{+}$ channels at creation, this difference becomes smaller at $1\text{\,}\mathrm{AU}$. This is primarily due to the oscillations between $\nu_{\tau}$ and $\nu_{\mu}$, which means that the $\nu_{\mu}+\overline{\nu}_{\mu}$ flux at $1\text{\,}\mathrm{AU}$ is approximately the average between the two neutrino flavours at creation, and this average is similar for the two polarisation channels (and closer to the unpolarised case).

Interestingly, we can also see that for the muon neutrinos, the $\tau_{L}^{-}\tau_{R}^{+}$ flux is lower than the $\tau_{R}^{-}\tau_{L}^{+}$ flux at creation, but larger at $1\text{\,}\mathrm{AU}$. This can be explained by the fact that the $\nu_{\tau}+\overline{\nu}_{\tau}$ flux is instead higher for the $\tau_{L}^{-}\tau_{R}^{+}$ channel than for the $\tau_{R}^{-}\tau_{L}^{+}$ channel, and sufficiently much larger to end up with a larger $\nu_{\mu}+\overline{\nu}_{\mu}$ flux for $\tau_{L}^{-}\tau_{R}^{+}$ at $1\text{\,}\mathrm{AU}$.

For the $W$ boson channels, it is interesting to see that the propagated flux at $1\text{\,}\mathrm{AU}$ is not very different from the production flux, apart from a slight shift towards lower energies, likely caused by interactions in the Sun. Thus a difference between the polarisation states remains present at the highest neutrino energies for the $W$ channels in the $1\text{\,}\mathrm{AU}$ flux.

Apart from oscillations, propagation effects also include neutrino interactions in the Sun. For the $m_{\chi}=$1000\text{,}\mathrm{GeV}$$ fluxes, we clearly see how interactions act to attenuate the fluxes and shift them towards lower energy. This effect also decreases the differences between the polarisation states. The high energy peak for the $W_{T}^{-}W_{T}^{+}$ channel in particular is damped by the interactions, making it less distinguishable from the $W_{L}^{-}W_{L}^{+}$ channel, which is less affected by this attenuation effect. As the probability to detect a neutrino in a neutrino telescope goes like $E_{\nu}^{2}$ one would naively have believed that the $W_{T}^{-}W_{T}^{+}$ channel would give a higher event rate. However, for large masses it actually goes the other way due to neutrino interactions in the Sun.

We can quantify this by looking at expected number of events in a neutrino telescope for these spectra. We do this by calculating the number of events in IceCube using their latest 3 year data analysis for dark matter searches from the Sun [54]. We use their given effective area (we will use the maximum of DeepCore and IceCube effective areas for each given energy) and calculate the expected number of events (per annihilation) for the lifetime of the detector in these three years for our different spectra (including both muon neutrinos and antineutrinos). Our results are shown in table 2. As one can see, the shape of the $W_{T}^{-}W_{T}^{+}$ spectrum gives a higher event rate than the $W_{L}^{-}W_{L}^{+}$ spectrum for a 100 GeV WIMP mass. However, for large masses, neutrino interactions in the Sun cause this harder spectrum to be absorbed more, and the end result is that the event rate for a 1000 GeV WIMP is actually lower than for the $W_{L}^{-}W_{L}^{+}$ spectrum. The same trend can be seen for the $\tau_{R}^{-}\tau_{L}^{+}$ spectra compared to the $\tau_{L}^{-}\tau_{R}^{+}$ spectra, even though in this case, the $\tau_{R}^{-}\tau_{L}^{+}$ spectrum still gives a higher event rate even for a 1000 GeV WIMP. The unpolarised Pythia6 results are typically in between our new polarised results, except for $W^{-}W^{+}$ at 1000 GeV, where the Pythia6 results is comparable to the $W_{T}^{-}W_{T}^{+}$ results (which can also be seen from the spectra in figure 8).

We end this section by noting that the differences we see for the $W^{-}W^{+}$ final state polarisations are slightly more than the 5–10% quoted by [24]. For the $\tau^{-}\tau^{+}$ channel, the effect can be even larger than this. However, for more massive dark matter, interactions tend to wash out some of the differences, even though some differences remain. This effect of washing out differences for more massive dark matter has also been observed in ref. [55, 7]

8 Summary and conclusions

In this study we have investigated the fluxes of secondary particles from dark matter annihilations, such as neutrinos, charged antiparticles and gamma-rays. A proper knowledge of these fluxes is necessary in order to analyse results in indirect searches for dark matter. We have simulated the dark matter annihilations with two of the most commonly used event generators, Pythia8 and Herwig7 (in combination with MadGraph5), in order to understand the impact that the differences between the modelling in the generators can have. For some of the steps in the simulation chain, the event generators rely on phenomenological modelling of the underlying physics, and this modelling in general differs between generators. In particular this concerns the QCD-related parts of the process, which is most relevant when the dark matter annihilation proceeds into colour charged final state particles such as quarks.

We have also looked at the effect that polarisation of the final state can have on the fluxes of the secondary particles. Polarisation of the final state is expected in many concrete dark matter models, for example if the annihilation proceeds through a chiral coupling or if there is a separation between the couplings to the longitudinal and transverse degrees of freedom of gauge bosons.

Concerning the comparison between Pythia8 and Herwig7, we find the smallest differences in the fluxes of neutrinos and positrons and the largest differences for antiprotons and, in some cases, gamma-rays, with the spectra varying by a factor of 2 or more for some energies. Differences in the leptonic spectra are typically of the order of 10-20% or lower. Apart from the gamma-ray fluxes, where a difference in the modelling of QED bremsstrahlung is a likely cause of the differences, it is expected that antiproton fluxes will show larger discrepancies than leptonic fluxes, since the antiprotons emerge from the hadronisation part of the event generation, which the event generators model differently. While uncertainties in the yields from event generators are often ignored in indirect detection analyses, these results illustrate the need for these uncertainties to be included, particularly for the case of hadronic processes.

The polarisation of the final state instead primarily affects the leptonic spectra. This is because the most energetic leptons are the ones coming from decays early in the chain, close to the hard process, whereas in the showering and hadronisation processes, producing the antiprotons and many of the gamma-rays, polarisation effects are washed out. Differences in the leptonic spectra for the different polarisations of the final state in many cases reaches a factor of ten, demonstrating the importance of taking into account the polarisations predicted by a model in determining these yields. The gamma-rays and antiprotons are not produced directly as decay products of the final states, and differences in these spectra are most likely traced to spectral changes for the other particles involved.

We have also looked more specifically at the neutrino fluxes from dark matter annihilations in the Sun. For this purpose, it is important to include the interactions of neutrinos with the solar material and neutrino oscillations when propagating the neutrinos from the annihilation point in the solar core to a neutrino telescope on Earth. We use the code WimpSim to account for these effects.

Interestingly, although the differences in the neutrino spectra for different polarisation states are relatively large at production in the annihilation process, these differences are to an extent reduced by the propagation effects. For example, for annihilations into $\tau$ leptons, in the polarisation state where the muon neutrino flux is larger, the tau neutrino flux is lower, and vice versa. Since these neutrino flavours mix significantly in the oscillations, the general size of the resulting muon neutrino flux at Earth ends up in between the fluxes of the two flavours at production, and this flux is similar for the two polarisation states. Hence, the polarisation differences are ultimately averaged out by the oscillations.

Neutrino interactions in the Sun also act to reduce the differences between the polarisation states. We have shown that for larger values of the dark matter mass, the attenuation of the fluxes at high energy, stemming from charged current interactions of the neutrinos with the solar material, dampens the high energy peak in the neutrino spectrum from transversely polarised $W$ bosons in the final state. This dampening reduces the differences in the neutrino flux from the $W_{T}^{-}W_{T}^{+}$ and $W_{L}^{-}W_{L}^{+}$ channels. At lower dark matter masses, where interactions are less important, there is however still a difference between the polarisations remaining in the propagated neutrino fluxes.

Finally, to quantify the above, we have looked at the differences in the event rates in IceCube. We find that the different polarisation states have an effect on the event rates which is typically around 10% in the cases considered but reaches over 50% in the $\tau$ annihilation channel at $100\text{\,}\mathrm{GeV}$.

Acknowledgments

The authors wish to thank Philip Ilten for assistance with $\tau$ lepton decays in Pythia and Marco Cirelli and Jared Evans for helpful discussions and correspondence. JMC gratefully acknowledges partial support from US Department of Energy (DOE) Grant No. DE-SC0020047 and expresses a special thanks to the Mainz Institute of Theoretical Physics (MITP), which is part of the DFG Cluster of Excellence PRISMA* (Project ID 39083149), for its hospitality and support. JE thanks the Swedish Research Council for support (contract 621-2014-5772).

Appendix A Yields for $m_{\chi}=$100\text{,}\mathrm{GeV}$$

In this section we show figures similar to figures 1-4, but for a dark matter mass of $100\text{\,}\mathrm{GeV}$.