Radiative corrections to vectorlike portal dark matter

A massive real scalar dark matter particle $S$ can couple to Standard Model leptons or quarks through a vector-like fermionic mediator $\psi$, a scenario known as the Vector-like portal. Due to helicity suppression of the annihilation cross section into a pair of SM fermions, it has been shown in previous works that radiative corrections, either at one-loop or through radiation of gauge bosons, may play a significant role both in determining the relic abundance and for indirect detection. All previous works considered the limit of massless final state quarks or leptons. In this work, we focus on a technical issue, which is to reliably determine the annihilation cross sections taking into account finite fermion masses. Following previous works in the framework of simplified supersymmetric dark matter scenarios, and building on an analogy with Higgs decay into fermions, we address the issue of infrared and collinear divergences that plagues the cross section by adopting an effective operator description, which captures most of the relevant physics and give explicit expressions for the annihilation cross sections. We then develop several approximations for the differential and total cross sections, which simplify greatly their expressions, and which can then be used in various phenomenological studies of similar models. Finally, we describe our method to compute the final gamma-ray spectrum, including hadronisation of the heavy fermions, and provide some illustrative spectra for specific dark matter candidates.


I. INTRODUCTION
Dark matter (DM) amounts to about 27% of the energy budget of our universe and, yet, little is known about its precise nature. A much studied possibility is that dark matter is made of a weakly interacting massive particle (WIMP). In this scenario the observed relic abundance Ω DM h 2 = 0.1186 ± 0.0020 [1] more or less naturally results from chemical freeze-out of a massive particle, provided its annihilation cross section is σv 3 · 10 −26 cm 3 · s −1 . Many such candidates have been proposed in the literature.
In the present note, we study further a real scalar particle coupled to Standard Model (SM) fermions through a vector-like fermion. This scenario has been dubbed the Vector-like Portal (VLP) in [2] following [3]. The simplest realization of the VLP is given by where S is a singlet real scalar (the DM candidate), f R is an SU (2) L singlet SM fermion (lepton or quark) and ψ a vector-like massive fermion. Clearly other combinations of SM multiplets could be considered, but this is not our purpose here. Stability of DM is ensured by imposing a discrete Z 2 symmetry, under which both S and ψ are taken to be odd and SM particles even. For simplicity and to avoid addressing flavour physics aspects, it is assumed that S couples dominantly to a single SM flavour. Also, in the sequel we will assume that the possible quartic coupling of S with the SM Higgs is small and may be neglected. As such, (1) falls in the category of so-called simplified DM models with a t-channel mediator, see e.g. [4][5][6][7]. This model may be considered as the scalar version of a bino-like Majorana DM candidate, with which it shares some basic properties, the first being that their s-wave annihilation is helicity suppressed. They differ in the fact that annihilation of a bino-like candidate is pwave in the chiral limit [8] while that of the real scalar S is d-wave [2,9]. However, in both cases, the helicity suppression is lifted by radiative corrections [10,11]. As discussed in several works, this has interesting phenomenological implications. In the case of coupling to leptons, radiative processes, either in the form of internal bremsstrahlung or annihilation at one-loop into, say, two gamma-rays, may lead to striking spectral features. Such spectral features are of interest for indirect searches for WIMPs (see, e.g. [12][13][14][15][16] for the Majorana case and specifically [2,9,17,18] for the scalar case). In the case of coupling to (light) quarks, radiative processes involving gluons on top of gammas may be relevant at the time of thermal freeze-out, thus impacting both the effective annihilation cross section and indirect signatures, see e.g. [19,20].
In previous works, annihilation of the particle S through internal bremsstrahlung was only considered in the chiral limit, neglecting the mass of the final state fermions. This was motivated by simplicity, but also by physics. Indeed, it is for light leptons and quarks that radiative corrections are more important, by lifting the helicity suppression. Also, it is in this limit that they lead to most spectacular spectral signatures, with a sharp gamma-ray spectral feature around E γ ∼ m dm when bremsstrahlung is dominated by emission from the intermediate particle 1 , a process called virtual internal bremsstrahlung (VIB) [12]. Instead, in this work, we want to consider the possibility that S couples dominantly to heavy fermions. The prime application would be annihilation into a top-antitop pair, but we will also consider the indirect detection signatures from annihilation into bb and τ + τ − .
The detailed phenomenological analysis of a top-philic candidate, including searches at the LHC and constraints from direct detection searches, is the object of a separate article [21]. In the present work, we specifically focus on more technical aspects of determining the total cross section, taking into account radiative corrections, as well as the spectra into gamma-rays relevant for indirect searches. Concretely, our goal is to keep track of the non-zero quark mass effects, the most important being that the s-wave part of the annihilation cross section into quark-antiquark is helicity suppressed. The issue we will have to face is that the total annihilation cross section is plagued by infrared (IR) divergences, associated to final state radiation (FSR) of soft gluons or gammas. According to the Kinoshita-Bloch-Nordsieck theorem [22,23] the full cross section is free of IR divergence. This involves properly taking into account radiative corrections at a given order in the gauge coupling. For the case at hand, this requires calculating the one-loop corrections to the annihilation cross section SS → ff . Although in principle straightforward, the calculations are involved.
In this note, we give a calculation of the annihilation cross section at next-to-leading order (NLO) following an effective approach advocated in [20], but which they applied to the case of bino-like DM (see also [24] for more general cases ). The expression of the cross section is free of IR divergences and may be applied to the case of annihilation of S in heavy quarks. The main idea behind [20] is to consider separately the emission of FSR and VIB gluons or gammas. The former is dominated by emission from final state fermions, which is the source for infrared divergences of the total cross section. Following the Weizsäcker-Williams approximation, the amplitude for emission of soft modes is obtained by multiplying the leading order (LO), treelevel amplitude by a universal factor (for fermions in the final state). For non-relativistic DM in an s-wave, this tree-level amplitude can be equivalently obtained from an effective contact interaction, which, in the scalar case, is given by the following 5-dimensional operator, where m q is the quark mass, Consequently, the IR divergence can be tackled by taking into account the one-loop correction to the effective interaction of (2). The rest, that is the emission of hard modes, is IR safe, and can be obtained by considering the full, underlying theory. This effective approach simplifies very much the calculations, while capturing the underlying physics, i.e. with limited error compared to a full NLO calculation [20]. The separation between soft and hard modes is implemented by a cut-off on the energy of the emitted gamma or gluon. Such strategy has been used for first calculations of NLO QCD corrections to the Higgs decay, in particular in [25] (see also [26]) and will follow this reference closely. Incidentally, for the soft part, the calculations are precisely the same as in the case of Higgs decay. They differ for the emission of hard modes, for which we will give complete expressions in the case of annihilating scalar dark matter.
The manuscript is organized as follows. In section II, we provide the calculation of the annihilation cross section of a real scalar DM particle into SM fermions through t-channel exchange of a vector-like fermion. We introduce the effective approach of [20] and give explicit expressions from the decomposition of the cross section for soft emission, the associated one-loop corrections, and hard emission, the latter including virtual internal bremsstrahlung. In section III, we study the differential cross sections (with gluon and gamma emission) and implications for indirect detection, in particular for the gamma-ray spectra. We draw our conclusions in section IV. Some lengthy expressions are relegated to the Appendices A and B.

II. TOTAL ANNIHILATION CROSS SECTION
In this section, we first revisit the basics of the model, and the reason why internal bremsstrahlung may be relevant. We then go on with the main steps of the calculation of the cross section for massive final state SM fermions.

A. Leading order annihilation cross section
We consider the amplitudes depicted by the Feynman diagrams in Fig. 1. The annihilation cross section of non-relativistic DM particles is then where v is their relative velocity and N c the number of colours. The helicity suppression (∼ m 2 q ) of the s-wave part of the cross section stems from the fact that, from Eq. (1), the S coupling to SM fermions is chiral while the quark-antiquark pair must have zero total helicity; matching the two requires a chirality flip. Note that the p-wave threshold behavior, ∝ (1 − m 2 q /M 2 S ) 3/2 , of this cross section is similar to that of Higgs decay into SM fermions, as both the scalar DM pair in an s-wave and the Higgs are J P C = 0 ++ . Incidentally, the s-wave part of the cross section can be derived from the low-energy effective interaction in Eq. (2) with as in Fig. 1. We keep the quark mass in the denominator, but assume that the DM particles interact at rest. To compare with [20], we mention that a similar result holds for s-wave annihilation of a pair of Majorana DM. The difference is that the initial state is instead in a J P C = 0 −+ , equivalent to a pseudo-scalar particle. Another, but related difference is that the cross section for scalar DM is d-wave suppressed in the chiral limit m q → 0, while it is p-wave in case of Majorana [2,9]. A well-know consequence of the above is that the cross sections relevant for thermal freeze-out and for indirect detection will differ in general. In particular, in the chiral limit, the LO cross section is suppressed if v 1. This suppression may however be alleviated taking into account radiative corrections [10,11].

B. First look at internal bremsstrahlung
For concreteness, we consider QCD corrections, and so leading emission of a gluon. Provided C F α s → Q 2 α where C F = 4/3 and Q is the SM fermion electric charge, (most of) our results can be applied to radiation of a gamma instead of a gluon. Now, a pair of S in an s-wave can annihilate into a pair of gluons at one-loop or through internal bremsstrahlung, a 3-body process shown in Fig. 2. Although suppressed by powers of α s or phase-space, these radiative processes may play an important role, both for indirect detection and for setting the relic abundance [19,20]. Annihilation into two gluons has been studied in details in [17,18], and this for an arbitrary quark mass. Here we focus on internal bremsstrahlung, taking into account quark mass effects.
The relevant amplitudes are depicted in Figs. 2a to 2c. We will refer loosely to Fig. 2a and 2b as final state radiation (FSR) and to Fig. 2c as virtual internal bremsstrahlung (VIB) respectively. Then, taking the S particles to be at rest, the amplitude for S(k 1 )S(k 2 ) → q(p 1 )q(p 2 )g(k) associated to the VIB diagram is given by and K = k 1 = k 2 ≡ (M S , 0). In this expression, * is a shorthand for the polarization vector * (k) and t a are the representation matrices for the fundamental of SU (3). The FSR amplitudes in Fig. 2a and 2b read altogether The last term in Eq. (6) is perhaps surprising, as one could have expected the combination of propagators D 1 D 2 to arise only from the VIB amplitude, see Eq. (4). Concretely, this term comes from the combination together with D −1 ik = 2i p i ·k. Actually, this term (which, incidentally, does not vanish in the limit m q → 0) is gauge dependent and so must compensate terms from Eq. (4) (notice that this implies that our distinction between FSR and VIB is not clear-cut). The total amplitude Using one verifies that the total amplitude is gauge invariant, k µ M µ IB = 0. To gain further insight, it may be instructive to look at Eq. (9) from an effective interaction perspective. To do so, we consider an expansion of Eq. (9) in r −1 = (M S /M ψ ) 2 assuming M ψ M S . Keeping only the dominant contributions, we get the following three terms, each of which is gauge invariant, The first two terms are ∝ m q . While they cannot be written in terms of local effective operators, they have a simple structure. The first term contains the familiar Weizsäcker-Williams eikonal factor which multiplies the LO amplitude for SS → qq and captures the IR divergences of the total annihilation cross section. The second term is IR finite and its numerator has the structure of a dipole interaction, with F µν = t a F a,µν and σ µν = i[γ µ , γ ν ]/2. Due to cancellations between contributions from both M VIB and M FSR amplitudes, it can be seen that this term comes entirely from the M FSR amplitude in Eq. (6). Both contributions involve a chirality flip, ∝ m q , like the leading order s-wave annihilation amplitude. These contributions are collectively depicted by the diagrams (a) and (b) in Fig. 3 and we refer to them as FSR amplitudes. Finally, the third term in Eq. (11) is local and can be derived from the following dimension eight operator (see diagram (c) in Fig. 3) O already introduced in [27,28]. This term comes from both the M VIB and M FSR amplitudes, but we call it VIB for short, as it reduces to it in the chiral limit m q /M S → 0. Incidentally, as it has no helicity suppression, it may be the dominant contribution to SS annihilation if m q M S . In the limit M ψ M S the situation is clear and simple. Concretely, one could use the amplitude of Eq. (11) to compute the annihilation cross section. The first term leads to IR divergences, but these can be tamed in the usual way, as we will see below. However, here we would like to be more general, first because the large M ψ /M S expansion spoils the spectral feature of VIB, which are most prominent when M ψ and M S are almost degenerate, and second because we have in mind candidates that could annihilate into heavy quarks, in particular the top, so that neglecting m q may not be a good approximation.
Anticipating on the results of the next sections, these considerations are illustrated in Fig. 4 where we depict the typical gluon or gamma-ray spectrum (at the partonic level) ωdN/dω as function of χ = ω/M S for a DM candidate with a strong VIB feature, thus for almost degenerate masses M ψ M S , but also a substantial contribution from FSR. The full spectrum, obtained from the amplitudes of Fig. 2, is shown as the solid (blue) line. The VIB feature is the peak near ω M S . Emission of soft bosons, corresponding to the Weizsäcker-Williams approximation, is shown as a dotted (yellow) line. As expected, it captures the behavior of the full spectrum at low energies ω M S . The two other curves correspond to spectra obtained by using the amplitudes of the effective theory. The dashed (green) curve is obtained from the amplitudes of Figs. 3a and 3b. Compared to the Weizsäcker-Williams approximation, it includes the emission of hard photons or gluons. The spectrum has a sharp edge feature, which is characteristic of FSR [29]. Finally, the dot-dashed (red) curve encompasses all the effective amplitudes of Eqs. (11)- (13). While the Weizsäcker-Williams approximation reproduces well the emission of soft gluons or gammas, the effective operator in Eq. (15), corresponding to the amplitude of Eq. (13), fails to fully reproduce the VIB spectral feature.

C. Radiative corrections
For the purpose of probing DM through indirect detection, we aim at determining the spectrum of quark and gluons emitted when DM annihilates through internal bremsstrahlung, dσv qqg /dω where ω is the gluon energy. The integrated cross section is also relevant for determining the relic abundance of the DM particle [19]. However, for finite quark mass, its expression suffers from IR and collinear divergences. The recipe to address these divergences is standard and involves computing not only the 3-body process, but also the one-loop corrections to the 2-body annihilation. Then, the Kinoshita-Lee-Nauenberg theorem (or Bloch-Nordsieck for QED) states that, order by order in the gauge couplings, IR divergences from phase-space integration, which in our case are O(α s ), are cancelled by those from loop corrections. Thus, in principle, we should compute the one-loop corrections depicted in Fig. 5. While such calculations have been performed for minimal supersymmetric candidates (see e.g. [30]), it is another matter to do so for a simplified model. Instead, inspired by the strategy of [20], we will separate the problem into the emission of soft and hard gluons. For emission of soft gluons, we will use an effective interaction to control and cancel the IR divergences that affect the cross section for soft modes, while keeping as much as possible the full, UV complete amplitudes to capture the VIB spectral features. The matching between the two regimes will be controlled by a cut-off on the energy of the emitted gluon, ω 0 . In equation, the total, next-to-leading order cross section we aim at will read where σv LO ≡ σv qq . Both ∆σv| soft and ∆σv| hard depend on the matching energy ω 0 but their sum does not. In this expression where λ is a fictitious mass of the gluon, introduced to regularize the cross section obtained by integrating over soft modes. The tilde onσ is there to mean that this cross section is unphysical, as it diverges for λ → 0. As usual, the λ dependence requires to take into account one-loop corrections to the LO cross section, which are computed using the effective theory.
The λ dependence will cancel in the sum of the two contributions. Clearly, the main advantage of this down-to-earth approach is that we will only need to calculate the one-loop corrections depicted in Fig. 6 to cancel infrared divergences. Incidentally, as this coupling has precisely the same structure as the Higgs coupling to SM fermions, much of the underlying physics is the same as that discussed in [25,26]. What is specific to the DM scenario is the emission of gluons by the vector-like mediator.

Soft gluon emission
We first consider the annihilation into a pair of massive SM quarks with the emission of a soft gluon. By this, we mean a real gluon with energy ω = | k| ≤ ω 0 , where ω 0 is a cut-off energy, which we take to be small compared to other characteristic mass scale in the theory, ω 0 {M S , M ψ , m q }, but larger than Λ QCD . In that limit, we describe the emission of a soft gluon using the Weizsäcker-Williams approximation, This differs from the first term in Eq. (11) by the factor 1/∆ with ∆ = 1 + M 2 ψ /M 2 S − m 2 q /M 2 S ≡ 1 + r − z, which stems from neglecting the soft gluon 4-momentum in the propagator of the mediator.
Integrating over phase space for final state fermions, we get the following differential cross section, valid for emission of a soft gluon of energy ω, where C F = 4/3 and χ = ω/M S . This expression involves the velocity of the final state quarks in the rest frame of the qq system (see e.g. [31]) To regulate the IR divergence that will arise when integrating Eq. (19) over the gluon energy, we have introduced a fictitious mass λ ω 0 for the soft gluon, which appears through µ = λ 2 /4M 2 S . Keeping only the leading terms in the limit λ → 0, the integrated cross section for emission of a soft gluon in the energy range λ ≤ ω ≤ ω 0 is given by By construction, this cross section is proportional to the leading order cross section σv qq , which corresponds here to the s-wave part of Eq. (3). This expression can be compared (and agrees) with that of [25] for decay of the Higgs, with which it shares the IR divergence term ∝ log(ω 2 0 /λ 2 ).

Virtual one-loop corrections
Physical cross sections should be free of IR divergences as λ → 0. To obtain a IR finite result, we need take into account the contributions of O(α s ) virtual one-loop corrections to the leading order cross section into qq. This stems from the interference term between LO and the one-loop corrections At one-loop, the IR divergent contributions come from the vertex correction and final state fermion wave-function corrections, depicted by the diagrams of Fig. 6. Using dimensional regularisation in D = 4 − 2 , the virtual correction to the effective vertex is given by (see [25] for comparison) which is both UV and IR divergent.
According to the LSZ reduction formula [32], we must also take into account the O(α s ) correction from the one-shell wave-function of the final state quark and anti-quark, with 2 Thus, the one-loop unsubstracted correction to annihilation into a quark-antiquark pair is Comparing the second term of this expression to the first term in Eq. (21), we see that, adding the O(α s ) one-loop corrections to the tree level cross section to the cross section for emission of soft gluons, the dependence on the fictitious gluon mass (i.e. the terms in log(λ 2 )) disappears [25], leaving only the dependence on the cut-off on the energy of the emitted gluon ∝ log(ω 2 0 ). The resulting expression has still a UV divergence, which must be appropriately cancelled. The renormalization prescription used in [20] is the same as the one advocated in [25] in the case of QCD corrections to Higgs decay into quarks. In this case, since the current quark mass stems from Yukawa coupling to the Higgs, the counter-term is that for quark mass renormalization, This term clearly cancels the UV divergent part of Eq. (24) but any other prescription would only differ from this choice by a constant term. Which choice one makes does not matter. Indeed, for fixed particle masses, the only free parameters in the annihilation cross section (meaning here at O(α s )) is the Yukawa coupling y f . Its value is fixed by matching to the cosmic relic abundance. All other parameters being kept fixed, a different renormalization prescription just amounts to fixing y f to a (slightly) distinct value. For definiteness, here we use the same prescription of [20] to renormalize our effective theory. 3 Doing so, we get the one-loop correction to the LO cross section 2 Intuitively, the need for wave-function renormalization stems from the fact that IR divergences come both from squaring the (a) or (b) amplitudes of Fig.2 and from their interference, the loop correction to the vertex cancelling the IR divergence from the latter. 3 The correspondence with the problem of QCD corrections to SM Higgs decay into quarks rests on the use of an effective vertex. In principle, a procedure we could follow is to match our effective theory with the more complete theory at the scale(s) at which one integrates out the heavy degree(s) of freedom. For instance, the one-loop corrections include the box diagram depicted in Fig. 5, which has a better UV behavior than in the effective theory. With the mass of vector-like quark, M 2 ψ acting as a cut-off, divergent terms 1/ could actually correspond to ∝ log(M 2 ψ /M 2 S ) contributions (see e.g. [32]). Using the matching procedure would only introduce minor corrections (at least compared to the major impact of taking into account bremsstrahlung). See [20] for further considerations on errors from using the effective approach.
Adding this contribution to Eq. (21) gives ∆σv| eff soft , which depends on the cut-off energy ω 0 but not on λ, i.e.
where the dots correspond to terms that are O(ω 0 0 ).

Hard gluon emission
It remains to determine the spectrum of hard gluons and their contribution to the total NLO cross section. For this, we use the full theory, including the effects of the vector-like particle, from the amplitudes of Fig. 2. Since we will put a cut-off on the energy of the gluon, no gluon mass term is required. The calculations, although cumbersome, are straightforward. The differential cross section can be written as In this expression, we have separated the terms that are divergent in the limit χ = ω/M S → 0 from those that are regular; the latter are collectively expressed as the function S 0 (χ), whose expression is extraordinarily long and not particularly illuminating; its full expression is given for the sake of reference in Appendix A. It contains in particular contributions that reduce to the known expression of virtual internal bremsstrahlung in the limit of massless quarks, see Eq. (B3). For finite quark masses, it also includes hard emission from final state quarks and interference terms between the latter and VIB. Integrating Eq. (28) over ω ≥ ω 0 , we get The first term in this expression involves the cut-off energy ω 0 . Adding ∆σv| full hard to the soft contribution gives a result that is independent of ω 0 . Our final expression for the cross section for s-wave annihilation is then which is one of our main results.

D. Discussion
The expression of Eq. (30) is free of infrared divergences and thus is a priori useful to determine the relic abundance of S particles and its indirect signatures. It is however too complex to be practical. Furthermore, and despite its complexity, it still has some limitations. First, the expression in Eq. (30) diverges close to threshold for quark-antiquark pair production. Second, due to collinear divergences, it is pathological in the opposite limit, m q M S , or β 0 → 1. The same issues appear in the decay of the Higgs, and the cure are, for the essential, the same, so we only discuss briefly how they should be addressed in principle, and then what we do in practice (see also [20] for a similar discussion). By the same token, we discuss in this section approximations that can be used to take into account the relevant aspects of radiative corrections to the DM model without resorting to the full complexities of the O(α s ) cross section.
We first dispose of the problem posed by the divergent behavior of the NLO cross section close to threshold for fermion-antifermion production. For M S m q , corresponding to β 0 = 1 − m 2 q /M 2 S → 0, the annihilation cross section behaves as The O(α s ) terms arise from expanding near zero velocity β 0 → 0 the two Spence functions in Eq. (30). As shown in [25], this singular behavior is spurious, as the cross section should be in a p-wave quark-antiquark final state, ∝ β 3 0 . It can be traced entirely to the virtual correction to the effective vertex in Fig. 6 or, in other words, from Eq. (26). Physically it signals the tendency to form a bound state, so in principle one should sum an infinite number of diagrams or, below threshold, take into account a possible quarkonium bound state [33]. We consider this to be beyond our scope. Now, for bb (or τ + τ − ) the mass of the dark matter and its charged partners is too small to provide viable DM scenario, i.e. we are always far from threshold. A case of potential interest is a top-philic scenario, with annihilation of DM into top-antitop pairs (gg) above (resp. below) threshold [21]. As the threshold for top-antitop corresponds to a very specific and narrow region of the parameter space, we may take the LO annihilation cross section as a proxy for the true cross section, σv qq (1 + K), with a K-factor that we estimate empirically (see [21]).
The cross section is also pathological in the opposite limit M S m q or β 0 → 1. In particular, the part that is proportional to σv LO receives a large negative logarithmic contribution [25,26], The dots represent the terms that are regular for small z = m 2 q /M 2 S . Importantly, they reduce to the cross section for the pure VIB process in the limit z → 0, see Eq. (B3), a regime in which σv qq → 0 so, a priori, the new logarithmic term is harmless for most of the parameter space. If not, one must in principle re-sum large logarithmic contributions to the cross section. A simple recipe to address this is to notice that, as in the case of the Higgs decay, the leading log term in Eq. (32) is precisely the one-loop O(α s ) correction to the quark mass operator [20]. A convenient way to have a regularized expression for the NLO cross section is thus to subtract from Eq. (30) the logarithmically divergent term in Eq. (30) and replacing in the expression for σv qq the quark mass parameter m q by the running mass [25]  Now that we have a complete and reliable description of the NLO effects, we go on discussing controllable approximations that can be made to obtain simple and practical expressions for the annihilation cross section with FSR and VIB. In Fig. 7a  VIB . For lower DM masses, the cross section is dominated by the chirally suppressed component ∝ σv qq . For increasing r = (M ψ /M S ) 2 , the VIB contribution is relatively suppressed, as it scales like σv VIB ∝ r −4 while σv qq ∝ r −2 , see Eq. (11). The cross-over between the two regimes may be read from Fig. 7a where the dotted lines correspond to the leading order σv qq cross section. Crossing of σv qq and the VIB cross section occurs roughly for M S /m q ≈ r 2π/0.21C F α s with m q ≡ m t in the figure and the factor 0.21 ≡ 7/2 − π 2 /6, see Eq. (B4).
The main lesson is that the NLO cross section σv NLO is reasonably approximated by the following simple expression, in which the leading 2-body cross section σv qq is added to the VIB cross section in the massless limit, σv NLO ≈ σv qq + σv This is little surprise, as our NLO calculation is built upon the effective operator in Eq.
(2), which should lead to the dominant contribution to the cross section when VIB emission may be neglected. Our calculations show that interference terms play little role, even when VIB is relevant. More difficult to assess is the error made by using the effective theory instead of the full one-loop amplitudes depicted in Fig. 5, but it should not be more than a few percent, based on the experienced gained in the case of Majorana dark matter [20]. In [21] we have used the approximation of Eq. (33), which is easy to implement in numerical codes for DM abundance calculations, like MicrOMEGAs [34], just adding the massless 3-body process associated to VIB, as in Eq. (B4), when it is relevant.

A. Differential cross sections
The total annihilation cross section is not only relevant for determining the relic abundance; it also sets the scale for indirect searches. But for the latter purpose, we need an handle on the differential cross sections, both in gamma-rays and in gluons, so as to build their spectra, where here ω stands for the energy of the emitted gamma or gluon and σ is the total annihilation cross section into this channel. This is a priori straightforward but in practice, things are more complicated, as to assess the indirect signature from DM annihilation, say into gamma-rays, what we need is to take into account both the contribution of gamma-ray produced directly by the annihilation process (i.e. prompt photons or gluons produced at the partonic level) and those that will emerge from the process of fragmentation into hadrons from both the final state quark-antiquark and the gluon from bremsstrahlung. This requires to resort to Monte-Carlo simulation tools, like Pythia [35]. The way we handle this is discussed in the next section. Here we focus on the differential cross section at the partonic level and on the possible simplifications one may use to get approximate results. Eq. (28) is in principle all we need to determine the spectrum of prompt gluons or gammas. Its expression is however not very convenient, as it involves all the terms given in the Appendix A. In practice, we may rely on rather simple approximations. The first is to replace the full expression of Eq. (28) by In this expression, the first term is the differential cross section for emission of a gamma or gluon using the effective theory, that is the amplitudes in Figs. 3a and 3b. Basically, this contribution is equivalent to the process of Higgs decay studied in e.g. [25]. The second term is the differential cross section for VIB calculated from the amplitude of Fig. 2c, where the subscript refers to the limit of a massless quark, see Eq. (B2) in Appendix B. The rational is that VIB is mostly relevant in the limit m q M S (and provided the mediator is not much heavier than the DM particle). Less obvious is that this expression works pretty well for intermediate regimes. This is illustrated in Fig. 8 where, for the specific choice of M S = 2 TeV, M ψ = 2.4 TeV and m q = m t , we show the differential cross section for prompt gamma-ray emission based on the full calculation (solid black line) compared to the one obtained from the expression in Eq. (35) (dotted blue line). The two curves are clearly very close to each other. We should emphasize that to get a better matching we have shifted the energy ω in the differential cross section for the VIB contribution to ω + m 2 q /M S , so as to take into account the finite quark mass effect on the end-point of the gamma-ray spectrum. Modulo this, the correspondence between the two expressions is surprisingly good. This simple decomposition into FSR plus VIB suggests a further approximation, which is to express FSR in terms of the LO 2-body annihilation cross section time a factor that takes into account the emission of a gamma-ray or gluon by the final state quarks, In this expression, the factor F is the standard splitting function for emission of gamma by a final state fermion, In case of emission of a gluon, one must of course replace the factor αQ 2 by α s C F in Eq. (36). The splitting function captures the collinear divergences that arise in the limit m q ω ∼ M S , hence for hard emission [29]. Integrating over ω down to a cut-off energy ω 0 leads to the characteristic Sudakov double logarithmic divergence ∝ log(m 2 q /M 2 S ) log(m 2 q /ω 2 0 ) which one can read in the expression of Eq. (29). We have checked explicitly that the expression of dσ FSR /dω reduces to the expression of Eq. (36) in the limit m q ω ∼ M S . The last term in Eq. (36) is non-universal, but keeping it gives a better matching to the exact result. This is shown in Fig. 8 where the dashed red line correspond to the differential cross section obtained by summing the contribution from Eq. (36) to VIB (in the massless limit). The matching is not as good as in the previous approximation, Eq. (35), but for practical applications, it is much simpler and also more physically transparent. It also paves the way to the determination of the final gamma-ray spectrum from hadronisation.

B. Gamma-ray spectra
Our final goal is to obtain the gamma-ray spectrum from annihilation of a S dark matter candidate, taking into account FSR and VIB emission both of gluons and gammas. One obvious way to proceed is to implement directly our differential 3-body process into Pythia, by first building a Monte-Carlo distribution using e.g. CalcHEP [36] which then can be hadronised using Pythia. This is the strategy we have used in the past for the case of coupling of the vector-like portal to light quarks [19], in which case we could altogether neglect the quark mass and so the 2-body annihilation process. Clearly this strategy also applies to heavy quarks in the limit m q M S → 0, provided the mediator is not much heavier than the DM candidate, see Fig.7a. In general, however, the fermion mass and the associated FSR may not be neglected. The problem is that infrared and collinear singularities associated to FSR lead to sharp peaks in the Monte-Carlo distribution, which, for numerical convergence purpose, require to introduce a cut-off on the energy and, a priori on the emission angle of emitted gluon or gamma to obtain a reliable numerical output. Doing so for each DM candidate is cumbersome and CPU time consuming. Now, as discussed at length in the previous sections, the total cross section(s) can be decomposed into soft and hard gluon or gamma emission. So, the next idea is to implement separately soft and hard emissions. The latter is free of infrared and collinear divergences and so can be reliably calculated using first CalcHEP to simulate the 3-body process and then Pythia for hadronisation. One still has to deal with divergences from soft emission. This however, is a radiative correction to a 2-body process that can be handled directly by Pythia, which has built-in simulation of FSR emission, including splitting of gluons, through Sudakov factors [35].
The separation between soft and hard modes however amounts to consider non-inclusive cross sections, which typically have Sudakov double logarithm divergences. So a question is where to put the cut-off between soft and hard modes. Indeed, Eq. (27) has a Sudakov double logarithm ∝ log(m 2 q /M 2 S ) log(m 2 q /ω 2 0 ) typical of non-inclusive cross sections and which can be traced to the behavior of the differential cross section in the regime of collinear emission, so we should avoid taking the cut-off at low energies where the log diverges. Instead we take it such that the FSR and VIB match smoothly. This is depicted in Fig. 9 where for concreteness we show the spectrum of gamma-rays (at the partonic level) for a candidate with a strong VIB feature. The blue, solid line, is the spectrum obtained using the full theory, and the red, dashed line the one from final state radiation using the 2-body effective theory. The cut-off ω c that separates high and low energy emission may be chosen as the energy where the differential cross section in the effective theory departs from the one in the full theory. Once ω c is determined, we can generate hard events (with CalcHEP and then Pythia for hadronisation) using the full theory on one hand, and soft events coming from qq with FSR events (using Pythia only) on the other hand. Using the theoretical expression for the differential cross section in the full theory, it is easy to generate (using CalcHEP) gluons or gammas with an energy larger than the cut-off ω c . However, within Pythia8 we did not find a simple way to extract partonic level events with gluons or photons with energy ω < ω c . Instead, we have first generated with Pythia8 the complete gamma-ray spectrum from the 2body process, including FSR. Next, we have simulated with CalcHEP the distribution of hard gluons and gammas events with ω > ω c using the analytical cross section for the 2-body process with FSR. After hadronisation of this part with Pythia8, we have subtracted the resulting gamma-ray spectrum from the one obtained in the first step, to get only the soft part of the gamma-ray spectrum. Finally we have added back the hard part from the full theory, which includes VIB effects, to get the final gamma-ray spectrum. As S annihilation proceeds through several final states (qq + γ + g, γγ, gg), the total photon spectrum, dN tot γ /dE γ , is finally given by the sum of the photon spectrum originating from each final state (dN tt+γ+g γ /dE γ , dN γγ γ /dE γ and dN gg γ /dE γ respectively), weighted with their respective branching ratios. To distinguish prompt emission gluons and gammas of energy ω, we use E γ for the energy of the gamma-rays in the final spectrum. The normalization is with respect to the total full inclusive annihilation cross section.  We may now consider, for the sake of illustration, different SM fermionic final states for which the fermion mass may be a relevant parameter. Concretely, we consider the τ + τ − as well as the bb (Fig. 10) and tt (Fig. 11) channels. We show spectra for fixed (m q /M S ) 2 ≡ z ratios, z = 0.1 2 , r = 1.2 2 , adjusting both the dark matter mass M S and mediator mass M ψ . Each spectrum is generated following the procedure depicted above, with a specific cut-off separating emission of soft and hard gluons or gammas. We have included for completeness the contribution from annihilation at one-loop into two gluons and into γγ, assuming a gamma-ray detector with resolution ∆E γ /E γ = 10%. The parameters of the DM model are chosen so as to illustrate the possible presence of a feature in the final gamma-ray spectrum, not taking into account other possible constraints (relic abundance, direct, indirect and collider constraints).
In that respect, most relevant are the two plots of Fig.11 with coupling of DM to the top quark, which corresponds to actual DM candidates. The phenomenology of such DM candidates are discussed in details in [21].

IV. CONCLUSIONS
We have studied further radiative corrections to a simple DM scenario, in the form of a real scalar particle annihilating into SM fermions through a heavy vector-like fermion. This topic has been already covered in several phenomenological studies [2,9,[17][18][19] that focused on coupling to light quarks or leptons. Of particular relevance in this regime is the helicity suppression of the 2-body annihilation cross section, akin to p-wave suppression of Majorana dark matter annihilating into SM chiral fermions [8,10,11]. This implies in particular that radiative corrections, in the form of one-loop annihilation into two photons (or gluons) and so-called virtual internal bremsstrahlung may play a significant role both in determining the relic abundance and for indirect searches. Due to infrared and collinear divergences that affect bremsstrahlung of massless gauge bosons, the extension of these results to heavy quarks (or annihilation into a τ + τ − pair) poses some technical problems, which we have tried to overcome in the present work. The motivation was manifold, but the main aim was to try and test simple approximations that could be then applied for phenomenological studies, in particular to the case of top-philic coupling, a topic of much interest, in particular from the perspective of simplified DM searches at the LHC. Such study has been the object of a separate publication [21] (see also [37][38][39][40]). Here, we focused on technical aspects of the calculations. In particular, following a proposal of [20], we have adapted an effective approach suited for emission of soft gamma or gluons and that circumvent several unnecessary steps in the regularization of infrared divergences. From an effective approach perspective, much of the calculations map to the equivalent problem of radiative corrections to Higgs decay into SM fermions [25,26]. This approach, which is quite systematic and pedestrian, has also a pedagogical value. In particular, it illustrates in a simple framework how the cancellation of infrared and collinear divergence takes place in the calculation of a total cross section, in agreement with standard quantum field theory theorems. At the end of the day, our main results include a full, explicit, albeit cumbersome expression for the total annihilation cross section for S DM into SM fermions at NLO in α s (or α). Building on this, we have studied simple approximations, both to the total and differential cross sections. The take home lesson is that the simple approximations discussed in the bulk of this article are well suited for phenomenological studies, as discussed in [21].
In the limit r → 1, In the opposite limit, r = M ψ /M S 1, VIB may be neglected. With a factor of σv qq instead of the tree level decay rate of the Higgs, we recover the expression of [25] (taking into account their errata) where (B6)