Bound-state formation, dissociation and decays of darkonium with potential non-relativistic Yukawa theory for scalar and pseudoscalar mediators

Dark matter models with light mediators featuring sizable interactions among dark particles enjoy an increasing attention in the model building community due to the elegance with which they can potentially explain the scaling relations governing galactic halos and clusters of galaxies. In the present work we continue our study of such models using non-relativistic and potential non-relativistic effective field theories (NREFTs and pNREFTs) and explore the properties of a Yukawa-type model with scalar and pseudoscalar interactions between a low-energetic scalar mediator and heavy dark matter fermions. In particular, we make first steps towards the formulation of such theories at finite temperature by providing the thermal bound-state formation rate and the thermal break-up of bound states from the self-energies of the dark-pair fields, that interact with the thermal environment. We estimate numerically bound-state effects on the dark matter energy density, that provide up to a 35% correction depending on the relative size of the model couplings.


Introduction and motivation
The quest for uncovering the true nature of dark matter (DM) comes as a compelling and challenging endeavor for both particle physics and cosmology communities.On the one hand, gravitational effects provide a strong evidence of the existence of DM, while possible alternative explanations seem inconclusive.On the other hand, we still lack the much anticipated discovery of signals of dark particles at colliders or experimental facilities dedicated to direct as well as indirect DM searches.
The non-observation of any DM signal has resulted in cornering many minimal DM models.In particular, models comprising weakly interacting massive particles (WIMPs) are under strong tension with current experimental observations [1].This naturally prompts us to consider the existence of richer dark sectors, containing multiple dark particles, cf.e.g.refs.[2][3][4].The dark particles can experience their own hidden forces, which are partially or entirely secluded from the Standard Model (SM) sector.Furthermore, the existence of light force mediators with masses much smaller than that of the actual DM particles may affect the DM dynamics in many ways.A straightforward consequence is that DM may experience sizable self-interactions that depend on the relative velocity among dark particles and can provide a dynamical explanation for the scaling relations governing galactic halos and clusters of galaxies [5][6][7][8][9][10][11][12][13].In this sense, from a phenomenological point of view a next-to-minimal dark sector is more than welcome.
Another remarkable property of DM models with light mediators is the possibility to have bound states within the dark sector [14][15][16].Depending on the details of the model, both stable and unstable bounds state may form.Under the conditions, that the dark particles have reached a thermal equilibrium in the early universe and the DM relic density is generated through the well-known freeze-out mechanism, bound-state formation and dissociation can become important for an accurate determination of the present-day energy density.Indeed, whenever bound states are formed, and not effectively dissociated or melted away in the thermal plasma, they provide an additional process for the depletion of DM particles in the early universe. 1 quantitative estimation of bound-state effects on the DM relic density is model dependent.Most of the past and recent literature is focused on vector mediators, that resulted in complementary and diversified approaches to calculate formation cross sections at zero [14,16,[20][21][22][23][24] and at finite temperature [25][26][27][28][29][30][31][32][33].A more careful look at scalar mediators has been initiated only recently [24,[34][35][36][37][38][39][40][41].Even without the addition of scalars in the dark sector, it has been shown that if the dark particles couple to the SM Higgs boson, and DM is sufficiently heavier than the SM scalar, Sommerfeld enhancement and bound-state formation are relevant [42,43].The in-vacuum bound-state formation has been scrutinized for both neutral and charged scalar mediators and, in the latter case, it was found that fast monopole interactions can dramatically change the relic density estimation even for small values of the couplings [39][40][41].
As for the treatment of the bound-state dynamics, one of the main approaches adopted so far relies on the Bethe-Salpeter [44] equation for the wave functions, and the reduction thereof to a Schrödinger equation using the instantaneous approximation in the nonrelativistic regime [35].The building blocks of this framework are relativistic four-point Green functions of single DM particles at zero temperature, whereas [29] features a finitetemperature version of this formalism.A complementary approach, that we adopt in this work, is to exploit the technology of non-relativistic effective field theories (NREFTs).Our main goal is to provide quantum field theoretic tools for computing thermal processes involving bound states, or in general heavy pairs interacting through a potential.It is easy to see that the problem at hand comes as a multi-scale system.One finds three typical and presumably well-separated scales relevant for the non-relativistic dynamics, which are M M v M v 2 , where M is the DM particle mass, while v stands for its typical velocity in a bound state.A Coulombic bound state satisfies v ∼ α, with α being the relevant coupling constant.Furthermore, our system comprises the mediator mass m, which is assumed to be much smaller than M , and the thermal scales, most notably the temperature of the early universe.Whenever the scalar mediator can interact with other light degrees of freedom in the dark sector, or with itself via a self-coupling, a thermal mass m T can be generated.In this case it must be added to the hierarchy of scales mentioned beforehand.
Bound states in general exhibit a strong sensitivity to thermal scales, and it is rather useful to make an analogy between a well-studied system that arises in QCD, namely heavy quarkonium in a quark-gluon plasma, and heavy DM pairs in the early universe.In the following we choose to call a bound state made of two DM particles a darkonium.Two processes are at play for the dissociation of quark-antiquark bound states in the medium, which are the gluo-dissociation [45] and the dissociation by inelastic parton scattering [46].The former describes a situation, where a thermal gluon hits the quark-antiquark pair in a singlet configuration and, if sufficient energy is available, breaks it into an unbound color-octet state.In this work we will address the analogue process induced by a thermal scalar being responsible for the breaking of a darkonium into a scattering state.The second dissociation process comes as a 2 → 2 scattering reaction, where the thermal particles in the bath exchange energy with the heavy quarkonium through a gluon exchange, turning a bound-state into an unbound one.A distinctive feature of an EFT approach is the availability of rigorous power counting rules that can be used to show that the gluodissociation process is dominant over the inelastic parton scattering for E > m D [47][48][49], where m D ∼ g s T is a thermal mass for the gluon called Debye mass and E ∼ M v 2 defines the ultra-soft scale.
It is now clear that in order to follow the intricate dynamics of bound states in a thermal medium, an EFT approach constitutes a powerful tool that (i) allows scrutinizing the various arrangements of the energy scales and (ii) can be employed to devise the proper field theory to carry out the derivation of cross sections and decay widths by using the relevant degrees of freedom at the scale of interest.We shall attack this problem by using a recently proposed EFT for non-relativistic fermion-antifermion pairs interacting via a scalar mediator.Building upon the well-known NREFTs and potential NREFTs (pNREFTs) of QED and QCD [50][51][52][53][54] as well as scalar Yukawa theory [55,56], we developed potential non-relativistic Yukawa theory (pNRY) [57] that is well suited to address the questions at hand.The relevant degrees of freedom of pNRY are heavy fermion-antifermion pairs and ultra-soft scalars.In the first paper on the subject, we have applied pNRY to compute the darkonium spectrum and the bound-state formation cross section at zero temperature.Here we make a first step towards the generalization of pNRY at finite temperature.We shall use the so-obtained EFT to derive the thermal cross section for bound-state formation and the reverse process, namely the dissociation rate via the absorption of a thermal scalar from the medium.The latter process is a genuine finite temperature effect, and we will show how it can be recast as a convolution of the thermal distribution of the scalar particles and an invacuum dissociation cross section, without the need of Boltzmann-like prescriptions.Our pNREFT can be also used to treat pair annihilations, which are known to be key ingredients for DM freeze-out calculations.Here, the Sommerfeld-enhanced cross section and boundstate decay width will emerge naturally.Our program closely follows the successful strategy outlined in the case of heavy quarkonium at finite temperature [48,58,59].It is worth mentioning that an EFT approach to search for DM in atomic-spectroscopy measurements has been recently developed for different types of force mediators, including scalar and pseudoscalar ones, in ref. [60].
In this work we focus on a particular hierarchy of scales, that reads and enables us to use pNRY as a starting point, as we explain in section 2. This hierarchy of scales is motivated by two main aspects.First, temperatures of order of or smaller than the binding energy are interesting for DM phenomenology, since a depletion of DM pairs in the form of bound states is more efficient in such a temperature window.Indeed, due to a suppression of the dissociation process, bound-state formation becomes more relevant and opens up an additional channel for the pair annihilation.Second, as guided by the potential non-relativistic QCD (pNRQCD) at finite temperature and the inherent power counting, it is expected that the bound-state formation (dissociation) via scalar emission (absorption) dominates over 2 → 2 scattering processes with plasma constituents.This latter effect will be addressed in a future study on the subject, as it requires the derivation of thermal self-energies for the scalar field at NLO. 2 A similar investigation for the same hierarchy of scales (1.1), however in the case of an Abelian DM model with a vector mediator, is in preparation [62].
The structure of the paper is as follows.In section 2 we introduce pNRY γ 5 , a variety of pNRY featuring scalar and pseudoscalar Yukawa interactions, and discuss the thermal propagators of the heavy pair and the scalar mediator as the main ingredients for the finite temperature treatment.Then, in section 3, the bound-state formation and dissociation is computed in pNRY γ 5 starting from the self-energies of the heavy-pair field.The thermal rates are naturally obtained in terms of quantum mechanical matrix elements, that are analytically simplified in the Coulombic regime for the first nS states.Annihilations of heavy pairs in pNRY γ 5 are discussed in section 4. Next, section 5 is devoted to a phenomenological study of the DM energy density, where the bound-state formation, bound-state decay width and the thermal width from scalar dissociation enter as key ingredients.Conclusions are offered in section 6, while some technical details underlying our results are collected in the appendices.

Non-relativistic EFTs for a scalar mediator
In this section we introduce the DM model featuring a light scalar force mediator between heavier DM particles, in particular fermions.Such model Lagrangian will be our fundamental theory, from which we can construct towers of low-energy theories by integrating out energy scales.Next, we proceed with the discussion of the pNRY Lagrangian, and we introduce the thermal propagators of the heavy-pair field and the scalar particle.

Dark matter model
We assume DM to be a Dirac fermion singlet under the SM gauge group that couples to a scalar particle via Yukawa-type interactions.The Lagrangian density of the model reads [63,64] where X is the DM Dirac field and φ is a real scalar field.The scalar self-coupling is denoted with λ, whereas the scalar and pseudo-scalar couplings are g and g 5 respectively.
The mass of the scalar mediator m is assumed to be much smaller than the DM particle mass M , m M .In our work, we adopt a simplified model realization, where the issues of the fermion mass generation and of the gauge group governing the dark sector are not addressed (cf.e.g.[65,66] for a simplified model with two mediators, scalar and vector, where the gauge invariance and spontaneous symmetry breaking in the dark sector is fully accounted for).Our aim is to consider the Lagrangian given in eq. ( 2.1) as one of the simplest representatives for minimal DM models [3,64] with a light scalar mediating interactions between DM particles.It is worth mentioning that the model can be much more involved, and can be extended to have a richer set of interaction terms [34,38,67].For example, an interaction of the form ρ λ φ 3 can be foreseen, that is responsible for an additional bound state formation process [38,39].
Then, L portal comprises the interactions between the scalar φ and other degrees of freedom that can belong to the SM sector or to the dark sector.One of the most common realizations of such a portal involves interactions with the SM Higgs boson.Portal interactions are welcome in order to introduce a mechanism that allows φ particles to decay and deplete their population.Indeed, the light scalar particles φ are abundant in the early universe and a substantial population survives after the freeze-out of the dark fermion [64,68]. 3A richer portal sector also allows alleviating and often removing the 3 As for the Higgs portal, the interaction terms with the scalar mediator φ read, before electroweak symmetry breaking, as L portal = −aφH † H − bφ 2 H † H; H is the SM Higgs doublet.The smallest portal couplings that ensure a thermalization between the SM and dark sector, and allow for the decay of the scale φ before BBN, are typically much smaller than the electroweak gauge couplings, see e.g.[64,69].We have explicitly verified this for the assumed range of numerical values for our model parameters g, g5 and M .This means that the portal interactions can be safely neglected in our EFT construction and the corresponding matching calculations.
tension with the experiments for the model, if one considers only the interactions of the scalar φ with the Higgs boson [64,70].
However, when dealing with a thermal environment, many interactions between the scalar mediator and other plasma constituents may endanger the assumed hierarchy of scales.In particular, a thermal mass of the form m T ∼ g T is generated, that can become larger than the in-vacuum mass or the typical ultrasoft scale M α 2 .Moreover, in the case of the bound-state formation process, when the scalar mass exceed the difference between an above-threshold scattering state and a bound-state, cf.eq.(3.4), the 1 → 2 radiative formation process becomes kinematically forbidden.For example, if one only considers the scalar self interaction, a thermal scalar mass m T = λ/12T is generated, that translates into the condition T M α 2 12/λ to attain the hierarchy of scales given in eq.(1.1).We defer the scrutiny of the impact of thermal masses, and a thermal width, of the scalar propagator to future works on the subject.Since we included a fermion-pseudoscalar interaction, we denote the following towers of EFTs with NRY γ 5 and pNRY γ 5 in order to make clear that the corresponding low-energy theories differ from the ones obtained in the case of the sole scalar interaction [57].
To ensure the correctness of the results presented in the next sections we make use of computer tools that were already employed in our previous work on the subject [57].Calculations within the fully relativistic DM model given in eq.(2.1) can be automatized using FeynArts [71], FeynRules [72] and FeynCalc [73][74][75].The matching between the full theory and the corresponding non-relativistic EFTs is done using FeynOnium [76], while we employ QGRAF [77] interfaced to FeynCalc via FeynHelpers [78] to generate the corresponding Feynman diagrams.Together with extensive pen and paper cross checks, this approach greatly facilitates the task of avoiding typos and unintentional mistakes.

NRY γ 5
We now want to proceed to the construction of the low-energy theory relying on the hierarchy of scales in eq.(1.1).As highlighted and extensively discussed in the existing studies of heavy quarkonia as well as hydrogen and muonic atoms at finite temperature [48,49,79], one can first integrate out the in-vacuum scales M and M α.Subsequently, all the smaller scales, including the thermal ones, can be set to zero.Let us stress that in practice the matching and the derivation of the low-energy theories is insensitive to the thermal scales.Therefore, we can first integrate out the scale M to obtain NRY γ 5 and then do the same for the scale M α thus arriving at pNRY γ 5 .We refer to ref. [57] for further details on the construction of these two EFTs.
NRY γ 5 is well suited to describe non-relativistic fermions and antifermions interacting with a soft scalar.In particular, fermion-antifermion annihilations into light mediators are accounted for by four-fermion operators.Schematically, the NRY γ 5 Lagrangian reads where ψ (χ) is the Pauli field that annihilates (creates) a heavy fermion, while all the scalar particles have energy and momenta much smaller than M .The inclusion of the pseudoscalar interaction as in eq.(2.1) introduces new operators in the bilinear sector with respect to NRY, whereas the set of four-fermion operators and L scalar remain the same.As for the modification of the bilinear sector, we find that the leading fermion-pseudoscalar interaction is suppressed with respect to the leading fermion-scalar interaction by k/M , where k is the soft or ultra-soft momentum carried by the field φ.The fermion bilinear at order 1/M with the matching coefficients set to their tree-level values reads As for the antifermion fields one finds (2.4) Some details on the derivation of the above fermion-bilinear Lagrangians are given in appendix A. The operators proportional to one power of g 5 are parity violating, and the notation [∇φ] stands for the derivative acting on the scalar field only.The operators proportional to g 2 5 and involving two powers of the scalar fields were absent in NRY.Our result for the bilinear sector with a pseudoscalar interaction agrees with the findings in ref. [80].
Most notably, the presence of the pseudoscalar coupling generates non-vanishing matching coefficients for the leading dimension-6 four-fermion operators.In this paper we are mostly interested in such leading order modification to the four-fermion sector, whereas the correction to the bilinear sector, and the corresponding ones in the pNRY γ 5 , will not be carried out in full details.The two independent dimension-6 operators read [51] (L 4-fermions and the corresponding matching coefficients are found to be The spectroscopy notation is borrowed from NRQED/NRQCD, so that one can classify the annihilations in terms of the total spin S of the pair, the relative angular momentum L and the total angular momentum J, by writing 2S+1 L J .The pseudoscalar interaction modifies the matching coefficients of the velocity suppressed operators as well.These dimension-8 operators are where we refer to ref. [51] for their explicit definitions.The matching coefficients, that generalize our result for the sole scalar interaction in ref. [57], read It is worth to mention that the vanishing of f ( 3 S 1 ) and also f ( 1 P 1 ) can be understood using symmetry arguments.Although the pseudoscalar Yukawa interaction violates parity, it still preserves the charge conjugation symmetry.The conservation of the charge conjugation selects particular combinations of spin and angular momentum of the annihilating fermionantifermion pair into two scalars. 4e conclude this section by calculating the non-relativistic annihilation cross section for the process X X → φφ without any kind of resummation of the soft physics, which implies that we are solely sensitive to the hard energy modes at the scale M .In order to compare to the results in the literature, we average over spin polarizations of the incoming fermion and antifermion.The cross section then reads where we choose the heavy fermion fields to be normalized non-relativistically as is customary in NRQCD [51].The relative velocity in the center-of-mass frame is given by The imaginary part of the non-relativistic amplitude in the numerator of eq.(2.9) can be easily computed using the four-fermion Lagrangians given in eqs.(2.5) and (2.7), and the matching coefficients in eqs.(2.6) and (2.8a)-(2.8c).We obtain (2.10) The result agrees with the literature [34,36] as for the terms proportional to α 2 v 2 rel and α 2  5 v 2 rel and the velocity independent term, whereas the term proportional to the product of the couplings at order v 2 rel is new.To the best of our knowledge, the αα 5 v 2 rel -piece was not included in former derivations.An analogous recasting of DM annihilations, as in the first line of eq.(2.10), in the context of supersymmetric models can be found in ref. [81].

pNRY γ 5
Having obtained NRY γ 5 , the next natural step is to derive the corresponding EFT at the ultra-soft scale in order to calculate processes relevant to this work, namely the bound-state formation cross section, the bound-state dissociation width and the bound-state decays.To this aim, we work with the pNRY Lagrangian [57] where the degrees of freedom are heavy fermion-antifermion pairs and low-energetic scalars.In the following we briefly summarize its form and elucidate on the modifications that occur when the pseudo-scalar interaction is considered.
In order to obtain NRY γ 5 , we integrated out the hard scale M .The next scale one has to integrate out according to the scale hierarchy in eq.(1.1) is M α, or the inverse of the typical distance between X and X in a pair. 5As long as we assume that the mass of the scalar satisfies m M α, we find ourselves in a Coulomb-like regime, where the velocity of the particles in a bound state scales as v ∼ α.As explained in section 2.2, we can set the thermal scales to zero in the matching between NRY γ 5 and pNRY γ 5 .In practice, one projects the NRY γ 5 Hamiltonian onto the particle-antiparticle sector via where the square brackets in the second line of (2.12) indicate that the spatial derivatives act on the scalar field only, which has to be understood as multipole expanded in the last line of eq.(2.12) as well.The typical size of the two coordinates is given by r ∼ 1/(M α) and R ∼ 1/(M α 2 ).The multipole expansion for the scalar fields in r R ensures that they carry only ultra-soft energies and momenta.To avoid cluttering the notation we suppress the spin indices of the bilocal fields that are contracted with each other.Each term in the pNRY γ 5 Lagrangian has a well-defined scaling.The time derivative scales as ∂ 0 ∼ M α 2 ∼ T , the inverse relative distance and the corresponding derivative obey r −1 , ∇ r ∼ M α, whereas the scalar field and the center-of-mass derivative satisfy gφ, ∇ R ∼ M α 2 ∼ T .Indeed, the dynamical scales active in the so-obtained EFT are the ultrasoft scale and the temperature.
The potential is understood as a matching coefficient and is organized as an expansion in α(M ) and λ(M ), as well as 1/M , the coupling α(1/r) and the relative distance r.The imaginary part of the potential comprises local terms of the form ϕ † δ 3 (r)ϕ as inherited from the four-fermion operators of NRY γ 5 accounting for fermion-antifermion annihilations [53,54,82,83], that we shall address explicitly in section 4. 6 In the second line of eq. ( 2.12), we see the appearance of a monopole and a quadrupole interactions as well as interactions involving the derivative in the relative distance.Such structures arise from the multipole expansion.In contrast to pNRQED and pNRQCD, the dipole term is absent.The last line of eq.(2.12) accounts for the scalar sector, where the scalar field φ should be equally understood as being multipole expanded.The matching coefficients d 1 , ..., d 5 are inherited from NRY and they read, at tree level, as follows [57] (2.13) Let us discuss the modifications of pNRY in the presence of the pseudo-scalar coupling.There are two main aspects to be addressed.The first one concerns the modification of the potential V (p, r, σ 1 , σ 2 ), or simply V in the following, from the inclusion of the pseudoscalar coupling.We have where V (p, r, σ 1 , σ 2 ) g denotes the part of the potential induced by pure scalar interactions, while V (p, r, σ 1 , σ 2 ) g 5 stems from pseudoscalar and mixed scalar-pseudoscalar contributions.
The potential is obtained by matching four-fermion off-shell Green's functions of NRY γ 5 onto pNRY γ 5 .The insertions of the scalar operator with the matching coefficient d 4 and the contributions from the four-fermion operators are suppressed in the power counting.This is why at the accuracy we are aiming at it is sufficient to consider tree level diagrams with two vertices involving the leading fermion-scalar and antifermion-scalar interactions shown in figure 1.The first diagram in figure 1 is identical to the one found in NRY, meaning that the leading contribution to V (p, r, σ 1 , σ 2 ) is the Coulomb potential V (0) = −α/r.In addition to that, one has to consider diagrams that involve pseudoscalar interactions with g 5 , which are a distinctive feature of NRY γ 5 .Here one finds a mixed scalar-pseudoscalar contribution and a pure pseudoscalar term.Doing the matching in the momentum space and Fourier transforming [84] the results to the position space we obtain where r = r/r with σ 1 (σ 2 ) being the spin matrix of the two-component particle (antiparticle) field in the fermion bilinear.In order to estimate the relative importance of such corrections, we set 1/r ∼ M α.Therefore, the mixed and pure-pseudoscalar terms are of order M α 3 (g 5 /g) and M α 4 (g 2 5 /g 2 ) respectively.As for the pure pseudoscalar potential, our result agrees with earlier derivations [67,85].In the case of g 5 ≈ g, the pseudo-scalar exchange induces an α-and α 2 -suppressed contributions as compared to the leading scalar Coulomb potential.In the following we assume the pseduoscalar coupling g 5 to be smaller than g, which renders the contributions from eq. (2.15) even more suppressed.Hence, we work with the leading Coulombic energy levels and the Bohr radius given by The second aspect to be discussed is the modification of the second line in eq.(2.12), namely the interaction between fermion-antifermion pairs and ultra-soft pseudoscalars.One may compute the corresponding vertices by multipole expanding terms containing g 5 in eq. ( 2.3) and eq.(2.4), where we refer to appendix B for further details.The main outcome is that the pseudoscalar induced ultra-soft transitions are suppressed by powers of α and g 5 /g with respect to the ultra-soft contributions displayed in the second line of eq.(2.12).In summary, upon neglecting the suppressed corrections from the pseudo-scalar coupling and assuming the hierarchy of scales given in eq.(1.1), we find that pNRY γ 5 is formally identical to the in-vacuum case already studied in ref. [57].However, there is a crucial difference: the vertices and propagators now have to be understood in a finite temperature field theory, since the EFT in eq.(2.12) comprises the temperature as a dynamical scale.
We conclude the section by stressing that a different choice of the coupling arrangement, namely g 5 > g, would change the low-energy theory eq.(2.12) in a number of ways.First, one would need to include the pseudoscalar and mixed scalar-pseudoscalar contributions to the potential given in eq.(2.15).In this case the leading order potential V is not Coulombic anymore, and nice analytic results for the quantum mechanical matrix elements should be superseded by a purely numerical evaluation.Second, one has to include the ultra-soft contribution originated by the pseudoscalar vertex as well.Phenomenologically, we expect such a regime to be important when g 5 /g ∼ 1/α, that lifts the suppression we rely on.

Thermal propagators
Let us introduce the finite temperature formalism needed for the calculation of the boundstate formation and dissociation rate in a thermal environment.In particular, we shall provide the Feynman rules for heavy fermion-antifermion pairs and light scalars of the Yukawa model as given in eq.(2.1).We work within the real-time formalism (RTF) of thermal field theory.Real-time expectation values depend on how the contour of the time integration in the partition function is deformed to include real times.The modified contour has two lines stretching along the real-time axis.The main practical consequence of this is that the number of degrees of freedom doubles.Fields on the upper branch enjoy the usual time ordering, whereas for the fields living on the lower branch the time ordering is reversed.The physical degrees of freedom describing initial and final states, are of type 1 (upper branch), while fields of type 2 stem from the lower branch.Propagators are represented by 2 × 2 matrices, and they can mix fields of type 1 with fields of type 2. The vertices, however, do not couple fields of different types.For further details, we refer to the standard textbook treatment [86,87] of the subject.
In order to introduce the finite temperature propagator of a darkonium system, we can build upon its QCD counterpart, the heavy quarkonium [48,58].The two systems share many similarities, since both are made of a heavy particle-antiparticle pair with a mass much larger than the temperature of the thermal bath.The non-relativistic propagator of a fermion-antifermion pair interacting through a potential V (r) reads [48] 7 where the Hamiltonian h at leading order comprises M α 2 terms, namely h (0) = k 2 /M + V (r).As first noted in ref. [48], since the [∆ ϕ (k 0 , k)] 12 component vanishes, the fermionantifermion field of type 2 never appears in amplitudes with final and initial states being physical fields of type 1.Hence, it comes as a great simplification in the non-relativistic theories at finite temperature to discard the type-2 fields when considering physical amplitudes.Notice that even though in RTF the heavy-pair propagator acquires a matrix form, it does not depend on the temperature.
Next, one needs the thermal propagator of the scalar mediator.In this work we assume the in-vacuum and thermal masses of the scalar to be negligible with respect to the binding energy.As discussed earlier, it suffices to assume a very small scalar self-coupling λ in order not to generate a sizable m T .The 2 × 2 free scalar propagator at finite temperature reads where n B (x) = 1/(e x/T − 1) is the Bose-Einstein equilibrium distribution.Due to the decoupling of the type-2 heavy pair fields, in the following we will only need the 11component of the heavy fermion-antifermion propagators.Accordingly, the relevant scalar propagator is given by [∆ φ (k 0 , k)] 11 , whereas the other entries of [∆ φ (k 0 , k)] are irrelevant at the order we are working in this paper.One-loop corrections to the scalar propagator, where all entries of [∆ φ (k 0 , k)] contribute, as well as the connection to additional physical processes between the heavy pairs and the thermal environment are deferred to another study in preparation [88].

Thermal cross section and thermal width
We are now in the position to carry out the calculation of the bound-state formation cross section in pNRY γ 5 at finite temperature, thus extending the result presented in ref. [57] for the in-vacuum case.Moreover, we compute the inverse process, namely bound-state dissociation, which is a genuine thermal effect.Let us stress that as far as the extraction of the bound-state formation and dissociation are concerned, when working at the accuracy we are aiming at in this work, pNRY and pNRY γ 5 are equivalent.This is not true when dealing with heavy-pair annihilations and bound-state decays, as will be shown in section 4.

Bound-state formation via scalar emission
The main advantage in addressing bound-state calculations within pNRY γ 5 is that we can describe the system using appropriate degrees of freedom at the energy scale of interest, which are heavy fermion-antifermion pairs as well as ultrasoft and thermal scalars.The bound-state formation cross section can be extracted from the imaginary part of the selfenergy diagram of the pair, that reads where |p denotes a scattering state.The self-energy diagrams contributing to Σ s are shown in figure 2.
Out of the three vertices in the second line of eq.(2.12),only the quadrupole and the derivative interactions contribute.As has already been observed in the literature [34,89], the monopole contribution is zero due to the vanishing overlap between the scattering and bound state wave functions.This still holds at finite temperature.The momentum region in the loop diagrams are those still dynamical in pNRY γ 5 according to the assumed hierarchy of scales, namely the ultrasoft scale and the temperature scale, which are taken to be of the same order.
As an example, let us explicitly show the contribution to the self-energy as originating from the quadrupole interaction, and elucidate on the main steps towards the thermal cross section.We adopt dimensional regularization (DR) with D = 4 − 2 .The corresponding self-energy, leftmost diagram figure 2, reads where one may notice the appearance of powers of the scalar three momenta in (3.2) as induced by the action of the derivative operator ∇ R on the scalar propagator.Notice that here r is a quantum mechanical operator that generates expectation values when acting on the scattering states.The scalar propagator contains both an in-vacuum and a finite temperature contribution, that are manifestly disentangled in RTF and at the order we perform the calculation.Next, we insert a complete set of bound states, thus ensuring that the internal propagator in the loop diagram describes the propagation of the discrete states of the particle-antiparticle spectrum at leading order, where p = M v rel /2 is the relative momentum of the pair in a scattering state.Finally, one can extract the imaginary part of the in-vacuum contribution with the standard cutting rules, whereas for the term involving the finite temperature contribution to the scalar propagator, which already gives an on-shell thermal distribution of propagating scalar fields, we simply select the imaginary part of the bound-state propagators with the relation 1/(x±iη) = P.V.(1/x)∓iδ(x) (cf.[48] for a similar derivation in pNRQCD), where P.V. stands for the Cauchy principal value.The result for the quadrupole-induced cross section reads This expression contains non-trivial quantum mechanical expectation values that must be explicitly evaluated in order to arrive at a final result.The thermal cross section (3.5) is a generalization of the in-vacuum counterpart derived in ref. [57] and constitutes an original result of the present work.It is easy to see that for T → 0 we readily recover the in-vacuum result.
The contributions from the derivative vertex (middle diagram in figure 2) and the mixed diagrams (rightmost diagram in figure 2), that comprises one quadrupole and one derivative vertices, have to be added as well.These contribute at the same order in the power counting [57].The finite temperature cross section, upon including all the relevant diagrams reads where the in-vacuum cross section is As expected from the quadrupole-induced thermal cross section in eq.(3.5), [1 + n B (∆E p n )] factors out in each term.The second line in eq.(3.7) corresponds to the derivative contribution, whereas the third line accounts for the mixed contribution from the quadrupolederivative interaction.Let us stress that the result given in eq.(3.6) is obtained without making use of any semi-classical construction from Boltzmann equations to describe particle collisions in the plasma.The thermal character of the cross section is rigorously derived from a quantum field theory at finite temperature, here pNRY γ 5 .One may recognize the typical factor that one would introduce according to the Boltzmann prescription to accompany the produced scalar particle φ in the final state.There is no such factor for the heavy darkonium, since in our setting M T .This means that the corresponding Bose distribution is exponentially suppressed and cannot show up in our EFT.Another relevant derivation that exploits a pNREFT can be found in ref. [33], where a carbon copy of QED for the dark sector is employed to derive a thermal cross section.However, the identification of the thermal cross sections appears to be partly based on the network of Boltzmann equations.
The cross section (3.6) can give us a first indication of the importance of thermal effects in the temperature region allowed for the hierarchy of scales of interest.The factorization of the in-vacuum and thermal parts of the cross section allows to rewrite the ratio state, while E 1 = −M α 2 /4 stands for the ground-state energy.In figure 3, we show the contour lines, solid and dashed respectively, for a 50% and 100% increase of the in-vacuum cross section due to thermal effects.One may notice that for sufficiently large value of ζ or sufficiently small relative velocities for a fixed α, the argument of the Bose distribution becomes independent of ζ, and this is reflected in the solid and dashed curves reaching an asymptotic behavior.As we adopt the variable |E 1 |/T for the different bound states, the same increase of the in-vacuum cross sections for the excited states occurs at smaller temperature.This essentially implies that |E n | < |E n | for n > n .A larger increase of the cross section requires larger temperatures, which corresponds to smaller values of The result in terms of pNRY matrix elements allows us to address the calculation of excited states as well.We provide explicit expressions for the ground state and the excited 2S and 3S states (cf.appendix C for details).The corresponding in-vacuum cross sections read 2 6  15 2 5  15 where the Coulombic S-wave Sommerfeld factor is defined as The result for the ground state in eq.(3.8a) has been adopted from ref. [57], whereas the expressions for the excited states are an original contribution of this work.The in-vacuum cross sections are consistent with the unitarity bound on the annihilation cross section in the non-relativistic regime [90], and have been lately revisited and extended to the boundstate formation process [16,35,91].The velocity scaling at ζ 1 obeys σ bsf ∝ 1/v 2 rel , which equally holds for our findings in eqs.(3.8a)-(3.8c).
Figure 4 shows the in-vacuum (solid lines) and thermal (dotted and dashed lines) cross sections for the 1S and 2S states, plotted as functions of ζ for two different values of |E 1 |/T .Some comments are in order.First, the in-vacuum and thermal cross sections are virtually identical for sufficiently small values of ζ, whereas the thermal cross sections get larger than the vacuum ones for ζ > 1.Second, the larger the temperature, i.e. the smaller the values of |E 1 |/T , the sooner σ bsf v rel | T gets enhanced by the finite-temperature nature of the emitted scalar particle.For |E 1 |/T = 3, the in-vacuum and thermal cross section for the ground-state are the same, whereas for ζ > ∼ 10 there is still a factor of 2 or 4 in the case of the excited state 2S or 3S respectively.
In addition to the bound-state formation in nS states, we can easily compute the formation in P -waves.The result for the |21m and |31m states, where we sum over the magnetic component of the angular momentum read 2 9  15 (3.10) We show the corresponding cross sections in figure 5, where they are plotted together with the corresponding bound-state formation cross section for 2S and 3S states.One may notice that for ζ > ∼ 1 the cross section is larger for the nP states.Moreover, we have We close the section by inspecting more closely the thermal effects in the ratios of the bound-state formation cross section σ bsf v rel for the different nS states.In particular, we consider σ nS bsf v rel /σ 1S bsf v rel for n = 2, 3.In figure 6, we show the ratio of the in-vacuum bound state formation cross sections (gray dotted line), as well as the corresponding ratio when thermal cross sections are considered, the latter for some benchmark values of |E 1 |/T .For ζ > ∼ 1 the ratio of the thermal bound-state formation cross section deviates from the invacuum case.As for the comparison between the 1S and 2S states (cf.left panel of figure 6), one may see how the ratio is almost constant in the whole range of ζ for |E 1 |/T = 1, and there is no sudden decrease for large ζ-values (dotted yellow line contrasted with dotted gray).In the case of the 3S state for the same choice |E 1 |/T = 1, the ratio apparently increases with ζ in the region ζ > ∼ 1.We note in passing, that within pNRY γ 5 the relative velocity of the scattering states is typically v rel < ∼ α (hence ζ > ∼ 1).Indeed, in this regime the relevant degrees of freedom are non-relativistic fermions with momenta of order M α and low-energetic scalars.

Bound-state dissociation
In the following we shall take pNRY γ 5 as a starting point for the computation of the dissociation rate of a bound-state, as triggered by the absorption of a scalar from the thermal medium.As it was done for the bound-state formation, we again consider the selfenergy of the heavy pair and exploit the optical theorem.The relevant diagrams are shown in figure 7, where one may see that this time the external wave function field describes a bound state, rather than a scattering state.The optical theorem provides the master formula that we use to extract the thermal width as follows where the subscript "bsd" stands for bound-state dissociation.In order to illustrate the main steps as it was done for the bound-state formation, we again consider the contribution that originates from the leftmost diagram of figure 7, which contains two quadrupole vertices.Formally, the self-energy of the bound state is identical to Σ Q in eq.(3.2).However, the incoming state has an overall negative energy, and one has to insert a complete set of scattering states with positive energy as in This readily explains the difference in the normalization and mass dimension of the bound and scattering state wave functions, that results in obtaining a width for the former and a cross section for the latter from eq. (3.11) and eq.(3.1) respectively.Then, the second and most relevant difference arises when extracting the imaginary part.As the sign of the energy difference has changed, the real process at T = 0 (first term in the scalar propagator [∆ φ (k 0 , k)] 11 ) becomes kinematically forbidden, and only the finite temperature contribution is left.This conforms with the scalar-induced dissociation as a genuine inmedium effect, as in the case of gluo-dissociation for a heavy quarkonium state [48,59].The thermal width from the pure quadrupole contribution reads where ∆E n p = E n − E p = −∆E p n < 0, and in the last equation we have used the energy conservation |k|+E n = p 2 /M to trade the relative momentum of the pair for the momentum of the incoming scalar particle.We notice that the matrix elements of pNRY have to be evaluated accordingly for |p| = M (|k| + E n ) fixed by the momentum conservation.Furthermore, one can see that the scalar needs to have a threshold momentum to trigger the thermal break-up of the bound state.
Our result in eq.(3.13) resembles expressions obtained for the gluo-dissociation of a color-singlet quark-antiquark in the same hierarchy of scales [59], and for the hydrogen atom in QED [49].In the form of the second expression, the thermal width (3.13) can be interpreted as a convolution of the scalar thermal distribution and an in-vacuum function of the scalar momentum.The latter can be taken to be the in-vacuum scalar-dissociation cross section, and we identify it with σ n bsd (|k|) as appearing in the factorized form9 Let us again emphasize that the dissociation rate is obtained from pNRY γ 5 at finite temperature, by calculating the thermal self-energy of the bound-state.The thermal average of the incoming scalar particle naturally arises from a quantum field theoretical derivation, rather than being inferred from a Boltzmann-like treatment.Indeed, we have never introduced rate equations to formulate thermal rates in the first place.In the context of DM, the dissociation rate from a vector boson (a dark photon for a U(1) Abelian DM model) has been given in ref. [16], where the convolution of the in-vacuum dissociation cross section with the thermal distribution of the vector boson has been introduced by relying on a network of Boltzmann-like equations for the heavy pairs.In order to provide the full thermal width at the leading order, and the corresponding bound-state dissociation cross section, all the diagrams in figure 7 have to be considered.The calculation is very similar to the case just shown, and from the total thermal width, Γ bsd = Γ Q bsd + Γ ∇ bsd + Γ mix bsd , we write the total dissociation cross section as follows The result is organized in a transparent way in terms of pNRY matrix elements, that can be evaluated in a quantum mechanical calculation for the given bound state.As a reference, in the following we list the result for the ground state where w 1 (|k|) ≡ |k|/|E 1 | − 1. Explicit expression for the 2S, 3S, 2P and 3P states can be found in appendix D. In figure 8 we show the dissociation rate for the first nS states.One may notice how for T < ∼ 0.3|E 1 |, the bound-state dissociation is larger for the shallower states.
By adopting the results available in the literature, we can write the bound-state dissociation cross section induced by a dark photon [16] for an Abelian gauge group.It reads . (3.17) For the sake of comparison, we also add the result for a heavy quarkonium [59], that illustrates the case of a non-Abelian vector boson, where the scattering state is in an octet configuration.Then, the unbounded pair in the final state experiences a repulsive interaction .18) where ρ ≡ 1/(N 2 c − 1), with N c being the number of colors.By comparing the powers of the corresponding fine structure constants, one may see how the bound-state dissociation induced by the thermal scalar is α 2 -suppressed with respect to the photo-and gluo-dissociation in eqs.(3.17) and (3.18) respectively.This can be traced back to the different vertices inducing the ultrasoft or thermal transitions (quadrupole and derivative for the scalar case and dipole interactions for the vector mediators).Upon plugging the dissociation cross sections into eq.(3.14), the scalar, photo-and gluo-dissociation rates are given in figure 8 (right panel) for the ground state.

Pair annihilations in pNRY γ 5
As a preparation for the next section where we address the extraction of the DM energy density, let us first discuss the annihilations of heavy pairs.This process is responsible for the depletion of the DM particles into lighter degrees of freedom.Annihilations can equally occur for a particle-antiparticle pair in a scattering state or in a bound state.At variance with the annihilation cross section presented in eq.(2.10) for the scattering states, where the heavy fermions were taken as free particles, here we shall include the effect of soft momentum exchange as mediated by the scalar particle φ.In order to derive pNRY γ 5 , we integrate out the momentum scale associated to the soft scalar exchange, thus obtaining the leading order Coulomb potential.Accordingly, the bilocal field ϕ in pNRY γ 5 satisfies the Schrödinger equation with the potential induced by the scalar exchange.In the evaluation of the annihilation diagrams we will also encounter the wave-function of the interacting pair that will lead us to quantum mechanical matrix elements entering our quantum field theoretical predictions.
In NRY γ 5 heavy pair annihilations are accounted for by local four-fermion operators.The same process can be also described in pNRY γ 5 , where the four-fermion operators of NRY γ 5 generate local terms in the imaginary part of potential V (p, r, σ 1 , σ 2 ) in eq.(2.12).As for the NRY γ 5 , the relevant ingredient is the imaginary part of the matching coefficients as given in eqs.(2.6) and (2.8a)-(2.8c).As mentioned earlier, the presence of the pseudoscalar coupling allows for velocity-independent pair annihilations (cf.eq.(2.10)) which are typically more relevant for the freeze-out dynamics that fixes the DM energy density.Nevertheless, a non-trivial dependence on the relative magnitude between g and g 5 can make velocity-dependent annihilations equally relevant, which is why in the following we explicitly include the corresponding operators.We write the annihilation term from the Lagrangian density of pNRY γ 5 as follows [54,92] L ann where S is the total spin of the pair (S 2 = 0 for spin singlets and S 2 = 2 for spin triplets), while T ij SJ and Ω ij SJ are spin projector operators (cf.e.g.[54,92]).We did not write the R and t dependence in the argument of the field ϕ to avoid cluttering the notation.Some comments are in order.First, the Lagrangian term given in eq.(4.1) is suited to provide the annihilation of the heavy fermion-antifermion pair both in a scattering or a bound state.The different normalization of the states, as outlined in section 3.1 and section 3.2, will automatically determine the corresponding observable: on the one hand a cross section for the annihilation of scattering states, on the other hand a bound-state decay width.Second, the operators are contact terms, as inherited from NRY γ 5 , and the wave functions and their derivatives contribute only at the origin (r = 0), as originally showed for NRQCD [82].Rigorous formulations for quarkonium wave functions have been developed in the pNRQCD formalism [52][53][54].
In pNRY γ 5 we can express the decay width from the self-energy of the field ϕ in the same fashion as we did it for the bound-state dissociation/formation at one loop.This amounts to a simple vertex insertion at tree-level as shown in figure 9. Then the decay width for nS bound states can be written as where the subscript "ann" stands for annihilation.At lowest order in the couplings, as one may see from eq. (2.6), only bound-states in a spin-singlet configuration can annihilate into a pair of scalars.Moreover, the operator in the second line of eq.(4.1) does not contribute because the derivative of nS wave functions vanishes at the origin.The corresponding decay width reads where in the last step we used |R (0) nS (0)| 2 = 4/(n 3 a 3 0 ) from the wave function given in appendix C, the Bohr radius and the binding energy at leading order from eq. (2.16).Let us also comment on the second term in the curly brackets in eq. ( 4.3).This term is generated by the operator in the third line of eq.(4.1), and corresponds to the wave function combination Re(R * nS ∇ 2 r R nS )(r), that diverges at r = 0.As noted in the context of NRQCD, and subsequently reinterpreted in pNRQCD, some perturbative matrix elements are indeed UV divergent and require regularization and renormalization.Here, we employ the relation Re [92], which holds in dimensional regularization, and leads to the expression in eq. ( 4.3).Alternative possibilities to regularize this divergence employ a hard cut-off [82], or well dimensional regularization in conjunction with the MS scheme [93].
Let us remark that, in principle, these annihilation rates can be also obtained directly in NRY γ 5 by relating the bound-state-to-vacuum matrix elements to the wave functions at the origin.The corresponding discussion in the context of NRQCD and heavy quarkonia can be found e.g. in section III C of [51].In this respect the situation in NRY γ 5 is even simpler as compared to NRQCD since our matrix elements are perturbative (at least in the Coulomb-limit) and do not contain covariant derivatives.The decays of heavy fermions into two scalars in NRY γ 5 correspond to the electromagnetic decays of η q or χ qJ with q = c, b into 2 photons in NRQCD.
For the ground-state and the first two excited 2S and 3S states, one finds from eq. ( 4.3), and the matching coefficients in eqs.(2.6) and (2.8a)-(2.8c) (4.4) The same derivation holds for nP J states, where the decay width reads It is easy to find the result for the two combinations of the total angular momentum J = 0, 2 by using the matching coefficients in eq.(2.8a)-(2.8c)and |R (0) nP (0)| 2 = 4(n 2 − 1)/(9n 5 a 5 0 ).In order to compute the annihilation cross section, we instead evaluate the annihilation vertex on a scatterings state, where we include the spin average factor 1/4 as it was done in eq.(2.10) where in the second line used the expressions for the radial wave function at the origin for a DM pair in a scattering state (cf.eq.(C.4)), while the S-and P -wave Sommerfeld factors are defined as follows with p = M v rel /2.We remark that also the quantity Re(R * 0 ∇ 2 r R 0 ) r=0 appearing in eq.(4.6) is divergent.We employ the analog expression already exploited for the bound states for the positive part of the spectrum, namely the scattering states.Hence, we trade the divergent combination for a finite expression using Re( We would like to highlight that in pNRY γ 5 , the annihilation cross sections automatically include the Sommerfeld enhancement originating from the attractive Coulomb potential, namely S(ζ) and S p (ζ) that agree with the literature (cf.e.g.[94][95][96]).In other words, when computing the annihilation cross section for scattering states, the resummation of multiple soft-scalar exchanges (ladder diagrams) is already taken care of, since pNRY is a quantum field theory of interacting pairs.In the limit α 5 → 0 our eq.(4.6) reproduces the result from ref. [91] for the same model.

Phenomenological discussion on the energy density
The depletion of DM through bound states in the early universe depends on many aspects: which darkonium states are formed, how fast they decay and how efficiently they are dissociated by the interactions with the medium.The thermal rates we have computed in section 3.1 and section 3.2 can serve as ingredients when plugged into some rate equations that govern the time evolution of bound and scattering states.Our derivation of thermal cross sections and dissociation rates is based on a quantum field theoretical treatment at finite temperature, meaning that the evolution equations should be worked out in the same fashion.A promising approach that brings together pNREFTs and the evolution equations of heavy pairs in a thermal environment has been recently put forward in the case of heavy quarkonium [97][98][99].However, a derivation of quantum evolution equations for DM pairs in a thermal environment is beyond the scope of the present work. 10evertheless, there already exists a treatment of bound states in terms of semi-classical Boltzmann equations, which allows estimating the bound-state effects on the DM energy density [16,101]. 11In the most general case, the situation is rather intricate since there is an equation for the single DM particle number density, denoted by n X , and an equation for the number density of each bound state, denoted by n ϕ b .All possible reactions have to be taken into account.Yet, under the assumption that the bound states are kept close to the chemical equilibrium through their direct and inverse decays, that is ensured by having Γ ann H, where H is the Hubble rate, one can obtain algebraic equations for each n ϕ in terms of its equilibrium counterpart n eq ϕ b and n X , and reuse these results in the Boltzmann equation for n X .Such an approach was originally suggested in [101], and it is often adopted in the literature on the subject.Here, the annihilations of DM particles are treated in an effective way with a single Boltzmann equation for the DM number density n X which reads dn where the Hubble rate can be expressed in terms of the energy density H = 8πe/3/M Pl , where M Pl is the Planck mass with M Pl = 1.22 × 10 19 GeV, while e = π 2 T 4 g eff /30 and g eff denotes the effective number of relativistic degrees of freedom.The effective thermally averaged cross section, when neglecting bound-to-bound transitions, is where the sum runs over bound states.The quantities involved are the thermally averaged annihilation cross section for the pair in a scattering state σ ann v rel , the thermally averaged bound-state formation cross section σ n bsf v rel , the state bound-state decay width Γ n ann and the bound-state thermal width Γ n bsd , the latter accounting for the dissociation process.As for the cross sections, we are going to use the standard definition of the thermal average (cf.e.g.[16,103]).Then, the combination of the decay width and dissociation width, that we label with I n ≡ Γ n ann /(Γ n ann + Γ n bsd ), determines to what extent DM annihilations via bound states are efficient.One typically has to wait until the temperature, that sets the scale for the energy of the light particles that hit the bound states, is of order of the binding energy of the bound states or smaller.A small population of ionized bound states corresponds to Γ bsd Γ ann and therefore to I n 1.A more quantitative assessment is model dependent, and one has to look at the rates that appear in eq.(5.2).In the following, we shall consider values of the pseudoscalar coupling that satisfy g 5 < g according to our discussion in section 2.3.
For the model under consideration benchmark values for I n are summarized in figure 10 for the ground state.Three benchmark values of the coupling constant are chosen as α = 0.05, 0.1, 0.3 for a fixed combination g 5 = g/5.In addition to that, in the right panel, we fix α = 0.1 and vary the ratio of the coupling g 5 /g for the ground state.Let us observe that in the scalar case, there is a visible dependence on α.This is different with respect to the case of the vector mediator model in ref. [16] (black-dashed line in the left panel of figure 10).In this latter case, the combination of the dissociation rate and bound-state decay is such, that the resulting ionization factor is α-independent. 12Also, for the model under study, we obtain larger values of I 1S with respect to the Abelian vector mediator, pointing to an earlier contribution of bound-state effects.Nevertheless, one has to bear in mind that the bound-state formation process mediated by a scalar mediator features an α 2 suppression [35], and so the overall impact of bound-state effects on the DM energy density is a non-trivial combination of I n and σ bsf v rel for different models.
In the right panel of figure 10, one may see how decreasing values of α 5 , that make the bound-state decay width in eq.(4.4) smaller, but at the same time leaving unaffected the dissociation thermal width, determines a larger ionized population of bound states in the same temperature window.The dashed lines correspond to different states (1S, 2S and 2P 0 ) for the same value of g 5 = g/5.Excited states remain ionized for longer times.
As a final summarizing result, we give the solution of the effective Boltzmann equation (5.1).As usual, we solve the evolution equation for the DM yield Y X = n X /s, where s = 2π 2 h eff T 3 /45 is the entropy density; we assume that the dark states are at the same temperature as the plasma of SM particles.The thermalization condition can be easily satisfied via interactions between the scalar mediator and the Higgs boson.The typical size of the couplings that ensure thermalization at high temperatures are typically much smaller than the couplings g and g 5 that we consider in this work.Therefore, the corresponding contribution to dark matter annihilations and bound-state decays can be safely neglected.
As for h eff and g eff , we take their temperature-dependent values from ref. [104].In the left plot of figure 11, the curves reproduce the observed DM energy density Ω DM h 2 = 0.1200 ± 0.0012 [105] in the parameter space (M, α).Both the Sommerfeld enhancement and bound-state effects are more prominent for large values of α as expected.One may notice the large impact upon including the Sommerfeld enhancement, both for the S-and P -wave annihilations (cf.eq.(4.6)), with respect to the free annihilation cross section (cf.eq.(2.10)).In contrast to that, the contributions of the bound-state formation and decays are moderately relevant, albeit still non-negligible.
In order to single out the bound-state effects, we consider the ratio between the DM energy density as obtained by including or omitting the bound-state formation and decays (namely the second term on the right-hand side of eq.(5.1)), in addition to the Sommerfeld enhancement for the pair annihilations.This is shown in the right plot of figure 11.The ratio is smaller than unity because, when neglecting bound-state formation, a smaller population of DM particles is annihilated away, and a larger energy density is found.The formation of bound-states and their decays into light scalars act as an additional channel for depleting DM.We notice that a smaller pseudoscalar coupling makes the bound-state formation term more relevant, and this traces back to the different powers of α and α 5 entering various cross sections. 13For the smaller ratio between the pseudoscalar and scalar couplings we consider, g 5 = g/10, and the largest mass compatible with the observed relic density for α max = 0.3, namely M = 13 TeV, we find a 30%-36% effect on the predicted energy density if bound-state formation is not included.The dotted lines correspond to the inclusion of the ground state only, whereas upon adding the excited states (nS and nP states with n = 2, 3) we obtain the dashed lines.Bound states with n > 4 induce a contribution smaller than the uncertainty in the observed DM energy density in our numerical implementation.
Upon solving the Boltzmann equation (5.1), we start the integration at temperatures larger than those strictly admitted by the hierarchy of scales (1.1), which implies that T M α is not always satisfied.This would amount to deriving the bound-state formation and dissociation rate at such temperatures, in the so-called screening regime, cf.refs.[48,49] for QCD and QED, and refs.[27,29] for a DM application. 14However, due to our assumption for the self-coupling of the scalar, and upon neglecting the explicit interactions in L portal , a thermal mass for the scalar mediator does not appear and, moreover, 2 → 2 scattering processes with the medium constituents mediated by a soft scalar are absent at least in the Hard Thermal Loop limit [107].Despite of the fact that bound-state formation and corresponding decays are more effective at later temperature stages, further investigations on this aspect can be worthwhile, and are left for future research on the subject.

Conclusions
In this paper we have taken the first step towards the generalization of a potential nonrelativistic effective theory for scalar mediators at finite temperature.In addition to the dynamically generated energy scales typical of bound states, there are also thermodynamic scales at play.Most notably, these are the plasma temperature and thermal masses.Our main goal is to contribute to the recent effort and the development of theoretical tools for estimating the near-threshold and bound-state effects on DM freeze-out calculations.In this scenario DM particles are indeed non-relativistic and slowly moving in the early universe plasma in the relevant temperature window.This means that suitable NREFTs can help in inspecting the corresponding dynamics and calculating thermal cross sections and widths.
We considered a simplified DM model where the dark particles are Dirac fermions and antifermions interacting via a light scalar mediator.With respect to our previous work, we have included the pseudo-scalar coupling between the mediator and the dark fermions.This 13 For the case when only the ground state is considered, and looking at the S-wave annihilations only, the effective cross section is σ eff v rel ≈ 2παα5 S(ζ) /M 2 + πα 4 S bsf (ζ) R1S/M 2 .Here, S bsf (ζ) can be read off from eqs. (3.6) and (3.8a).
14 Pair annihilations would be affected in terms of the thermal potential at such temperatures, and the corresponding modifications to the wave functions from a plasma-modified Schrödinger equation, cf.refs.[25-27, 29, 106].From an EFT perspective, one has to integrate out the temperature scale before obtaining a pNREFT, the latter finally attained when removing potential scalars.
very fact has important consequences for the processes that drive the thermal freeze-out of the dark fermions, and ultimately on the corresponding phenomenology.In particular, velocity independent pair annihilations are possible.We recast the relevant processes in terms of low-energy theories, where the degrees of freedom are heavy DM pairs, both in a scattering and in a bound states, as well as ultrasoft and thermal scalars.In our framework, the bound-state formation, bound-state dissociation, annihilations of the unbound pairs and decays of bound states are inferred from (thermal) self-energies of the bilocal fields in pNRY γ 5 .We investigated a specific hierarchy of scales and set of choices for the couplings, that allow us to neglect thermal masses in our present study.According to our assumed hierarchy of scales, the bound-state formation cross section factorizes into an in-vacuum contribution and a thermal factor.The latter arises from the thermal character of the scalar mediator (cf.eq.(3.6) and eq.(3.7)).The thermal factor implements and resembles the Bose enhancement of the emitted scalar in the formation of a bound-state through radiative emission of the mediator.
A special effort was dedicated to obtain explicit analytic results for the matrix elements of the corresponding pNREFT under the assumption of Coulombic wave functions.We have derived the bound state formation for the ground state, and the excited states in nSand nP -wave configurations for the principal quantum numbers n = 2, 3.
The thermal width for the bound state, namely its thermal break-up as induced by a thermal scalar particle, was extracted from the thermal self-energy of the corresponding wave function field.It is worth noting that this approach features a manifest factorization between an in-vacuum dissociation cross section and the thermal distribution of the scalar particles in the medium (cf.eq.(3.14)).Most importantly, the in-vacuum dissociation cross section can be written in terms of generic pNRY γ 5 matrix elements.These quantities represent genuine quantum mechanical expectation values that, depending on the states under consideration and the assumed nature of the potential, can be evaluated by solving the corresponding Schrödinger equation analytically or numerically.The scalar-induced thermal break up of a bound state is the analog of the photo-and gluo-dissociation, the latter being particularly relevant for heavy quarkonium physics in heavy ion collisions.
As third and last piece of information, we considered scattering state annihilations and bound state decays, that were written in terms of local operators in pNRY γ 5 .Here the matrix elements arising in the course of the calculation can be naturally understood as expectation values of quantum mechanical operators, where the wave functions follow from solving Schrödinger equation with the corresponding pNRY γ 5 potential.This is an obvious advantage as compared to NRY γ 5 , where the connection between the given matrix element and the corresponding quantum mechanical expression is not always straightforward.In this respect the situation in pNRY γ 5 is much more clear.First, in the case of an unbound above-threshold pair, the Sommerfeld enhancement shows up when computing the annihilation cross section.Second, for the negative part of the spectrum, namely bound states, the bound-state decay width is recovered.As an original result of our work, we provided the Sommerfeld enhanced annihilation cross section for the model in eq.(2.1) at O(α 2 v 2 rel , αα 5 v 2 rel , α 2 5 v 2 rel ).We exploit the so-obtained rates in an effective Boltzmann equation, routinely used in phenomenological studies, in order to estimate the effect of the Sommerfeld enhancement, bound-state formation and decays on the DM energy density.We find that a large effect is played by the Sommerfeld enhanced annihilations for the scattering states.In this model, that features a non-vanishing pseudoscalar coupling, S-wave velocity-independent annihilations are possible, in addition to the P -wave annihilations of the sole scalar interaction case.The former boost DM annihilations in a prominent way, and the DM mass compatible with the observed relic density is pushed to much larger values with respect to those obtained with free annihilations for the same values of the couplings.Moreover, despite the bound-state effects being more moderate, we find them to be non-negligible.Their impact depends on the ratio g 5 /g between the couplings in the Lagrangian model (2.1).In particular, decreasing the ratio g 5 /g makes bound-state formation more prominent.For g 5 /g = 0.25 (g 5 /g = 0.1) we find a 16% (35%) effect on the DM energy density with respect to the case one includes only the Sommerfeld enhancement for the scattering states.The effect of the excited state decays accounts for few per-cents depending again on the particular values of the couplings.A highlight of our study is the observation, that even when the bound-state formation is driven by a scalar mediator, bound-state effects may affect the determination of the DM energy for the simplified model (2.1).The main reason for this is rooted in the presence of a pseudoscalar interaction, and in a non-trivial interplay between bound-state decay widths, bound-state formation and dissociation, and Sommerfeld enhancements.
We conclude by remarking that, in order to make connections with experimental constraints on the parameter space of the model, one has to specify the portal interactions and inspect the corresponding possible modifications on the various ingredients presented in our study.the so-obtained operator basis may not be the most useful one and that by construction we miss operators whose Wilson coefficients happen to be loop-induced.Here we would like to refer to appendix A.3 of ref. [57] for a detailed treatment of the purely scalar Yukawa theory.In the following we will make use of the concepts and the notation introduce there.
The leading order ansatz Ŝ for the mixed case at hand reads while the Hamiltonian is given by It is worth noting that the ansatz given in eq.(A.1) is identical to Ŝ required for the pure pseudoscalar Yukawa theory.When doing the decoupling of even and odd operators we need to keep in mind that βγ 5 and α i 1 . . .α i 2n γ 5 with n ∈ N are odd, while α i 1 . . .α i 2n+1 γ 5 is even.
At O(1/M ) the bilinear Lagrangian is just a direct sum of the non-relativistic Lagrangians for the pure scalar [55,56] and pseudoscalar [80] Yukawa theories.Operators with tree-level matching coefficients induced by the mixing of the two interactions start occurring at O(1/M 2 ).For the sake of completeness, let us explicitly provide the O(1/M 2 )piece of the bilinear part of the NRY γ 5 Lagrangian with tree-level matching coefficients obtained using FWT transformations where in (∂ 0 φ) the time derivative acts only on the φ field, while ∇ always acts on everything to its right.

B Derivation of pNRY γ 5 from NRY γ 5
In section 2.3 it was stated that pNRY γ 5 can be obtained from NRY γ 5 by projecting the Hamiltonian of NRY γ 5 onto the particle-antiparticle sector.Since the technicalities behind this approach are rarely discussed in the literature, here we would like to provide some more details on the corresponding calculation for the bilinear part of the NRY γ 5 Lagrangian.Switching to H γ 5 effectively means multiplying the expressions given in eqs.(2.3) and (2.4) by −1 and removing the kinetic terms with time derivatives ψ † i∂ 0 ψ and χ † i∂ 0 χ.Then, we must sandwich the Hamiltonian between the following ket and bra states In order to evaluate the resulting expression we need to move all ψ-and χ † -fields to the very right and all χ-and ψ † -fields to the very left, since To this end we make use of the equal-time anticommutation relations obeyed by the Pauli fields In the case of spatial derivatives acting on one of the fields, it is convenient to let the derivative formally act on the Dirac delta e.g. as in When evaluating the integral over x we can always integrate by parts to move ∇ x away from δ (3) (x − x ).
The bilocal field ϕ ij carries two Pauli indices and behaves as Its Hermitian conjugate is given by15 and consequently When projecting the Hamiltonian we encounter following Pauli structures involving the bilocal field Sometimes NREFT practitioners prefer to trade the field χ for χ c defined as The reason for this is that while the basis (ψ, χ) is best suited for studying annihilations and formations of heavy fermions, the set (ψ, χ c ) often turns out to be more convenient when dealing with scattering processes.Switching from χ to χ c on the level of the pN-REFT Lagrangian can be easily achieved via a suitable field redefinition, which amounts to introducing where in the first step we transposed the initial expression with respect to its Pauli indices and the employed eqs.(B.13).Furthermore, we used that In this context it is worth noting that, in general, one can eliminate χ in favor of χ c already at the level of the NREFT Lagrangian and then directly project onto In this case one would first apply eq.(B.15) to the bilinear part of the NREFT Lagrangian and then make use of the Fierz transformations in the Pauli space when rewriting the four-fermion operators accordingly.Hence, even though physical results clearly do not depend on whether we use χ or χ c , the former turns out to be more useful to study production or annihilation of heavy particles (as is often done in NREFTs), while the latter is very convenient for scattering processes typically arising in pNREFTs.
Finally, for the sake of completeness, let us provide the explicit O(1/M ) result for the projection of H γ 5 onto the subspace spanned by eqs.(B.1) before the multipole expansion Then we only need to switch to the center-of-mass coordinates via and multipole expand the scalar fields in r up to the desired order using to arrive at the final bilinear piece of the pNRY γ 5 Lagrangian at O(1/M ).

C Matrix elements
The matrix elements presented in the main body of the paper can be evaluated analytically in the Coulomb limit, where we essentially employ analytic solutions for the bound and scattering state wave functions of the Hydrogen atom.The bound state wave function reads where Y m (θ, φ) is the spherical harmonic, while 1 F 1 (a; b; z) denotes the confluent hypergeometric or Kummer's function.The scattering wave functions arise as solutions of the Schrödinger equation with Z = −1 for the Hydrogen atom.Solving this equation in the parabolic coordinates and choosing the solution that describes the wave before approaching the source at the origin we can expand the resulting wave function in plane waves.This yields Notice that the scattering state wave function is decomposed in terms of partial waves with definite angular momentum and that we choose the relative momentum p to be aligned along the z-axis.This agrees with the conventions commonly used in quantum mechanics textbooks [108,109] and DM literature [34,38].

C.1 Expectation values of r 2
We start with the discussion of matrix elements involving the position operator squared This integral obviously factorizes into a radial and an angular part.The integration over the solid angle can be straightforwardly evaluated using the following identity 1 F 1 ( + 1 + iζ; 2 + 2; +2ipr) 1 F 1 + 1 − n; 2 + 2; 2r na 0 (C.9) for generic quantum numbers n and can be calculated with the aid of the relation (cf.section 7.622 in [110]) (C.12) Notice that the arising phase cancels out when squaring the above quantity. 17Let us also remark that if we are interested only in the expectation values for fixed quantum numbers (e.g.n = 1, = 0, m = 0 etc.), then the radial integral in eq.(C.9) simplifies and can be evaluated using a much simpler identity (cf.section 7.522 in [110]) that reads with Re(b) > 0 and Re(z) > max(Re(k), 0).Finally, let us provide explicit results for the phenomenologically relevant lowest-lying S-and P -wave matrix elements.The matrix elements squared for the first three nS states 16 Notice that 1F * 1 ( + 1 − iζ; 2 + 2; −2ipr) = 1F1( + 1 + iζ; 2 + 2; 2ipr) for p, ζ ∈ R. 17 It should be noted, that for nS states, the remaining phase e −2i(n−1) arccot(ζ/n) is exactly canceled by the factor that appears in 2F1(1 − n, 1 + iζ, 2, 4inζ/(n + iζ) 2 ) when setting n to specific integer values.The matrix element is then purely real.We have explicitly verified these results by calculating the matrix elements either via the master formula eq.(C.12) or on a case-by-case basis using eq.(C.13).Furthermore, we conducted a numerical check that gave a perfect agreement over a large range of ζ values.

C.2 Expectation values of ∇ 2
As far as the evaluation of the matrix element where we would like to remark that matrix elements for the P -wave states |21 ± 1 and |31 ± 1 vanish.The so-obtained results are also in perfect agreement with our numerical checks.

C.3 Expectation values of r i r j
Last but not least, we also need to deal with the tensor matrix elements given by p|r i r j |n m = 1 M 2 d 3 r r i r j ψ * (r)ψ n m (r).(C.18) Here it is useful to observe that the first 9 spherical harmonics Y lm with ≤ 2 and the corresponding m-values can be written as linear combinations of x, y, z, x 2 , y 2 , z 2 , xy, xz and yz.Inverting this linear system we can express each of these 9 quantities as a linear combination of the spherical harmonics.More explicitly, introducing T ij = r i r j we have Notice that while the above results can be found in many quantum mechanics textbooks (cf.e.g.[108]), the algorithm itself naturally works also for higher rank tensors beyond

D Bound-state dissociation cross sections
In section 3.2 we provided only one explicit result (|1S state) for the bound-state dissociation cross section.For the sake of completeness, here we list the results for the remaining Sand P -wave states that can be calculated using the matrix elements given in appendix C.

Figure 1 .
Figure 1.Leading order NRY γ5 diagrams for the potential matching with scalar and pseudo-scalar interactions required to calculate V (p, r, σ 1 , σ 2 ).

Figure 2 .
Figure 2.Self-energy of a heavy pair in a scattering state (solid double line) induced by the quadrupole (crossed vertex) interaction, derivative interaction (triangle vertex) and their mixing.The internal single line stands for a pair in a bound state, the dashed line is the scalar mediator.

Figure 3 .
Figure 3. Contour levels for the ratio of the thermal and in-vacuum cross sections 1 + n B (∆E p n ) in the (|E 1 |/T, ζ)-plane.The solid (dashed) lines stand for a 50% (100%) increase of the in-vacuum cross section due to the Bose enhancement.

4 /M 2 )Figure 5 .
Figure 5. Bound-state formation cross section for the |2P and |3P states as a function of ζ.The corresponding bound-state formation for the |2S and |3S states are shown.The in-vacuum and thermal cases, for |E 1 |/T = 1 are considered.

Figure 7 .
Figure 7.Self-energy of a heavy pair in a bound state (solid line) induced by the quadrupole (crossed vertex) interaction, derivative interaction (triangle vertex) and their mixing.The internal single line stands for a heavy pair in a scattering state.

Figure 9 .
Figure 9. Annihilation diagrams in pNRY γ5 for a scattering state (left) and a bound state (right).

Ω 5 g 5 =g/ 4 g 5 =g/ 5 g 5 Figure 11 .
Figure 11.Left: parameter space in the (M, α)-plane that reproduces the observed DM energy density for g 5 = g/5.The purple dotted curve accounts for the ground state only, whereas the purple-dashed curve corresponds to the inclusion of excited states.Right: ratio between the DM energy density as obtained with and without bound-state effects.The parameter α is fixed as a function of the mass M to account for the observed DM energy density when bound-state effects are considered.Dashed (dotted) lines do (not) include excited states.
The ratio σ nS bsf v rel /σ 1S bsf v rel for n = 2, 3 is shown for both the in-vacuum (gray-dotted line) and thermal cross sections (colored curves).