Squark-pair annihilation into quarks at next-to-leading order
S. Schmiemann J. Harz B. Herrmann M. Klasen K. Kovařík


LAPTH-008/19
MS-TP-19-07
TUM-HEP-1194-19





 saskia.schmiemann@uni-muenster.de
 
	Institut für Theoretische Physik, Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Straße 9, D-48149 Münster, Germany
  
  

 julia.harz@tum.de
 
	Physik Department T70, James-Franck-Straße, Technische Universität München, D-85748 Garching, Germany
  
  

 herrmann@lapth.cnrs.fr
 
	Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LAPTh, F-74000 Annecy, France
  


 michael.klasen@uni-muenster.de
 
	Institut für Theoretische Physik, Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Straße 9, D-48149 Münster, Germany
  


 karol.kovarik@uni-muenster.de
 
	Institut für Theoretische Physik, Westfälische Wilhelms-Universität Münster, Wilhelm-Klemm-Straße 9, D-48149 Münster, Germany
  



 
              


The Minimal Supersymmetric Standard Model (MSSM) is under intense scrutiny at the LHC and in dark matter searches. Interestingly, scenarios with light squarks of the third generation remain not only viable, but also well motivated by the observed Standard-Model-like Higgs boson mass and dark matter relic density. The latter often requires important contributions from squark pair annihilation. Following up on previous work, we present in this paper a precision analysis of squark pair annihilation into quarks at next-to-leading order of QCD including Sommerfeld enhancement effects. We discuss all technical details of our one-loop, real emission and resummation calculations, their implementation in the precision tool DM@NLO, as well as the numerical impact on the annihilation cross section and cosmological relic density in phenomenological MSSM scenarios respecting in particular current LHC constraints. We demonstrate that including these radiative corrections leads to substantial shifts in the preferred parameter regions by up to 20 GeV.


12.38.Bx,12.60.Jv,95.30.Cq,95.35.+d
                             

                              
Squark-pair annihilation into quarks at next-to-leading order
    K. Kovařík
    December 30, 2023
=============================================================












§ INTRODUCTION



Strong evidence for the existence of dark matter, along with the fact that neutrinos are massive, are compelling signs for the need of physics beyond the Standard Model. Even though dark matter still evades direct detection in Earth-based experiments such as LUX <cit.> and XENON1T <cit.>, there is overwhelming evidence from cosmological data such as the cosmic microwave background that dark matter exists in the Universe. Moreover, the relic density of cold dark matter (CDM) has been determined to the unprecedented precision of

    Ω_ CDMh^2  =  0.1200 ± 0.0012

as measured by the Planck satellite and interpreted within the ΛCDM cosmological model <cit.>. The indicated uncertainty corresponds to the 68% confidence level, and h stands for the present Hubble expansion rate H_0 in units of 100  km  s^-1  Mpc^-1.

The Standard Model does not contain any suitable candidate for dark matter with the required properties.
The leading candidate therefore remains a weakly interacting massive particle (WIMP), which leads to the correct relic density via the freeze-out mechanism. However, alternative candidates and mechanisms do exist, e.g. in the form of the freeze-in mechanism <cit.>).
The Standard Model is therefore extended to include new particles which provide the required dark matter candidate. The new particles are usually protected from decaying by introducing an ad-hoc Z_2-symmetry, where all new particles are Z_2-odd and the Standard Model particles are Z_2-even. One such model, which was actually not introduced to address the existence of dark matter, is the Minimal Supersymmetric Standard Model (MSSM), where the conserved Z_2-symmetry is the R-parity. In most MSSM scenarios, the lightest supersymmetric particle (LSP) is the lightest neutralino χ̃^0_1, which is stable and interacts only weakly. The lightest neutralino is an extremely well-studied candidate for cold dark matter.

The theory prediction for its relic abundance, Ω_χ̃^0_1h^2, is related to the number density n_χ of the neutralino, which can be computed solving the Boltzmann equation <cit.>

    n_χ/ t =  
            -3 H n_χ - ⟨σ_ annv ⟩[ n^2_χ - (n_χ^ eq)^2 ]  ,

where H denotes the (time-dependent) Hubble parameter, n_χ^ eq the number density in thermal equilibrium, and v the Møller velocity of the annihilating particles <cit.>. All specifics about the interaction of dark matter with other particles in the chosen particle physics model is contained in the annihilation cross section σ_ ann, which accounts for all possible annihilation and co-annihilation processes. Its thermal average can be expressed as

    ⟨σ_ ann v ⟩ = ∑_i,j⟨σ_ij v_ij⟩n_i^ eq/n_χ^ eqn_j^ eq/n_χ^ eq ,

where the double sum runs over all Z_2-odd particles of the theory. The ratios of the equilibrium number densities are proportional to

    n_i^ eq/n_χ^ eq ∼ exp{ -m_i - m_χ/T} ,

T being the temperature. From this last equation, it becomes obvious that annihilations involving particles other than the lightest neutralino are suppressed if these particles are heavy compared to the neutralino. On the other hand, co-annihilations with the next-to-lightest supersymmetric particle (NLSP) will be important or even dominant if the mass difference is rather small. Typical examples in the MSSM are co-annihilations of the neutralino with a scalar top quark or a scalar tau lepton. For smaller mass differences between the LSP and NLSP, even pair-annihilations of the next-to-lightest particle contribute in a sizeable manner, and can even become dominant in the total annihilation cross section. In case there are more than two almost mass degenerate particles, (co-)annihilations between all particles have to be taken into account.

In the present paper, we focus on the case, where the masses of one or two squarks, the lightest stop and/or the lightest sbottom, are close to the neutralino mass. A case of a light scalar top quark is very well motivated. First, a light scalar top is necessary to achieve the electroweak baryogenesis in the MSSM <cit.>. Second, scenarios with light scalar tops satisfy the experimental constraints from LHC searches and can also contribute to a successful prediction of the mass of the lightest Higgs boson in the MSSM <cit.>.

The relic density of dark matter in scenarios with a light stop which is almost mass degenerate with the lightest neutralino is very sensitive to the mass difference of the two particles. Any small shift in the predicted relic density can cause a large shift of the parameter region where the relic density is compatible with the experimental limits given by Eq. (<ref>). In this analysis we focus on next-to-leading supersymmetric QCD (SUSY-QCD) corrections to the corresponding annihilation and co-annihilation cross sections in scenarios with a light scalar quark. These corrections have the potential to significantly modify the annihilation cross section and thereby also the relic density.

The impact of such radiative corrections of order α_s on the relic density has been demonstrated for gaugino pair annihilation into quarks <cit.>, gaugino-squark co-annihilation into final states containing a quark <cit.>, and squark-antisquark annihilation into electroweak final states <cit.>. Moreover, electroweak corrections to neutralino annihilation have been presented <cit.>, leading to similar conclusions concerning their impact on the relic density.

Including radiative corrections to the total annihilation cross section not only shifts the parameter regions corresponding to the correct relic density, but it also reduces the theoretical uncertainty of the relic density prediction. The theoretical uncertainty from scale and scheme variations on the annihilation cross section and the neutralino relic density has been evaluated for specific subclasses of processes in Ref. <cit.>. 

After the work presented in Refs. <cit.>, with the present paper, we make a first step towards completing the missing processes sensitive to radiative corrections of order α_s. More precisely, we present such corrections for squark-squark annihilation into quark-quark pairs. The discussion of squark-antisquark annihilation into quarks and gluons is left for forthcoming publications. 



In the following, in Sec. <ref>, we start by discussing the phenomenological importance of the processes under consideration in this work. We also present two reference scenarios featuring important contributions of the processes of our interest. In Sec. <ref>, we then detail the analytical calculation of the radiative corrections. We discuss in particular points that are beyond the discussion presented in Refs. <cit.> and analyze the impact that the radiative corrections have on the corresponding cross sections. The impact of the corrections on the relic density in the two reference scenarios is presented in Sec. <ref>. Our conclusions are given in Sec. <ref>.







§ PHENOMENOLOGY OF SQUARK ANNIHILATION



The analysis presented in this paper concentrates on the contributions from squark-pair annihilation to the total annihilation cross section σ_ ann of neutralino dark matter. We investigate scenarios in the phenomenological MSSM (pMSSM), where the processes

    t̃_1 t̃_1     ⟶  t t  , 
    b̃_1 b̃_1     ⟶  b b  , 
    t̃_1 b̃_1     ⟶  t b  ,

play an important role. Supersymmetry and the MSSM in particular have been extensively tested by searches at the Large Hadron Collider (LHC) and at experiments aiming at the detection of direct signals from elastic collisions of dark matter with heavy nuclei such as XENON1T. In order to take into account the most important experimental constraints from the searches for supersymmetry, we use the results of an analysis performed by the ATLAS collaboration in the light of recent searches at the LHC <cit.>.[A similar study of the pMSSM has been conducted by the CMS collaboration <cit.>.] The ATLAS analysis is performed in the pMSSM with 19 parameters (defined at the SUSY scale) and is based on a sample of 5· 10^8 parameter points. Applying constraints from ATLAS SUSY searches, electroweak precision observables such as Δρ and (g-2)_μ, flavor observables such as b→ sγ, and requiring that the neutralino is the LSP and a dark matter candidate with the relic density less than 0.1208 (for details on further contraints see Ref. <cit.>) leads to a subset of about 300.000 viable points. We have analyzed this subset in order to examine in which regions of parameter space the above processes contribute significantly.




In order for the contribution from the annihilation of third-generation squarks to the total dark matter cross section to be significant, one (or more) scalar quarks have to be almost mass degenerate with the lightest neutralino. This is not an unnatural requirement because a light scalar top quark is necessary to explain the measured mass of the Standard Model Higgs boson within the MSSM. Moreover, scenarios where a scalar top is almost mass degenerate with the lightest neutralino are quite frequent, as the mass degeneracy gives rise to different topologies in collider searches making their testing more challenging and their exclusion less likely. Another aspect of scenarios where scalar quarks are the next-to-lightest supersymmetric particle (NLSP) is that the lightest neutralino is mostly bino-like. Higgsino-like and wino-like lightest neutralinos mostly lead to scenarios with other gauginos being the NLSP. A consequence of the lightest neutralino being bino-like is that the annihilation of dark matter is typically not efficient enough for the relic density to reach the value determined by the Planck collaboration given in Eq. (<ref>). Therefore in scenarios with bino-like neutralino dark matter, some enhancement mechanism is needed for them to be consistent with the relic density measurements. As we will discuss below, in the scenarios analyzed here, the enhancement comes from the presence of LSP-NLSP co-annhilations as well as NLSP annihilations. 







 §.§ Reference scenarios



As mentioned above, the numerical part of the present study will be based on two reference scenarios inspired by the findings presented in Ref. <cit.>. More precisely, we will focus on two pMSSM scenarios, whose most relevant soft-breaking parameters and particle masses are presented in Tab. <ref>. It is to be noted that, although the input soft mass parameters of the two scenarios are identical to those of two actual scenarios given in Ref. <cit.>, the resulting physical masses slightly differ from those associated with the ATLAS study due to the fact that we are using a different computational setup. The actual shift in the physical masses is small so that all experimental constraints are still satisfied and the phenomenology is not altered. 

Both scenarios feature bino-like neutralinos, the bino mass parameter M_1 being smaller than the wino and higgsino mass parameters M_2 and |μ|. The key parameters of the third-generation squarks of our interest are the “left-handed” stop and sbottom mass parameter M_q̃_L, and the “right-handed” stop and sbottom mass parameters M_t̃_R and M_b̃_R. In both scenarios, squarks of the first and second generation, the sleptons, and other electroweak gauginos are heavier such that they do not influence the phenomenology discussed here. 

In our setup, starting from the soft-breaking terms defined at the scale Q_ SUSY indicated in Tab. <ref>, we obtain the physical mass spectrum using the spectrum generator SPheno 3.3.3 <cit.>. The mass spectrum is then handed over to micrOMEGAS 2.4.1 <cit.> making use of the SUSY Les Houches Accord 2 <cit.>. In addition to the actual value of the relic density, micrOMEGAs also provides the contributions of all individual channels contributing to σ_ ann given in Tab. <ref> for the two chosen reference scenarios.

As can be seen, the processes given in Eqs. (<ref>) to (<ref>) contribute in a significant manner for both scenarios. More precisely, in Scenario I, the scalar top pair annihilation is the second most important process and together with the processes previously analysed in Refs. <cit.> makes up more than 45% of the total annihilation cross section. In this scenario, the mostly “right-handed” scalar top t̃_1 is the NLSP, and the mass difference between the lightest neutralino and the NLSP is about 20 GeV. Moreover, the process t̃_1 t̃_1 → t t is enhanced by the relatively low gluino mass. Scalar bottom quarks are heavy in this scenario, such that the corresponding annihilation channels are negligible. The process t̃_1 t̃_1 → tt amounts to about 30% of the annihilation cross section, while neutralino-stop co-annihilation <cit.> accounts for about 15%, such that our next-to-leading SUSY-QCD corrections affect almost 50% of the total annihilation cross section. 

The importance of different relevant contributions to the total annihilation cross section in and around this scenario is shown in the first four plots of Fig. <ref>. The part of the parameter space where the lightest neutralino is not the LSP and hence also not the dark matter candidate is indicated in grey. Different shades of green indicate the relative importance of several NLSP annihilation and LSP-NLSP co-annihilation processes. We see that in general the co-annihilations are most important when the mass splitting between the LSP and the NLSP (in our case the lightest neutralino and the lightest scalar top) is larger, i.e. about 150 GeV, and the dominant contribution shifts to the NLSP annihilations as the mass splitting gets smaller. The parameter region where the dark matter relic density is within 2σ of the experimental value given in Eq. (<ref>) is highlighted in all plots in orange. The neutralino relic density is computed using micrOMEGAs. For scenarios which contain lighter neutralinos and lighter stops than the reference Scenario I, the total annihilation cross section is dominated by the same processes as in the reference Scenario I, and the region where the relic density agrees with the experimental measurement follows an almost straight line where the mass splitting is constant. As we move along the region with the correct relic density towards scenarios with heavier LSP and NLSP, we reach a point where the LSP and NLSP are similar in mass to the light gluino (m_g̃=1495.5 GeV). In these scenarios, we have three particles with almost degenerate masses and gluino annihilations and co-annihilations with the stop dominate the total annihilation cross section.    


The situation is different for Scenario II. Here, the “left-handed” mass parameter M_q̃_L is much smaller than the “right-handed” masses M_t̃_R and M_b̃_R, such that the relevant physical states t̃_1 and b̃_1 are mainly “left-handed” with almost degenerate masses. The mass difference between them and the lightest neutralino is about 30 GeV. As a consequence, processes containing both t̃_1 and b̃_1 contribute to the annihilation cross section σ_ ann, as can be seen in Tab. <ref>. The three processes of our interest contribute to more than 50% of the total annihilation cross section. As we shall discuss later in Sec. <ref>, the mixed annihilation t̃_1 b̃_1 → t b dominates as compared to stop-pair or sbottom-pair annihilation. In the last four plots of Fig. <ref>, we show the relative importance of the channels of our interest in the vicinity of Scenario II. Again, the viable region of parameter space where the relic density is within 2σ of the experimental value determined by the Planck satellite closely follows the border between the neutralino and stop LSP regions. In most scenarios along this border the mixture of the contributing processes is similar to the one presented in Tab. <ref>. However, around M_q̃_L∼ 1600 GeV, where the scalar top mass reaches about half of the heavy Higgs mass (m_H^0 = 3451.2 GeV), the composition of the contributing processes changes. The stop-anti-stop annihilation processes enhanced by the Higgs exchange grow in importance. In contrast to the situation around Scenario I, for large masses of the neutralino dark matter and the scalar top NLSP, the annihilation and co-annihilation processes are not efficient enough to produce the required observed relic density which can be partly compensated by lowering the mass difference. For even larger masses the region where the relic density is compatible with the Planck measurement features a stop LSP, such that neutralino dark matter would be excluded for M_1 ≳ 1800 GeV.  




 §.§ Leading order




Having shown that the processes in Eqs. (<ref>), (<ref>) and (<ref>) are important in large regions around the two scenarios introduced in the previous section, we now turn to review important features of the leading-order cross sections of these processes. The Feynman diagrams for the processes in question are shown in Fig. <ref>. The matrix elements of all three processes considered here have contributions from t-channel or u-channel exchanges of strongly interacting gluinos as well as from electroweak gauginos. Therefore the cross sections can be symbolically written as

    σ = σ_s(α_s^2) + σ_se(α_s α_e) + σ_e(α^2_e) ,

where σ_s is the cross section proportional to the square of the strong coupling constant α_s^2, σ_se is the cross section originating from the interference of the strong and electroweakly interacting parts of the scattering amplitude and σ_e is the purely electroweak cross section proportional to the square of the electromagnetic coupling constant α_e^2.

The decomposition of the total cross section into contributions from different channels and interferences of the three processes under consideration here, is shown in Fig. <ref>. The cross sections for t̃_1 t̃_1 → t t and b̃_1 b̃_1 → b b in the top left, top right and bottom right panels in Fig. <ref> show the expected hierarchy, in which the gluino t-channel and u-channel exchanges dominate the cross section and are about an order of magnitude larger than the next largest contribution which is the interference of the gluino exchange with the electroweak t- and u-channels. The contribution from the interference between the gluino exchange diagrams and the gaugino exchange diagrams is yet another order of magnitude larger than the purely electroweak contribution. As argued before, in scenarios where the processes in Eqs. (<ref>) to (<ref>) are important, the lightest neutralino is bino-like and the gluino mass is relatively small. These facts imply that the neutralino-squark-quark coupling and the gluino-squark-quark coupling differ mainly by the coupling constant. Therefore, the hierarchy observed in Fig. <ref> is simply due to the ratio of the different coupling constants α_s(√(m_t̃_1m_t̃_2))/α_e(m_Z).  




The only process where this hierarchy is not present is the annihilation t̃_1 b̃_1 → t b. Here the hierarchy observed in the other two processes is modified due to a few factors. First, there is no gluino u-channel exchange. Then, the gluino mass in this scenario is larger, and this processes proceeds also through a chargino u-channel exchange. The larger gluino mass together with the missing u-channel suppresses the gluino contribution compared to the other two processes. Moreover, in the case of the higgsino-like chargino exchange, the Yukawa component of the chargino-squark-quark coupling is not suppressed as in the case of bino-like neutralino. The combination of these effects results in the interference between the gluino and the chargino exchange being suppressed with respect to the pure gluino contribution only by a factor of about two. On top of that, the electroweak contribution is comparable with the gluino-chargino interference.




 §.§ Color decomposition




Another important aspect of the processes we investigate is the fact that both initial and final state particles carry color. The color structure of the initial and final state will be extremely relevant later in the discussion of the next-to-leading SUSY-QCD corrections and their resummation. Both scalar quarks in the initial state (and also the quarks in the final state) transform under the fundamental representation of the SU(3) group (denoted here as 3 due to the dimensionality of the representation). The two particle system however transforms under a tensor product of the corresponding representations 3⊗3 which can be decomposed via a Clebsch-Gordan decomposition into SU(3)-invariant subspaces as 

    3⊗3 = 3⊕6 .

In order to construct a color basis adapted to our matrix element, we can use the Clebsch-Gordan coefficients of the decomposition

    C^{3}_α  a_1a_2   =   1/√(2)ϵ_α a_1 a_2 ,    α =1,2,3 ,
    
        C^{6}_α  a_1a_2   =   1/2(δ_α_1 a_1δ_α_2 a_2+δ_α_1 a_2δ_α_2 a_1) ,   α_i = 1,…,6 ,

where the indices a_1,2 can take the values 1 to 3 (for details see Ref. <cit.>). The basis is constructed by considering that SU(3) color symmetry is an exact symmetry of the theory and so the color is conserved between the initial and final states. That means if a pair of initial state particles transforms in an irreducible representation of the SU(3) group, the pair of final state particles must transform in the same representation. After proper normalization, we can combine the Clebsch-Gordan coefficients into the following basis relevant for our processes

    C^{3,3}_a_1a_2a_3a_4=1/√(2N_c(N_c-1))(δ_a_1 a_3δ_a_2 a_4-δ_a_1 a_4δ_a_2 a_3) ,

and

    C^{6,6}_a_1a_2a_3a_4=1/√(2N_c(N_c+1))(δ_a_1 a_3δ_a_2 a_4+δ_a_1 a_4δ_a_2 a_3) .

The matrix element can be expanded in this basis as

    M_stij = M_3  C^{3,3}_stij + M_6  C^{6,6}_stij ,

where s,t,i and j are the color indices of the incoming and the outgoing particles. Given the orthonormality of the basis, the triplet and sextet parts of the amplitude can be determined as

    M_3 = M_stij C^{3,3}_stij, 
    	     
    	M_6 = M_stij C^{6,6}_stij .

In the case of the annihilation process t̃_1 t̃_1 → t t or b̃_1 b̃_1 → b b, the triplet and the sextet matrix elements are a linear combination of the gluino and gaugino t-channel and u-channel exchanges. At tree-level the explicit expression for the triplet part of the matrix element is 

    M_3 = (N_c^2-1)/2√(2N_c(N_c-1))    (-M_g̃^t + M_g̃^u)
        + N_c (N_c-1)/√(2N_c(N_c-1)) (M_χ̃^t - M_χ̃^u)  .

Analogously, the sextet part of the matrix element is

    M_6 = (N_c^2-1)/2√(2N_c(N_c+1))    (M_g̃^t + M_g̃^u)
       + N_c(N_c+1)/√(2N_c(N_c+1))(M_χ̃^t + M_χ̃^u) .

The same decomposition can be performed for the process t̃_1 b̃_1 → t b and the explicit results given in Eqs. (<ref>) and (<ref>) can be used after setting M_g̃^u = 0 and interpreting M_χ̃^u as the u-channel chargino exchange. The squared amplitude is then in all cases given simply by

    |M|^2 = |M_3|^2+|M_6|^2 ,

where due to the orthonormality of the color basis, there is no interference between the triplet and the sextet matrix elements.

The leading order triplet and sextet cross sections for the relevant processes are shown in Fig. <ref>. The general behavior of the color decomposed cross sections for the processes t̃_1 t̃_1 → t t and b̃_1 b̃_1 → b b is very similar. Both processes contain identical particles in the initial state and are symmetric with respect to their interchange. Given that the color basis vector C^{3,3} is anti-symmetric with respect to the same interchange, the partial wave of the the triplet cross section is a p-wave, making its contribution to the relic density subdominant. For these two processes only the sextet color combination contributes. In the case of the last process t̃_1 b̃_1 → t b, the symmetry argument does not apply and both color combinations contain an s-wave and contribute equally to the relic density.






§ NEXT-TO-LEADING ORDER



In this section we will discuss the details of our analytical calculation of the full SUSY-QCD corrections to squark pair annihilation into a pair of quarks. We first concentrate on the virtual corrections and treatment of the UV divergencies in the case of squark pair annihilation. We continue with the discussion of the treatment of IR divergencies. Finally, we address the Sommerfeld enhancement, its treatment and impact on the full SUSY-QCD correction to the squark pair annihilation.

Before we discuss specific details of the next-to-leading order calculation, we will address the systematics of SUSY-QCD corrections to processes which at leading order have both strong and electroweak contributions (see Eq. (<ref>)). If we consider any radiative corrections (SUSY-QCD or electroweak) to the processes in question, the cross section including the next-to-leading order corrections can be symbolically written as

    σ^ NLO = σ^ Tree   + Δσ_s^ NLO(α_s^3) + Δσ_se^ NLO(α_s^2α_e)
       + Δσ_e^ NLO(α_sα_e^2) + Δσ_ee^ NLO(α_e^3) .

The SUSY-QCD corrections contribute to the Δσ_s^ NLO, Δσ_se^ NLO and Δσ_e^ NLO parts of the NLO cross section whereas the electroweak corrections would contribute to the Δσ_se^ NLO, Δσ_e^ NLO and Δσ_ee^ NLO parts. Both classes of corrections, the SUSY-QCD and the electroweak, are ultraviolet and infrared finite and gauge independent by themselves making them formally consistent. 

The first and leading term in the NLO correction is Δσ_s^ NLO(α_s^3) which receives contributions only from SUSY-QCD corrections. In particular these are the SUSY-QCD corrections to the gluino exchange diagrams interfered with the gluino tree-level contribution. These corrections are the main result of this analysis.

The following term Δσ_se^ NLO(α_s^2α_e) receives contributions from three sources - from the interference of the SUSY-QCD corrected gluino exchange with the electroweak gaugino exchange, from the interference of the SUSY-QCD corrected electroweak gaugino exchange with the gluino diagrams and the last contribution would come from electroweak corrections to the gluino exchange interfered with the gluino tree-level. The last contribution is not included in this analysis and even though it is formally of the same order, due to the small size of the electroweak corrections which are typically a factor 10 smaller than SUSY-QCD ones this last contribution is the smallest of the three. This way our analysis provides also the leading corrections in the term Δσ_se^ NLO(α_s^2α_e).  

The third term Δσ_e^ NLO(α_sα_e^2) contains the interference of the SUSY-QCD corrected electroweak gaugino exchange with the leading order electroweak gaugino diagrams as well as electroweak corrections to both parts of the interference between the gluino and the electroweak gaugino exchange.

The last term Δσ_ee^ NLO(α_e^3) is not considered here as it contains only electroweak corrections to the electroweak parts of the cross section.

The analysis presented here does not consider electroweak corrections as they are for the most part subleading and contribute about 1% to 3% correction <cit.>. In some instances however, the electroweak corrections and specifically the Yukawa corrections can become important <cit.>. Even though we do not calculate electroweak corrections in this analysis, the leading effects of the enhanced Yukawa corrections are taken into account as described in <cit.>. In particular, these become relevant in the case of chargino exchange in the Scenario II (neutralino exchanges are not as enhanced due to the lightest neutralino being a pure bino in both scenarios).

As the discussion below shows, the SUSY-QCD corrections presented here are the dominant corrections even in scenarios with large tanβ and are even more dominant owing to the presence of the Sommerfeld enhancement.






 §.§ Virtual corrections and renormalization



The class of processes considered here – the squark pair annihilation to a pair of quarks – include strongly interacting particles in the initial state, in the final state, and even in the intermediate state (in the case of the gluino t- and u-channel exchanges). As a consequence, the next-to-leading order SUSY-QCD corrections include contributions from vertex corrections, propagator corrections, and box corrections. The corresponding diagrams are displayed in Figs. <ref>, <ref> and <ref>, respectively. The next-to-leading order corrections to the squark-pair annihilation contain one-loop diagrams which are ultraviolet (UV) and infrared (IR) divergent. The UV divergencies are cancelled by renormalization of the parameters of the theory and the fields. In order to cancel the IR divergences, one has to properly define an infrared safe cross section, which is done by including also 2→ 3 processes with an additional gluon being radiated (see Fig. <ref>). 

All one-loop diagrams have been calculated in the SUSY-invariant dimensional reduction scheme () <cit.> where, similar to the minimal subtraction scheme (), the number of space-time dimensions is set to D=4-2ε in order to regularize otherwise divergent loop integrals. We have used the standard Passarino-Veltman reduction <cit.> in order to reduce the tensor loop integrals in the one-loop amplitudes to only a few scalar integrals which have then been evaluated using well-known results, e.g., Refs. <cit.>. Our analytical calculations were performed and verified with the help of publicly available tools  <cit.>,  <cit.>, and  <cit.>.

In order to cancel UV divergencies of the one-loop amplitude, we renormalize the MSSM by introducing counterterms to the relevant parameters and fields. Because we consider SUSY-QCD corrections to a process involving scalar quarks, quarks and intermediate gluinos, the relevant parameters are the ones that receive corrections proportional to the strong coupling constant α_s. Every renormalization scheme is characterized by a careful selection and definition of its input parameters. In a series of previous analyses <cit.>, we have put forward a renormalization scheme which combines the advantages of both on-shell and  renormalization schemes and treats consistently the renormalization of the quark and the squark sector. In these sectors the input parameters are chosen to be the on-shell masses m_t, m_t̃_1, m_b̃_1, m_b̃_2, the mass of the bottom quark m_b, and trilinear couplings of the third generation A_t and A_b. The last three parameters are defined in the  renormalization scheme. We define our renormalization scale as μ_R=Q_SUSY=√(m_t̃_1 m_t̃_2). In the following, we will comment on new aspects of the renormalization relevant for our analysis presented here such as the renormalization of the mass of the gluino and its wave function. 




  §.§.§ Gluino mass and wave-function renormalization



In all our analyses we adopt a convention where not only the complete next-to-leading order corrections to the cross section should be rendered UV finite, but also all building blocks such as the n-particle irreducible Green's functions should be UV finite as well. This choice requires to introduce wave-function renormalization constants not only to the fields that correspond to the initial and final state particles but also to fields that give rise to internal propagators. In our case, the only strongly interacting particle that appears in a propagator in our amplitude is the gluino which has not yet been treated within the  analysis.

In the case of the gluino, both wave function and mass have to renormalized in order for the vertex corrections and propagator corrections to be separately UV finite. To this end, we introduce counterterms to the gluino wave function δ Z^L,R_g̃ and the gluino mass δ m_g̃ as

    ψ_g̃   →   (1+ 12 δ Z^L_g̃P_L+ 12 δ Z^R_g̃P_R) ψ_g̃ ,
    
    	m_g̃   →    m_g̃ + δ m_g̃ .

All gluino counterterms are determined by considering the gluino two-point Green's function. The one-loop contribution to the two-point Green's function is given by the gluino self-energy diagrams shown in Fig. <ref>. This contribution can be parametrized as

    Π(k) = k (P_L Π^L(k^2)    + P_R Π^R(k^2))
       +  m_g̃(P_L Π^SL(k^2) +  P_R Π^SR(k^2)) ,

where Π^L,R(k^2) and Π^SL,SR(k^2) are form factors which receive contributions from the corresponding self-energy diagrams. 

Even though the gluino is not an external particle in the processes considered in this analysis, we still require that the residue of the propagator at one-loop order is set to unity. This condition fixes the wave-function renormalization constant using the form-factors as

    δ Z^L/R_g̃   = -Π^L/R(m_g̃^2) + 12 (Π^SL/SR(m_g̃^2)
       -Π^SR/SL(m_g̃^2))
    	- m_g̃^2  (Π̇^L/R(m_g̃^2)+Π̇^R/L(m_g̃^2)
       + 
    	Π̇^SL/SR(m_g̃^2)+Π̇^SR/SL(m_g̃^2)) ,

where Π̇^i(m_g̃^2) = ∂/∂ k^2.Π^i(k^2)|_k^2=m_g̃^2. Using the gluino wave-function counterterm renders both the propagator and vertex corrections separately UV finite. Moreover, given that the gluino is not an external particle, renormalization of its wave-function is not necessary for UV finiteness of the full next-to-leading order amplitude and so the full amplitude is independent of the gluino wave-function counterterm. This constitutes another consistency check of our analytical calculation.

The mass counterterm is determined from the on-shell condition which requires that the gluino mass m_g̃, which is an input parameter, is identical with the position of the pole of the gluino propagator. It is given as

    δ m_g̃ = 12 m_g̃   (Π^L(m_g̃^2)+Π^R(m_g̃^2). 
       .+Π^SL(m_g̃^2)+Π^SR(m_g̃^2)) .





  §.§.§ Some remarks on the renormalization of the squark sector



Even though we have discussed the details of the renormalization of squark parameters in Ref. <cit.>, we would like to remark here on one feature of the renormalization scheme relevant for evaluating the results of this analysis. As discussed in Ref. <cit.>, we use the relation between the non-diagonal squark mass matrices for up-type and down-type squarks

    U^q̃( [ m^2_LL m^2_LR; m^2_RL m^2_RR ]) 
    	(U^q̃)^† = ( [ m^2_q̃_1        0;        0 m^2_q̃_2 ])	 ,

where

    m^2_LL   =    M_Q̃^2 + (I^3L_q - e_q s_W^2)cos2β m_Z^ 2 + m_q^2 ,
    
    	m^2_RR   =    M_{Ũ, D̃}^2 + e_q s_W^2 cos2β m_Z^ 2 + m_q^2 ,
    
    	m^2_LR   =    m^2_RL = m_q( A_q - μ (tanβ)^-2 I^3L_q) ,

to relate the input parameters in the whole squark sector, which are defined in different renormalization schemes. In the next step, we determine the dependence of the soft supersymmetry-breaking squark mass parameters M_Q̃^2 and M_{Ũ, D̃}^2 of the three on-shell masses m_b̃_1, m_b̃_2, m_t̃_1 and the input parameters contained in m^2_LR (A_q, m_q). In the sbottom sector the parameters M_Q̃^2 and M_D̃^2 are contained in the matrix elements m^2_b̃,LL and m^2_b̃,RR, which are given by the on-shell masses as

    m^2_b̃,LL   =   12 (m^2_b̃_1+m^2_b̃_2) ±12 √((m^2_b̃_1-m^2_b̃_2)^2 - 4 m^4_b̃,LR) ,
    
        m^2_b̃,RR   =   12 (m^2_b̃_1+m^2_b̃_2) ∓12 √((m^2_b̃_1-m^2_b̃_2)^2 - 4 m^4_b̃,LR) .

One notices that there are two possible values for the parameters M_Q̃^2 and M_D̃^2 and consequently also for the third parameter M_Ũ^2 which can be found in one of the diagonal elements of the non-diagonal scalar top mass matrix and is related to the first two parameters
through

    m^2_t̃,RR   = 1/m^2_t̃,LL-m^2_t̃_1
       ×(m^2_t̃_1 m^2_t̃,LL-m^4_t̃_1 + m^4_t̃,LR) .

The parameter M_Q̃^2 is common to both elements m^2_b̃,LL and m^2_t̃,LL. Given the freedom to choose from two possible solutions for the squark soft supersymmetry-breaking mass parameter, we can end up with two possibly very different mass matrices and two different sets of mixing matrices (three out of four masses of the squarks would be the same in both cases as they are used as input). In order to ensure a naturally small correction to the mixing matrices when changing between our and the  renormalization scheme, we always select the solution which preserves the hierarchy between the mass matrix elements in the scalar top quark sector m^2_t̃,LL and m^2_t̃,RR which was present in the pure  scheme.




 §.§ Real corrections



After treating the ultraviolet divergencies and removing them by renormalization, we are left with the infrared (IR) divergencies in the one-loop amplitude. The IR divergencies are also regularized by performing calculations in a general dimension D=4-2ε and are subsequently removed by considering an IR-safe observable, which in our case is a cross section with one additional gluon in the final state. The corresponding diagrams are depicted in Fig. <ref>. 

The cross section for the radiation of an additional gluon cannot be calculated analytically. On the other hand, the cancellation of the IR divergencies between the virtual corrections and the real radiation cross section has to be performed analytically. One of the methods how to extract the IR divergencies out of the real radiation cross section is the phase-space slicing method <cit.>. This method is based on the fact that the infrared divergence is connected to a specific configuration of the momentum of the gluon. The soft infrared divergence arises when the additional gluon's energy vanishes whereas the collinear infrared divergence arises when the additional gluon is radiated collinearly to the momentum of an external massless particle. The phase-space slicing method uses kinematical cuts to divide the three-particle phase-space into regions where either one or both of the aforementioned configurations of the gluon 4-momentum occur and the remainder where the cross section is IR finite. In the singular regions the full matrix element is replaced by an approximation which can be integrated analytically making the IR divergence explicit. In the non-singular region the full matrix element can be integrated numerically without any obstacles.

In our case, all external particles are massive so the 2→ 3 cross section contains only a soft IR divergence. In the singular region where the gluon's energy is smaller than an arbitrary small cutoff Δ E, we use the soft gluon approximation which factorizes the differential cross section as

    ( σ/Ω)_soft = F ×( σ/Ω)_Tree,

where F contains the integral over the phase-space of the gluon

    I_ab =  μ^4-D∫_|k⃗| ≤Δ E^D-1 k/(2π)^D-41/k^0(a.b)/(k.a)(k.b),

and therefore also the divergence. Here, k is the 4-momentum of the gluon and a and b are 4-momenta of two external particles which can emit a gluon. These integrals are given in Refs. <cit.>.

The cutoff Δ E introduced to separate the singular from the non-singular region in the three particle phase-space enters into the calculation of the soft-gluon radiation as well as into the integration of the 2→ 3 cross section over the non-singular phase-space. In principle the dependence on this cutoff should disappear in the sum of the contributions from both phase-space regions, but in practice the independence on the cutoff is limited by the numerical stability of the integration over the non-singular region. We have investigated the dependence on the cutoff and found that the integration in our case is stable and independent of the cutoff for a relatively large interval of cutoffs around Δ E = 10^-5√(s).

The complete result after we have included all virtual corrections, the counterterms and the real radiation 2→ 3 cross section is UV and IR finite. In contrast to the leading order result which consists of a cross section with two particles in the final state and is implemented in , the complete result is a consistent combination of a one-loop corrected cross section with two particles in the final state and a leading order cross section with three particles in the final state. In , the complete result replaces the leading order result of .




 §.§ Sommerfeld resummation




When calculating the relic density in our case, an important contribution comes from the annihilation of squarks moving with non-relativistic velocities. If annihilating, non-relativistic particles couple to much lighter force mediators which in our case are the gluons, the annihilation cross section is modified due to the well-known Sommerfeld effect <cit.>. The reason for this modification is that the exchange of n gluons between the initial state squarks (see Fig. <ref>) contains a correction proportional to (α_s/ v_rel)^n. This correction becomes significant and can spoil the perturbative expansion when the relative velocity of the squark pair v_rel is comparable to the strong coupling constant α_s. In such a case these contributions have to be resummed to all orders leading to the Sommerfeld effect. 

Small relative velocities occur naturally in the freeze-out regime, E_kin∼ T_FO∼ m_χ̃^0_1/25 and therefore the Sommerfeld resummation is expected to be relevant in the case of dark matter annihilation in general and in our case in particular. 



As for our processes of interest, q̃_i q̃'_j → qq', the cross section is dominated by the s-wave component (see Fig. <ref>). We can factorize the resummed cross section as

    (σ v)_resum = S_0,[3] (σ v)^Tree_q̃q̃→ qq, [3] + S_0,[6] (σ v)^Tree_q̃q̃→ qq, [6] ,

where we have split the leading-order cross section according to its color contribution to the triplet and sextet configurations (see Sec. <ref>). S_0,{[3],[6]} indicate the corresponding s-wave Sommerfeld factors, whose evaluation we discuss in the following. In the non-relativistic limit, the resummation of the gluon exchange diagrams as shown in  Fig. <ref> amounts to solving the Schrödinger equation with the corresponding Coulomb potential. The Coulomb potential including gluon loops at next-to-next-to-leading order was evaluated in <cit.> and extended by fermion loops in Ref. <cit.>. In the  scheme, the Coulomb potential reads <cit.>


    V^[R](r)    = C^[R] α_s(μ^[R]_C)/r
       ×{ 1 + α_s(μ^[R]_C)/4π[2 b_0 (ln(μ^[R]_C r) + γ_E )+a_1 ]  } ,


with γ_E=0.5772 being the Euler-Mascheroni constant. Furthermore, we have defined 


    b_0    = 11/3C_A - 4/3T_f n_f , 
    
    	a_1    = 31/9C_A - 20/9T_f n_f ,


where b_0 corresponds to the one-loop β-function coefficient with C_A=3 and T_F=1/2. We treat the top as the only massive quark, such that we set the number of massless quarks to five (n_f=5). The Coulomb potential, given in Eq. (<ref>), describes the interaction of any non-relativistic colored particles transforming in general SU(3)-representations R_1 and R_2. The color structure of such a scattering process can be decomposed as

    R_1⊗R_2 = R^'⊕R^'' .

The color factor C^[R] is given in terms of the quadratic Casimir operators of the relevant SU(3)-representations as

    C^[R^i] = T_1^aT_2^a = 1/2 [(T_1^a + T_2^a)^2 - (T_1^a)^2 - (T_2^a)^2]
       = 1/2 (C_2^R^i - C_2^ R_1 - C_2^ R_2 )      R^i = R^',R^'' .

In the case considered here, the two squarks in the initial state both transform under the fundamental representation of SU(3) and the color decomposition is 3⊗ 3 = 3⊕ 6.
Using the quadratic Casimir operators for the fundamental and the sextet representations, we obtain <cit.>

    C^3   = -1/2(1 + 1/N_c) =  -2/3 ,         
    
    C^ 6   =  1/2(1 - 1/N_c)= 1/3 .


The Sommerfeld factors are then obtained by solving the Schrödinger equation


    [ -2/m_red∇^2 + V^[ R](𝐫) - ( √(s)+iΓ_t̃_1) ]
    		   𝒢^[ R]( 𝐫;√(s)+iΓ_t̃_1)
       = δ^(3)(𝐫)


with the reduced mass m_red=(m_q̃  m_q̃^')/(m_q̃ + m_q̃^') of the two annihilating particles q̃ and q̃^'. The solution of the Schrödinger equation with the NLO Coulomb potential defined in Eq. (<ref>) is given by the Green's function 𝒢^[R](𝐫,E+iΓ_t̃_1) = 𝒢^[R](𝐫, 𝐫^'=0,E+iΓ_t̃_1). 
The Sommerfeld factor which is used to correct the cross section in Eq. (<ref>) is given by a ratio of two Green's functions at the origin (r = 0) <cit.> 


    S_0,[R] = [𝒢^[R](0,E+iΓ_t̃_1)]/[𝒢_0(0,E+iΓ_t̃_1)] ,

where the Green's function 𝒢_0(0,E+iΓ_t̃_1) stands for the solution of the Schrödinger equation without any Coulomb potential.

The solution to Eq. (<ref>) at the origin is well known <cit.>, and we consider here all terms up to NLO,


    G^[R](0;√(s)+iΓ_t̃_1)   = i m^2_red v_ s/π + C^[R]α_s(μ^[R]_C) m^2_red/π
       ×[ g_LO + α_s(μ^[R]_C)/4πg_NLO] ,


where the LO and NLO contributions are given by


    g_LO = L - ψ^(0), 
       g_NLO= 	β_0 [ L^2 - 2L(ψ^(0) - κψ^(1)) + κψ^(2) + 
    		(ψ^(0))^2  
       -3ψ^(1) - 2κψ^(0)ψ^(1) 
    		+ 4 _4F_3(1,1,1,1;2,2,1-κ;1)] 
        + a_1 [L-ψ^(0)+κψ^(1)] .


In Eqs. (<ref>)-(<ref>) and in the following, we use 
the short-hand notation


    κ   =i C^[R]α_s(μ^[R]_C)/2v, 
    
    	v_ s    =  √(√(s)+iΓ_t̃_1-2m_avg/2 m_red),
    
    	L     =  lni μ^[R]_C/4 m_redv_ s .


Moreover, ψ^(n)=ψ^(n)(1-κ) is the n-th derivative of ψ(z) = γ_E + d/dzlnΓ(z) with the argument (1 - κ), _4F_3(1,1,1,1;2,2,1-κ;1) a hypergeometric function, and m_avg=(m_q̃+ m_q̃^')/2 the average mass of the two incoming particles. Note that in the case of identical initial state particles the parameter v_ s in Eq. (<ref>) is the non-relativistic velocity of one of the incoming particles, and should not be confused with v_rel=2v_ s, the relativistic, relative velocity of the two annihilating particles. In order to calculate the Sommerfeld factor in Eq. (<ref>), we also need the Green's function that solves the system without any potential term in Eq. (<ref>), which is given by


    [𝒢_0(0,E+iΓ_t̃_1)] = m^2_red v_ s/π .


Finally, we need to fix the scale μ_C that appears in the potential and has an impact on the evaluation of α_s in the Sommerfeld factor. We follow here the treatment presented in Ref. <cit.> and set 

    μ^[R]_C = max{ 4 m_red v_ s, μ^[R]_B } ,

where 4 m_red v_ s is motivated by the typical momentum exchange of the gluons in the ladder diagram and the scale μ_B^[R] corresponds to twice the inverse Bohr radius r_B. It is defined via 

    μ_B^[R]≡ 2/r_B = 2 C^[R] m_redα_s(μ_B^[R]) .

In order to obtain μ_B^[R], we solve Eq. (<ref>) iteratively. 

As the box diagrams in the full NLO calculation also contain the velocity-enhanced part of the one-gluon exchange, which is at the same time already included in the Sommerfeld resummation, we have to subtract this contribution in order to avoid any double counting.  

To isolate the velocity-enhanced term from the box contribution, we expand the box contribution in the relative velocity (for details see App. <ref>). We then construct the subtracted cross section (σ v)^sub_NLO based on the expanded matrix element of the box diagrams given in Eq. (<ref>). 
The leading velocity-enhanced term of the subtracted cross section 
is

    (σ v)^sub_NLO∼∑_ R(α_s(μ_R) π/v_rel)   C^[𝐑]_box  (σ v)^ R_Tree .

Comparing Eq. (<ref>) with the next-to-leading order part of the Sommerfeld resummation which arises from the imaginary part of g_LO, namely g^sub_LO = i π/2, and reads 

    (σ v)^NLO_resum=  ∑_ R(C^[R]α_s(μ^[R]_C) π/2 v_ s)   (σ v)^ R_Tree ,

we see, that given C^[R] =  C^[𝐑]_box and v_ s=v_rel/2, the two expressions differ in the scale at which the strong coupling constant is being evaluated. While in the perturbative NLO calculation α_s is evaluated at the renormalization scale μ_R = √(m_q̃m_q̃^'), in Sommerfeld resummation the characteristic scale μ_C is used. By choosing to use Eq. (<ref>) to avoid the double counting, we make use of the fact that the natural scale used in the description of the interaction between incoming particles at small velocities is μ_C. This is consistently used to all orders in the resummed cross section (σ v)_resum given by Eq. (<ref>).








The full next-to-leading order cross section including consistently also the Sommerfeld resummation reads


    (σ v)^full= (σ v)_NLO + (σ v)_resum - (σ v)^sub_NLO .


Given the large trilinear couplings in both scenarios, it might be interesting to study the Sommerfeld enhancement coming from the exchange of Higgs bosons <cit.>. In the regime of Sommerfeld enhancement, also bound state formation can potentially occur, giving rise to new annhilation channels and thus altering the relic density prediction. This has been previously studied for stop-antistop annihilation for both gluon <cit.> and Higgs exchange <cit.>. Such studies, however, are far beyond the scope of this work and are left for future analyses. 












 §.§ The NLO cross section results




In this section we present the first result of our analysis, which is the impact of SUSY-QCD next-to-leading order corrections on the annihilation cross sections of scalar top or bottom pairs. Apart from the cross section (or more precisely σ v) we also show in arbitrary units the Boltzmann distribution function which is involved in the calculation of the thermal average ⟨σ v ⟩ at the freeze-out temperature (grey shaded area). It should serve as a reminder that the cross section contributes to the determination of the relic density only in a limited range in the center-of-mass momentum p_ cm. 




  §.§.§ Scenario I



In the first scenario introduced in Sec. <ref>, the mass splitting between the lightest neutralino and the lightest stop quark is relatively
large. As a result the dark matter annihilation cross section receives important contributions not only from the stop pair
annihilation into top quarks but also from the neutralino-stop co-annihilation into a top quark and a gluon and other final states. 
The results of SUSY-QCD corrections for these processes are shown in Fig. <ref>. As described in the previous section, the 
next-to-leading order cross section consists of vertex corrections, propagator corrections, box corrections, counterterm contributions and the real radiation
cross section, which has to be added to render the prediction infrared finite. The corrections from each contribution are not shown in 
Fig. <ref> due to the cancellations of the ultraviolet and infrared divergencies between the contributions, which make each contribution on its own ill defined. In addition to the NLO cross section, we include also the enhanced higher-order contributions stemming from the non-relativistic Coulomb correction.



In case of the annihilation of a pair of scalar top quarks, both initial particles are colored and in the limit of vanishing relative velocity of the squark pair, the Coulomb corrections dominate the full corrected cross section. The origin of these corrections is the exchange of multiple gluons between a pair of slowly moving squarks in the initial state and the details were discussed in Sec. <ref>. The effect of the Coulomb corrections strongly depend on the color multiplet, in which the pair of squarks transform. Based on the color decomposition presented in Sec. <ref>, the annihilation cross section t̃_1  t̃_1  →  t   t in Scenario I is dominated by the contribution where the squark pair forms a SU(3)-sextet (see Fig. <ref>). In this representation the multiple exchange of the gluons can be described by a repulsive non-relativistic QCD potential, as discussed in detail in Sec. <ref>. That is why the Coulomb corrections in this case cause a reduction of the cross section. As already discussed in the previous section, the next-to-leading order cross section contains the one-loop contribution included also in the Coulomb enhancement. This one-loop contribution can be traced back to all box diagrams in Fig. <ref>, where one gluon is exchanged between the incoming squarks. The contribution from this class of diagrams dominates the one-loop cross section for small velocities and is so large that it causes the cross section with one-loop corrections to be negative (see the green dashed line in Fig. <ref>). As discussed in Sec. <ref>, in order to prevent double counting, we remove the part of the box contribution which is already included in the Coulomb resummation. This allows us to quantify the pure one-loop correction to the annihilation cross section without any enhancement (red dash-dotted line in Fig. <ref>). We see that the one-loop correction without the enhancement is a large positive correction of about 30-40% over a large range of p_ cm.

Comparing the result for the non-enhanced NLO cross section with the full result, which is the sum of the non-enhanced NLO cross section and the Sommerfeld corrections shows that the latter are important for all relevant values of p_ cm. Starting at the largest value of p_ cm∼ 600  GeV which is still relevant for the determination of the relic density, we observe that the Coulomb corrections already reduce the constant 30% NLO correction by a few percent. For smaller relative velocities corresponding to p_ cm∼ 150  GeV, the NLO correction is fully cancelled by the Coulomb corrections, and for very slow velocities the Coulomb corrections take over and the overall correction is large and negative. For almost vanishing velocities the total cross section becomes negative due to this large negative correction. Even though a negative cross section is unphysical, the fact that the cross section vanishes should not be very surprising. The dynamics of the squark pair in the regime when Coulomb corrections are very large (meaning for vanishing velocities) correspond to a motion of the pair in a highly repulsive QCD potential. This in turn means that large repulsive forces repel one squark from the other reducing the probability of annihilation and thereby reducing the cross section. The cross section becomes negative only for p_ cm < 10  GeV, which is irrelevant for the relic density determination, as for such small momenta it is multiplied by an almost vanishing Boltzmann distribution function.

In summary, we can conclude that SUSY-QCD corrections to t̃_1  t̃_1  →  t   t are sizeable either through the one-loop corrections for large p_ cm or the enhanced Coulomb corrections for small p_ cm.

The co-annihilation processes important in this scenario were discussed in detail in Refs. <cit.>. In Fig. <ref> we also show the effect of the SUSY-QCD corrections on the co-annihilation cross sections in Scenario I. We see that the next-to-leading order corrections in the case of co-annihilations are substantial ranging from -30% in the case of the co-annihilation into top quark and Higgs-boson final state to +50% in the case of the top gluon final state. There are a few substantial differences such as the fact that the corrections are negative in the case of co-annihilations with electroweak bosons or Higgs bosons in the final state or that there is a large difference between our leading order result and the  result for the co-annihilations with electroweak and Higgs bosons, which can be traced to a different definition of underlying parameters in our renormalization scheme. The next-to-leading order correction with respect to the  result is largely reduced to at most -10%. Given that our leading order prediction for the co-annihilations into top quark and a gluon coincides with  prediction, the large next-to-leading order correction gives directly also the correction with respect to the  result. 




  §.§.§ Scenario II



In the second scenario, the choice of parameters such as tanβ and the gaugino and squark mass parameters cause the masses of the lightest neutralino, the lightest scalar top and bottom quarks to be almost degenerate. This leads to different processes contributing significantly to the total dark matter annihilation cross section. The smaller mass difference renders co-annihilations ineffective and the fact that three particles are mass degenerate leads to a larger number of annihilations. Moreover, the large value of tanβ enhances the gaugino exchange in the case of the stop-sbottom annihilation, and this together with a different color structure makes this annihilation dominant in the case of Scenario II. 

The full next-to-leading order results for three dominant processes are shown in Fig. <ref>. The processes t̃_1  t̃_1  →  t   t and b̃_1  b̃_1  →  b   b have very similar features to the annihilation of a pair of scalar top quarks in Scenario I. The main difference in this scenario is the process t̃_1  b̃_1  →  t   b, which has an entirely different decomposition of the leading order cross section in terms of the t- and u-channel exchanges combined with a different color decomposition, which is essential in explaining the behaviour of the NLO cross section. Similar to the already discussed case of t̃_1  t̃_1  →  t   t, the NLO correction contains a velocity enhanced term, which is already resummed in the Sommerfeld correction. In order to avoid double counting, we define again the non-enhanced NLO correction σ_NLO_B where we subtract the term which is already accounted for by the Sommerfeld resummation (red dash-dotted curve in Fig. <ref>). As one can see in Fig. <ref>, the non-enhanced NLO correction is substantial in all processes in Scenario II. In the case of stop pair or sbottom pair annihilations, this NLO correction is compensated by a large and negative Sommerfeld correction which is here derived from a repulsive QCD potential. The color decomposition of the t̃_1  b̃_1  →  t   b shows (see Fig. <ref>) that in contrast to the other processes, the cross section is here dominated by the part where the initial stop and sbottom quarks transform as a SU(3)-triplet. In this color configuration a pair of slowly moving squarks experiences an attractive strong force, which leads to a large enhancement of the annihilation cross section. Comparing the full result (solid blue line in Fig. <ref>) with the result containing just the non-enhanced NLO corrections, we see that the Sommerfeld enhancement is important over the whole region in p_cm that is relevant for the calculation of the relic density. 

It is worth mentioning that due to the large value of tanβ in Scenario II and to the fact that the chargino u-channel exchange gives an important contribution to the cross section t̃_1  b̃_1  →  t   b, the Yukawa corrections to the chargino exchange can give a non-negligible contribution. We have included the tanβ dependent Yukawa corrections even beyond next-to-leading order in the full result and show their effect separately in Fig. <ref> (yellow dashed line). Even though the Yukawa corrections are non-negligible, they are small (about 3%) compared to the remaining SUSY-QCD corrections or to the Sommerfeld enhancement. 

The full correction to t̃_1  b̃_1  →  t   b is larger than 50% over the whole relevant range of p_cm and can even exceed 100% (without threatening perturbativity as this correction originates from a resummation). We will show in the next section that using the fully corrected annihilation cross sections in Scenario II has a large impact on the relic density.  









§ IMPACT ON THE NEUTRALINO RELIC DENSITY



We finally come to the discussion of the impact of the corrections presented in Sec. <ref> on the neutralino relic density Ω_χ̃^0_1h^2. To this end, we have implemented the corrections into the numerical code DM@NLO, which is used as an extension to micrOMEGAs. In practice, this means that the Boltzmann equation is still integrated using the latter, while the cross section calculation of the relevant processes (see Tab. <ref>) is performed by DM@NLO instead of CalcHEP. 

In the following, we will illustrate the impact of the corrections by comparing the relic density obtained using the full DM@NLO NLO calculation to the values obtained using the tree-level calculation of micrOMEGAs / CalcHEP.




 §.§ Scenario I



We start by examining the impact of NLO corrections in the vicinity of Scenario I, where we compare the relic density obtained from the micrOMEGAs calculation to the one obtained using our full NLO result as presented in Sec. <ref>. The impact is illustrated in Fig. <ref>, where we show the corresponding viable regions of parameter space in the M_1-M_t̃_R plane. As can be seen, the favoured parameter region where the calculated relic density satisfies the experimental constraint, Eq. (<ref>), is shifted towards smaller mass parameters in order to compensate the increased annihilation cross section. It is important to note that this shift is larger than the width of the band which corresponds to the Planck 2σ-uncertainties. 

The situation changes for higher masses, where the processes discussed in this work and corrected at the NLO level, are less relevant. The correction of the remaining processes relevant in this part of parameter space are left for future work.

In Fig. <ref>, we show the same results in the vicinity of Scenario I, but this time projected onto the plane of the physical neutralino and stop masses. Note that, here, the variation of the physical masses solely stems from varying the parameters M_1 and M_t̃_R, respectively, while all other soft parameters, including those that in general may influence the neutralino and stop masses, are kept fixed to the values of Table <ref>. From Fig. <ref>, we see that the cosmologically favoured region of parameter space is shifted by about 7 GeV in both the neutralino and the lighter stop mass. 

These results lead to the conclusion that the corrections presented in this work are relevant when performing an extraction of either physical masses or fundamental model parameters from cosmological data. Let us note that this conclusion is the same as for the previous analyses of other processes entering the calculation of Ω_χ̃_1^0h^2 <cit.>.

In order to get a better understanding of the impact of the different contributions to the annihilation cross section σ_ ann, and consequently to the relic density Ω_χ̃^0_1 h^2, we have performed a one-dimensional scan along the region where  predicts the correct relic density varying simultaneously the parameters M_1 and M_t̃_R. The result of this scan is shown in Fig. <ref> as a function of both parameters while all other parameters were fixed to the values given in Tab. <ref>. 

First, it can be noticed that our tree-level prediction differs from the micrOMEGAs result. This is a direct consequence of the corresponding difference in the annihilation cross sections, as discussed in Sec. <ref> and shown, e.g., in Fig. <ref>. Taking into account the corrections discussed in Sec. <ref>, it can be seen that the total correction is split into two parts associated with the relevant classes of processes, namely t̃_1 t̃_1 → tt and χ̃^0_1 t̃_1 → qg,qV,qϕ. Fig. <ref> shows that the correction to the relic density in Scenario I is dominated by the corrections to the co-annihilation processes even though at leading order these processes contribute about a factor two less than t̃_1 t̃_1 → tt. This is a consequence of the Sommerfeld supression of the annihilation cross section as discussed in Sec. <ref>. Moreover, we see that for lower bino/squark mass parameters M_1/M_t̃_R the correction to the co-annihilation processes is numerically more important than for large mass parameters. This is explained by the fact that the relative importance of the co-annihilation processes is higher in this region of parameter space (see, e.g., Fig. <ref>). Moving towards higher values of M_1, the relative importance of the stop-pair annihilation increases and, consequently, the associated correction to the relic density becomes more important. 




Overall, the relic density obtained using our full (i.e. NLO including resummation) caluclation is about 6% smaller than the one obtained by using micrOMEGAs. Again, we emphasize that this shift is more important than the uncertainty given by the Planck measurement, which is, at the 2σ confidence level of about 2%.




 §.§ Scenario II



Let us now focus on Scenario II, where not only t̃_1 t̃_1 → tt is relevant, but also the related processes b̃_1 b̃_1 → bb and t̃_1 b̃_1 → tb give sizeable contributions to the annihilation cross section σ_ ann. Therefore they need as well to be corrected at the NLO level including the resummation, as discussed in Sec. <ref>.













Again, we start by depicting the parameter region compatible with the measured value for the relic density (see the second row of plots in Fig. <ref>) in the vicinity of Scenario II, in this case in the M_1-M_q̃_L plane, which are the relevant neutralino and squark mass parameters. Similar to Scenario I, as discussed above, the viable regions with respect to the relic density is shifted towards lower masses. Here, the shift is again more important than the uncertainty and is much larger than in Scenario I. It corresponds to a shift of about 17 GeV in the bino mass parameter M_1 and about 15 GeV in the left-handed squark mass parameter M_q̃_L. 

In terms of physical masses, shown in Fig. <ref>, this corresponds to a shift of about 17 GeV for the neutralino mass, and of about 15 GeV in the lighter stop mass. Once more, these findings underline the importance of the presented corrections in the light of precision cosmology. It is to be noted that for this part of the analysis only M_1 and M_q̃_L have been varied and the results have been projected on the so obtained plane of the physical neutralino and stop masses, while all other parameters have remained fixed to the values given in Tab. <ref>.

In order to decompose our full NLO prediction for the relic density into the contributions from different processes, we show in Fig. <ref> the NLO corrected contributions to the relic abundance from individual processes along the region where  predicts the correct relic density varying simultaneously the parameters M_1 and M_q̃_L. We see that already using our tree-level annihilation cross section with differently defined input parameters shifts the relic density by a few percent (black dashed line). In all remaining contributions we use our tree-level annihilation cross section for all relevant processes. Starting from our tree-level, a very small correction of about 1% in the whole region comes from including NLO corrections only to the t̃_1 t̃_1 → tt and b̃_1 b̃_1 → bb. It is important to point out that the reason for this extremely small correction is the Sommerfeld enhancement, which in this case in fact suppresses the cross section. This is due to the repulsive nature of the dominant SU(3)-sextet contribution. Even if the full correction to the cross section was large and negative for small p_ cm (see Fig. <ref>), interestingly the correction to the thermal averaged cross section is still positive leading to a drop in relic density. The largest contribution comes from the annihilation cross section of t̃_1 b̃_1 → tb. The first reason is the large contribution of this process to the total annihilation cross section already at tree-level (see Tab. <ref>). The second reason is the large Sommerfeld enhancement emerging from the attractive potential of the dominant SU(3)-triplet contribution,  which in this case makes the full correction to the cross section extremely large. We see that depending on the dominant contribution of the color decomposition, SU(3)-sextet or SU(3)-triplet, the Sommerfeld corrections either suppresses or enhances the cross section such that the total SUSY-QCD correction to the relic density over the whole range is about 25%. This results in the visible shift of the preferred parameter region, which is much larger than the experimental uncertainty given in Eq. (<ref>).  







§ CONCLUSION



Scenarios in the MSSM with light stops are still very appealing due to their potential to address many problems that the MSSM with heavy particles might have. In such scenarios, squark pair annihilations into quarks are often very important processes that govern the annihilation of dark matter, which is typically the lightest neutralino. 

We have analyzed two such example scenarios, which pass all current experimental constraints. We focused on the SUSY-QCD corrections to squark pair annihilations into quarks. We have reviewed the details of the one-loop calculation and of the Sommerfeld enhancement. We have shown that the one-loop corrections of the cross sections are sizeable even without the Sommerfeld enhancement. The Sommerfeld corrections are shown to cause two different effects depending on the nature of the strong force between the pair of incoming scalar quarks. In the case of annihilations between the same type of squarks, the enhancement turns into a reduction of the cross section as here the squarks experience a strong repulsive force. If the scalar quarks are different, however, the cross section is strongly enhanced due to the attractive strong force. Finally, we have investigated the impact of these corrections on the predicted relic density. We have demonstrated in our typical scenarios that even with the Sommerfeld reduction, the corrections are larger than the experimental uncertainty. In case of an enhancement, the corrections cause a 25% shift in the preferred parameter region where the relic density satisfies the experimental constraints. 



The authors would like to thank A. Pukhov for providing us with the possibility to interface our DM@NLO code with micrOMEGAs. J.H. was  supported  by  the  DFG  Emmy  Noether  Grant  No.  HA  8555/1-1, and earlier by the Labex ILP (reference ANR-10-LABX-63) part of the Idex SUPER, and received financial state aid managed by the Agence Nationale de la Recherche (ANR), as part of the programme Investissements d'avenir under the reference ANR-11-IDEX-0004-02.  The work of B. H. is supported by Investissements d'avenir, Labex ENIGMASS, contrat ANR-11-LABX-0012. The work of S. S. was supported by the DFG through RTG 2149 Strong and Weak Interactions - from Hadrons to Dark Matter. The figures presented in this paper have been generated using MatPlotLib <cit.>.




§ SMALL VELOCITY EXPANSION OF THE BOX CONTRIBUTION





In order to subtract the velocity enhanced part of the NLO contribution that is already included in the Sommerfeld resummation, we expand the corresponding contribution of the box diagrams in the relativistic relative velocity <cit.>

    v_rel = √(λ(s,m_1^2,m_2^2))/s-m_1^2-m_2^2 ,

with λ(a,b,c) = a^2 + b^2 + c^2 - 2ab - 2 ac - 2bc.

All box diagrams that contain the velocity enhanced contribution feature an exchange of a massless gluon between the incoming pair of scalar quarks. A generic diagram showing the masses of internal and loop particles is shown in Fig. <ref>. The matrix element contains tensor coefficients 

    D_i ( p_10^2,p_21^2,   p_32^2,p_30^2,p_20^2,p_31^2,M_0^2,M_1^2,M_2^2,M_3^2 )  = 
        D_i ( m_1^2,m_3^2,m_4^2,m_2^2,t_i,s,0,m_1^2,M_2^2,m_2^2 )  ,

where t_i are the Mandelstam variables t or u depending on the box diagram. The velocity enhanced terms in the box contribution are contained in the scalar 4-point integrals D_0 ( m_1^2,m_3^2,m_4^2,m_2^2,t_i,s,0,m_1^2,M_2^2,m_2^2 ) <cit.>. The full enhanced box matrix element consisting of box diagrams, where different gauginos with mass m_χ are exchanged, can be written using the corresponding color decomposed tree-level matrix element (see Eq. (<ref>)) as

    M_box   = ∑_t_i∑_χ M_box^χ, t_i = 
       (C^[3]_box M^Tree,χ, t_i_3  C^{3,3}_stij + C^[6]_box M^Tree,χ, t_i_6  C^{6,6}_stij)×
        2α_s/4π(s-m_1^2-m_2^2)(t_i-m_χ^2) D_0 ,

where the tensor integral has the arguments as in Eq. (<ref>) with M_2=m_χ and the color factors C^[𝐑]_box are given as

    C^[3]_box = -N_c+1/2N_c ,         C^[6]_box = N_c-1/2N_c .

The scalar integral for the specific arguments from Eq. (<ref>) can be written as

    D_0    = x_13/m_1m_2 (t_i-m_χ^2)(1-x_13^2)×
       {2ln(x_13)[-c_ϵ/ϵ-ln(μ m_χ/m_χ^2-t_i)+ln(1-x_13^2)].
       .+ln^2(x_12)+ln^2(x_23)+Li_2(x_13^2). 
       .+∑_k,l=±1ℒi_2(x_13,x_12^k,x_23^l)-π^2/6} .

The generalized polylogarithm in Eq. (<ref>) is defined as <cit.>


    ℒi_2 (x_1,…,x_n) =    Li_2(1-∏_i=1^n x_i) + [ln(∏_i=1^n x_i)-∑_i=1^nln(x_i)][ln(1-∏_i=1^n x_i)
        -θ(|∏_i=1^n x_i|-1)×.
        
       .(ln(-∏_i=1^n x_i)-1/2ln(∏_i=1^n x_i)-1/2∑_i=1^nln(x_i))
        ] .


The variables x_ij are defined using the loop masses M_i and M_j as well as the invariant combinations of 4-momenta p^2_ij, as given in Eq. (<ref>), as

    x_ij= √(1-4M_iM_j/p_ij^2-(M_i-M_j)^2)-1/√(1-4M_iM_j/p_ij^2-(M_i-M_j)^2)+1 .

In our case the only velocity dependent x_ij is x_13 which for M_1=m_1, M_2=m_χ and M_3=m_2 gives

    x_13   = -2(s-m_1^2-m_2^2)-2√(λ(s,m_1^2,m_2^2))/4m_1m_2
       = -(1-v_rel)/√(1-v_rel^2) ,

having used the relative velocity as defined in Eq. (<ref>). Given the color factors C^[𝐑]_box are independent of the exchanged gaugino, we can simplify the enhanced matrix element as

    2(M_boxM^†_Tree   ) =  ∑_t_i,t_j∑_χ,χ' 2(M_box^χ,t_i(M^χ', t_j_Tree)^†)  
     = ∑_t_i,t_j∑_χ,χ'(C^[3]_box    (M^χ,t_i_Tree,3(M^χ',t_j_Tree,3)^†) 
        + C^[6]_box (M^χ,t_i_Tree,6(M^χ',t_j_Tree,6)^†)) 2 (F^χ,t_i_box) ,

where F^χ,t_i_box is given by

    F^χ,t_i_box   =   α_s/π(s-m_1^2-m_2^2)/2m_1m_2 x_13/(1-x_13^2)×
          {2ln(x_13)[-c_ϵ/ϵ-ln(μ m_χ/m_χ^2-t_i)+ln(1-x_13^2)].
          .+ln^2(x_12)+ln^2(x_23)+Li_2(x_13^2). 
          .+∑_k,l=±1ℒi_2(x_13,x_12^k,x_23^l)-π^2/6} .

We first expand the expression in the relative velocity v_rel retaining just the leading term. The pre-factor can be expressed in terms of the relative velocity using

    (s-m_1^2-m_2^2)/2m_1m_2 x_13/(1-x_13^2) = -1/2v_rel .

Taking the real part of the F^χ,t_i_box factor results in

    2 (F^χ,t_i_box)    = -α_s/π1/v_rel {ln^2(x_12)+ln^2(x_23)+Li_2(1)
       +∑_k,l=±1ℒi_2(-1,x_12^k,x_23^l)-π^2/6} .
 
This expression seems to be implicitly dependent on the mass of the gaugino m_χ through the variables x_12 and x_23. However, this dependence vanishes after a more careful analysis, making the factor 2 (F^χ,t_i_box) generic for all underlying hard processes.

We will show the universality explicitly for a simple case where m_1=m_2 (i.e. x_12 = x_23) and where all x_ij are purely real. In such a case ℒi_2(-1,x_12^k,x_23^l) reduces to a simple polylogarithm Li_2(1+x_12^k+l). The factor then reduces to 

    2 (F^χ,t_i_box)    = -α_s/π1/v_rel {2ln^2(x_12)+2 Li_2(2)
       + Li_2(1+x_12^2)+Li_2(1+x_12^-2)} ,
 
which can be simplified using <cit.>

    Li_2(z) = -Li_2(z/z-1)-1/2ln^2(1-z)   z<1  .

The use of this identity eliminates all dependence on x_12 and the factor greatly simplifies to

    2 (F^χ,t_i_box) = -α_s/π1/v_relπ^2 .

Even though we have derived this particularly simple result in a special case, the same can be obtained in the most general case as well. With such a universal factor
the one-loop contribution to the squared matrix element from the enhanced box contribution is just

    2(M_boxM^†_Tree) = (-α_s/v_relπ) (C^[3]_box    |M^Tree_3|^2 
       + C^[6]_box |M^Tree_6|^2) .

This expression is compatible with the next-to-leading part of the Sommerfeld enhancement after we realize that in the non-relativistic case the relative velocity can be easily related to the velocity used in the Sommerfeld enhancement for identical incoming particles (see Eq. (<ref>)) as v_rel = 2v_ s.



99
LUX:2016
  D. S. Akerib et al. [LUX Collaboration],
  
  Phys. Rev. Lett.  118, no. 2, 021303 (2017)
  
  [arXiv:1608.07648 [astro-ph.CO]].
  

XENON:2019
    E. Aprile et al. [XENON Collaboration],
     
    arXiv:1902.03234 [astro-ph.CO]
    

Aghanim:2018eyx 
  N. Aghanim et al. [Planck Collaboration],
  
  arXiv:1807.06209 [astro-ph.CO].
  

Hall:2009bx
  L. J. Hall, K. Jedamzik, J. March-Russell and S. M. West,
  
  JHEP 1003 (2010) 080
 
  [arXiv:0911.1120 [hep-ph]].
  
  


Klasen:2013ypa 
  M. Klasen and C. E. Yaguna,
  
  JCAP 1311, 039 (2013)
  
  [arXiv:1309.2777 [hep-ph]].
  
  
Bernal:2017kxu
  N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen,
  
  Int. J. Mod. Phys. A 32 (2017) no.27,  1730023
 
  [arXiv:1706.07442 [hep-ph]].
  
  
  
Belanger:2018sti
  G. Bélanger et al.,
  
  JHEP 1902 (2019) 186
 
  [arXiv:1811.05478 [hep-ph]].
  
  


Gondolo:1990dk 
  P. Gondolo and G. Gelmini,
  
  Nucl. Phys. B 360, 145 (1991).
  
  

Griest:1990kh 
  K. Griest and D. Seckel,
  
  Phys. Rev. D 43, 3191 (1991).
  
  

Edsjo:1997bg 
  J. Edsjo and P. Gondolo,
  
  Phys. Rev. D 56, 1879 (1997)
  
  [hep-ph/9704361].
  

Cannoni:2013bza
  M. Cannoni,
  
  Phys. Rev. D 89, 103533 (2014)
  
  [arXiv:1311.4508 [hep-ph]].
  


Delepine:1996vn
  D. Delepine, J. M. Gerard, R. Gonzalez Felipe and J. Weyers,
  
  Phys. Lett. B 386 (1996) 183
  
  [hep-ph/9604440].
  
  
Liebler:2015ddv
    S. Liebler, S. Profumo and T. Stefaniak,
    
    JHEP 1604 (2016) 143
    
    [arXiv:1512.09172 [hep-ph]].
    
    
Arbey:2012bp
  A. Arbey, M. Battaglia, A. Djouadi and F. Mahmoudi,
  
  Phys. Lett. B 720 (2013) 153
  
  [arXiv:1211.4004 [hep-ph]].
  
  
Herrmann:2007ku 
  B. Herrmann and M. Klasen,
  
  Phys. Rev. D 76, 117704 (2007)
  
  [arXiv:0709.0043 [hep-ph]].
  
  
Herrmann:2009wk 
  B. Herrmann, M. Klasen and K. Kovařík,
  
  Phys. Rev. D 79, 061701 (2009)
  
  [arXiv:0901.0481 [hep-ph]].
  
  
Herrmann:2009mp 
  B. Herrmann, M. Klasen and K. Kovařík,
  
  Phys. Rev. D 80, 085025 (2009)
  
  [arXiv:0907.0030 [hep-ph]].
  

Herrmann:2014kma 
  B. Herrmann, M. Klasen, K. Kovařík, M. Meinecke and P. Steppeler,
  
  Phys. Rev. D 89, no. 11, 114012 (2014)
  
  [arXiv:1404.2931 [hep-ph]].
  

Freitas:2007sa 
  A. Freitas,
  
  Phys. Lett. B 652, 280 (2007)
  
  [arXiv:0705.4027 [hep-ph]].
  
  
Harz:2012fz 
  J. Harz, B. Herrmann, M. Klasen, K. Kovařík and Q. L. Boulc'h,
  
  Phys. Rev. D 87, no. 5, 054031 (2013)
  
  [arXiv:1212.5241 [hep-ph]].
  

Harz:2014tma 
  J. Harz, B. Herrmann, M. Klasen and K. Kovařík,
  
  Phys. Rev. D 91, no. 3, 034028 (2015)
  
  [arXiv:1409.2898 [hep-ph]].
  

Harz:2014gaa 
  J. Harz, B. Herrmann, M. Klasen, K. Kovařík and M. Meinecke,
  
  Phys. Rev. D 91, no. 3, 034012 (2015)
  
  [arXiv:1410.8063 [hep-ph]].
  
  
Boudjema:2005hb 
  F. Boudjema, A. Semenov and D. Temes,
  
  Phys. Rev. D 72, 055024 (2005)
  
  [hep-ph/0507127].
  
  
Baro:2007em 
  N. Baro, F. Boudjema and A. Semenov,
  
  Phys. Lett. B 660, 550 (2008)
  
  [arXiv:0710.1821 [hep-ph]].
  

Baro:2009na 
  N. Baro, F. Boudjema, G. Chalons and S. Hao,
  
  Phys. Rev. D 81, 015005 (2010)
  
  [arXiv:0910.3293 [hep-ph]].
  
  
Harz:2016dql 
  J. Harz, B. Herrmann, M. Klasen, K. Kovařík and P. Steppeler,
  
  Phys. Rev. D 93, no. 11, 114023 (2016)
  
  [arXiv:1602.08103 [hep-ph]].
  

ATLAS2015pMSSM 
  G. Aad et al. [ATLAS Collaboration],
  
  JHEP 1510, 134 (2015)
  
  [arXiv:1508.06608 [hep-ex]].
  

CMS2016pMSSM 
  V. Khachatryan et al. [CMS Collaboration],
  
  JHEP 1610, 129 (2016)
  
  [arXiv:1606.03577 [hep-ex]].
  

Porod:2003um
  W. Porod,
  
  Comput. Phys. Commun.  153 (2003) 275
  
  [hep-ph/0301101].
  

Porod:2011nf
  W. Porod and F. Staub,
  
  Comput. Phys. Commun.  183 (2012) 2458
  
  [arXiv:1104.1573 [hep-ph]].  
  

Belanger:2001fz
  G. Bélanger, F. Boudjema, A. Pukhov and A. Semenov,
  
  Comput. Phys. Commun.  149 (2002) 103
  
  [hep-ph/0112278].
  
  
Belanger:2006is
  G. Bélanger, F. Boudjema, A. Pukhov and A. Semenov,
  
  Comput. Phys. Commun.  176 (2007) 367
  
  [hep-ph/0607059]. 
  

Allanach:2008qq
  B. C. Allanach et al.,
  
  Comput. Phys. Commun.  180 (2009) 8
  
  [arXiv:0801.0045 [hep-ph]].
  

Beneke:2009rj
  M. Beneke, P. Falgari and C. Schwinn,
  
  Nucl. Phys. B 828 (2010) 69
  
  [arXiv:0907.1443 [hep-ph]].
  
  
Carena:1999py
  M. Carena, D. Garcia, U. Nierste and C. E. M. Wagner,
  
  Nucl. Phys. B 577 (2000) 88
  
  [hep-ph/9912516].
  

Siegel:1979wq
  W. Siegel,
  
  Phys. Lett.  84B (1979) 193.
  
  
  
Stockinger:2005gx
  D. Stöckinger,
  
  JHEP 0503 (2005) 076
  
  [arXiv:hep-ph/0503129].
  
  
Passarino:1978jh
  G. Passarino and M. J. G. Veltman,
  
  Nucl. Phys. B 160 (1979) 151.
  

Denner:1991kt
  A. Denner,
  
  Fortsch. Phys.  41 (1993) 307
  [arXiv:0709.1075 [hep-ph]].
  

Denner:2010tr
  A. Denner and S. Dittmaier,
  
  Nucl. Phys. B 844 (2011) 199
  [arXiv:1005.2076 [hep-ph]].
  

Dittmaier:2003bc
  S. Dittmaier,
  
  Nucl. Phys. B 675 (2003) 447
  [arXiv:hep-ph/0308246].
  

FeynArts
      T. Hahn,
      
      Comput. Phys. Commun. 140 (2001) 418
      [arXiv:hep-ph/0012260].  
      

FeynCalc
    R. Mertig, M. Bohm and A. Denner,
    
    Comput. Phys. Commun.  64, 345 (1991).

    

FORM
      J. A. M. Vermaseren,
      
      arXiv:math-ph/0010025.
      
      
GieleGlover
      W. T. Giele and E. W. N. Glover,
      
      Phys. Rev. D 46 (1992) 1980.
      
      
HarrisOwens
  B. W. Harris and J. F. Owens,
  
  Phys. Rev. D 65 (2002) 094032
  [arXiv:hep-ph/0102128].
  

tHooft:1978xw
  G. 't Hooft and M. J. G. Veltman,
  
  Nucl. Phys. B 153 (1979) 365.
  

Sommerfeld:1931
    A. Sommerfeld, 
    
    Annalen der Physik, Vol. 403, Issue 3, pp.257-330.


Fischler:1977yf
  W. Fischler,
  
  Nucl. Phys. B 129 (1977) 157.
  
  

Billoire:1979ih
  A. Billoire,
  
  Phys. Lett.  92B (1980) 343.
  
  

Kauth:2011bz
  M. R. Kauth, A. Kress and J. H. Kuhn,
  
  JHEP 1112 (2011) 104
  
  [arXiv:1108.0542 [hep-ph]].
  

Beneke:2010da
  M. Beneke, P. Falgari and C. Schwinn,
  
  Nucl. Phys. B 842 (2011) 414
  
  [arXiv:1007.5414 [hep-ph]].
  

Cassel:2009wt
  S. Cassel,
  
  J. Phys. G 37 (2010) 105009
  
  [arXiv:0903.5307 [hep-ph]].
  

Hagiwara:2008df
  K. Hagiwara, Y. Sumino and H. Yokoya,
  
  Phys. Lett. B 666 (2008) 71
  
  [arXiv:0804.1014 [hep-ph]].
  

Kiyo:2008bv
  Y. Kiyo, J. H. Kuhn, S. Moch, M. Steinhauser and P. Uwer,
  
  Eur. Phys. J. C 60 (2009) 375
  
  [arXiv:0812.0919 [hep-ph]].
  

Harz:2017dlj
  J. Harz and K. Petraki,
  
  Phys. Rev. D 97 (2018) no.7,  075041
  
  [arXiv:1711.03552 [hep-ph]].
  

Harz:2018csl
  J. Harz and K. Petraki,
  
  JHEP 1807 (2018) 096
  
  [arXiv:1805.01200 [hep-ph]].
  

Harz:2019rro
  J. Harz and K. Petraki,
  
  arXiv:1901.10030 [hep-ph].
  

Matplotlib
  J. D. Hunter, Matplotlib: A 2D graphics environment, Computing In Science & Engineering 9, 3 (2007).
  

Cannoni:2016hro
  M. Cannoni,
  
  Int. J. Mod. Phys. A 32 (2017) no.02n03,  1730002
  
  [arXiv:1605.00569 [hep-ph]].

Lewin:1981  
    L. Lewin, 
    Polylogarithms and Associated Functions, New York: North-Holland (1981).





