Finite-temperature modification of heavy particle decay and dark matter annihilation

We apply the operator product expansion (OPE) technique to the decay and annihilation of heavy particles in a thermal medium with temperature below the heavy particle mass, mχ. This allows us to explain two interesting observations made before: a) that the leading thermal correction to the decay width of a charged particle is the same multiplicative factor of the zero-temperature width for a two-body decay and muon decay, and b) that the leading thermal correction to fermionic dark matter annihilation arises only at order T4/mχ4. The OPE further considerably simplifies the computation and factorizes it into model-independent matrix elements in the thermal background, and short-distance coefficients to be computed in zero-temperature field theory.


Introduction
In the early Universe all interactions between particles take place in a thermal plasma, leading to modifications of the cross sections and corresponding rates with respect to their zero-temperature values. Many interesting phenomena occur in the regime where the temperature T is much smaller than the energy scale M of the hard process. Examples are the freeze-out of weakly interacting massive dark matter particles and the decay of heavy particles during big-bang nucleosynthesis. In these cases the finite-temperature modification is expected to be of order g 2 × (T /M ) a , where g represents the interaction strength and a is some number.
Modifications to particle decay widths due to a thermal bath of photons have been studied in different cases and remarkably simple results were found. In particular, the leading O(T 2 ) correction vanishes for the decay of the neutral Higgs boson [1] and heavy Majorana neutrino [2], while for a charged particle two-body decay it is simply proportional to the tree level width [3]. Interestingly, the same result was found for the three-body decay of the muon. The more involved case of annihilation of two neutral fermions χ into two charged fermions via t-channel exchange of a massive scalar has been computed in [4]. The primary purpose of that work was to demonstrate the infrared divergence cancellation in relic density computations, but the leading temperature-dependent correction was also considered and found to appear only at O(T 4 ).
The simplicity of these results calls for an explanation, which we provide here in the framework of the operator product expansion (OPE). In the case of particle decay this treatment is inspired by, and in many ways very similar to, the OPE of inclusive heavy quark decay in QCD [5]. The thermal background plasma takes the role of the soft QCD vacuum quantum fluctuations. For dark matter annihilation we apply a generalized operator production expansion to the (imaginary part of the) χχ → χχ forward scatteringamplitude. In addition to providing a physical understanding of the observed universality JHEP09(2016)031 and simplicity of the O(T 2 ) corrections, the OPE offers a considerable simplification of the calculations. The approach presented here can be straightforwardly applied to other physically motivated situations such as the co-annihilation of charged states during freezeout. We note that OPE methods have been developed systematically before for thermal field theory in the QCD context for SVZ sum rules at finite temperature [6,7], and for the study of more general spectral functions in the low-energy QCD plasma [8]. Further, in [9] the calculation of the thermal production rate of a right-handed neutrino has been represented in an OPE form.
2 Charged particle decay at finite temperature For definiteness we consider the spin-averaged total width of a fermion ψ with electric charge q into another fermion χ of the same electric charge and other neutral particles in an unpolarized thermal bath of photons and SM fermions f at a temperature T . Specific examples are those of [3], ψ → χφ, where φ is a neutral scalar, and the more realistic muon decay µ → eν µνe . We assume that the temperature of the bath is small compared to the ψ mass, m ψ ≫ T , but can be of the same order or even much larger than mass of the other charged particle. For simplicity, we further assume (as in [3]) that the decay occurs at rest with respect to the thermal bath. The decay rate will be modified compared to the zero-temperature value by interactions with the plasma.
The decay of ψ is caused by some weak interaction where J µ is the fermion current and O µ represents the neutral fields. For the toy situation ψ → χφ, we follow [3] and adopt J µ =χP L ψ (P L = 1−γ 5 2 ) and O µ = φ, while for muon By the optical theorem the decay width can be expressed as to lowest order in the weak coupling λ, but to all orders in the electromagnetic interaction.
Here Im {T µν } refers to the imaginary part of and L µν to the neutral particles, which are unaffected by the plasma, integrated over phase space. |ψ; T denotes the ψ one-particle state with momentum p = m ψ;T v and non-relativistic normalization in the thermal bath. The scale hierarchy T ≪ m ψ allows us to separate the hard decay process from the effects of the thermal bath by performing the OPE of the correlation function (2.3). The short-distance physics is encoded in the Wilson coefficients, which can be computed at zero temperature, while the thermal modifications all reside in the matrix elements of local operators computed in the thermal bath. The situation is analogous to the calculation of the (zero-temperature) semi-leptonic decay width [10][11][12] of a heavy b-hadron H b . The JHEP09(2016)031 ψ particle plays the role of the b quark, while the thermal bath substitutes the hadronic medium of soft light quarks and gluons in H b . The finite-temperature ψ mass m ψ;T is the analogue of the B-hadron mass, the zero-temperature mass corresponds to the quark mass.
For the construction of the OPE we can decompose p µ = m ψ v µ + k µ , where v is the four-velocity of the plasma and m ψ the zero-temperature mass. Since we assume that the decaying particle is at rest with respect to the plasma, and T ≪ m ψ , k is a soft momentum with scaling k ∼ T . The relevant OPE of the time-ordered product in (2.3) is where the Wilson coefficients C µν 0 and C µν 2 are specific to the particular decay process. However, they depend only on m ψ v, and are to be computed by matching at T = 0. Hence only the zero-temperature mass m ψ enters. In general, the background plasma breaks Lorentz invariance and it might be necessary to keep additional operators, which are not scalars, and have non-vanishing matrix elements in the thermal bath. However, since we assume that the particle decays at rest with respect to the plasma, the vector v coincides with the one already introduced by the particle state itself, and (2.4) for muon decay is the same as appears in [13]. The neglected term can contribute at most at O(T 3 /m 3 ψ ), which is smaller than the putative leading thermal correction.
As the second (magnetic) operator does not contribute in an unpolarized medium, it remains to evaluate the matrix element ψ; T |ψψ |ψ; T of the leading operator in the thermal background. Since this is independent of the short-distance decay, it already follows here that the thermal correction must be a universal modification of the zero-temperature decay width. The matrix elements of the OPE in heavy particle states have themselves a non-trivial 1/m ψ expansion. Using the equation of motion, we can write [11] ψψ =ψ/ vψ The usefulness of this equation stems from the fact that the first term is related to a conserved current, and hence its matrix element is known exactly, from which it follows that Therefore the leading thermal correction is a direct effect of the matrix element of the kinetic energy operator O k ≡ −ψ (iD ⊥ ) 2 2m 2 ψ ψ evaluated in the thermal background. We may interpret K ψ = ψ; T |O k |ψ; T as the average kinetic energy of the ψ particle acquired by the interactions with the photons in the plasma. By dimensional analysis The diagrams defining the one-loop contributions K 1 , K 2 and K 3 , respectively, to the kinetic operator matrix element. The filled square represents the operator insertion and diagrams 1 and 2 involve the usual QED vertices aside from the operator vertex.
This derivation explains in a rather straightforward manner three observations originally made in [3]: 1) that soft and collinear divergences cancel in the sum of virtual corrections and emission and absorption processes, 2) the leading finite-temperature correction is O(T 2 /m 2 ψ ), and 3) a universal factor multiplying the tree-level decay width. 1 In contrast to the QCD analogue of the semi-leptonic decay of H b , where the soft physics is non-perturbative, the matrix element K ψ of the kinetic operator in the thermal plasma can easily be computed perturbatively. The one-loop diagrams are depicted in figure 1, in terms of which the matrix element is given by the sum Since we are interested in the temperature-dependent correction, we retain only the thermal part of the equilibrium photon propagator where f B (k 0 ) = (e |k 0 |/T − 1) −1 is the Bose-Einstein distribution of the photons in the rest frame of the plasma. Writing down the diagrams explicitly, we find the spin-averaged matrix elements (2.12) Note that the particle mass in the thermal background, m ψ,T , rather than the Lagrangian mass m ψ should be used in the evaluation of the low-energy matrix elements, but since the difference of the mass-squares is O(T 2 ), it contributes to a higher-order correction, and both can be identified. The above expressions can be simplified by using p = m ψ v with v 2 = 1, and v · k ⊥ = 0. The straightforward evaluation then results in where we defined τ = T /m ψ . This has to be multiplied by q 2 for a particle with electric charge q in units of the positron charge. Inserting this into (2.8) gives, explicitly, in complete agreement with [3].
Comparison with the derivation of this result in [3] highlights the power of the OPE approach. It also provides a physical interpretation of the correction as a time dilatation effect due to the average kinetic energy of the particles due to collisions with the photons of the plasma.
The finite-temperature modification of the decay width of a Majorana fermion [2] mentioned in the introduction was obtained through effective field theory methods, which are also based on systematic scale separation. The OPE nevertheless provides a more direct approach to the inclusive decay width in the same way as the full development of heavy quark effective theory is not required to compute the 1/M expansion of the inclusive or semi-leptonic b-hadron decay width in QCD.

Dark matter annihilation
Effective field theory and the OPE can also be applied to two-particle annihilation in the thermal medium, provided the temperature is small compared to the annihilating particles' mass. This considerably simplifies the diagrammatic analysis of [4] and provides an understanding of the temperature scaling of the leading thermal correction. Although the method is general, we discuss below the annihilation of a heavy, electrically neutral Dirac fermion into a pair of light charged fermions (m f ≪ T ), and refer to the heavy fermion as the "dark matter" particle to establish contact with [4].
We assume at first that the annihilation process χχ → ff occurs through the local four-fermion operator where 1/Λ 2 is an unspecified coefficient of mass dimension −2, and Γ A , Γ ′ A are Dirac matrices, which may be multiplied by up to one covariant derivative.
The total spin-averaged annihilation cross section in the thermal background follows from the optical theorem, where the state |χχ; T represents the annihilating pair in the thermal photon background, p ≡ p 1 + p 2 is the total incoming momentum and s ≡ p 2 the center-of-mass energy squared. Once again, we assume that the center-of-mass frame of the annihilation is at rest with respect to the plasma, in which case p = √ s v with v defining the plasma frame.
Since the annihilating particles do not couple to the thermal bath, the χ field part of (3.2) is readily done by contracting the fields with the external state, which results in JHEP09(2016)031 some tensor L AB . The non-trivial part is the time-ordered product of the fermion current J A =f Γ ′ A f . Since the final state particles are very energetic relative to the soft degrees of freedom of the plasma, we perform the OPE of which the matrix element within the thermal vacuum |Ω T needs to be taken. Up to dimension 4, the operators O i are constructed from contractions with the metric tensor and the plasma velocity v of with Γ a general Dirac matrix in spinor space. Apart from the unit operator there is no operator of dimension lower than 4. This allows us to deduce immediately that the leading order thermal correction is at least of the order O(T 4 ) or O(m 2 f T 2 ). The thermal matrix elements are easily computed, see appendix A. The photon "condensate" is given by The light fermion operator m ff Γf is generated only with Γ = 1 1 due to parity invariance and helicity conservation and is power-suppressed for m f ≪ T (and exponentially suppressed for m f ≫ T , but we do not consider this limit): The addition of a covariant derivative alleviates the m f suppression. In the limit m f → 0 we obtain The short-distance coefficients from contracting the χ fields and the OPE depend on the momenta p 1 and p 2 of the two annihilating particles. Since p 1 + p 2 = p = √ sv, for n ≡ i + j = 2, 3, 4 (i, j = 1, 2) form a complete set of scalar operators. Using (3.5), (3.7) their matrix elements in the state |Ω T are given by with e χ = E χ /m χ = √ s/(2m χ ) the rescaled energy of the two annihilating particles in their center-of-mass frame. The leading 2 finite-temperature correction of order T 4 /m 4 χ to the annihilation cross section is obtained by multiplying the scalar operator matrix elements with their Wilson coefficients (including the χ-field part of the operator (3.1)), and hence takes the form (C f 1 ≡ 0) The short-distance coefficients can be obtained by a standard zero-temperature calculation but depend on the specific model. To validate the OPE approach, we reconsider the model of [4], where the dark matter particle annihilates into a pair of fermions through the t-channel exchange of a scalar. In contrast to [4], we now assume that χ is a Dirac fermion, which avoids the P-wave suppression of the zero-temperature annihilation cross section to a pair of fermions. More precisely, we consider the process χχ → f + f − in an extension of the Standard Model by an SU(2) × U(1) singlet Dirac fermion χ and a scalar doublet φ = (φ + , φ 0 ) T . The relevant terms in the Lagrangian read where the SM fermions form a left-handed doublet f = (f 0 , f − ) T . The t-channel exchange of the heavy scalar φ occurs through a chiral Yukawa-type interaction.
We further assume that m φ ≫ m χ and integrate out the scalar mediator field in a first step. This leads to local annihilation operators O ann = J µ t J µ of the general form of (3.1), where the "current" operators are given by at tree level. Since the χ field is electrically neutral we have already simplified the expression by eliminating derivatives on χ (χ) in favour of the momentum p 1 (p 2 ) of the particle, which literally arises only after taking the matrix element.
In the second step we construct the OPE (3.3) of the above currents. We shall provide results explicitly only for the leading term in the expansion in 1/m 2 φ in the first step, which allows us to simplify J t → λ 2 2m 2 φ J 0 . In this case, the calculation of the Wilson coefficients of 2 We recall that we assume m f ≪ T , hence the neglected O(m 2 f T 2 /m 4 χ ) correction is sub-leading.

JHEP09(2016)031
the photon and fermion bilinear operators is very similar to the corresponding calculation in QCD, except that one cannot assume Lorentz invariance of the vacuum state, in which the operators are eventually evaluated. The computation and the relevant diagrams are sketched in appendix B. Defining C i = 2 Im (C i ) for the Wilson coefficients appearing in (3.13), we find for the coefficients of the photon operators (using the notation of [4], ξ = m φ /m χ , ǫ = m f /2m χ and τ = T /m χ ): 20) where the proportionality factor reads αλ 4 /(48m 6 χ e 5 χ ξ 4 e 2 χ − 4ǫ 2 5/2 ). We have kept here the full dependence on the fermion mass m f , since for the contribution from the thermal photons our calculation is accurate also for m f = O(T ), and this provides a further check as discussed below. However, within the assumption m f ≪ T ≪ m χ adopted for the thermal fermion contribution, these terms can be dropped, in which case the coefficients simplify to (3.23) In case of the Wilson coefficients from the fermion operators, we set m f = 0 from the start (since we also neglected the condensates proportional to m 2 f T 2 from the beginning), and obtain Assembling the coefficient functions and matrix elements according to (3.13) we obtain the following two results: whose magnitude is 3.5 times larger than the thermal photon correction.

JHEP09(2016)031
As a further check of the procedure we computed the tree-level, zero-temperature cross section, expanded in the heavy mediator mass, which follows from the coefficient function of the unit operator in the OPE, as well as the next ξ −6 term in the expansion by keeping the 1/m 2 φ suppressed currents in the expression (3.16) for J t . In [4] the corresponding and further results were obtained in the same model except for the Majorana nature of the dark matter particle, which involves further diagrams and leads to a suppressed zero-temperature cross section in the S-wave limit e χ = 1. We have computed the Dirac case in the diagrammatic approach adopted in [4] and confirmed the results for the zero-temperature cross section and thermal corrections given above.
As was the case for charged particle decay, the OPE computation results in a significant conceptual and technical simplification. The infrared and collinear divergences present in individual diagrams in the diagrammatic approach are absent from the beginning. The OPE also provides direct insight into the temperature dependence of the thermal correction, and explains the absence of a O(T 2 ) correction for the neutral fermion dark matter annihilation observed in [4] by the non-existence of a gauge-invariant operator of dimension 2. In general it provides a transparent interpretation of the thermal correction in terms of model-independent condensates in the thermal plasma and factorizes the model dependence in short-distance coefficients, which can be computed in conventional zero-temperature field theory.

Conclusions
In this paper we applied the operator product expansion technique to the decay and annihilation of heavy particles in a thermal medium with temperature below the heavy particle mass, m χ . This allows us to explain two interesting observations made before: a) that the leading thermal correction to the decay width of a charged particle is the same multiplicative factor of the zero-temperature width for a two-body decay and muon decay [3], and b) that the leading thermal correction to fermionic dark matter annihilation arises only at order T 4 /m 4 χ [4]. The OPE further considerably simplifies the computation and factorizes it into model-independent matrix elements in the thermal background, and short-distance coefficients to be computed in zero-temperature field theory.
While we considered here several specific examples, the method itself is general. An interesting extension, which combines the features of charged particle decay and neutral particle annihilation, concerns the annihilation of charged particles. In the context of dark matter relic density calculations this situation may arise if co-annihilations of dark matter particles with charged states are important. In this case it is expected from the above that the leading thermal correction appears already at O(T 2 /m 2 χ ) as in charged particle decay. However, since the correction is a combination of a typically weak interaction with the plasma and of T /m χ suppression factors, it is generally still rather small.

JHEP09(2016)031
with f F (ω) the Fermi-Dirac distribution. Decomposing the general Dirac matrix Γ into the basis {1 1, γ 5 , γ µ , γ µ γ 5 , σ µν } we find that there are only two non-vanishing matrix elements of dimension four operators. There first involves Γ = 1 1, in which case it follows immediately that Ω T | m ff f |Ω T = O(m 2 f T 2 ) is suppressed for m f ≪ T . On the other hand, where the relative coefficient of the two Lorentz tensors follows from the equation of motion i / Df = 0 in the m f → 0 limit. The coefficient d 1 is found by multiplying with v µ v α , which yields

B Calculation of the Wilson coefficients
We briefly sketch the calculation of the coefficient functions of the operators (3.8) in the operator production expansion (3.3). As discussed in the main text, at O(ξ −4 ) we only need to consider the OPE of the product of two currents J µ 0 =f P R γ µ f . Referring to (3.2) and (3.13), the matching coefficients of the photon operators can be obtained by taking the one-photon matrix element. Denoting C i = 2 Im (C i ), the matching relation reads Here tr χ µν = 1 4 spin χχ| (χP L γ µ χ)(χγ ν P R χ) |χχ = 1 4 Tr ( / p 2 −m χ )P L γ µ ( / p 1 +m χ )γ ν P R (B.2) denotes the part that comes from contracting the fields in the dark matter current J µ = χP L γ µ χ with the external state, while O µν (p) refers to the Fourier-transformed operator product on the left-hand side of (3.  Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.