Gauge-Independent Approach to Resonant Dark Matter Annihilation

In spontaneously broken gauge theories, transition amplitudes describing dark-matter (DM) annihilation processes through a resonance may become highly inaccurate close to a production threshold, if a Breit-Wigner (BW) ansatz with a constant width is used. To partially overcome this problem, the BW propagator needs to be modified by including a momentum dependent decay width. However, such an approach to resonant transition amplitudes generically suffers from gauge artefacts that may also give rise to a bad or ambiguous high-energy behaviour for such amplitudes. We address the two problems of gauge dependence and high-energy unitarity within a gauge-independent framework of resummation implemented by the so-called Pinch Technique. We study DM annihilation via scalar resonances in a gauged U(1)$_X$ complex-scalar extension of the Standard Model that features a massive stable gauge field which can play the role of the DM. We find that the predictions for the DM abundance may vary significantly from previous studies based on the naive BW ansatz and propose an alternative simple approximation which leads to the correct DM phenomenology. In particular, our results do not depend on the gauge-fixing parameter and are consistent with considerations from high-energy unitarity.


Introduction
Unveiling the nature of Dark Matter (DM) constitutes one of the biggest challenges in Cosmology and possibly in Particle Physics. Even though several pieces of evidence coming from both cosmological and astrophysical scales confirm its presence, the actual composition of the DM itself remains elusive to us thus far. In numerous models, the DM is assumed to be a new kind of massive particles that were in thermal equilibrium with the Standard Model (SM) particles in the early Universe. In such thermal DM models, the DM relic abundance crucially depends on the rate at which DM particles annihilate into thermal bath states. This rate is proportional to the thermally averaged DM-pairannihilation cross-section summed over all possible final states. Most remarkably, if this cross-section has a value which is typical to the one governing SM weak interactions, then the predicted DM density turns out to be in the right ballpark in agreement with observations. Note that the most accurate determination of the DM relic density comes from measurements of the Cosmic Microwave Background (CMB) by the Planck satellite in conjunction with a number of other experiments and astrophysical data (for a recent analysis, see [1]).
A popular DM scenario [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] is based on the working hypothesis that there exists a mediator with couplings to both DM and SM particles, and whose mass happens to be approximately twice that of the DM. In such scenarios, one has to consider DM annihilation processes at energies where the tree-level propagator of the mediator becomes singular. A standard way to avoid such singularities is to employ a Breit-Wigner (BW) ansatz for the propagator with a constant decay width [17]. Then, the amplitude gets regularized by the non-zero width of the mediator and can thus get drastically enhanced by many orders of magnitude, when compared to other non-resonant contributions. In Quantum Field Theory (QFT), the BW form of the propagator usually arises from a Dyson series summation of self-energy graphs of the mediator. In the so-called on-mass-shell (OS) scheme of renormalization [18], the dispersive parts of the self-energies renormalize the masses, whilst their absorptive part is related to the decay width of the mediator.
From the phenomenological point of view, the channel of resonant DM annihilation turns out to be an attractive option, as it leads to suppressed DM-nucleon crosssections thereby avoiding the tight constraints emanating from the null results of direct DM searches. Therefore, it should not be too surprising that this resonant region in question may become the only viable region in the parameter space of a given model with DM mass in the GeV range that survives after all direct detection limits on the DM-nucleon cross-section were imposed [19,20].
As was first noted in [14], the BW approximation with a constant decay width in the propagator of the mediator can become very inaccurate close to the DM production threshold, especially when the respective DM channel contributes significantly to the decay width of the mediator. To partially remedy this problem, one is compelled to use an effective BW propagator with momentum-dependent or running width for the exchanged particle in the s-channel. In this article, we will go beyond the previously used non-relativistic approach [14]. In general, an s-dependent width results from the imaginary (absorptive) part of the self-energy of the mediator. This quantity is only gauge-independent when evaluated at the pole of the propagator, but becomes gauge-variant in the off-shell region. This was a well-known problem in QFT and pertains to the question whether a consistent gauge-independent definition of off-shell Green's function for unstable particles exists in spontaneously broken gauge theories [21]. To deal with this issue, a number of recipes and methods have been put forward by several authors, such as the Laurent series expansion [22,23], the complex mass scheme [24,25], the fermion loop scheme [26], and the effective theory approach [27].
In this paper we discuss the problems that arise in the relativistic treatment of resonant DM annihilation processes in spontaneously broken gauge theories and show how these can be avoided in a resummation approach implemented by the PT. As an archetypal model, we consider a gauged U(1) X complex-scalar extension of the SM that includes a massive stable gauge boson X as a candidate particle for a Vector DM (VDM). We explicitly demon-strate how the transition amplitude for DM-pair annihilation has the proper high-energy unitarity limit in compliance with the Equivalence Theorem [43,44] and its generalized version (GET) [45,46]. By virtue of the PT resummation method adopted in this paper, we obtain predictions for the yield of DM abundance that may vary significantly from previous considerations based on the naive BW approximation. In particular, we illustrate the gauge independence of our results and their consistency with high-energy unitarity. In the same context, we present an alternative simple approximation which leads to the correct DM phenomenology.
The paper is organized as follows. In sec. 2 we discuss the process of resonant DM-pair annihilation and the problems that arise in the relativistic treatment of a VDM in the vicinity of a resonance. In sec. 3 we describe the PT resummation approach and apply it to transition amplitudes that are typical to DM annihilation processes, such as XX →f f , where f is a SM fermion. By means of this approach, we obtain a Born-improved amplitude which is gauge independent and possesses the proper asymptotics in the high-energy limit. In sec. 4 we use the Born-improved amplitude derived in the previous section to compute the annihilation cross-sections and DM relic density. This enables us to assess the significance of our results by comparing them with those derived with other methods. In sec. 5 the key points of our analysis were summarised. Technical details of this study were relegated to a number of appendices. Specifically, appendix A provides further details of the VDM model under study, appendix B presents the Feynman rules for this model in the R ξ and background field gauges (BFG), and appendix C contains analytical expressions of one-loop vertex corrections in the PT. Finally, appendix D complements our proof of the GET for the DM annihilation process XX →f f .

Resonant Dark Matter Annihilation
Annihilation of DM in the vicinity of a resonance has certain features that renders it distinct from a generic scenario of a thermal DM. First, the thermally averaged crosssection displays strong dependence on the temperature T that can lead to a prolonged period of effective DM annihilation which substantially changes the comoving DM density, even though the DM itself may have already decoupled from the thermal bath. Second, the cross-section may be significantly enhanced at small velocities leading to a strengthening of the signal probed by DM indirect searches [5,8]. Finally, the coupling of the DM to SM particles can become rather suppressed. This last property can be quite challenging for collider or direct detection experiments as it usually gives rise to weaker constraints for this area of parameter space. It also raises the temperature of kinetic decoupling which, in connection with the high T -variability of the annihilation cross-section, may affect the predictions of relic density [14,47,48].
In the Born approximation, the transition amplitude for a resonant DM annihilation process does not describe consistently the dynamics of a possible unstable mediator with mass M , because its tree-level propagator ∆ 0 (s) ≡ (s − M 2 ) −1 becomes singular at s = M 2 . The standard way to treat this singularity is to perform a Dyson series summation of the mediator's self-energy Π(s) which results in a replacement of the tree-level propagator ∆ 0 (s) with a resummed one, ∆(s) ≡ (s − M 2 + Π(s)) −1 . In the OS scheme of renormalization [18], this modification amounts to the well-known BW approximation [17], with Π(s) ≈ m Π(s = M 2 ) = iZ −1 M Γ, where Z is the wave-function renormalization (set to 1 at the one-loop level), and M and Γ are the renormalized mass and width of the resonant particle in this scheme. In this way, one obtains a finite analytic expression for the amplitude in the resonance region. For instance, the cross-section for a 2 → 2 annihilation process, which proceeds via the s-channel and has two identical DM particles of mass m i in the initial state and two particles of equal mass m f in the final state, is approximately given by (2.1) In the above, we followed the notation and conventions of [5], after taking into account an extra factor of 2 because of identical particles in the initial state. Hence, B i and B f in (2.1) are the branching ratios for the mediator state to decay into initial and final states, denoted as i and f , respectively. Nevertheless, when the dominant contribution to the mediator's self-energy Π(s) comes from particles which are also the initial states of the transition amplitude, then the BW approximation is in general not applicable [14]. In this case, one finds high inaccuracies in the transition amplitudes that result from the phase space factor 1 − s/(4m 2 ) present in Π(s), which varies substantially for s > ∼ 4m 2 . To deal with this problem, one has to consider the energy-dependence of the imaginary part of the self-energy, which is sometimes described by the running width of the mediator Γ(s) ≡ m Π(s)/M . However, this approach encounters other serious theoretical problems. As the resummation process relies on taking a subset of graphs from each order of the perturbative expansion, the subtle cancellations which guarantee the gauge-independence of the amplitudes at each order of the expansion are spoiled and the resulting expression has explicit dependence on the gauge fixing parameter(s) chosen to calculate Π(s). The gauge-dependent amplitude results in an ambiguous annihilation cross-section, and as such its physical relevance becomes rather obscured.
In order to address the above issue, we adopt the PT framework which enables us to calculate such resonant DM annihilation amplitudes following a gauge-independent and self-consistent approach. As mentioned in the introduction, the so-derived Born-improved amplitudes will satisfy several desirable field-theoretic properties as a consequence of the analyticity, unitarity and gauge-invariance of the S-matrix. Our approach implemented by the PT will be illustrated within a VDM model [49][50][51][52][53][54]. This is an extension of the SM augmented by local U(1) X symmetry and by a complex scalar field S. The field S has a portal interaction with the SM Higgs doublet, and also develops a non-zero vacuum expectation value (VEV) that breaks U(1) X spontaneously. The resulting gauge boson X acquires a mass due to the Higgs mechanism and can be made stable, provided its kinetic mixing with the hypercharge gauge boson is forbidden by a Z 2 symmetry. In this U(1) X extension of the SM, the vector boson X can play the role of a viable DM. This VDM model contains two scalar mass eigenstates, denoted as h 1 and h 2 , where the first scalar Figure 1. The tree-level (left) and the Born-improved (right) Feynman diagram for the process XX → ff .
represents the SM-like Higgs boson, while the second one will serve as a resonant mediator. The couplings of the h 1 and h 2 scalars to SM matter and VDM are governed by the interaction Lagrangian, where α is the scalar mixing angle and g x is the U(1) X gauge coupling. Our primary focus in this analysis will be on the resonance region, where m 2 ≈ 2M X with a scalar mixing angle α 1. Further details of the VDM model under consideration may be found in appendix A.

Gauge Independence and High Energy Unitarity
In this section, we will describe a gauge-invariant resummation approach to resonant transition amplitudes, within the Pinch Technique (PT) framework. As an illustrative example, we will be studying the annihilation process XX →f f , where X is a U(1) X gauge field which assumes the role of the DM and f collectively stands for a SM fermion. In the U(1) X scenario under study, the DM annihilation process, XX →f f , proceeds via the exchange of a coupled system of two mixed scalar resonances h 1 and h 2 . As illustrated in the left diagram of fig. 1, the tree-level amplitude of such a reaction reads Here, summation over the repeated indices i, j = 1, 2 is implied, m 1,2 are the masses of the Higgs scalars h 1,2 , and V XXh i µν and V h jf f denote the tree-level expressions for the XXh i and h jf f vertices, respectively. Moreover, µ (p 1 ) and ν (p 2 ) denote the polarizations of the X-vector bosons with four-momenta p 1 and p 2 , and u(r 2 ) and v(r 1 ) are the usual Dirac spinors for the fermion-antifermion pair f andf with momenta r 2 and r 1 , respectively.
As stated in sec. 2, the tree-level amplitude (3.1) becomes singular near the resonance regions s ≈ m 2 1,2 , requiring resummation of an infinite series of self-energies Π ij (s). Our aim is to calculate a Born-improved amplitude for the DM annihilation process XX → f f , which does not depend on the gauge-fixing parameter ξ X used to parameterize the unphysical degrees of freedom of the gauge field X. In this context, we will adopt a gaugeindependent resummation approach (illustrated in the right diagram of fig. 1) which is implemented by the PT in order to properly treat the coupled system of the two resonances, h 1 and h 2 , in a fashion analogous to [33]. In this framework, in sec. 3.1, we obtain a resummed propagator which is a 2 × 2 matrix and is denoted as ∆ ij . Subsequently, in sec. 3.2 we discuss the necessity of including one-loop corrections to the XXh i vertices and present the tree-level Ward Identities (WIs) satisfied by the PT self-energies and PT vertices which ensure the validity of the ET [43,44], as well as of the GET [45,46]. Finally, we show in sec. 3.3 that the so-derived Born-improved amplitude XX →f f exhibits proper high-energy asymptotic behaviour as expected from unitarity considerations.

Pinch Technique Resummation
To start with, we first calculate the X-field contribution to the h i h j self-energies Π ij (s) in the VDM model of interest to us. In the renormalizable class of R ξ gauges, upon inclusion of the would-be Goldstone bosons G X and the pertinent ghost fields, these are given by where ξ X is the gauge-fixing parameter associated with the U(1) X gauge field. Resummation of this self-energy inevitably leads to a gauge-dependent resonant amplitude and the appearance of unphysical thresholds at s = 4 ξ X M 2 X , when ξ X = 1. We may now calculate X, Z, W , t, h k , h l contributions to the h i h j self-energies Π ij (s) employing the PT, as done in [37,38]. For the case of the VDM, these are given by where is the 't Hooft-Veltman function [55] defined in n = 4−2 dimensions and the coupling V h ikl is specified in table 1. Note that the above results for the PT self-energies Π ij (s) are also in agreement with that obtained in [39], upon appropriately replacing the couplings of that model. Moreover, the same results can be derived by making use of the established equivalence between the PT and the covariant Background Field Gauge (BFG) for ξ Q = 1 [56][57][58][59]. For reader's convenience, the Feynman rules in the BFG are presented in Appendix B. Evidently, the h i h j self-energies Π ij (s) do not depend on the gauge-fixing parameters ξ X , ξ W , ξ Z and all display an absorptive part at the physical thresholds s = 4M 2 The full PT self-energies Π ij (s) can then be systematically resummed, yielding the PT resummed propagator In the above, we have suppressed the indices i, j = 1, 2 labelling the Higgs scalars h 1,2 and the s-dependence of the propagators and self-energies.
2 ) −1 is the diagonal tree-level propagator matrix. From (3.9), one then gets . As we will see in the next section, in addition to its gauge independence, a Born-improved amplitude for XX →f f must have the proper high-energy asymptotic behaviour, as dictated by the ET.

The Generalized Equivalence Theorem
Replacing naively the tree-level propagator matrix ∆ 0 (s) with the resummed one ∆(s) in (3.1) modifies the transition amplitude A XX→f f (s) not only in the vicinity of the resonance region, but also changes drastically its high-energy limit as s/(4M 2 X ) → ∞. This is a generic problem for most forms of BW propagators [60,61] and is due to the presence of non-zero self-energies Π ij (s) or decay widths in ∆(s) that may posses non-trivial s-dependence. The latter distort subtle cancellations which are triggered by the WIs that result from the gauge invariance of the classical action, thus leading to a different highenergy limit from the one expected in the Born approximation.
To ensure the proper high-energy asymptotic behaviour of a scattering process involving massive gauge fields in the initial or final state, the amplitude of such process must obey the GET [45,46], which is a consequence of the classical WIs mentioned above. Applying the GET to the amplitude of the process XX →f f gives Equation (3.11) establishes a relation between the amplitude with initial longitudinally polarized X-bosons, denoted here as X µ L (p 1 ) = µ L (p 1 ) and X ν L (p 2 ) = ν L (p 2 ), and the amplitudes with corresponding would-be Goldstone bosons G X (p 1,2 ) or the energetically suppressed remainders, x µ (p 1 ) and x ν (p 2 ). Specifically, x µ (p 1 ) is defined as which satisfies the identities Note that an analogous definition and set of identities hold for x ν (p 2 ) as well.
Proceeding as in [37,38], we may reinforce the validity of the GET stated in (3.11) by considering one-loop h i XX, h i XG X and h i G X G X vertices within the PT framework. More explicitly, the following substitutions for the tree-level h i -couplings need to be considered: where the relevant PT one-loop vertices, were calculated by means of the BFG method and presented in Appendix C. One can show that when PT resummed vertices as given in (3.14)-(3.16) are employed, the relation (3.11) is satisfied in terms of the PT resummed (Born-improved) amplitudes.
The proof of (3.11) for the Born-improved amplitudes relies on the fact that the PT one-loop vertices V h i XX µν , V h i XG X µ and V h i G X G X , as well as the PT self-energies Π ij , Π XG X µ and Π G X G X , satisfy tree-like WIs that are identical to those derived from the classical action. In detail, these tree-like PT WIs read It is not difficult to verify that the PT WIs (3.17)- (3.20) are indeed satisfied by the tree-level couplings of the theory, and V h i G X G X , after making the replacements: Further details pertinent to the proof of (3.11) are given in Appendix D.

High-Energy Asymptotics of the Born-Improved Amplitude
We will now study the high-energy behaviour of the PT resummed amplitude for the process XX →f f . The high-energy asymptotics of such an amplitude is similar to the one of the SM processf f → ZZ mediated by the Higgs boson, which was considered in [37,38] and studied in more detail in [62].
As we are interested in deriving a minimal Born-improved amplitude, we will ignore the real (dispersive) parts of the self-energies and vertex corrections. In the OS scheme, these are either subdominant or suppressed by higher orders in the resonance region [33]. Hence, we will only consider the imaginary (absorptive) parts of the PT resummed propagators ∆ ij (s) and the PT vertices Γ h i XX µν . In this minimal self-consistent PT framework, the Born-improved amplitude takes on the form p 1 , p 2 ) are the PT one-loop vertices, whose absorptive parts are given in Appendix C.
To analyse the high-energy asymptotics of A XX→f f , we consider the longitudinal X boson polarizations which in the centre-of-mass (CoM) frame can be expanded as follows (3.22) For instance, taking into account the absorptive effects due to the opening of the XX threshold only, the amplitude exhibits the following asymptotic behaviour: where g is SU(2) L gauge coupling and M W is the W ± -boson mass. This should be contrasted with the corresponding tree-level result, It is not difficult to see that up to higher order terms O(gg 3 x ), the PT resummed amplitude given in (3.23) approaches the tree-level result (3.24) in the high-energy limit. Evidently, the Born-improved amplitude (3.21) displays the expected high-energy asymptotics.
It is instructive to see how the energetically constant terms ∝ s 0 cancel in the highenergy limit of both the tree-level and Born-improved amplitude by virtue of the PT WIs. In the Born approximation, such cancellation is a consequence of the orthogonality of the mixing matrix R, since Beyond the Born approximation, we may employ the PT WI (3.19) to derive an equivalent WI for the PT resummed vertex Γ h i XX µν in the high-energy limit, where we have ignored contributions from V h i G X G X (q, p 1 , p 2 ) and Π G X G X (p 1,2 ) that grow as ln(s/M 2 X ) or go to s 0 , for p 2 1 = p 2 2 = M 2 X [62]. With the help of (3.26), the high-energy limit of the Born-improved amplitude for the longitudinal X bosons may then easily be evaluated as follows: Consequently, the leading high-energy asymptotics of A X L (p 1 )X L (p 2 )→f f is expected to be proportional to s −1 ln(s/M 2 X ) or s −1 , as given in (3.23). The energetically subleading terms are not arbitrary either, but obey the GET stated in (3.11) [cf. Appendix D].
In addition to the resonant process XX → ff , there are also the vector and scalar annihilation channels, XX → (ZZ, W + W − ) and XX → h k h l , in the VDM model, as shown in fig. 2. In order to obtain the Born-improved amplitudes for these processes, we proceed as before. We replace the tree-level propagator ∆ 0 (s) in the s-channel graphs with the PT resummed propagator ∆(s) given in (3.10). Likewise, we replace the tree-level couplings with their PT resummed counterparts [38,62]: Note that the PT vertex corrections to h i ZZ and h i W + W − satisfy tree-like Ward identities analogous to those given in (3.17)-(3.20). As can be seen from fig. 2, besides the resonant s-channel graphs, there exist also nonresonant t-and u-channel diagrams and four-point vertices contributing to the processes XX → h k h l which do not require any improvement. The high-energy behaviour of these amplitudes can be studied using the GET. Because the structure of the vertex V h jkl contains a piece proportional to R 2j , there is no cancellation of the leading terms in energy for both tree-level and Born-improved s-channel diagrams as opposed to what happens for the process with ff in the final state. Therefore, we expect that such graphs and the total amplitude will go asymptotically to a constant (∝ s 0 ) at high energies. One the other hand, if we replace the initial states X with their respective Goldstone bosons G X , the presence of the propagator suppresses the s-channel amplitude and makes the latter fall as s −1 . Consequently, thanks to the ET, the high energy behaviour of the total amplitude is dictated by the non-resonant part of the amplitude with Goldstone bosons in the initial state, i.e. G X G X → h k h l . This part is not affected by resummation and coincides with the tree-level result, as can easily be inferred from the Feynman diagrams depicted in the second line of fig. 2. This implies that for s/(4M 2 X ) → ∞, the s-channel contribution to the total Born-improved amplitude, which we denote here as A s X L X L →h k h l , has to behave exactly as its tree-level counterpart, i.e. A s X L X L →h k h l ∝ s 0 . Indeed, this is the case, as the PT ensures that the classical WI (3.26) be satisfied at both tree and Born-improved levels. Therefore, we obtain the same high-energy limit for the Born-improved amplitude A s X L X L →h k h l as the one given by the tree level amplitude, i.e.
A s Observe that according to the ET, the high-energy limit of the total Born-improved amplitude A X L X L →h k h l will be given by the contact diagram G X G X h k h l which tends to a constant as s/(4M 2 X ) → ∞.

Annihilation Cross-Sections and Evolution of Dark Matter Density
In this section, we apply the PT resummation approach presented in the previous section in order to compute the DM annihilation cross-sections, as well as the evolution of the DM density within the VDM model. To this end, we adopt the Born-improved amplitude (3.21) with the PT resummed propagator and one-loop dressed vertices for all the annihilation channels, viz.
where w specifies the SM final state, i.e. w = Z, W, h i , f , and Γ h jw w represents the proper h jw w-vertex upon appropriate contraction with polarization vectors and spinors 1 . For the annihilation channel into scalars (when w = h i ), we also include the non-resonant contributions shown in fig. 2. The total annihilation cross-section σ is obtained as usual, by averaging the squared amplitude over initial X-boson polarizations and summing over all annihilation channels w in the final state, their polarizations and other possible internal degrees of freedom. Furthermore, we have introduced in our analysis the parameter which provides a measure of proximity of the h 2 -boson mass to the DM XX threshold.
1 Strictly speaking, one should include one-loop corrections to hiZZ-and hiW + W − -vertices within the PT framework similarly to what was done in sec. 3. However, these effects are found to be numerically negligible. Specifically, the one-loop corrections to hiXX-vertex change the cross-sections by as much as 2 %, when the dark gauge coupling gx approaches its perturbativity limit √ 2π. The respective one-loop absorptive corrections to hi vertices with the SM Z and W ± vector bosons scale as g 4 , so simple estimates indicate that these are at most at per mille level. Consequently, for most applications of phenomenological interest, the absorptive loop corrections to hiZZ-and hiW + W − -vertices may safely be neglected. 2π. Note that the inset plots display the respective cross-sections at the close vicinity of the DM threshold region. Figure 3 shows the product σv of the annihilation cross-section σ with the relative velocity v in the CoM frame for the processes XX → tt and XX → h 1 h 1 , as functions of the normalized energy squared parameter (4. 3) The results are obtained by utilising the PT resummed propagator, a resummed BW propagator with self-energies calculated in the R ξ class of gauges with ξ = 0, 1, 5 and in the unitary gauge (UG). In addition, fig. 3 displays the quantity σv computed at the tree-level and in the BW approximation assuming a constant or an s-dependent decay width (SDW).
In the BW case, we set Γ i = m Π ii (m 2 i )/m i . In the SDW approximation, we include in the partial decay width Γ h 2 →XX the exact form of the phase space factor corresponding to the XX channel. This gives the leading energy-dependent correction to Γ 2 near the XX threshold, i.e.
For δ > 0, the above expression has to be analytically continued. This SDW approximation works exceptionally well near the Higgs pole, where s ≈ m 2 2 . Moreover, it does not modify the high-energy behaviour of the tree-level amplitudes, because Γ 2 (s) approaches a constant value when s/(4M 2 X ) → ∞. From fig. 3, we can see that in the region where s/(4M 2 X ) 1, the PT result is close to the one that has been evaluated in the R ξ gauges with ξ = 0, 1, 5. However, at large s, the Born-improved (PT) cross-section (times the relative velocity v) for the process XX → tt varies between the one computed in the R ξ gauge and that in the tree-level approximation, whereas the cross-section calculated in the unitary gauge is suppressed. The latter is a result of the Z-and W ± -boson contributions to Π ij (s) which display a scaling behaviour ∝ s 2 /(16M 4 X ) when the energy is well above the corresponding W + W − and ZZ thresholds. For the annihilation channel XX → h 1 h 1 , the differences between the various results are smaller, because the leading (in s) parts of the propagator do not cancel asymptotically. Moreover, the PT result approaches exactly the tree-level one as discussed earlier. Observe that near the threshold the naive BW and tree-level results deviate dramatically from those found in the PT, in contrast to the approximation with s-dependent XX phase space factor (4.4) which leads to a very good agreement in this region. For δ < 0, the large width Γ 2 used in the BW approximation results in an underestimated cross-section right at the threshold. On the other hand, for δ > 0, the width Γ 2 contains only SM contributions which are suppressed by the small mixing angle α. For this reason, the BW approximation is close to the tree-level result. Finally, we note the presence of the fictitious threshold for ξ = 5. This is an example of a gauge artifact that originates from a naive resummation of self-energies in the R ξ class of gauges.
The evolution of DM density is mainly governed by the size of the thermally averaged annihilation cross-section σv rel , which may be computed as a function of x ≡ M X /T by  Figure 4. Thermally averaged total annihilation cross-section σv rel (x) as a function of x ≡ M X /T , for g x = 0.3, 1 and √ 2π and the remaining model parameters same as in fig. 3. The solid line corresponds to the cross-section (4.5) obtained using the PT with resummed propagator and corrected h i XX vertex while the dashed lines are for the self-energies calculated using R ξ gauges for several values of the gauge parameter ξ, the unitary gauge (UG), standard BW approximation with constant (BW) or s-dependent widths (SDW) given by (4.4).
employing the standard formula [63]: where the dimensionless parameters is defined in (4.3) and K n (x) denotes the modified Bessel function of order n. In fig. 4, we present the values of σv rel for g x = 0.3, 1 and √ 2π 2 . The results were calculated in the PT, the Feynman gauge (ξ = 1), the unitary gauge (UG) and in the BW approximation with the usual constant decay widths and sdependent widths (SDW) as stated in (4.4). We see that as the dark gauge coupling g x decreases, the standard BW approximation starts to describe more accurately the resonant cross-section. This follows from the fact that XX contribution to the self-energy Π ij (s), which strongly depends on the energy close to the threshold, ceases to dominate over the nearly constant contributions from the SM states. Then, as presented in the lower panel of fig. 4 for g x = 0.3, the thermally averaged cross-section σv rel forms a hilltop located at x = 1/|δ|. This is the value of the temperature for which the thermal distribution (4.5) is centered at the resonance peak s = m 2 2 . In fig. 5, we display numerical comparisons of thermally averaged cross-sections computed by means of the PT and the standard BW approximation on the plane defined by the parameters g x and δ. We observe that the BW approximation is applicable only if the parameters g x and |δ| are sufficiently small. Otherwise, its simplistic use for relative velocities v = 100 km/s (typical in astrophysical searches of DM annihilation signal) can lead to annihilation cross-sections that are much smaller (larger) than the PT result, for negative (positive) values of δ. On the other hand, for thermally averaged cross-sections calculated at temperatures close to that of chemical decoupling (x = 20), the BW approximation gives results that are comparable to or larger than those derived by the PT resummation method.
In the following, we will analyze the evolution of the relic DM density. Following the recent studies [14,48], we take into consideration the effect of early kinetic decoupling, which turns out to be an important phenomenon for scenarios with resonant DM annihilation. This effect appears, because the resonantly enhanced annihilation rate corresponds to scattering processes of suppressed strength, whose role is to maintain the kinetic equilibrium between the DM and the thermal bath in the early Universe. Moreover, as readily observed from fig. 4, the thermally averaged cross-section σv rel grows by many orders of magnitude when x varies from the chemical decoupling temperature of x ≈ 20 until x ∼ O(10 4 ). Therefore, the effective DM annihilation can be prolonged to a period when the DM does no longer have the same temperature as that of the SM thermal bath.
In order to determine the evolution of DM density, we assume that the DM distribution does not deviate significantly from the thermal one. Therefore, we may describe it with the DM temperature T X . To compute the DM relic abundance, we solve numerically a set of 2 The value gx = √ 2π corresponds to the saturation of the perturbativity bound for the quartic coupling λH,S ≤ 4π for the resonant region of the VDM model, for which 2MX ≈ m2.
In the above, m Pl is the Planck mass, n X denotes the number density of DM, s is the entropy density, H is the Hubble parameter, g X = 3, whereas g and h are respectively the effective numbers of relativistic degrees of freedom for energy density and entropy. Besides the collision terms ∝ σv rel , the coupled equations (4.6) and (4.7) contain also the following thermal averages: where E and p are the energy and momentum of the DM and z ≡ (s − 1)( 2 + − 1). Finally, the scattering momentum exchange rate γ(T ) may be expressed as where m SM is the mass of the SM state upon which the DM scatters, g ± (ω) = 1/[exp(ω/T )± 1] denotes its phase-space distribution and σ T = dΩ (1 − cos θ dσ/dΩ) is the standard transfer cross-section for elastic scattering. The sum runs over all relativistic degrees of freedom present in the SM thermal bath. The different predictions for the thermally averaged cross-sections σv rel calculated using various methods manifest themselves in different evolutions of the DM yield Y (x). As can be seen from fig. 6, the result obtained in the unitary gauge or in the Landau gauge ξ = 0 deviates significantly from the PT evaluation. Likewise, the standard BW approximation completely fails to describe the underlying dynamics if g x = √ 2π. Instead, an evaluation in the usual Feynman gauge ξ = 1 (where no unphysical thresholds occur) or in an s-dependent width approximation for the BW propagator of the resonant mediator turns out to be both very good approximations of the PT resummation approach.

Conclusions
In processes of resonant DM annihilation, one central difficulty arises from the use of a BW approximation with a constant width for the propagator of the exchanged particle in the s-channel. In particular, if the respective DM channel contributes significantly to the decay width of this mediator particle, the evaluation of the resonant transition amplitude can become very inaccurate close to the DM production threshold. This problem may be partially circumvented by utilising a running s-dependent width for the mediator [14], which results from the imaginary (absorptive) part of the mediator's self-energy. However, in spontaneously broken gauge theories, such a running width becomes gauge-dependent in the off-shell region.
In this paper, we addressed this issue within the gauge-independent framework of the Pinch Technique (PT). The PT has many desirable field-theoretic properties, including analyticity, unitarity and the gauge invariance of the classical action. As an illustrative scenario, we considered the SM extended by a local U(1) X symmetry which may lead to a massive stable gauge field X that could successfully play the role of a Vector DM. By adopting a resummation approach implemented by the PT, we calculated the PT resummed propagator for the scalar h i h j system (with i, j = 1, 2). In addition, we calculated the relevant one-loop corrections to the h i XX-coupling in the PT. As we have shown in this work, the latter are necessary in order to obtain a transition amplitude for DM-pair annihilation which has the proper high-energy limit in agreement with both the Equivalence Theorem and the Generalized Equivalence Theorem.
The PT resummation method enabled us to derive a Born-improved amplitude which was used to obtain self-consistent and accurate predictions for the DM abundance of the X vector boson. We have illustrated how these predictions differ from those that one finds using other methods of resummation. In our comparative analysis, we have taken into account the effect of early kinetic decoupling of the dark-sector particles from the SM thermal bath. We have shown that the standard BW approximation is not applicable to large part of the parameter space in which the DM contribution to the mediator's selfenergy dominates over the SM one. Its mere use results not only in distorted predictions for the DM relic abundance, but also in wrong annihilation cross-sections for indirect detection experiments. To deal with this problem, we devised an alternative simple approximation that utilises an energy-dependent width which yields the correct DM annihilation rates and relic abundance. In this study, we presented results for resonant DM annihilation processes only. Nevertheless, it is straightforward to extend our considerations to DM elastic scatterings. We postpone such an analysis to a future communication.
In this VDM model, the VEVs of the scalar fields are given by ( H , S ) = 1 √ 2 (v, v x ) which generate the gauge-boson masses where g, g and g x are the SU(2) L , U(1) Y and U(1) X gauge coupling constants, respectively. The scalar potential of the VDM model reads with H = H + , H 0 T . To determine the physical scalar spectrum, we expand the scalar fields linearly about their respective VEVs as follows: Our next step is to diagonalise the square mass matrix M 2 for the fluctuations (φ H , φ S ). This can be done by carrying out an orthogonal O(2) rotation R which is usually parameterised by a mixing angle α, such that M 2 diag = R −1 M 2 R. In this way, we obtain where h 1 is the observed Higgs particle with a mass m 1 = 125 GeV. All scalar vertices of the theory may be described in terms of the four extra parameters: M X , m 2 , g x and α, in addition to the SM ones. Specifically, the quartic couplings of the VDM potential may be expressed as follows: with c α ≡ cos α and s α ≡ sin α. p 2

D Generalized Equivalence Theorem for Resummed Amplitudes
The proof of the Generalized Equivalence Theorem (GET) for the amplitude A XX→f f goes along the lines of the one given in [38] for the SM processf f → ZZ. It is based on the classical WIs that are satisfied by the PT Higgs self-energies and the PT hZZ vertex.
In order to prove the GET for the process XX →f f , we have to check that the following identity holds true (summation over the repeated indices i, j, k implied):

(D.2)
Writing µ (p 1,2 ) = p µ 1,2 /M X + x µ (p 1,2 ) and using (3.18) and (3.19), one arrives from (D.2) at an equivalent expression This last expression can be shown to be valid for both the tree-level and its PT extension.
To do so, we first observe that (s − m 2 k )δ ki + Π ki = ∆ −1 ki . Then, using (3.13) and (3.20), one finds that the main bulk of the terms in (D.3) cancel. Thus, one is only left to prove that The latter is true due to the orthogonality of the mixing matrix R. This concludes our proof of the GET for the amplitude A XX→f f .