Non-relativistic and potential non-relativistic effective field theories for scalar mediators

Yukawa-type interactions between heavy Dirac fermions and a scalar field are a common ingredient in various extensions of the Standard Model. Despite of that, the non-relativistic limit of the scalar Yukawa theory has not yet been studied in full generality in a rigorous and model-independent way. In this paper we intend to fill this gap by initiating a series of investigations that make use of modern effective field theory (EFT) techniques. In particular, we aim at constructing suitable non-relativistic and potential non-relativistic EFTs of Yukawa interactions (denoted as NRY and pNRY respectively) in close analogy to the well known and phenomenologically successful non-relativistic QCD (NRQCD) and potential non-relativistic QCD (pNRQCD). The phenomenological motivation for our study lies in the possibility to explain the existing cosmological observations by introducing heavy fermionic dark matter particles that interact with each other by exchanging a light scalar mediator. A systematic study of this compelling scenario in the framework of non-relativistic EFTs (NREFTs) constitutes the main novelty of our approach as compared to the existing studies.

On top of being desirable from the phenomenological and observational points of views, the possibility of a richer dark sector, that comprises more than one particle, is fairly common in many DM models, cf. e.g. [21][22][23]. The dark particles can enjoy their own hidden forces, which are far less constrained than the interactions between DM and Standard Model (SM) degrees of freedom. Furthermore, the existence of light (i.e. with masses much smaller than that of the actual DM particles) mediators may affect the DM dynamics in multiple ways. Most notably, whenever DM particles are slowly moving with non-relativistic velocities, light mediators can induce bound states in the dark sector in the early universe and/or in the dense environment of present-day haloes [18,[24][25][26]. As for the above-threshold states, the effect of repeated mediator exchange manifests itself in the so-called Sommerfeld enhancement for an attractive potential [27,28]. In this context the role of a light mediator can also be played by SM particles. For a sufficiently heavy DM these may even be weak gauge bosons [24,25,[29][30][31][32] or the Higgs boson [30,33]. This latter option is becoming increasingly relevant as null searches for new physics at the LHC are pushing the scale of possible novel particles, including many thermally produced DM candidates, into the multi-TeV region. 1 Depending on the model at hand, one may find unstable bound states, that usually appear in symmetric DM models, as well as stable bound states (the latter are part of the present-day DM energy density). Typically, the annihilating particle-antiparticle pairs feel an attractive potential that can not only drastically change the annihilation cross section via Sommerfeld enhancement but also induce bound-state formation [18,26]. Once bound states are formed, and not effectively dissociated in the thermal plasma, they provide an additional channel for the depletion of DM particles in the early universe. The relic density determination has to be adjusted accordingly, since substantial annihilations may still occur after the chemical decoupling. This typically results in (i) mapping out different combinations of DM masses and couplings that reproduce the observed DM cosmological abundance Ω DM h 2 = 0.1200 ± 0.0012 [37]; (ii) a reinvestigation of DM phenomenology due to the interplay between the model parameters that fix the relic density and guide the experimental strategies. The stable bound states that often arise in asymmetric DM models affect the detection strategy and experimental searches for both indirect [38][39][40][41][42] and language and methods of NREFTs. Apart from the NRY we also consider the pNREFT version of the scalar Yukawa theory, which we call potential non-relativistic scalar Yukawa theory (pNRY). It is worth noting that the effect of adding interactions between heavy fermions and the Higgs to the conventional pNRQCD (which naturally leads to Yukawa potentials) has been considered e. g. in [73,74] when studying tt-production near threshold.
At variance with the previous works, here we are interested in the pure Yukawa theory that lacks any interactions with gauge bosons such as photons, gluons, W or Z. Furthermore, we would like to abstain from introducing any additional symmetries apart from what is already present in the scalar Yukawa theory. In our view, this approach allows us to investigate and highlight the essential features of non-relativistic scalar Yukawa interactions in a clear and transparent fashion without making any assumptions on the nature of the underlying higherenergy theory. The aim of the present work is to revisit the construction of the NRY by extending the treatment of [68,69] and to explore the consequences of the resulting NREFT and pNREFT for the DM phenomenology, where we are interested in describing the interactions of heavy Dirac fermions X with a much lighter scalar field φ. To the best of our knowledge, pNRY as a pNRQED-like theory that contains solely Yukawa interactions is presented in this work for the first time. In both cases we explore possible hierarchies of scales and discuss the appropriate power-counting rules. In this paper, we shall focus on the zero temperature case, and only marginally comment on the finite temperature generalization.
It is worth noting that the DM model under consideration has some intriguing properties that are unique to heavy fermions exchanging a scalar. First, as opposed to the vector mediator case, the annihilation of heavy particle-antiparticle pairs at leading order in the velocity and 1/M expansion proceeds via a P -wave process. More explicitly, one finds that the matching coefficients of 4-fermion dimension-6 operators vanish at O(α 2 v 0 ), whereas the first non-vanishing contributions show up in the velocity suppressed dimension-8 operators. Second, the pNRY exhibits, already at the Lagrangian level, the absence of electric-dipole transitions and the presence of monopole and quadrupole interactions between a heavy pair and the scalar mediator. In the context of pNREFTs, monopole interactions were discussed for super-symmetric Yang-Mills theories at weak coupling in [75]. Finally, in the case of vector mediators, pNREFTs have been already fully, or at least to some extent, exploited in the context of DM with and without co-annihilating partners [49,51,52,[76][77][78].
The structure of the paper is as follows. In section 2 we briefly introduce the simplified model that we take as our high-energy (in the EFT sense) theory. Then, in section 3 we address the construction of the low-energy NREFT (denoted as NRY) for non-relativistic fermions and antifermions exchanging a scalar. Here we shall give the set of operators as an expansion in 1/M , v and coupling constants, and discuss the symmetries and power counting rules of the low-energy theory. In section 4 we apply the NRY formalism to describe DM interactions and provide the results for the matching coefficients. As far as the fermion bilinears are concerned, we shall be content with tree-level matching coefficients. The matching for 4-fermion dimension-6 and dimension-8 operators will be carried out at O(α 2 ). These operators encode the hard contribution to the annihilation cross section for the process XX → φφ. In section 5 we proceed to the derivation of the pNREFT (denoted as pNRY), whose degrees of freedom are bound states, scattering states with kinetic energy of order M v 2 and ultrasoft scalar particles. We perform the potential matching at O(M α 4 ) and then provide an application of pNRY to the derivation of the discrete spectrum and the calculation of the bound-state formation cross section. Conclusions and outlook are offered in section 6.

Dark matter model
In this section we briefly introduce the DM model under consideration and discuss the relevant degrees of freedom. We assume DM to be a Dirac fermion singlet under the SM gauge group that it is coupled to a scalar particle with a Yukawa-type interaction. The Lagrangian density of the model reads [79,80] where X is the DM Dirac field and φ is a real scalar field. The scalar self-coupling and the Yukawa coupling between the fermion and the scalar fields are denoted as λ and g respectively. The mass of the scalar mediator m is assumed to be much smaller than the DM particle mass M , m M . Here we adopt a simplified model realization, where the question of the fermion mass generation and of the gauge group governing the dark sector are ignored. 2 Our aim is to consider the Lagrangian given in eq. (2.1) as one of the simplest representatives for a family of minimal DM models [22,80] with a light scalar mediating interactions between DM particles. It goes without saying that such a scenario admits different realizations that can be much more involved than a single Yukawa interaction (cf. e.g. [55,59,83]).
Next, L portal accounts for the interactions between the scalar φ and other degrees of freedom that can be either in the dark sector (e.g. all particles lighter than φ), and/or in the SM sector. The most common realization of such a portal involves interactions with the SM Higgs boson. In general, portal interactions are needed because the light scalar particles φ are abundant in the early universe and a substantial population is still present after the freeze-out of the dark fermion. Hence, there has to be a mechanism that allows φ particles to decay and deplete their population so that the scalar does not happen to dominate the energy density of the Universe [80,84]. The minimal model of eq. (2.1) is moderately under tension if one considers the interactions of the scalar φ with the Higgs boson and hence SM fermions. Especially the interactions with quarks severely constrain the model via direct detection experiments [80,85]. However, these tensions can be removed in a number of different ways [83,86]. Since a detailed phenomenological analysis is beyond the scope of this work, we do not specify L portal further and merely focus on the complementary terms in eq. (2.1) to derive the low-energy field theories relevant for the bound-state dynamics. This sets the stage for our NREFT and pNREFT formulations and paves the way for more thorough investigations (also with respect to the DM phenomenology) in future works.

Non-relativistic Yukawa theory
In the following we would like to discuss the procedure of constructing a tower of non-relativistic EFTs for a heavy Dirac fermion X that interacts with a light scalar field φ via a scalar Yukawa interaction. Our main motivation is to investigate the properties of XX bound states such as spectra, production and decays in a rigorous and model-independent way. In order to proceed systematically, it is useful to disentangle low-energy modes relevant for the bound-state formation from high-energy modes that are naturally present in the UV-complete theory described by eq. (2.1). Nevertheless, the contributions from large energies and momenta are not simply discarded: their effects will be incorporated into Wilson coefficients multiplying the operators that appear in the EFT Lagrangians. The process of determining these coefficients by comparing Green's functions of two theories at low-energies is called matching. The EFT description can be systematically improved by including higher-order operators compatible with the symmetries of the underlying theory. The effects of these operators can be quantified using EFT power-counting rules, so that at each order in the relevant expansion parameters only a finite number of operators must be taken into account. This leads to a comprehensive description of the low-energy physics, that allows us to make predictions for the physical observables of interest (e.g. cross sections or decay rates) in a simple and straightforward fashion.
Obviously, we need to assume a certain hierarchy between the scales relevant for the nonrelativistic bound states (see figure 1). The largest of these relevant scales is the heavy fermion mass M . An important scale below M is the typical size of the relative momentum between the fermions in a bound state, |p| ∼ M v, where v is the relative velocity of the particles. Notice that this scale is also related to the typical bound state size r, where 1/r ∼ M v (one can use the Bohr radius a 0 for Coulombic bound states for the size estimate). As our fermions are heavy and non-relativistic, we have v c. We assume that v is sufficiently small with at least v 2 ≤ 0.3. In nature v 2 ∼ 0.3 is found e.g. in heavy quarkonia made of a charm and an anti-charm quarks. The non-relativistic description is still applicable to such systems, but the velocity expansion converges rather slowly. On the other hand, for XX bound states with v 2 ∼ 0.1 (as in bbquarkonia), corrections of O(v 2 ) should be sufficient for a reliable phenomenological analysis. The typical bound state energy of an XX system scales as M v 2 . In the following we denote the scales M , M v and M v 2 are hard, soft and ultrasoft respectively.
For simplicity, we would like to consider the situation where the mass of the scalar m is of the same order, or smaller, than the ultrasoft scale m < ∼ M v 2 . In practice, this corresponds to considering Coulombic states induced by the scalar mediator, which is the regime typically studied in the existing literature (however see e.g. [57,74] for numerical studies with finite m). Furthermore, should the full theory feature a scale Λ below which perturbation theory ceases to be applicable (such as Λ QCD in strong interactions), this scale should be much smaller 3 than M v 2 . This ensures that the scales M and M v can be integrated out perturbatively.
Integrating out all degrees of freedom with energies and momenta of order M and above we obtain an EFT known as the Non-relativistic Yukawa Theory (NRY) [68,69]. The degrees of freedom of NRY are Pauli spinor fields ψ and χ describing a particle and an antiparticle respectively 4 as well as soft and ultrasoft scalar fields φ. The Lagrangian of NRY is a double expansion in 1/M and v. While M explicitly appears as a parameter in L NRY , this is not the case for the velocity v. Therefore, to determine the velocity scaling of the given operator it is necessary to work out power-counting rules that assign powers of velocity to the typical operator building blocks, i.e. couplings, fields and derivatives.
The fact that the energies and momenta of the φ fields can be soft or ultrasoft leads to additional complications in the power-counting. In particular, the scaling of scalar mediators involved in potential exchanges between the heavy fermions will, in general, differ from that of the on-shell φ fields in the external states. In other words, the power counting of NRY is not 3 In principle, it would be sufficient to demand only Λ M , which would allow us to integrate out the scale M perturbatively. The procedure of integrating out the scale M v without relying on the perturbative expansion in a small coupling has been discussed in [87,88]. However, to keep the present discussion as simple as possible, we assume perturbativity at least up to scales much smaller than the bound state energy. 4 To be more precise, ψ annihilates a fermion, while χ creates an antifermion. This property is most easily seen in the operator approach, where the free-field Fourier decomposition of ψ contains only a single particle annihilation operatorâ(p, s), while that of χ is proportional to the antiparticle creation operatorb † (p, s).
homogeneous, as it has been already discussed in [68,69]. This is less of a problem for production and decay calculations, but turns out to be rather inconvenient when looking at the bound state properties. In section 5 we will show how to circumvent this problem by devising yet another EFT (pNRY) that works at energies much smaller than M v. By construction, NRY is valid only at scales of order M v and below. In this energy region it must reproduce the full theory, so that both theories have identical infrared (IR) behavior. Formally, the Lagrangian of NRY contains infinitely many operators suppressed by increasing powers of M . This corresponds to the statement that both theories coincide in the limit M → ∞. However, one can always employ the velocity scaling rules to determine which operators contribute at the given order in v. This is why in practice we will only need to consider a small set of relevant operators.

Symmetries and NRY Lagrangian
A crucial property of an EFT is that it must encompass the symmetries of the underlying full theory. Therefore, to construct the Lagrangian of NRY we must write down all possible operators compatible with the symmetries present in the scalar Yukawa theory. For example, each operator must be invariant under charge conjugation, parity and time reversal. Lorentz symmetry is still present in NRY, but it is not manifest. 5 One of the implications thereof is the invariance under rotations in the 3-dimensional space. In addition to that, we will also encounter some symmetries that manifest themselves only when particles and antiparticles are treated as separate degrees of freedom and are not obvious when looking at the relativistic full theory Lagrangian.
The procedure of enumerating all operators that may appear in the given NREFT order by order in 1/M can be found e.g. in [93,94]. This problem can be also approached using the Hilbert series framework adapted to non-relativistic theories [95]. A more explicit way to obtain the fermion-bilinear piece of L NRY is to subject the full theory Lagrangian given in eq. (2.1) to a sequence of Foldy-Wouthuysen-Tani (FWT) transformations [96,97] or to use the equations of motion (EOM) method as in the Heavy Quark Effective Theory (HQET) [98][99][100][101] (cf. [102] for a pedagogical introduction to EOM). Both approaches can be iterated order by order in 1/M and lead to effective Lagrangians that incorporate relevant operators together with their tree-level matching coefficients. At this point it is important to stress that these techniques should not be employed mindlessly for a number of reasons. First of all, it is well-known (cf. e.g. [103]) that FWT and EOM by construction miss all operators that are allowed by symmetries but happen to have vanishing tree-level matching coefficients. 6 This should not come as a surprise, since both procedures essentially correspond to the tree-level matching. Second, the so-obtained 5 A thorough discussion of the Poincaré invariance in NREFTs such as NRQCD and pNRQCD can be found in [89][90][91][92]. 6 It is clear that such operators can still become relevant at higher loop orders and hence must be included in the Lagrangian.
operator basis is not guaranteed to be the most useful one and may contain redundancies. Field redefinitions can be used either to completely eliminate some of the appearing operators or to trade them for other operators. For example, in the case of NRQCD, NRQED or NRY one can get rid of operators with time derivative acting on the heavy fermions by introducing suitable redefinitions of these fields. Notice that field redefinitions leave only on-shell Green functions unchanged but alter the offshell ones. This is why the matching between the full theory and the NRY should be performed for on-shell Green functions. Nonetheless, as long as one keeps in mind the above facts, FWT and EOM can be regarded as a useful aid when working out a new NREFT containing heavy fermions. We demonstrate an explicit application of these tools to the scalar Yukawa theory up to O(1/M 2 ) in appendix A.
At O(1/M 2 ) the most general Lagrangian compatible with the symmetries of the scalar Yukawa theory can be written as where ψ (χ) is the Pauli field that annihilates (creates) a heavy fermion, while φ is the light scalar mediator. The anticommutators are defined as {a, b} = ab + ba. Furthermore, σ stands for Pauli matrices and we have ∂ i = ∇ i . Notice that the derivatives in the bilinear fermion and antifermion sector act on all the fields (scalar and spinors) on the right. The c i and c i are the matching coefficients of fermion and antifermion bilinears respectively, while d i belong to the scalar sector. The L 4-fermions part of the Lagrangian contains 4-fermion contact interactions that describe annihilations/decays of XX pairs 7 . These operators are necessary, since a heavy-fermion annihilation process such as XX → φφ cannot be described via the fermion-bilinear part of the NRY Lagrangian. In this case the scalar fields must carry energies of O(M ), yet these modes have been integrated out when constructing the NRY. This is why such processes must be described via 4-fermion interactions, where the effects of the high energy modes are incorporated in the imaginary parts of the Wilson coefficients multiplying these operators [62,104].
The NRY Lagrangian enjoys a heavy fermion spin symmetry (HFSS) up to corrections of 7 XX production can be described by vacuum expectation values of 4-fermion operators containing XX-Fock states between the fermion bilinears. Such objects are therefore not included in L 4-fermions but will enter the corresponding production cross sections.
O(1/M 2 ), where the first spin-flipping operator shows up. It is interesting to observe that in the case of NRQCD or HQET the heavy quark spin symmetry is broken already at O(1/M ). However, since NRY has no gauge symmetry and an operator proportional toXφ∇ · σX is forbidden by parity, the spin flip may occur only through an operator involving at least two spatial derivatives. The validity of the HFSS up to O(1/M 2 ) implies particularly small splittings in the spin-symmetry multiplets of XX bound states, which is an intriguing feature of the NRY phenomenology. Another symmetry of L NRY that should be familiar to NREFT practitioners is the heavy fermion phase symmetry ψ → e iα ψ, χ → e iβ χ, α, β ∈ R, (3.2) which implies separate conservation of the number of particles and antiparticles.

Power counting
To derive the power-counting rules of the theory we can make use of the standard arguments 8 used in NRQCD [62]. To this end it is useful to employ a quantum mechanical perspective before the second quantization, where we can interpret ψ as a wave function interacting with an external potential φ. The wave function normalization condition together with our previous estimate of the typical bound state radius r ∼ 1/M v readily suggests that d 3 x ∼ 1/(M v) 3 and therefore ψ ∼ (M v) 3/2 . A spatial derivative acting on ψ probes its typical 3-momentum, so that ∇ψ ∼ M v ψ. The equation satisfied by ψ at the lowest order in the 1/M expansion reads where we have anticipated the tree-level results for the matching coefficients c 1 and c 2 (cf. eq. (4.1)).
Here gφ plays the role of the leading-order contribution to the interacting part of the quantum mechanical Hamiltonian. Using the virial theorem for bound states we can estimate that ∂ 0 ψ ∼ M v 2 ψ and gφ ∼ M v 2 . The same argument applies also to the scaling of the χ field. In a similar manner, we may also pass to a picture in which the wave function φ satisfies the following Klein-Gordon-Schrödinger equation (again at lowest order in 1/M ) where the last term scales as gM 3 v 3 . If we assume that the typical momentum of φ scales as M v, then the virial theorem implies that gM v ∼ φ. Hence, g 2 ∼ v for a Coulombic state and φ ∼ M v 3/2 . Notice also that λφ 3 ∼ M 3 v 7/2 may seem much less suppressed than (gφ) 3 ∼ M 3 v 6 . However, since λ does not appear in the fermion-bilinear part of the NRY Lagrangian, a diagram involving λ must also contain at least one insertion of gφ that couples directly to the fermion current. This accounts for an extra suppression of processes involving the scalar self-coupling with λ.
Notice also that if the energy and momentum of φ scale as M v 2 , we find v 4 φ ∼ gM v 3 and consequently g 2 ∼ v 3 , upon using gφ ∼ M v 2 . In this case we would actually need less operators to describe the same observable at the given order in v as compared to the previous counting. Yet, to be on the safe side, in the following we will adopt the more conservative counting with ∇φ ∼ M vφ. We summarize the scaling rules in figure 1.

Applications of NRY to dark matter
In this section we adapt the general discussion of section 3 to the DM phenomenology, and derive the matching coefficients of the low-energy version of the model Lagrangian eq. (2.1), namely the parameters of the NRY (3.1). The effective Lagrangian comprises unknown coefficients that have to be fixed by the matching procedure. In practice, one computes on-shell Green's functions in the full theory in eq. (2.1) and in the effective theory and demands their equality at a matching scale µ match with m, M v 2 , M v µ match M . The relative size of the smaller scales is irrelevant here. Through the matching coefficients, which could be obtained at arbitrary loop order, the low-energy theory is also organized as an expansion in the couplings g and λ. As it is common in DM models, we assume the scalar self-coupling to satisfy λ ∼ g 2 , which facilitates the organization of the perturbation series. Furthermore, we define α ≡ g 2 /4π to organize the power counting of the low-energy theories. In this work λ will barely play any role.
A non-relativistic regime for dark particles is relevant both for annihilations during the thermal freeze-out, as well as in the present-day galactic halos. In the latter case, typical DM velocities are of order 10 −4 -10 −3 in units of c, cf. e.g. [105,106]. In the former case, DM particles are kept in chemical equilibrium through interactions with the thermal bath until T M and gradually freeze out at temperatures T ∼ M/25. 9 Annihilations continue even during later stages where the DM particles are still in kinetic equilibrium. In this situation most of the energy of a DM particle is sourced by its mass and, for non-relativistic species, the typical momentum is |p| = √ T M = M T /M . One usually identifies an average velocity v ≈ T /M , which is smaller than unity in the regime of interest. For above-threshold particle-antiparticle pairs feeling a Coulomb-like potential, the regime v ∼ α signals the potential energy E pot ∼ α/r ∼ M αv being the same as the kinetic energy M v 2 [46]. For XX in a bound state, the velocity estimate of the relative motion is fixed at v ∼ α. In a perturbative regime, this again gives a velocity smaller than unity. Since the temperature of the plasma is T M , the temperature scale is treated on the same footing with other smaller scales, and does not affect neither the matching and nor the form of the NRY (cf. discussions in [107,108] for NRQED and NRQCD at finite temperature and [109] for an explicit derivation of an NREFT for Majorana fermions in a thermal bath).
In summary, at energies much smaller than M , the degrees of freedom are non-relativistic Dirac fermions and antifermions, including bound states and near-threshold states, and scalars with energies and momenta much smaller than M . The NRY presented in eq. (3.1) is then a suitable field theory that describes non-relativistic DM particles and their dynamics. The first two lines of eq. (3.1) encode interactions between the non-relativistic fermion (and antifermion) and the light scalar mediator. An important difference with respect to NRQED/NRQCD is the lack of gauge symmetry. Hence, the effective Lagrangian eq. (3.1) contains no covariant derivatives and the form of effective operators containing scalars and fermions is not constrained accordingly. We discuss the matching of the bilinear sector in section 4.1.1. The last line of eq. (3.1) comprises 4-fermion operators, which account for DM pair annihilations in the low-energy theory. The corresponding cross section XX → φφ is a key ingredient for the determination of the relic density governed by the freeze-out as well as present day annihilations in the Milky-Way. In this work, we do not consider pair annihilations induced by interactions in L portal . We address the 4-fermion operator in section 4.1.2.

NRY Matching
We now discuss the derivation of the matching coefficients of the low-energy theory given in eq. (3.1). As already anticipated, this procedure amounts to enforcing the equality of on-shell scattering amplitudes in the full theory (2.1) with on-shell scattering amplitudes constructed with the general expressions of the NRY in terms of ψ, χ and φ and the unknown matching coefficients. The matching scale provides a UV cut-off for the low-energy theory, above which NRY is not reliable. Clearly, the fundamental theory and the NREFT have a different UV behavior, whereas the infrared (IR) properties are the very same. Only high-energy modes of order M (integrated out in NRY) contribute to the matching coefficients of eq. (3.1). In other words, when computing scattering amplitudes there can be residual IR contributions that are not included in the matching coefficients, because they appear on both sides of the matching condition.
Most of the calculations done in the course of this work (e.g. determination of the matching coefficients, derivation of the Feynman rules, manipulations of the EFT Lagrangians etc.) were carried out not only by pen and paper but also using software tools for automatic calculations. For the latter we employed Mathematica packages FeynArts [110], FeynRules [111] and FeynCalc [112][113][114]. The automation of non-relativistic calculations was significantly simplified by making use of interface to QGRAF [116] diagram generator was added to the development version of the Feyn-Helpers [117] extension. This allowed us to generate Feynman diagrams for non-relativistic EFTs in a straightforward fashion. All new functions that were developed while working on this project should be made publicly available and properly documented in the upcoming versions of FeynCalc, FeynOnium and FeynHelpers.

Fermion bilinear and scalar sector
Let us discuss the matching coefficients of the bilinear fermion and antifermion sectors, first two lines in eq. (3.1). This amounts to comparing scattering amplitudes with one incoming and one outgoing fermion and scalar mediators (one, two or three of the latter field). The diagrammatic representation of the matching for a fermion interacting with one single scalar field is given in figure 2. In this work, we consider the matching of the NREFT Lagrangian at tree-level as far as the fermion (antifermion) bilinear is concerned. For the trilinear coupling this means that it is sufficient to work at order g. However, we remark that this procedure is general and applicable to the matching at any loop order. We collect some details in the appendix A, whereas here we list the results for the matching coefficients that read For c 3 , c 4 , c 3 and c 4 , we have considered the matching of diagrams with two and three external scalars respectively. Consistently with the findings from the FWT and EOM methods (cf. appendix A), these matching coefficients are found to vanish at tree-level. The matching coefficients c 1 , c D , c S may receive O(g 2 ), O(λ) corrections (not addressed in this work 10 ), whereas c 3 , c 4 may start getting non-trivial contributions at one-loop level. The Figure 3: Diagrammatic matching between the relativistic theory (diagrams on the l.h.s) and the corresponding four-particle local interactions in the NREFT (diagrams on the r.h.s). The latter correspond to the dimension-6 and dimension-8 operators of the four-fermion sector (3.1) respectively.
coefficients of the kinetic terms c 2 and c 2 are fixed to unity to all orders in perturbation theory owing to the reparametrization invariance.
There is an important aspect we want to highlight. As one may read off from eq. (4.1), there is a relative sign difference between the particle and antiparticle interactions with the scalar field. At order O(1/M 0 ) this is in contrast to the situation in NRQED and NRQCD, where the signs are the same. This very difference will be the reason behind the appearance of monopole and quadrupole interactions in the lower energy EFT that we will derive in section 5, instead of typical dipole interactions of pNRQED and pNRQCD.
As for the scalar sector described in third line of eq. (3.1), we equally perform the matching at tree-level only. Our guidance here is again the power counting of the pNRY that will be given in section 5. Postponing the one-loop matching of the NRY to a future work on the subject, one can simply obtain the matching coefficients at tree-level to be

Four-fermion operators and annihilation cross section
As anticipated, the NRY can readily describe heavy pair annihilations in terms of local 4-fermion operators in eq. (3.1). The inclusive annihilation rate can be recast in terms of an amplitude that conserves the number of the heavy particles by means of the optical theorem: the imaginary part of the loop amplitude with four external heavy fermion legs is related to the cross section of the process XX → φφ [62,104], cf. figure 3. For the model at hand, it is known that the annihilation cross section is velocity suppressed [79]. This will be reflected in a vanishing contribution from the velocity-(or derivative-) independent operators, that are of dimension-6. They read [62] ( 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 . Then, we consider dimension-8 operators, which comprise higher powers in 1/M . These are compensated by derivatives acting on the fermion (antifermion) fields, that induce velocity suppressed contributions due to ∇ψ ∼ M v (cf. section 3). As we shall see, they provide the leading contribution to the annihilation process XX → φφ. Of course, higher dimensional operators (further suppressed in the velocity expansion) are allowed as well but will not be considered in this work. The explicit structure of the dimension-8 operators in eq. (3.1) reads [62] ( where the operators explicitly included are where ↔ ∂ is the difference between the derivative acting on the spinor to the right and on the spinor to the left, namely χ † ↔ ∂ ψ ≡ χ † (∂ψ) − (∂χ) † ψ. The notation T (ij) for a rank 2 tensor stands for its traceless symmetric components, T (ij) = (T ij + T ji )/2 − T kk δ ij /3. As pointed out in [62], we may also have operators with the derivative acting on the product of the spinor fields ψ † and χ or χ † and ψ. The matrix elements of such operators are proportional to the total momentum of the pair XX, that is zero in the rest frame of the particle-antiparticle pair.
The detailed derivation of the matching coefficients can be found in appendix A, whereas here we merely list the results Both matching coefficients of the dimension-6 operators in eq. (4.3) are zero at the order we are working. Accordingly, they do not contribute to the pair-annihialtions. We observe that the annihilating fermions are always in the spin triplet configuration with the orbital angular momentum L = 1 and definite total angular momentum values J = 0, 2. The allowed combinations of L, S, and then J, are constrained by the symmetry of the fundamental Lagrangian eq. (2.1), that are also inherited by the NRY eq. (3.1). Since the scalar does not carry any spin, the conservation of parity and the total angular momentum forbids S = 0, and impose ∆L = 0 or even. We conclude this section by reproducing the non-relativistic annihilation cross section for the process XX → φφ. 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 used the non-relativistic flux factor with the relative velocity in the center of mass frame, so that v rel = |v ψ − v χ | = 2v. The imaginary part of the non-relativistic amplitude M NR can be readily computed using the Lagrangian from eq. (4.4) and the matching coefficients eq. (4.15), and upon setting v = v in eq. (A.52) (since we consider the scattering amplitude of ψχ → ψχ) we obtain The result agrees with the literature, cf. e.g. [55,57].
As was pointed out in [56], and can be inferred from the benchmark point used in [57], large values of α in this model can be of particular phenomenological interest. It is in the reach of the NRY, and subject of another work [118], to derive such higher order corrections in the matching coefficients of the bilinear and 4-fermion sector, and to inspect their impact on, e.g. the annihilation cross section.

Potential non-relativistic Yukawa theory
In the previous section, we integrated out the hard degrees of freedom with energies of order M as well as fermion/antifermion fluctuations of the same order. Here we want to integrate out soft degrees of freedom with energies of the order of the relative momentum of the pair p ∼ M v. 11 The corresponding effective theory takes the form of a pNREFT and we can rely on the techniques employed in the derivation of pNRQED as long as we assume that the mass of the scalar satisfies m M v. This condition implements a Coulomb-like regime and implies the scaling v ∼ α for the velocity. Moreover, we are allowed to treat the scalar mediator as effectively massless in the matching between the NRY and the pNRY, upon relying on the scale hierarchy M α M α 2 , m, irrespective of the relative size of the smaller scales. We comment on the case m ∼ M α later in this section.
Let us come to the construction of the pNRY Lagrangian. First of all, as the two-point functions are not sensitive to the relative momentum of the pair, the fermion bilinears of the NRY from eq. (3.1) and the pNRY will look the same. That said, one has to keep in mind that only scalar fields with ultrasoft momenta are kept in the latter EFT. Conversely, diagrams with four-fermion external legs are sensitive to the relative momentum and non-trivial contributions will be generated: they are the potential terms in the pNREFT Lagrangian [121,122]. The important point to be stressed here is that the appearance of the potential terms can be seen as the effect of integrating out soft scalars, and hence the potential can be extracted by matching the NRY to the pNRY.
In order to elucidate on the distinction between soft and ultrasoft scalars, and to introduce the degrees of freedom of pNRY, we project the NRY onto the particle-antiparticle sector as follows where i, j are spin indices, while the state |φ US contains no heavy particles/antiparticles and an arbitrary number of scalars with energies much smaller than M α. Here, ϕ ij (t, x 1 , x 2 ) is a wave function representing the XX system. After the projection it will be eventually promoted to a bi-local field. As a next step, one recognizes the relative distance of the pair r = x 1 − x 2 to have typical size of the inverse of the soft scale M α, or the inverse Bohr radius of a Coulombic state, cf. eq. (5.9). This can be considered a small scale as compared to the typical wave-length of the ultrasoft scalars, which is of the order of the inverse of M α 2 . According to the projection (5.1), scalar fields now appear at points x 1 and x 2 , and we can ensure that they are ultrasoft by expanding them about the center of mass coordinate R = (x 1 + x 2 )/2 upon using r R. One has to evaluate the leading order interaction between particle-antiparticle pairs and the ultrasoft scalar field, namely the combination g(φ(x 1 ) + φ(x 2 )), around the relative coordinate 11 A distinction between potential photons, i.e. photons with k 0 ∼ M α 2 and k ∼ M α from soft photons, i.e. photons with k 0 ∼ k ∼ M α, can also be considered [119,120] and this could apply to the scalar mediator as well. This distinction is not that relevant in the formulation we are following since, as done for the QED case, both potential and soft photons are integrated out at the same time when matching NRQED to pNRQED. up to O(r 2 ) as follows As a consequence, the dipole terms, namely the ones linear in r, cancel exactly. This is a peculiar feature of the Yukawa-type theory eq. (2.1) and its low-energy version eq. (5.4) which distinguishes them from pNRQED (and pNRQCD) where dipole transitions naturally arise. In summary, we write the Lagrangian in terms of the wave function field ϕ(r, R, t) and the ultrasoft scalars φ(R, t) as follows where the square brackets in the second line of (5.4) 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. (5.4) as well. To avoid cluttering the notation we suppress the spin indices of the bilocal fields that are contracted with each other. At variance with NRY (3.1), each term in the pNRY has a well-defined scaling: ∂ 0 ∼ M α 2 , the inverse relative distance and the corresponding derivative r −1 , ∇ r ∼ M α, whereas the scalar field and the center-of-mass derivative gφ, ∇ R ∼ M α 2 . The potential is understood as a matching coefficient and is organized as an expansion in α(M ) and λ(M ), as well as 1/M (inherited from NRY), the coupling α(1/r) and the relative distance r (as proper of pNRY). In the following, we parametrize it as V = V (0) + δV . In the case of m M α the leading order potential reads V (0) = −α(1/r)/r, which is a Coulomb potential. It is important to specify the scale at which α is evaluated, since it helps to keep track of the matching between the full theory from eq. (2.1) and the NRY from eq. (3.1), and between the latter with pNRY given in eq. (5.4). The corrections to the potential δV are needed to compute observables, such as the energy spectrum, at next-to-leading order. Here, we aim at extracting the potential at the first non-trivial order M α 4 , as an application of pNRY and the corresponding power counting. We discuss the potential matching in the following section 5.1. As for the kinetic terms, we included contributions up to order M α 4 as well.
Next, in the second line of eq. (5.4), we see the appearance of a monopole and a quadrupole interactions as well as interactions involving the derivative in the relative distance. Such structures are found by performing the so-called multipole expansion of the ultrasoft fields, here the scalar mediator, and terms up to order M α 4 are retained. The absence of dipole transitions in this model was already pointed out in [55,57] when dealing with the calculation of the bound-state formation. In our approach, the absence of such terms is already manifest at the level of the Lagrangian, where the degrees of freedom and their interactions are spelled out at the energy scale of interest for bound-state calculations. These are the wave function field of the particle-antiparticle pair and ultrasoft scalars, that interact via monopole and quadrupole interactions. The term with the derivative ∇ r in the second line of eq. (5.4) arises from the 1/M 2 spin-independent operator in eq. (3.1), and is of the same order as the quadrupole term, namely M α 4 . Indeed, the contributions of spin-dependent operators are subleading in the power counting. The presence of the ultrasoft scalar-derivative coupling conforms with the findings of [55]. 12

Potential matching
In this section we address the matching between NREFT and pNREFT. This procedure brings us to the systematic derivation of the particle-antiparticle potential, which can be understood as a matching coefficient of the pNRY. This procedure is rather general, and it has been generalized at finite temperature in [107,123] for pNRQED, as well as for pNRQCD [108]: any scale larger than M α 2 may contribute to the potential. For the moment, let us assume that the scalar mass is much smaller than the soft scale M α and can, therefore, be treated on the same footing with the ultrasoft scale M α 2 in the matching. At variance with on-shell Green's functions exploited in the matching between the relativistic theory eq. (2.1) and NRY eq. (3.1), here we enforce the equality between off-shell four-fermion Green's functions. The reason is that this is the typical situation of particles in a bound state or near-threshold. We give a diagrammatic representation of the potential matching in figure 4. Then, external momenta of the fermions are soft ∼ M α, whereas the energies are at the ultrasoft scale M α 2 , so that one can expand in such a scale. For the sake of the matching performed here, the ultrasoft scale can be put to zero. This amounts to simplifying the scalar propagator to −i/k 2 and taking the fermion propagators in loop diagrams of the NREFT side of the matching as static [121,124]. Loop diagrams on the pNREFT side vanish in dimensional regularization, because they are scaleless when expanded in the ultrasoft scale.
Since the NRY is organized as an expansion in α and 1/M , we can readily understand which are the terms we need in the potential V up to some order M α n [121,122]: one has to carry out the matching by combining inverse powers of the hard scale and the couplings as (1/M ) a α b with a + b ≤ n − 1. In this work we aim at potential corrections up to order M α 4 . In practice, for a given diagram, one counts the powers of α, inverse powers of M , and then multiplies by the soft scale M α in order to obtain the dimension of an energy. One has to consider tree-level, one- 12 The authors of [55] do not use EFT methods and hence provide no power-counting rules. The coupling gφ(R, t)∇ 2 R /M 2 has been considered in the non-relativistic Hamiltonian, which is however α 2 suppressed with respect to gφ(R, t)∇ 2 r /M 2 . Figure 4: Green's function in NREFT and pNREFT, where the momentum assignment is explicitly shown. The momentum transfer is k = p − q, where k is the momentum carried by the scalar mediator. The potential depends on p and k and the spin of the particle and antiparticle. and two-loop diagrams for the matching similarly to what has been carried out in the pNRQED case [122]. As explained above, in the diagrams involving heavy fermions in loops, we employ static fermion propagators, which allows us to simplify the derivation [121,124]. We address here the tree-level diagram only, whereas a more detailed discussion on the loop diagrams is deferred to appendix B.1. We collect the relevant tree-level diagrams in figure 5. In the upper row, the leftmost diagram provides the leading contribution to the potential of order M α 2 . It is easy to see that the middle and rightmost diagrams give instead a contribution of order M α 4 : the two inverse powers of M in the vertex are compensated for the soft scale M α, namely α × 1/M 2 × (M α) 3 . A quantum correction to the matching coefficients c D and c S would give a contribution of order M α 5 . Next, let us discuss the diagrams in the lower row of figure 5, as an example of diagrams which are suppressed in the power counting. The leftmost diagram accounts for the insertion of the corrected scalar propagator, and the contribution scales as α × d 4 /M 2 × (M α) 3 ∼ d 4 M α 4 . Since the matching coefficient is d 4 = O(α, λ), i.e. it vanishes at tree-level, this diagram is beyond our desired accuracy. By applying similar power counting arguments, one sees that the 4-fermion dimension-6 operators would contribute at order f ( 1 S 0 )M α 3 and f ( 1 S 0 )M α 3 . However, as shown in section A.4, such matching coefficients vanish at O(α) and, therefore, they could contribute only beyond the required accuracy. Finally, dimension-8 operators need not be considered, as they are further suppressed. We find that one-and two-loop diagrams do not contribute at order M α 4 , but only beyond (see discussion in appendix B.1). Then, the potential matching receives contributions only from the upper diagrams in figure 5, and the corresponding ones with the vertices c D and c S on the antiparticle line (not displayed). The resulting potential reads where S = diag(σ 1 , σ 2 )/2 is the total spin matrix, with σ 1 (σ 2 ) being the spin matrix of the two-component particle (antiparticle) field in the fermion bilinear, and we used c D = −c D and c S = −c S . As expected, the first term is the same as in pNRQED. This finding is consistent with the assumption m M α and the contribution corresponds to a Coulomb potential upon performing the Fourier transform. The third term (up to a relative sign difference) equally appears in the pNRQED potential. Yet the second term is characteristic of the pNRY and differs from what one finds in pNRQED. This difference can be traced back to the different momentum combinations entering the non-relativistic Lagrangian in eq. (3.1) as compared to the NRQED case. Upon performing the Fourier transform (cf. e.g. [125] for a collection of useful formulas), we obtain the following potential in the position space with p = −i∇ r , L = r × p .
Solving this equation yields Coulombic energy levels E n and the Bohr radius a 0 given by Let us also remark that since eq. (5.8) describes a fermion-antifermion bound state, it features a factor of 2 in front of the −gφ term and a potential V (0) as compared to eq. (3.4) for a single fermion. Before moving to the applications of pNRY, one more comment is in order. It is well known that the Yukawa potential induced by a scalar mediator is universally attractive so that not only particle-antiparticle but also identical fermions can form bound states. This is very different from e.g. QED, where e + e + -or e − e − -interactions are repulsive. We have explicitly verified that the pNRY for particle-particle (antiparticle-antiparticle) pair admits the very same form as eq. (5.4) upon performing the replacements χ j (t, x 1 )). However, since identical Dirac fermions cannot annihilate into scalars, the bound-states XX andXX are completely stable in the context of the scalar Yukawa theory.

Scalar mass of order M v
The Yukawa potential is usually understood as a screened potential of the form −αe −mr /r. The mass of the mediator leads to a finite-range interaction with r ∼ 1/m, as opposed to the Coulomb case. Our calculation of the potential in eq. (5.5) assumed the scalar mass to be much smaller than the momentum transfer of order M v (we restore v instead of α in order to be generic and not necessarily in the Coulomb regime v ∼ α). Hence, we consistently neglected the mediator mass in the matching and it did not appear in the corresponding potential in eq. (5.6).
Let us briefly discuss how the matching calculation changes when m ∼ M v, so that the scale m cannot be neglected. The main difference resides in the scalar propagator that enters the where one recognizes the leading term to be a Yukawa screened potential. Moreover, in this case the pole of the bound-state propagator receives a finite mass shift [108,126]. 13 The corresponding one-loop diagrams are shown in figure 6, and we find Notice that eq. (5.10) reduces to eq. (5.6) in the limit m → 0. On the other hand, one can also study corrections to the Coulombic regime by expanding eq. (5.10) in mr 1. Let us also remark that in the case of a vanishing scalar mass, the loop integrals in figure 6 are scaleless, and vanish accordingly in dimensional regularization.

Bound-state spectrum at order M α 4
As a first non-trivial application of the pNRY, we carry out the derivation of the discrete spectrum at order M α 4 . We follow the setting of [121,122] put forward for pNRQED and consider the two-point function of the field ϕ. The potential and kinetic contributions are simple insertions into the ϕ propagator (see figure 7 leftmost and middle diagrams). This brings us to the evaluation of quantum mechanical expectation values when projected onto bound-state wave functions. In the limit m M α, we can approximate the state as Coulombic and compute the 13 Despite of the fact that these references discuss finite temperature calculations, the zero-temperature case follows the same pattern. The mediator mass m is integrated out together with the inverse distance between the pair 1/r ∼ M v. The potential and mass shifts are understood as matching coefficients of the pNREFT.
expectation values on such unperturbed states. Additionally, one has to consider the ultrasoft contributions to the binding energy, namely those originating from the loop corrections to the propagator of an ultrasoft scalar. Similarly to the case of pNRQED, the leading non-vanishing ultrasoft contributions arise from one-loop self-energy diagrams. We assume that the scalar mass can be as large as the binding energy M α 2 . Applying the power counting of the pNRY, we see that only the diagram with two monopole vertices contributes within our accuracy, displayed in figure 7, and it scales as M α 3 . The ultrasoft contribution is finite, and some details of the calculation are provided in appendix B.1. If the scalar mass is much smaller than the energy scale M α 2 , the scalar propagator can be expanded in m/M α 2 and the corresponding loop integral vanishes in dimensional regularization. Accordingly, the one-loop self-energy diagram yields no contribution to the spectrum.
Owing to the presence of the spin-orbit coupling in the Hamiltonian, the latter does not commute with L and S separately. However, the combination L · S can be rewritten in terms of the squared operators J 2 , L 2 and S 2 (and the Hamiltonian commutes with all of them), see e.g. [127]. In this case it is common to label the states with |n j , since n, and j are good quantum numbers. The corrections to the spectrum from the ultrasoft exchange δE US , kinetic energy δE kin , and the potential terms δE δV read , (5.12) and (5.14) A comment is in order on the form of the ultrasoft contribution δE US . One can check its appearance in a complementary way. It has been already stated that this term corresponds to the propagation of an ultrasoft scalar with momentum/energy of order M α 2 . Assuming its mass to be of the same order, as done here, this amounts to expanding the potential in eq. (5.10) for mr 1. The contribution at order mα from the potential expansion adds up with the one from δM in eq. (5.11), giving δE US , whereas the contribution at order αm 3 cancels against the one in δM . In the case of a massive vector boson as a force mediator, the monopole contribution from the potential completely cancels against the mass correction (see e.g. [108] for the QCD case) so that there is no analogue of δE US at order mα.

Bound-state formation cross section
Let us come to an application of the pNRY that establishes a connection to the recent developments in DM phenomenology. As noted in the original works [18,26], the formation of unstable bound states can trigger another channel for DM annihilations, and consequently affect the estimates of the present-day relic density. In general, bound-state formation and decay are not only relevant for the early universe but can also provide enhanced signals in the annihilations of DM in the galactic halos and affect the corresponding experimental signatures.
Here we deal with the bound-state formation via a radiative transition, where an abovethreshold scattering state emits a scalar particle and turns into a bound state. In terms of the model degrees of freedom one has (XX) open → (XX) bound + φ. This process occurs at the ultrasoft scale, so that the energy difference between the initial and the final state is of order M α 2 . Such interactions can be naturally accommodated in our pNRY, where the field ϕ accounts for both scattering states (with positive energies) and bound states (with negative energies). Formally, one can think of it as of splitting ϕ ≡ ϕ s + ϕ b in the particle-antiparticle Fock space.
As it was done for the hard annihilations into scalars within NRY, here we again make use of the optical theorem. The process of interest can be calculated from the self-energy of the pair in a scattering state by extracting its imaginary part. We show example diagrams in figure 8. Loop diagrams involve scales that are still dynamical in pNRY, namely the energy scale M α 2 and the mass of the scalar m, which we have assumed to be much smaller than M α for the derivation of the pNREFT. No specific relation between their relative importance was needed in the matching between NRY eq. (3.1) and pNRY eq. (5.4). However, at this stage it is important to clarify their relative size. In the following derivation, we consider the case m < ∼ M α 2 , so that the m must be retained in the scalar propagator. The method that we describe in the following is suitable for both in-vacuum and finite temperature calculations, provided that the thermal scales are smaller then the typical relative momentum of the pair. In this case, one can take the pNREFT Lagrangian from eq. (5.4) as a starting point [107,108,123,128] and incorporate T = 0 that would enter as a dynamical scale together with the in-vacuum parameters m, M α 2 . 14 We leave a comprehensive construction of EFTs for scalar mediators at finite temperature for future work on the subject. Let us come to the self-energy diagrams relevant for the derivation of the cross section. We show the ones that involve two monopole or two quadrupole vertices in figure 8. The calculation is done in dimensional regularization with D = 4 − 2 . Let us start with the left diagram in figure 8, where the corresponding self-energy reads with E φ = √ k 2 + m 2 being the energy of the scalar mediator. In order for the self-energy to acquire a full meaning, it has to be projected on the external scattering states, labeled with the relative momentum quantum number p = M v rel /2, so that P 0 = E p = p 2 /M . Next, we also insert a complete set of bound states so that the internal propagator in the loop diagram describes indeed the propagation of the discrete states of the spectrum. More explicitly, One can easily extract the imaginary part of the self-energy using Cutkosky cutting rules at zero temperature [129] that impose the kinematic condition 0 < k 0 < ∆E p n . The result for the inclusive cross section to produce all possible bound states reads which has the correct dimension of inverse energy squared. 15 However, one can readily see that the cross section in eq. (5.18) vanishes because of the orthogonality between the scattering and bound-state wave functions Ψ p (r) and Ψ n (r), that appear in the expectation value p|n = d 3 r Ψ * p (r) Ψ n (r) = 0. For the same reason mixed monopole-quadrupole diagrams give no contribution to the total cross section. These findings nicely agree with the results of [55,57], where the authors consider interactions of Dirac fermion DM with a scalar mediator. The same pattern is observed also in the case of DM being a non-relativistic scalar particle coupled to a scalar mediator [56,130].
14 The situation is of course different if the temperature and other thermal scales, for example thermal masses, are of comparable size or larger than M α. Then the pNREFT has to be derived accordingly and it would be qualitatively different from eq. (5.4). See [107,123] for QED and [108,128] for QCD. 15 One can simply see this by recalling the energy dimension of bound and scattering states kets, given by  Let us now consider the quadrupole contributions induced by the right diagram in figure 8. The corresponding self-energy reads where one may notice the appearance of powers of the scalar three momenta in eq. (5.19), that are induced by the action of the derivative operator ∇ R on the scalar propagator. This diagram can be evaluated in the same fashion as the monopole contribution. The result for the cross section reads which also has the correct mass dimension and features non-trivial expectation values. The form of the prefactor highlights the effect of the mediator mass being of the same order as the ultrasoft scale: this setting obviously features a strong suppression of the formation rate. On the other hand, in the case of m M α 2 one can simply expand the result in eq. (5.20) accordingly. The total cross section also receives contributions from the spin-independent relativistic correction operator −gφ∇ 2 r /M 2 in eq. (5.4). This amounts to additional diagrams involving two insertions of this operator, as well as diagrams with one such insertion and a quadrupole vertex, whereas the insertion of a monopole vertex yields a vanishing amplitude. Putting everything together, our final result for the total cross section reads When considering the ground-state formation cross section, namely the state |n = |100 , and neglecting the mediator mass, we obtain the following result using Coulomb bound state and scattering wave functions Our result in eq. (5.22) agrees with the earlier findings in the literature [55,57] in the limit ζ ≡ α/v rel 1 . The first and the second term in eq. (5.22) stem from the = 0 and = 2 partial waves in the scattering wave function respectively. 16 Their relative size as compared to the total cross section can be inferred from figure 9, where they are depicted as dashed orange and dot-dashed red curves. The brown curves correspond to the total cross section, the solid line is our result eq. (5.22), whereas the dotted one is the large ζ limit [55,57].
The recasting of the bound-state formation cross section in the language of pNRY offers a clear organization in terms of quantum mechanical expectation values. In particular, as for the ground-state formation, the = 2 contribution only comes from the p|r i r j |n matrix element (it develops a non-trivial angular dependence), whereas all the four matrix elements in eq. (5.21) contribute to the = 0 term.

Conclusions
Self-interacting dark matter is mostly welcome in the attempt of reproducing the observed galactic structures, and it appears to work better than collisionless dark matter. Typically, self-interactions between non-relativistic dark matter particles are induced by the exchange of a light mediator. In addition to the desired velocity-dependent interactions that accommodate the halos of different-sized objects, it may well be that such self-interactions induce DM bound states. Most notably, depending on the model at hand in terms of its field content, masses and couplings, the impact of bound-state formation can play a rather important role in the 16 As done in ref. [57], by choosing the coordinate system in such a way that p points into the z-direction, the scattering wave function can be expanded into partial waves as Ψp(r) = ∞ =0 Ψ p (r) with Ψ p (r) = r|p .
determination of the present-day DM energy density. This may lead to sizable changes in the parameter space compatible with the cosmological abundance, making it necessary to revisit the relevant experimental bounds. In this work we studied a model that represents a family of minimal DM models, where a light scalar mediator induces self-interactions between Dirac fermion DM via a Yukawa-type interaction. Making use of the assumed hierarchy of well separated dynamical scales M M v M v 2 , we employed EFT techniques to study the resulting bound-state dynamics. In particular, we carried out a rigorous derivation of NRY and pNRY for a scalar force carrier, in the spirit of NRQED and pNRQED and their QCD counterparts. These EFTs are known to be very useful and successful tools for investigating and calculating observables relevant to bound and near-threshold states. As for NRY we extended and generalized the formulation already available in the literature [68,69], whereas, to the best of our knowledge, the explicit construction of pNRY was carried out in this paper for the first time.
We started with the derivation of the NRY from the first principles, where we identified the relevant degrees of freedom (non-relativistic Pauli fields and a scalar) and worked out the powercounting. We explicitly included 1/M 2 -operators in the bilinear sector and 1/M 4 -operators in the four-fermion sector, which allowed us to reproduce the first non-trivial contribution to the annihilation cross section XX → φφ at leading order in the fermion-scalar coupling. In the bilinear sector the matching was performed at tree-level.
Then, we resolved the power-counting ambiguity in the NRY due to the soft and ultrasoft scales being still intertwined, by constructing the corresponding pNRY. The degrees of freedom of this low-energy theory were found to be particle-antiparticle pairs (represented by a bilocal field) interacting with ultrasoft scalars. The scalars were enforced to be ultrasoft by performing a multipole expansion in the relative center of mass coordinate r. This way the presence of characteristic monopole and quadrupole interactions at the level of the pNRY Lagrangian was made manifest. The same is true also for the arising spin-independent relativistic correction with derivatives acting on the heavy-pair field. Dipole interactions that typically appear in pNRQED, turned out to be absent in pNRY.
We explicitly computed the DM fermion-antifermion potential, which naturally arises at the level of the pNRY Lagrangian as a matching coefficient. This paves way for a systematic inclusion of quantum and relativistic corrections in future works on the topic. In the Coulombic m M α regime of pNRY the scalar-induced potential turned out to share some similarities with the one in pNRQED. However, we also found that the spin-orbit term comes with an overall opposite sing as compared to the pNRQED case, while the contribution induced by an operator proportional to r · (−i∇ r ) has no correspondence in the electromagnetic potential. We also performed the potential matching for the setting where the scalar mass is of order M v, thus recovering a Yukawa screened potential at leading order. Furthermore, we explained that pNRY can describe not only particle-antiparticle interactions but also bound states formed by identical Dirac fermions such as XX andXX.
As a first application of the pNRY, we computed the bound-state spectrum at the nextto-leading order, namely O(M α 4 ), which constitutes a new result presented in this work. In particular, we stressed the advantages of using a pNREFT for such bound state calculations as compared to non-EFT approaches. Our calculation was done in the Coulombic approximation, namely m M α, that still allows for the mediator mass to be as large as the ultrasoft scale M α 2 . The ultrasoft contribution to the spectrum, in the case of the scalar mass being not much smaller than the binding energy, was found to provide the leading contribution of order M α 3 .
A further application of the pNRY that was presented in this paper is the derivation of the bound-state formation cross section by taking the imaginary part of the heavy-pair field selfenergy diagram. In particular, the contributions from monopole-induced diagrams turned out to be vanishing due to the orthogonality of the wave function of the discrete and continuous spectrum, in full agreement with the previous findings in the literature. On the other hand, we identified the leading contribution to the cross section to be induced by quadrupole interactions and relativistic corrections. The final expression was written in terms of quantum mechanical expectation values that naturally arise in pNREFT calculations. We also performed an explicit analytic evaluation of these quantities for the Coulombic regime, thus agreeing with the earlier findings [55,57] in the limit of ζ 1. The fact that we were able to obtain the previously unknown full analytic result for arbitrary ζ-values using pNRY can be regarded as another highlight of this work.
To conclude, we would also like to provide a brief overview of future research directions in this field in conjunction with NRY and pNRY. The minimal model addressed in our work can be varied in different ways, such as a Majorana DM rather than Dirac DM, or a more general interaction with a pseudo-scalar force carrier and cubic self-interactions of the scalar/pseudo-scalar mediator. These equally compelling realizations can be handled within the EFT approach presented here. Moreover, an accurate derivation of the relic density requires calculation of various processes (e.g. bound-state formation, dissociation and Sommerfeld enhancement) to be done at finite temperature. We believe that the EFTs presented in this work can be regarded as a starting point for such finite temperature calculations, as it was the case with the corresponding generalizations of NRQED/NRQCD and pNRQED/pNRQCD. Especially in the heavy-ion phenomenology related to heavy quarkonia, pNRQCD has proven to be an extremely useful tool to scrutinize different hierarchies between in-vacuum and thermal scales, and to calculate relevant observables in a controlled and systematic way. It goes without saying that the presence of thermodynamical scales can significantly modify and affect the relevant cross sections also in the DM phenomenology. Therefore, the derivation of the finite temperature versions of NRY/pNRY constitutes a worthwhile and phenomenologically relevant task that we hope to address in the subsequent publications.

Acknowledgments
The work of S.B. is supported by the Swiss National Science Foundation under the Ambizione grant PZ00P2 185783. V. S. acknowledges the support from the DFG under grant 396021762 -TRR 257 "Particle Physics Phenomenology after the Higgs Discovery." The authors are grateful to Jacopo Ghiglieri for reading the manuscript and providing useful comments, and to Miguel Escobedo for stimulating discussions at early stages of our work. They would also like to thank Matthias Steinhauser for making them aware of [73,74] and Joan Soto for [88].

A Matching coefficients of NRY
In this appendix, we provide a detailed derivation of the matching coefficients that enter the NRY Lagrangian eq. (3.1). As for the bilinear fermion (antifermion) sector we work at leading order and we discuss the derivation in section A.1. Then, in section A.2 we derive the NRY by using the equation of motion method that allows us to (i) do a non-trivial check of the soobtained matching coefficients at tree-level; (ii) write the NRY in a covariant fashion, implement the reparametrization invariance, here at order 1/M , and consequently fix c 2 = 1 at all orders. Finally, in section A.4 we provide the derivation of the matching coefficients of the dimension-6 and dimension-8 operators.

A.1 Matching of the fermion bilinear with scattering amplitude
The derivation of the matching coefficients for the fermion bilinear at order 1/M 2 can be conducted in the following way. We write down the scattering amplitude for the process ψ(p) → ψ(q) + φ(k) in the fundamental theory eq. (2.1) and expand the resulting expression in powers of p/M , q/M . To this aim, we need to rewrite Dirac spinors in terms of the twocomponent Pauli spinors using non-relativistic normalization [62] with E p = p 2 + M 2 . Furthermore, we take the γ matrices in the Dirac basis and decompose them into Pauli matrices. By momentum conservation at the vertex we have k = p − q, which is the momentum carried by the scalar, so that the result for the diagram in figure 2 (left) reads for the particle interaction with φ, whereas for the antiparticle. Upon identifying the 2-spinors ξ and η with the ψ and χ fields respectively, we can compare the so-obtained expressions to the amplitudes induced by NRY eq. (3.1) and read off the values of the matching coefficients listed in eq. (4.1). It is interesting to remark that the Pauli structures appearing in eqs. (A.2) and eq. (A.3) differ by an overall minus sign, which is different from the situation that one finds in NRQED (no sign difference). This can be traced back to the vector structure of the electromagnetic interaction that features an additional γ matrix in the fermionic current, so that the couplings of the particle and the antiparticle to the temporal component of the photon field A 0 (k) have the same sign.

A.2 Equations of motions method
This method exploits the equations of motion of the high-and low-energy excitations of the relativistic field X [98][99][100][101][102] to derive the corresponding nonrelativistic EFT. The following derivation closely follows [102], where the same exercise is carried out for QCD. In practice, one starts from the full relativistic Lagrangian, where are we interested solely in the fermion-bilinear piece given by We then decompose the relativistic four-component spinor X as where (1 ± / v)/2 are velocity-dependent projectors with / v ≡ v µ γ µ . In the rest frame of the pair, where v µ = (1, 0), such operators reduce to (1 ± γ 0 )/2 and project onto the particle and antiparticle components of the Dirac field X.
Making use of the properties of the velocity projectors, we find where it is now clear that H v comprises the large energy modes of order M that we want to integrate out from the Lagrangian. To this aim, it is useful to introduce the derivative ∂ µ ⊥ = ∂ µ − v · ∂v µ and to rewrite the equation of motion for the field H v Substituting the expression for H v as given in (A.7) into (A.6) we find which is still exact. Now we can directly expand L heavy in 1/M and, up to order 1/M 2 , the Lagrangian reads In order to eliminate all terms containing (v · ∂)h v beyond O(1/M 0 ), we need to introduce a suitable field redefinition [124] given by Let us stress that the Lagrangian in eq. (A.11) may describe not only non-relativistic systems made of heavy Dirac fermions of the same mass, but also bound states formed out of a heavy and a light fermion, which might be another interesting DM scenario worth exploring in more details using our EFT framework. This statement is completely analogous to the well-known fact [124] that the HQET Lagrangian is equally suitable for studying properties of heavy-light mesons and heavy quarkonia: Both theories share the same Lagrangian but differ in their power-counting.
To complete our derivation for the case of the NRY, we switch to the rest frame with v µ = (1, 0), employ the relation σ i σ j = δ ij + i ijk σ k and identify h v with the particle component ψ of the X field. Thus, we obtain that agrees with the particle sector of eq.

A.3 Foldy-Wouthuysen-Tani method
The main idea behind the Foldy-Wouthuysen-Tani (FWT) [96,97] method is to introduce a sequence of unitary transformations that decouple the upper and lower components of the Dirac spinor order by order in 1/M . Consequently, in the non-relativistic limit the Dirac equation splits into two separate equations for Pauli fields describing particles and antiparticles respectively. The procedure of applying FWT transformations to QED can be found in various QFT textbooks (cf. e.g. [131][132][133] that we will partially follow here) and is often taught in advanced quantum mechanics courses. Therefore, we do not claim any originality for most of the material presented below. Once the technicalities behind the QED case are understood, it is a simple exercise to repeat the same procedure for the scalar Yukawa theory. The results for the pseudoscalar case can be found in [72]. First of all, let us introduce the concept of even and odd operators. Even operators are those that do not interchange upper and lower components of the Dirac spinor X, so that particles and antiparticles remain decoupled. Odd operators, on the contrary, are responsible for the mixing between particles and antiparticles. Schematically, we can writê whereÊ is an even andÔ is an odd operator. In the context of the Dirac Hamiltonian we have α i = γ 0 γ i and β = γ 0 , where the former is odd, while the latter is even. The Dirac spinor field X satisfies i∂ t X =ĤX, (A.14) which yields the familiar Dirac equation in the case of a non-interacting Hamiltonian.
For the sake of clarity, let us first discuss the generic case, without making an explicit reference to a particular theory. Our starting point for applying the FWT procedure is the unitary transformation X → X =Û X, (A.15) withÛ = e iŜ , whereŜ is some operator. Then the time evolution of the transformed field becomes i∂ t X = e iŜ (Ĥ − i∂ t )e −iŜ X . (A. 16) Here∂ t means that the partial derivative acts only on e −iŜ but not on X . Therefore, the transformed field satisfies i∂ t X =Ĥ X , where we used that∂ t is non-vanishing only when it acts on a time-dependent function.
In the case of an interacting theory (e. g. QED) it is usually not possible to choose anŜ such, that the upper and lower components of X decouple to all order in the 1/M expansion. Instead, one proceeds by starting with an ansatz that removes all odd terms at O(1/M 0 ) and then calculatesĤ to the desired order in 1/M , say 1

A.4 Matching of the dimension-6 and dimension-8 operators
Here we closely follow the tree-level matching between QCD and NRQCD in the 4-fermion sector described in [62]. We derive the contribution to the amplitude XX → XX in the center of mass reference frame. Hence, we take the incoming X andX to have momenta p and −p, whereas the outgoing X andX have momenta p and −p respectively. Due to energy conservation and same masses involved in the process for the DM states, one has |p| = |p | ≡ p. 17 The form of the non-relativistic Dirac spinors has been already given in eq. (A.1). The matching is performed by enforcing on-shell four-fermion Green's function in the full theory (2.1) and in the NRY (3.1). For completeness, let us briefly describe the matching of the four-fermion operators at order α. The relevant tree-level diagrams are shown in figure 10. It is clear that no imaginary part can arise at this order. Moreover, the diagram on the left can be precisely reproduced in the NRY because the scalar can carry a soft momentum (it is indeed the diagram appearing in the potential matching in figure 5). Only the diagram on the right contributes to the matching at this order, and provides a contribution to the real part of the matching coefficients. We find the only non-vanishing coefficient at order α to be f ( 3 P 0 ) = 3πα, while the matching coefficients of all dimension-6 operators vanish. Going to order α 2 , we find two one-loop diagrams that contribute to the process XX → φφ and we show them in figure 3. We are interested in their imaginary parts, which can be extracted by using the standard cutting rules, namely putting the internal scalar fields on shell. Expanding up to second order in the velocities v = p/E, v = p /E, and writing the final result as in [62], we find the annihilation contribution to the scattering amplitude to be where the subscript t + u indicates the sum of the t-and u-channel diagrams of figure 3. We remark that there is no term of order v 0 , implying that in the Yukawa theory (2.1) annihilations are velocity suppressed and start at order v 2 . The matching coefficients of the dimension-6 operators are zero, also at order α 2 (as expected by symmetry arguments). On the contrary, In this section we would like to discuss one-and two-loop diagrams that need to be analyzed for the potential matching of the pNRY. The systematic analysis is partly based on the pNRQED matching in the Feynman gauge [122], where the temporal component of the photon field has to be considered in loop diagrams (at variance with the Coulomb gauge). Then we have to consider (i) one loop diagrams as given in figure 11 (possible contribution at O(M α 3 )); (ii) the same diagrams with a kinetic insertion p 2 /2M in one of the fermion lines at a time, (possible contribution at O(M α 4 )); (iii) again the same diagrams with external energy insertions arising from the expansion of the propagators around zero external energy (possible contribution at O(M α 4 )); (iv) two-loop diagrams involving scalar propagators without kinetic/external energy insertion (possible contribution at O(M α 4 )). We checked explicitly that the same arguments put forward for the QED case holds here. The sum of the two diagrams in figure 11 indeed vanishes. Then, the same diagrams with a kinetic or an external energy insertion vanish individually, due to an odd number of the static propagators involved. The last set (iv) equally vanishes, since they are an iteration of the one-loop diagrams (i), as shown in [134]. In addition to the previous class of diagrams, we have to consider possible contributions arising from other topologies, namely those that are induced by the interactions between fermions (antifermions) and two or three scalars. Before discussing the diagrams in some detail, let us remind that the coefficients c 3 (c 3 ) and c 4 (c 4 ) vanish at tree-level, so that c 3 , c 4 = O(α), O(λ) at least. Actually, as far as c 3 (c 3 ) is concerned, one finds that the matching coefficients go like O(α 2 ), because the tree-level topology is reproduced in the NRY, and then there is no contribution to c 3 (c 3 ) at order α. The one-loop diagrams involving the vertices with two scalar fields are collected in figure 12 (upper row). They all contribute at order M α 5 or higher. Finally, two example diagrams with a three-scalar-fields vertex are given in figure 12 (lower row). By applying the power counting one sees that they all go beyond the accuracy of this work, M α 9/2 . All other diagrams involving c 4 and c 3 are further suppressed.

B.2 Master integrals
Here we provide explicit analytic results for some of the 1-loop integrals that we encountered in the course of calculations done in this work