\section{Introduction}

We consider the following initial value problem for an inviscid Burgers-Hilbert equation for $u(t,x)$:
\begin{align}
\label{bheq}
\begin{split}
&u_t+ u u_x = \mathbf{H}\left[ u\right] ,\\
&u(0,x) = u_0(x),
\end{split}
\end{align}
where $\mathbf{H}$ is the spatial Hilbert transform on $\mathbb{R}$, and
the initial data $u_0$ is sufficiently small
\begin{equation*}
\Vert u_{0}(x)\Vert _{H^2(\mathbb{R})}\leq \epsilon \ll 1.
\end{equation*}
Here, $H^k(\mathbb{R})$ denotes the standard Sobolev space of functions with $k$ weak $L^2$-derivatives.

Equation (\ref{bheq}) was proposed in \cite{biello} as a model equation for nonlinear waves with constant frequency.
It also provides a formal asymptotic approximation
for the small-amplitude motion of a planar vorticity
discontinuity located at $y = u(t,x)$ \cite{biello, marsden}.
Moreover, even though the equation is quadratically nonlinear, this approximation
is valid over cubically nonlinear time-scales.

Equation (\ref{bheq}) is nondispersive, and solutions of the linearized equation oscillate in time but do not
exhibit any dispersive decay. Nevertheless,
in comparison with the inviscid Burgers equation, small smooth solutions of (\ref{bheq}) have an enhanced, cubically nonlinear lifespan.
This enhanced lifespan was observed numerically in \cite{biello} and proved in \cite{hunter} by use of a change of the independent variable, which was suggested by a normal form transformation of the equation.

In this paper we give a different, and simpler, proof of the enhanced lifespan of smooth
solutions of (\ref{bheq}) from the one in \cite{hunter}.  Our main result is the following theorem:

\begin{theorem}
\label{th:main}
$a)$ Let $k \ge 2 $. Suppose that the initial data $u_0 \in H^k(\mathbb{R})$ for the equation \eqref{bheq} satisfies $$ \Vert u_{0}\Vert _{H^2(\mathbb{R})}\leq \epsilon \ll 1 .$$
Then  there exists a solution $
u \in C\left( I^\epsilon ;{H}^{k}\left( \mathbb{R}\right) \right) \cap C^1 \left(I^\epsilon;{H}^{k-1}\left(\mathbb{R}\right) \right)$
of (\ref{bheq})
defined on the time-interval $I^\epsilon=\left[-{\alpha}/{\epsilon^2},{\alpha}/{\epsilon^2}\right]$, where $\alpha >0$ is a universal constant, so that
\begin{equation*}
\Vert u\Vert _{L^{\infty}(I^{\epsilon}; H^{s}(\mathbb{R}))}\lesssim \Vert u(0)\Vert _{ H^s(\mathbb{R})}, \qquad 0\leq s\leq k.
\end{equation*}
\\
$b)$ Consider $u_1, u_2$ two solutions of \eqref{bheq} as in part $(a)$. Then
\begin{equation*}
\Vert u_1-u_{2}\Vert _{L^{\infty}(I^{\epsilon}; L^2(\mathbb{R}))}\lesssim \Vert u_1(0)-u_{2}(0)\Vert _{ L^2(\mathbb{R})}.
\end{equation*}
\end{theorem}
The quadratically nonlinear terms in (\ref{bheq}) are non-resonant, and a standard approach for problems of this type  is to use a normal form transformation (see for instance \cite{shatah}) to remove them. However, this is a quasi-linear problem, so the normal form transformation is unbounded, and this method cannot be applied directly. Our alternative approach is to use the normal form transformation in order to derive a modified energy functional for the original problem \eqref{bheq} which evolves according to a cubic law.

The same modified energy method also applies to the linearization of the Burgers-Hilbert equation, and allows us to control differences of solutions on a cubic time-scale.

We remark that with some additional work one can replace the $H^2$ norm for the data in part $(a)$ of the above theorem with $H^{s_0}$ for any $s_0> 3/2$. Similarly, the $L^2$ norm in part $(b)$ can be replaced by $H^s$ for a smaller range $0 \leq s \leq s_0 -1$.

\section{The Normal Form Transformation}

Straightforward energy estimates applied to (\ref{bheq}) yield
\[
\frac{d}{dt} \| \partial_x^k u\|_{L^2}^2 \lesssim \|u_x\|_{L^\infty}
\|u\|_{H^k}^2
\]
for $k\ge 2$, which gives a lifespan of smooth solutions of the order $\epsilon^{-1}$,
as in the standard local existence theory for quasi-linear hyperbolic PDEs \cite{kato, majda}.

The quadratically nonlinear terms in (\ref{bheq}) can be  removed by a normal form transformation \cite{biello,hunter}
\begin{equation}
\label{nft}
v = u + \mathbf{H}[\mathbf{H}u\cdot \mathbf{H}u_x ].
\end{equation}
The transformed equation has the form
\begin{equation}
v_{t}+ \mathcal{Q}(u)=\mathbf{H}\left[v \right]
\label{veq}
\end{equation}
where $\mathcal{Q}(u)$ is cubic in $u$ but involves two spatial derivatives.
Straightforward energy estimates for (\ref{veq}) yield
\begin{equation*}
\frac{d}{dt} \| \partial_x^k v\|_{L^2}^2 \lesssim \|u_x\|^2_{L^\infty}
\|u\|_{H^{k+1}}^2.
\end{equation*}
However, since the above normal form transformation is not invertible,
the energy estimates for $v$ do not close.
Thus, the $v$-equation is cubically nonlinear, but there is a loss of derivatives. On the other hand,
the $u$-equation is quadratically nonlinear, but there is no loss of derivatives.

The reason for this disparity in energy estimates is that the $H^k$-norms of $u$ and $v$ are not comparable.
From (\ref{nft}), we have
\begin{equation}
\|\partial_x^k v\|_{L^2}^2 = \|\partial_x^k u\|_{L^2}^2 + 2\langle \partial_x^k u, \partial_x^k \mathbf{H}[\mathbf{H}u \cdot \mathbf{H}u_x]\rangle
+ \|\partial_x^k\mathbf{H}[\mathbf{H}u \cdot \mathbf{H}u_x]\|_{L^2}^2.
\label{vkenergy}
\end{equation}
The second term on the right-hand side of (\ref{vkenergy}) is comparable to the $H^k$-norm of $u$ because
\[
\langle \partial_x^k u, \partial_x^k \mathbf{H}[\mathbf{H}u \cdot \mathbf{H}u_x]\rangle
= \left(k+\frac{1}{2}\right) \langle \mathbf{H} u_x, (\partial_x^k \mathbf{H} u)^2\rangle
+ l.o.t.
\]
where $l.o.t.$ stands for lower order terms. The third term on the right-hand side of (\ref{vkenergy}) is not comparable to the $H^k$-norm, but it is quartic and therefore irrelevant to the cubically nonlinear energy estimates.

This discussion suggests that we define a modified energy functional by dropping the higher-derivative quartic term
from the right-hand side of (\ref{vkenergy}) to get
\begin{equation}
\label{henergy}
E_k(u) = \frac{1}{2} \|\partial_x^k u\|_{L^2}^2 +  \langle \partial_x^k u,\partial_x^k
 \mathbf{H}[\mathbf{H}u \cdot \mathbf{H}u_x]\rangle.
\end{equation}
As we will show, this modified energy is equivalent to the standard $H^k$-energy and satisfies
cubically nonlinear estimates without a loss of derivatives.

\section{The Modified Energy Method}

In this section, we use the modified energy functional introduced in \eqref{henergy} to prove Theorem~\ref{th:main}. We first show that $E_k$ energy is equivalent to the standard $H^k$ energy.

\begin{lemma}
\label{lemma:lemaenergy}
Let $E_k$ be defined by (\ref{henergy}). Then
\begin{equation*}
E_k(u) = \frac{1}{2} \|\partial_x^k u\|_{L^2}^2 (1+ O(\| \mathbf{H}u_x\|_{L^\infty})).
\end{equation*}
\end{lemma}

\begin{proof}
Denote
\begin{equation*}
T_u f := \mathbf{H}\left[ \mathbf{H}u \cdot \mathbf{H}f_x\right] .
\end{equation*}
This operator is essentially skew-adjoint. Moreover, we have
\begin{equation*}
\begin{aligned}
\langle f, T_u f \rangle &= - \langle \mathbf{H}f, \mathbf{H}u \cdot \mathbf{H}f_x \rangle \\
&= \frac{1}{2} \int_{\mathbb{R}} \mathbf{H}u_x \vert \mathbf{H}f\vert ^2 dx = O(\| \mathbf{H}u_x\|_{L^\infty}) \|f\|_{L^2}^2.
\end{aligned}
\end{equation*}

We write
\begin{equation*}
\partial_x^k  \mathbf{H}\left[\mathbf{ H}u \cdot \mathbf{H} u_x\right] = T_u \partial_x^k u + l.o.t.
\end{equation*}
where $l.o.t.$ stands for terms of the form $\mathbf{H} [ \partial_x^j\mathbf{ H}u_x \cdot \partial_x^{k-1-j}\mathbf{H} u_x ]$
with $0 \leq j \leq k-1$. For the leading term we use the previous estimate for $T_u$, and for the lower order terms we use
interpolation,
\[
\| \partial_x^j\mathbf{ H}u_x \cdot \partial_x^{k-1-j}\mathbf{H} u_x\|_{L^2} \lesssim \|\mathbf{H}u_x\|_{L^\infty} \|\partial_x^{k-1} u_x\|_{L^2}
\]
Thus, we obtain
\begin{equation*}
\langle \partial_x^k u,\partial_x^k
 \mathbf{H}\left[ \mathbf{H}u \cdot \mathbf{H} u_x\right] \rangle =  O(\| \mathbf{H}u_x\|_{L^\infty}) \|\partial_x^k u\|_{L^2}
\end{equation*}
\end{proof}

The energy estimate is discussed in the following lemma:
\begin{lemma}
For the energy (\ref{henergy}) the estimate below holds: for any $\delta>0$,
\begin{equation*}
\frac{d}{dt} E_k(u) \leq  C_{k,\delta} \Vert u_x\Vert _{H^{\frac12+\delta}}^2\|u\|_{H^k}^2
\end{equation*}
where $C_{k,\delta}$ is a constant depending on $k$ and $\delta$ only.
\end{lemma}

\begin{proof}
A direct computation yields
\begin{equation*}
\begin{aligned}
-\frac{d}{dt} E_k(u) = &\int_{\mathbb{R}}  \partial_x^k (u u_x)  \partial_x^k
 \mathbf{H}\left[ \mathbf{H}u \cdot \mathbf{H}u_x\right] \, dx  \\
 +&\int_{\mathbb{R}} \partial_x^k u \left\lbrace  \partial_x^k
 \mathbf{H}\left[ \mathbf{H}[u u_x]   \cdot \mathbf{H} u_x\right]  +  \partial_x^k
 \mathbf{H}\left[ \mathbf{H}u  \cdot \mathbf{H} [u u_x]_x\right] \right\rbrace   dx .
 \end{aligned}
\end{equation*}
Note that all cubic terms cancel. This is because the energy was
determined using the normal form transformation, which at the level of
the energy amounts to eliminating the cubic terms.
 Using the antisymmetry of $H$ and integrating by parts
  we rewrite the above expression in the form
\[
-\frac{d}{dt} E_k(u) = \int_{\mathbb{R}}  - \partial_x^k \operatorname{\mathbf{H}} (u u_x)  \partial_x^k
 \left[ \mathbf{H}u \cdot \mathbf{H}u_x\right]  + \partial_x^k \operatorname{\mathbf{H}} u_x  \partial_x^k
 \left[ \mathbf{H}u \cdot \operatorname{\mathbf{H}} (u u_x) \right] dx.
\]
We observe that in each term in the expansion of the above integrand
has the form
\[
\operatorname{\mathbf{H}} \partial^{\alpha} u  \operatorname{\mathbf{H}} \partial^{\beta} u \operatorname{\mathbf{H}} (\partial^\gamma u\  \partial^\delta u)
\]
where
\[
\alpha + \beta + \gamma + \delta =  2k+2
\]
We observe that as long as
\[
1 \leq \alpha, \beta, \gamma, \delta \leq k, \qquad \alpha+\beta \geq 3,
\delta+ \gamma \geq 3
\]
such terms can be estimated directly in terms of the right hand side
of the energy relation in the Lemma. We will call such terms ``good''.

We return to the expression above, and observe that a cancellation occurs
if the second  $\partial_x^k$ operator in each of the two expressions applies
to the second factor. Using a full binomial expansion, we are left with the
expressions
\[
I_\alpha =  \int_{\mathbb{R}}  - \partial_x^k \operatorname{\mathbf{H}} (u u_x)  \cdot \partial_x^{\alpha}
  \mathbf{H}u \cdot  \partial_x^{k-\alpha}\mathbf{H}u_x  +
\partial_x^k \operatorname{\mathbf{H}} u_x \cdot
  \partial_x^{\alpha}\mathbf{H}u \cdot  \partial_x^{k-\alpha} \operatorname{\mathbf{H}} (u u_x)  dx
\]
where $1 \leq \alpha \leq k$.  Distributing the remaining derivatives,
we get good terms if any derivatives fall on the $u$ factor
(an additional integration by parts  is needed to see this for the second term).
Thus we have
\[
I_\alpha =  \int_{\mathbb{R}}  - \operatorname{\mathbf{H}} ( u  \partial_x^k u_x)  \cdot \partial_x^{\alpha}
  \mathbf{H}u \cdot  \partial_x^{k-\alpha}\mathbf{H}u_x  +
\partial_x^k \operatorname{\mathbf{H}} u_x \cdot
  \partial_x^{\alpha}\mathbf{H}u \cdot  \operatorname{\mathbf{H}} (u   \partial_x^{k-\alpha} u_x)  dx
+ good.
\]
If $u$ is commuted out the two expressions above cancel. Hence
\[
I_\alpha =  \int_{\mathbb{R}}  - [\operatorname{\mathbf{H}},u]  \partial_x^k u_x  \cdot \partial_x^{\alpha}
  \mathbf{H}u \cdot  \partial_x^{k-\alpha}\mathbf{H}u_x  +
\partial_x^k \operatorname{\mathbf{H}} u_x \cdot
  \partial_x^{\alpha}\mathbf{H}u \cdot  [\operatorname{\mathbf{H}} ,u ]  \partial_x^{k-\alpha} u_x  dx
+ good.
\]

The commutator $[\mathbf{H},u] f$ vanishes unless the frequency of $u$
is at least as large as the frequency of $f$.  Hence we can always
move derivatives from $f$ onto $u$.
Then for the first term in $I_{\alpha}$ we can use directly the commutator estimate below:
\begin{equation}\label{com}
\Vert \left[\mathbf{H},u \right]\partial_{x}^{k}u_x\Vert_{L^2}\lesssim\Vert  u_x \Vert _{L^{\infty}}\Vert \partial_{x}^ku  \Vert _{L^{2}}.
\end{equation}
This is a consequence of the commutator estimates result obtained by Dawson,
McGahagan, and Ponce in \cite{ponce}. We recall it below:
\begin{lemma} Let $\mathbf{H}$ denote the Hilbert transform. Then for
  any $p\in (1, \infty)$ and any $l, m\in \mathbb{Z}^+$ there exists
  $c=c(p,l,m) > 0$ such that
\begin{equation*}
\Vert \partial ^{l}_{x}[\mathbf{H}, a]\partial_{x}^{m}f\Vert_{L^{p}}\leq c\Vert \partial_{x}^{l+m}a\Vert _{L^{\infty}}\Vert f \Vert_{L^{p}}.
\end{equation*}
\end{lemma}

For the second in $I_{\alpha}$ we need to integrate once by parts since the first factor
has too many derivatives. If $\alpha < k$ this is done directly. If $\alpha = k$
then we move a derivative on the third factor. We are left with two terms which are nontrivial, namely
\[
J_k := \int |\partial_x^k \operatorname{\mathbf{H}} u|^2  [\operatorname{\mathbf{H}} ,u_x] u_x dx
\]
and
\[
J_{k-1} := \int |\partial_x^k \operatorname{\mathbf{H}} u|^2  [\operatorname{\mathbf{H}} ,u] u_{xx} dx.
\]
For these we need the pointwise bound
\begin{equation}\label{danny}
\|  [\operatorname{\mathbf{H}} ,u_x] u_x\|_{L^\infty} + \|[\operatorname{\mathbf{H}} ,u] u_{xx}\|_{L^\infty} \lesssim \|u_x\|_{H^{\frac12+\delta}}^2.
\end{equation}
Indeed, both terms on the left can be estimated in $H^{\frac12+\delta}$ using standard Littlewood-Paley
decompositions.
\end{proof}

By a Gronwall type argument, the two lemmas above lead to an $\epsilon^{-2}$ lifespan, proving part (a) of Theorem~\ref{th:main}.

Next, in order to prove part (b) of Theorem~\ref{th:main}, we consider the linearized Burgers-Hilbert equation:
\begin{equation}
\label{linearization}
w_t+wu_x+uw_x=\mathbf{H}w.
\end{equation}

Beginning with the normal form transformation  (\ref{nft}) for the Burgers-Hilbert equation, we obtain the normal form transformation  for the linearization (\ref{linearization}):
\begin{equation*}
q:=w+\vert \partial_x\vert (\mathbf{H}w\cdot\mathbf{H}u).
\end{equation*}
where $\vert \partial_x\vert = \mathbf{H}\,\circ \, \partial_{x}$.
Inspired by this we can define a modified linearized energy
\begin{equation}
\label{lin energy}
E_{lin}(w)=
\int w^2+2\vert \partial_x\vert w\cdot \mathbf{H}u\cdot\mathbf{H} w\,dx
\end{equation}
This is equivalent to the $L^2$ norm of $w$,
\[
E_{lin}(w) = (1+O(\|\operatorname{\mathbf{H}} u_x\|_{L^\infty}) \|w\|_{L^2}^2.
\]
This is proved using integration by parts in order to write
\begin{equation*}
E_{lin}(w)
=\int w^2-\vert \partial_{x}\vert u\cdot \left( \mathbf{H}w\right) ^2\,dx.
\end{equation*}

We want to prove that $E_{lin}$ satisfies good cubic bounds; a Gronwall type  argument implies part $(b)$ of Theorem~\ref{th:main}.
\begin{lemma}
For the energy (\ref{lin energy}) the estimate below holds: for any $\delta>0$,
\begin{equation*}
\frac{d}{dt} E_{lin}(w) \leq C_{k,\delta}  \Vert u_x\Vert _{H^{\frac12+\delta}}^2 \|w\|_{L^2}^2
\end{equation*}
where $C_{k,\delta}$ is a constant depending on $k$ and $\delta$ only.
\end{lemma}

\begin{proof}
A straightforward computation leads to
\begin{equation*}
\begin{aligned}
\frac{d}{dt}E_{lin}(w)&=\int 2\mathbf{H}w \cdot\mathbf{H}\left[wu_x \right]\cdot \mathbf{H}u_x+ 2\mathbf{H}w \cdot\mathbf{H}\left[uw_x \right]\cdot \mathbf{H}u_x
\\
&\hskip1in+ \vert \mathbf{H}w\vert^2 \cdot\vert \partial_{x}\vert (uu_x) \, dx\\
&=\int 2\mathbf{H}w \cdot\mathbf{H}\left[wu_x \right]\cdot \mathbf{H}u_x+2\mathbf{H}w\cdot \left[\mathbf{H},u \right]\partial_{x}w\cdot \mathbf{H}u_x
\\
&\hskip1in +
\vert \mathbf{H}w\vert^2 \partial_{x} \left( \left[ \mathbf{H}, u \right] \partial_{x}u \right)  \, dx.
\end{aligned}
\end{equation*}
The first term is estimated directly and for the second we use the commutator bound \eqref{com}. The third term can be controlled because we can rewrite
\[
\partial_{x} \left( \left[ \mathbf{H}, u \right] \partial_{x}u \right) = \left[ \mathbf{H}, u_x \right] \partial_{x}u + \left[ \mathbf{H}, u \right] \partial_{x}^2u,
\]
where the right hand side can be estimated by the pointwise bound \eqref{danny}.
\end{proof}


\section{Introduction}
The unprecedented control over ultracold atomic Fermi gases in optical lattices has made them prime candidates for studying and engineering novel quantum phases as well as probing fundamental theories such as the BEC-BCS crossover \cite{Greinter2003,Martin2003,Solomon2003,Hulet2003,Regal2004, Martin2004}. This progress has been facilitated by the tunability of the two-body interactions using Feshbach resonances or by changing the lattice confinement. If the interspecies two-body $s$-wave scattering length $a_s$ is much larger than the range $r_0$ of the underlying two-body potential, the few- and many-body behavior of equal-mass two-component Fermi gases is universal, i.e. completely
determined by $a_s$.

Presently significant efforts are directed at creating ultracold atomic Fermi gas mixtures composed of two chemically distinct species \cite{Gupta2011, Grimm2010, Martin2012}. This introduces a new parameter, the mass ratio $\kappa$ between the two species. This new parameter affects the many-body physics of the system, allowing one to realize novel quantum phases such as the interior gap superfluid \cite{Wilczek2003}. Here we show that
at the few-body level this additional degree of freedom leads to new three-body resonances which may destabilize the system, making it harder for experiments to explore novel quantum phases with unequal-mass mixtures. In particular, we study a system of two
identical heavy fermions with mass $M$, which interact with a light particle through a short-range potential with positive $s$-wave scattering length $a_s$.

Previous studies revealed two intriguing properties. First,  Kartavtsev and Malykh~\cite{KM} predicted the existence of a universal trimer state with energy $E_{\text{u},1}$ for $\kappa_1<\kappa<\kappa_2$ and the existence of two universal states with energies $E_{\text{u},1}$ and $E_{\text{u},2}$ for $\kappa_2<\kappa<13.606$; $\kappa_1$ and $\kappa_2$ were found to be $8.173$ and $12.917$, respectively. Second, Endo {\it et al.}\/~\cite{Endo} investigated how the universal trimer states, which are completely determined by the $s$-wave scattering length and the mass ratio $\kappa$, are connected to Efimov trimers, which have been predicted to exist for $\kappa \gtrsim 13.606$. By analyzing the trimer system within the framework of the Skorniakov-Ter-Martirosian equation with a momentum cutoff
$\Lambda_c$, Endo {\it et al.}\/ predicted the existence of a third class of trimer states, termed crossover trimers, which were shown to continuously connect the universal trimers described by Kartavtsev and Malykh and Efimov trimers.

Our study employs, as in Ref.~\cite{KM}, the hyperspherical coordinates.
However, while Ref.~\cite{KM} enforced that the hyperradial
wavefunction vanishes at hyperradius $R=0$, we explore the entirety of physically allowed boundary conditions by introducing a short-range three-body or hyperradial phase $\delta(R_0)$. We determine the eigenspectrum as a function
of the value of the
three-body phase $\delta$, the hyperradius $R_0$ at which the hyperradial boundary condition is imposed and the mass ratio $\kappa$.
The universal states of Kartavtsev and Malykh are recovered for $R_0 \to 0$ and $\delta(R_0)=\pi/2$. However, for other boundary conditions we find {\it deviations from universality}, which are linked to the {\it existence of a non-universal three-body state.} Analogous non-universal three-body states have previously been shown to exist \cite{Daily,DailyPRA} (see also Refs.~\cite{Petrov, Nishida, Werner}) in the $a_s\to \infty$ limit. The existence of the non-universal state for positive $a_s$ is described accurately within a quantum defect theory (QDT) framework. Moreover, within a two-state model, deviations from universality are explained as being due to the coupling between the non-universal state and the universal states.
Our work provides a simple intuitive Schr\"odinger equation
based description of the energy spectra of heavy-light trimers and
an alternative means to understanding the
connection between universal trimers and Efimov trimers.

The remainder of this paper is organized as follows:
Section~\ref{sec:framework} describes the hyperspherical framework.
Section~\ref{sec:numeric} determines the three-body energies as
functions of $\delta(R_0)$ and $\kappa$ by solving the hyperradial
Schr\"odinger equation numerically.
Section~\ref{sec:analytic} develops an analytical description, which accounts for the universal and non-universal states of the energy spectrum. Finally,
Sec.~\ref{sec:conclusion} concludes.

\section{System Hamiltonian}\label{sec:framework}
We model the interactions between the heavy and light particles by a zero-range two-body pseudopotential with $s$-wave scattering length $a_s$. If we denote the heavy fermions of mass $M$ as $1$ and $2$ and the light particle of mass $m$ as $3$, the Hamiltonian is given by
\begin{equation}
\label{eq:Htot}
H_{\text{tot}}=
-\frac{\hbar^2}{2M}(\vec \nabla_1^2+
\vec \nabla_2^2)-\frac{\hbar^2}{2m}\vec\nabla_3^2+V_{\text{int}},
\end{equation}
where $\vec r_i$ is the position vector of the $i$th particle, $\vec \nabla_i^2$ is the Laplacian of the $i$th particle and
$V_{\text{int}}=V_{\text{ps}}(\vec r_{13})+V_{\text{ps}}(\vec r_{23})$ with
\begin{equation}
\label{eq:Vij}
V_{\text{ps}}(r_{ij})=\frac{2 \pi \hbar^2 a_s}{\mu_{\text{2b}}} \delta^{(3)} (\vec r_{ij}) \frac{\partial}{\partial r_{ij}} r_{ij}.
\end{equation}
Here $\mu_{\text{2b}}$ is the reduced mass of the heavy-light pair, $\mu_{\text{2b}}=m\frac{\kappa}{\kappa+1}$, and $r_{ij}=\vert \vec r_{ij} \vert = \vert \vec r_i -\vec r_j \vert$. The pseudopotential
$V_{\text{ps}}$
imposes the Bethe-Peierls boundary condition on the three-body wave function in the limit $r_{ij}\to 0$.

To solve the Schr\"odinger equation for the Hamiltonian
$H_{\text{tot}}$, we separate off the center of mass motion and write the relative wave function $\Psi$ in terms of the hyperradius $R$ and five hyperangles, collectively denoted by $\Omega$~\cite{KM, Rittenhouse2010}. The hyperrradius $R$, which provides a measure of the overall size of the system, is defined by $\mu R^2=M [(\vec r_1-\vec R_{\text{cm}})^2+(\vec r_2-\vec R_{\text{cm}})^2]+m(\vec r_3-\vec R_{\text{cm}})^2$,
where $\mu$ is the three-body reduced mass associated with the hyperradius,
$\mu=m \kappa/\sqrt{2\kappa+1}$,
and $\vec{R}_{\text{cm}}$ denotes the center of mass vector.
We expand the relative wavefunction $\Psi$ in terms of a set of weight functions $F_n(R)$ and adiabatic channel functions $\Phi_n(R;\Omega)$,
which depend parametrically on the hyperradius $R$~\cite{Macek},
\begin{equation}
\label{eq:wf} \Psi(R,\Omega)=\sum_n R^{-5/2} F_n(R)\Phi_n(R;\Omega).
\end{equation}
The adiabatic channel functions $\Phi_n(R;\Omega)$ satisfy the hyperangular Schr\"odinger equation at fixed $R$,
\begin{equation}
\label{eq:haSE} \left[\Lambda^2+\frac{2\mu R^2}{\hbar^2}V_{\text{int}} (R, \Omega)\right]\Phi_n(R;\Omega)=\left[s_n^2(R)-4\right]\Phi_n(R;\Omega).
\end{equation}
In Eq.~\eqref{eq:haSE},
$\Lambda$ denotes the grand angular momentum operator,
which accounts for the kinetic energy
associated with the hyperangles $\Omega$ \cite{Rittenhouse2010}.
Inserting Eq.~\eqref{eq:wf} into the relative
Schr\"odinger equation yields a set of coupled
equations for the weight functions $F_n(R)$.

In the following, we employ the adiabatic approximation,
which neglects the coupling between the different
channels~\cite{Rittenhouse2010,Greene1980,Lin1995}.
This approximation has been shown to provide a
qualitatively correct description for a number of three-body
systems~\cite{KM,Esry1996,Esry1996He,Jonsell2002}.
In this approximation, the hyperradial Schr\"odinger
equation for the lowest adiabatic channel reads
\begin{equation}
\label{eq:Hhs}
-\frac{\hbar^2}{2\mu}
\left[\frac{d^2}{d R^2}-\frac{s_0^2(R)-1/4}{R^2}+Q_{00}(R)\right]
F_0(R)=EF_0(R),
\end{equation}
where $Q_{00}(R)$ is the diagonal correction
to the adiabatic energies,
\begin{equation}
Q_{00}(R)=
\langle \Phi_0(R;\Omega)\vert \frac{\partial^2}{\partial R^2}\Phi_0(R;\Omega)\rangle.
\end{equation}

We determine the three-body energies $E$ using a two step process. First we find the hyperangular eigenvalues $s_0(R)$ and coupling elements $Q_{00}(R)$. Then we solve the radial Schr\"odinger equation, Eq.~\eqref{eq:Hhs}. For states with $L^\pi=1^-$ symmetry and zero-range interactions, the scaled hyperangular eigenvalues $s_0(R)$ can be obtained semi-analytically by solving the
transcendental
equation~\cite{KM,Rittenhouse2010}
\begin{widetext}
\begin{equation}
\frac{R}{a_s}=\frac{(s_0^2-1)(1+2\kappa)^{1/4}\sec(\frac{\pi s_0}{2})\left[s_0 \kappa\,\, {}_2 F_1 \left(\frac{1}{2} (3- s_0), \frac{1}{2}(3+ s_0), \frac{5}{2}, \frac{\kappa^2}{(1+\kappa)^2}\right)-3 (1 + \kappa) \sin(\frac{\pi s_0}{2})\right]}{3 s_0 (1+\kappa)^{3/2}}.
\end{equation}
\end{widetext}
The hyperangular eigenvalues are completely determined by
$R/a_s$, $\kappa$ and $L^\pi$.
In the limit $R/a_s \to 0$,
$s_0$ goes to a constant.
For the purpose of this study we are only concerned
with positive values of $s_0^2(R)$ in the small
$R/a_s$ limit, i.e. we only consider mass ratios
$\kappa \lesssim 13.606$.
Figure~\ref{fig:s0vsk} shows the
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig1.pdf}
\caption{The hyperangular eigenvalue $s_0(R_0)$ as a function of the mass ratio $\kappa$ for $R_0/a_s=0.0001$. The dashed lines labeled I and II mark the mass ratios $\kappa_1$ and $\kappa_2$, at which the universal states with energies $E_{\text{u},1}$ and $E_{\text{u},2}$ first become bound. }
\label{fig:s0vsk}
\end{center}
\end{figure}
scaled eigenvalue $s_0(R_0)$ as a function of $\kappa$ for $R_0/a_s=0.0001$.
As we discuss in more detail below,
the value of $s_0(R)$ determines the allowed
short-range boundary condition of $F_0(R)$.
In the large $R/a_s$ limit, $Q_{00}$ vanishes and
the quantity $\frac{\hbar^2(s_0^2(R)-1/4)}{2\mu R^2}$
approaches $-E_{\text{2b}}$,
where  $E_{\text{2b}}$ denotes the dimer binding energy,
$E_{\text{2b}}=\hbar^2/(2\mu_{\text{2b}}a_s^2)$.

Equation~\eqref{eq:Hhs} has two linearly independent solutions,
which scale as $R^{s_0+1/2}$ and $R^{-s_0+1/2}$ in the small $R/a_s$
limit and are referred to as the regular and irregular solutions,
respectively.
For $s_0(R)\geq 1$ and $R_0/a_s\to0$,
the irregular solution is not normalizable and thus does not
contribute,
implying that $F_0(R)$ goes to zero as $R\to0$~\cite{Petrov,Nishida,Werner}.
On the other hand, in the regime $0<s_0<1$ both solutions
are well-behaved and must be included when constructing the general
solution~\cite{Petrov,Nishida,Werner}.
We parameterize the short-range boundary condition of the wavefunction $F_0(R)$ at $R_0$ using the logarithmic derivative $\mathcal{L}(F_0(R_0))$,
\begin{equation}
\label{eq:logd}R_0 \mathcal{L}(F_0(R_0))=
R_0\frac{\frac{d}{d R} F_0(R)}{F_0(R)}\vert_{R_0}=
\tan (\delta(R_0)),
\end{equation}
where $-\pi/2\leq\delta(R_0) \leq \pi/2$.
Using this parameterization, we cover all possible short-range phases.
In the special case of $\delta(R_0)=\pi/2$ and $R/a_s\to 0$,
this boundary condition and the resulting
three-body energies agree with those of Ref.~\cite{KM}.
We refer to this boundary condition as ``universal boundary condition'' since the three-body states are completely determined by the regular solution. We refer to the corresponding states as ``universal states'' with energies $E_{\text{u},1}$ for $\kappa_1\leq \kappa\leq \kappa_2$, and $E_{\text{u},1}$ and $E_{\text{u},2}$ for $\kappa_2\leq \kappa\leq 13.606$. In the next section we discuss deviations from universality that occur if $\delta(R_0)$ is not equal to $\pi/2$. These deviations increase with increasing $\kappa$ (for a fixed $R_0$) and increasing $R_0$ (for a fixed $\kappa$).

\section{Numerical Treatment} \label{sec:numeric}
In this section we examine the behavior of the three-body
bound states
with $L^{\Pi}=1^-$ symmetry,
obtained numerically using the shooting algorithm,
as a function of the logarithmic derivative boundary
condition for selected mass ratios.
The three-body energies $E$ for $\kappa = 8.6$ and $10$ are shown
in Figs.~\ref{fig:Ecomb}(a) and (b), respectively.
\begin{figure}
\begin{center}
\includegraphics[width=0.45\textwidth]{fig2.pdf}
\caption{(Color online)
The dotted, dashed, and dash-dotted lines show the scaled three-body
energy $E/E_{\text{2b}}$ as a function of $\delta$ using
$R_0/a_s=0.0001$, $0.0003$ and $0.001$, respectively,
for (a) $\kappa=8.6$ and (b) $\kappa=10$.
The solid horizontal lines show the energy $E_{\text{u},1}$
of the universal state and the solid vertical lines show
the energy $E_{\text{nu}}$ of the non-universal state
(see Sec.~\ref{sec:analytic}) for $R_0/a_s=0.001$.}
\label{fig:Ecomb}
\end{center}
\end{figure}
Dotted, dashed and dash-dotted lines show $E$ as a function of the three-body phase $\delta$ for $R_0/a_s=0.0001$, $0.0003$ and $0.001$, respectively.
In Fig.~\ref{fig:Ecomb} the three-body energies have been scaled by the zero-range two-body binding energy $E_{\text{2b}}$. The three-body state becomes unbound with respect to the break-up into a dimer and an atom for $E/E_{\text{2b}}>-1$. For $\kappa=8.6$ and $\delta=\pi/2$, the system supports one three-body bound state. For $R_0/a_s=0.0001$ and $\kappa=8.6$, $E/E_{\text{2b}}$ is nearly constant for $\delta\gtrsim -0.5$. At $\delta\approx -0.5$, referred to as the critical angle $\delta_c(R_0)$, the energy rapidly goes to a large negative value and a second bound state, whose energy is approximately equal to
$E_{\text{u},1}$, is supported for $\delta\lesssim\delta_c(R_0)$. We refer to the feature in the vicinity of $\delta_c$ as ``three-body resonance''. As $R_0/a_s$ increases  [see dashed and dash-dotted lines in Fig.~\ref{fig:Ecomb}(a)] the width of the three-body resonance increases. Note, however, that the deviations from $E_{\text{u},1}$ are small for all $R_0/a_s$ considered, except for three-body phases very close to $\delta_c(R_0)$. As $\kappa$ increases [see Fig.~\ref{fig:Ecomb}(b)], the overall behavior of the energy spectrum is unchanged. The key trends
with increasing $\kappa$ are that, at a fixed $R_0$,
the energy away from $\delta_c(R_0)$ becomes more
negative and both $\delta_c(R_0)$ and the
``width'' of the three-body resonance increase (see
also symbols in Figs.~\ref{fig:deltac} and \ref{fig:beta}).

For $\kappa\gtrsim \kappa_2$, the three-body system supports a second bound state.
As an example, Fig.~\ref{fig:Ek13} shows the three-body bound state
energy for $\kappa=13$ and $R_0/a_s=0.0001$ as a function of $\delta$.
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig3.pdf}
\caption{The solid lines show the scaled three-body energies $E/E_{\text{2b}}$ as a function of $\delta$ for $\kappa=13$ and $R_0/a_s=0.0001$. The three-body resonance is located at $\delta_c\approx0.15$. }
\label{fig:Ek13}
\end{center}
\end{figure}
Away from $\delta_c(R_0)\approx0.15$, there exist
two bound states whose energies depend weakly on $\delta$. For $\delta\gtrsim\delta_c$, the corresponding hyperradial functions $F_0(R)$ possess $0$
nodes and $1$ node, respectively.
For $\delta\lesssim\delta_c$,
the corresponding hyperradial functions $F_0(R)$
possess $1$ node and $2$ nodes, respectively. This reflects the fact that a new bound state is being pulled in at $\delta\approx \delta_c(R_0)$.

In general, we find that $\delta_c(R_0)$ depends relatively weakly on $R_0$ as long as $R_0/a_s\ll1$. Moreover,
the three-body energy depends relatively strongly on $R_0$ in the vicinity of $\delta_c$ but comparatively weakly on $R_0$ away from $\delta_c$. This suggests that the states near $\delta_c$ and away from $\delta_c$ can be classified as non-universal and universal, respectively. This interpretation is corroborated by our analysis of the hyperradial wavefunction $F_0(R)$ for $\kappa=10$, $R_0/a_s=0.0001$ and $\delta\approx \delta_c$. The main part of Fig.~\ref{fig:k10drop} shows that the wavefunction has an appreciable amplitude in the small $R/a_s$ region, signaling non-universal behavior.
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig4.pdf}
\caption{The solid line shows the wave function $F_0(R)$ for
$\kappa=10$,
$R_0/a_s=0.0001$, $E/E_{\text{2b}}=-1.4$, and $\delta=-0.33$.
As can be seen in Fig.~\ref{fig:Ecomb}(b),
this angle is close to $\delta_c$.
For comparison, the inset shows the universal wavefunction
for $\delta=1.53$ (i.e. away from the three-body resonance).}
\label{fig:k10drop}
\end{center}
\end{figure}
Away from $\delta_c$, in contrast, $F_0(R)$ has a vanishingly small amplitude in the
small $R/a_s$
region (see inset of Fig.~\ref{fig:k10drop}). We find that $F_0(R)$ depends fairly weakly on $R_0$ when $\delta$ is away from $\delta_c$, lending further support to our assertion that the three-body system behaves universally in this region. For larger $R_0/a_s$, we observe deviations from universality
for a larger range of three-body phases.  Similar non-universal behavior has been previously reported for the three-body system at unitarity~\cite{Daily,DailyPRA} and for the four-body system with positive $s$-wave scattering
length~\cite{Blume4b}.

\section{Analytical Treatment}\label{sec:analytic}
Section~\ref{sec:analytic}A applies a QDT framework
to predict the short-range phase $\delta_c(R_0)$ at
which the non-universal state first becomes bound.
Section~\ref{sec:analytic}B develops a two-state model to describe the behavior of the universal and non-universal states as a function of $\delta(R_0)$.
\subsection{QDT Treatment}
To predict the short-range phase $\delta_c(R_0)$ at which the non-universal state first becomes bound with respect to the break-up into a dimer and an atom we apply QDT~\cite{Daily,DailyPRA,Fano,Seaton}.
In the short-range limit, the wavefunction $F_0(R)$ can be approximated
by~\cite{Daily,DailyPRA}
\begin{eqnarray}
\label{eq:SR} F_{\text{SR}}(R)=\nonumber \\
\sqrt{R}[ J_{s_0(R)}(\sqrt{\frac{2\mu E}{\hbar^2}}R)-\tan (\pi\alpha) \,\, Y_{s_0(R)}(\sqrt{\frac{2\mu E}{\hbar^2}} R)],
\end{eqnarray}
where $J_{s_0(R)}$ and $Y_{s_0(R)}$ denote the Bessel functions of the first and second kind, respectively. The quantum defect $\alpha$ controls the relative contribution of the regular solution $J_{s_0(R)}$ and the irregular solution $Y_{s_0(R)}$. A new three-body state is expected to be pulled in when the hyperradial solution is dominated by the irregular solution, i.e. for
$\alpha=1/2$. The critical angle $\delta_c(R_0)$ is then given by
\begin{eqnarray}
\tan \delta_c(R_0)=\operatorname{Re}[R_0 \mathcal{L}(F_{\text{SR}}(R_0))],
\end{eqnarray}
where $\operatorname{Re}$ denotes the real part and $F_{\text{SR}}(R_0)$ is
evaluated for $E=-E_{\text{2b}}$ and $\alpha=1/2$.

The solid, dotted and dashed lines in Fig.~\ref{fig:deltac} show
$\delta_c(R_0)$ determined using the QDT framework for
$R_0/a_s=0.0001$, $0.0003$ and $0.001$, respectively.
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig5.pdf}
\caption{(Color online) The critical short-range phase $\delta_c(R_0)$ as a function of $\kappa$. The solid, dotted and dashed lines show the QDT prediction for $R_0/a_s=0.0001, 0.0003$ and $0.001$, respectively.  The circles show $\delta_c(R_0)$ obtained by analyzing the numerical solutions of the hyperradial Schr\"odinger equation, Eq.~\eqref{eq:Hhs}, for $R_0/a_s=0.0001$. The inset shows a blow up of the large $\kappa$ region.}
\label{fig:deltac}
\end{center}
\end{figure}
It is interesting to note that $\delta_c(R_0)$ is very weakly dependent on $R_0$ for small $\kappa$. The symbols in Fig.~\ref{fig:deltac} show the critical angle $\delta_c(R_0)$ obtained by analyzing the numerical solutions of the hyperradial Schr\"odinger equation, Eq.~\eqref{eq:Hhs}, for $R_0/a_s=0.0001$. The numerical results are in excellent agreement with the QDT prediction. For $\kappa\gtrsim13$ (see inset of Fig.~\ref{fig:deltac}), the dependence of $\delta_c$ on $R_0$ becomes more pronounced.
Figure~\ref{fig:deltac} shows that the three-body system supports a non-universal bound state not only for $\kappa\gtrsim 8.619$, but also for $\kappa\lesssim8.619$
[for $\kappa=8.619$ one has $s_0(0)=1$].
Thus, it may be surprising at first sight that the system supports
a non-universal state for $\kappa \lesssim 8.619$
since the irregular solution cannot be normalized if $s_0>1$ for $R_0/a_s\to 0$. However, since we impose the boundary condition at a finite $R_0/a_s$ and not at $R_0/a_s=0$, the resulting wavefunction can, even though it contains an admixture of the irregular solution, be normalized. Correspondingly, non-universal states can exist for $\kappa<8.619$.
In fact, we find that the system supports a non-universal bound state even for $\kappa<\kappa_1$, i.e. for mass ratios where universal states are not supported.

\subsection{Two-state model}
In this section we develop a two-state model that describes
the behavior of the universal and non-universal states
as a function of $\delta(R_0)$. In our model the three-body resonance is an avoided crossing at $\delta_c(R_0)$ between the universal state with energy
$E_{\text{u}}$ and a non-universal state with energy $E_{\text{nu}}$.
In Fig.~\ref{fig:scaledE} we plot the quantity
$(E-E_{\text{u}})/E_{\text{2b}}$ as a function of $\delta-\delta_c$.
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.5\textwidth]{fig6.pdf}
\caption{(Color online)
Thick solid, dotted and dashed lines show the
quantity $(E-E_{\text{u},1})/E_{\text{2b}}$  as a function
of $\delta-\delta_c$ for $R_0/a_s=0.0001$ and $\kappa=10,\, 12.4$, and $13$,
respectively.
The circles show the quantity $(E-E_{\text{u},2})/E_{\text{2b}}$
for $R_0/a_s=0.0001$ and $\kappa=13$.
The thin dashed lines show the results of the two-state model.}
\label{fig:scaledE}
\end{center}
\end{figure}
The thick solid, dotted and dashed lines correspond to
$\kappa=10, 12.4$ and $13$, respectively, with $R_0/a_s=0.0001$ and
$E_{\text{u}}=E_{\text{u},1}$. The circles correspond
to $\kappa=13$, $R_0/a_s=0.0001$ and $E_{\text{u}}=E_{\text{u},2}$.

Motivated by Fig.~\ref{fig:scaledE}, we write the energy of the non-universal state as $E_{\text{nu}} =E_{\text{u}}+\mathcal{A}\frac{\hbar^2}{2\mu R_0^2}(\delta-\delta_c(R_0))$, where $\mathcal{A}$ is a dimensionless scaling constant; in our analysis we use $\mathcal{A}=0.00007$. The term proportional to $\mathcal{A}$ determines the slope at which the non-universal state crosses the universal state. In our parameterization of $E_{\text{nu}}$, $\delta_c(R_0)$ should be interpreted as the three-body phase at which the universal and non-universal states cross. Our numerical results show that this crossing point is nearly identical to the QDT three-body phase at which the non-universal state first becomes bound. Hence we
use the QDT value in our two-state model. The dimensionless two-state Hamiltonian is then
\begin{equation}
\label{eq:2st} H= \left( \begin{array}{cc}
E_{\text{u}}/E_{\text{2b}} &\beta/E_{\text{2b}}\\
\beta^\star/E_{\text{2b}} &E_{\text{nu}}/E_{\text{2b}} \\
\end{array} \right),
\end{equation}
where $\beta=\beta(R_0, \kappa)$ is the coupling between the universal and non-universal states.
Within this model, $\beta$ is real and
the scaled three-body energies $E/E_{\text{2b}}$
are given by the eigenvalues of the two-state Hamiltonian.
We apply the two-state model separately to the universal states with energies $E_{\text{u},1}$ and $E_{\text{u},2}$, yielding separate $\beta$ values for the two universal states. The values of $\beta$ are obtained by fitting the energy spectrum predicted by the two-state Hamiltonian to the numerically determined energies. Thin dashed lines in Fig.~\ref{fig:scaledE} show the results for the two-state model.
It can be seen that the two-state model provides
a quantitatively correct description of the
three-body spectra. Moreover, it provides
an intuitive physical picture in which deviations
from universality arise due to the
coupling of the universal states to the non-universal state.

The circles, diamonds and squares in Fig.~\ref{fig:beta}
show $\beta$ for $E_{\text{u}}=E_{\text{u},1}$ and
$R_0/a_s=0.0001, 0.0003$ and $0.001$, respectively, as a function of $\kappa$.
\begin{figure}[h]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig7.pdf}
\caption{(Color online)
Circles, diamonds and squares show the off-diagonal coupling $\beta/E_{\text{2b}}$ used in the two-state model, Eq.~\eqref{eq:2st}, as a function of $\kappa$ for $R_0/a_s=0.0001,0.0003$ and $0.001$, respectively, on a logarithmic scale. We used $E_{\text{u}}=E_{\text{u},1}$ in the two-state model. The dependence of $\beta/E_{\text{2b}}$  on the value of $R_0$ increases with increasing $\kappa$. }
\label{fig:beta}
\end{center}
\end{figure}
The value of $\beta/E_{\text{2b}}$
increases with increasing $\kappa$ for fixed $R_0$.
Moreover, the dependence of $\beta/E_{\text{2b}}$ on the value of $R_0$ increases with increasing $\kappa$. Since $\beta$ determines the coupling between the universal state and the non-universal state, it can be used to quantify the deviations from universality. For mass ratios larger than those considered in this work, the heavy-light trimer system supports three-body bound states with Efimov character \cite{Efimov1970, Efimov1973, Braaten2006}. Within the zero-range framework employed here, the exact number and energy of the Efimov trimers supported depends on the short-range hyperradial boundary condition. If the hyperradial boundary condition is fixed and $\kappa$ is varied, the energy spectrum changes smoothly from (i) deviating from the universal spectrum only in a small region around $\delta_c$ to (ii) deviating from the universal spectrum for a fairly large range of $\delta$ to (iii) supporting three-body states with Efimov character for $\kappa\gtrsim13.606$. Thus, the deviations from universality discussed in this paper can be interpreted as connecting the universal states predicted by Kartavtsev and Malykh and the Efimov trimers.

Our findings are in qualitative agreement with the results reported in Ref.~\cite{Endo}. In that work, the authors imposed a momentum cutoff $\Lambda_c$ in the Skorniakov-Ter-Martirosian equation and used the quantity $\frac{E-E_{\text{2b}}}{E_{\text{u},i}-E_{\text{2b}}}$ ($i=1$ or $2$)
to determine the ``boundary'' between the universal trimers
predicted by Kartavtsev and Malykh and the crossover
trimers in the $(\kappa,(\Lambda_c\, a_s)^{-1})$ parameter space.
The momentum cutoff was introduced in two ways, using a sharp
and a Gaussian cutoff.
We speculate that in our formulation a change in the three-body phase and/or hyperradius $R_0$ corresponds to a change in $\Lambda_c$.  A crossover trimer, in turn, corresponds to a trimer whose energy deviates appreciably from $E_{\text{u},i}$. While the deviations from $E_{\text{u},i}$ are, in our formulation, linked to the coupling of the universal state to a non-universal state, the treatment by Endo {\it et al.}\/ does not seem to yield an analogous physical picture. Additionally we speculate that while Ref.~\cite{Endo} employs two different models for the momentum cutoff $\Lambda_c$, it does not explore the entirety of
the $(R_0, \delta(R_0))$ parameter space.
In the future it will be interesting to investigate the precise connection between the formulations in the coordinate and momentum spaces by, e.g., comparing the wavefunctions. Such a comparison is needed to check if the above correspondencies are correct.

\section{Conclusion}\label{sec:conclusion}
In this paper we studied a system of two identical heavy
fermions of mass $M$ and a light particle of mass $m$ with zero-range
two-body interspecies interactions.
In particular, we looked at deviations of the three-body bound state
energies from the universal energies $E_{\text{u},1}$
for $\kappa>\kappa_1$ and $E_{\text{u},1}$ and $E_{\text{u},2}$
for $\kappa_2<\kappa\lesssim 13.606$ as a function of the
hyperradial short-range boundary condition.
We imposed a short-range phase $\delta(R_0)$ using a
logarithmic derivative boundary condition at various hyperradii $R_0$.
This parameterization allowed us to explore the full range of
possible short-range boundary conditions.

We found that
(i) for $\delta(R_0)=\pi/2$ the universal states with
energies $E_{\text{u},1}$ and $E_{\text{u},2}$,
predicted by Kartavtsev and Malykh~\cite{KM}, are recovered;
(ii) the three-body
states deviate from universality in the vicinity of a
three-body resonance located at
the short-range phase $\delta_c(R_0)$,
at which the non-universal state is first bound;
(iii) the deviations from universality increase
with increasing mass ratio $\kappa$
(at fixed $R_0$) and with increasing $R_0$ (at fixed $\kappa$);
(iv) QDT accurately predicts the values of $\delta_c(R_0)$;
(v) a two-state model quantitatively describes the behavior of
the universal and non-universal states as a function of $\delta(R_0)$;
and, finally,
(vi) the non-universal bound state exists
for $\kappa < \kappa_1$ even though universal bound states
are not supported in this regime.

\section{Acknowledgements}
This work was supported by the National
Science Foundation through a grant for the Institute for
Theoretical Atomic, Molecular and Optical Physics at
Harvard University and Smithsonian Astrophysical Observatory.
DB acknowledges support by the NSF through grant PHY-1205443.


\section{Introduction}

Significant experimental activity is currently devoted to the search for Dark Matter (DM)
by looking at the excesses in cosmic ray production through DM annihilations (or decays) in the galactic halo.
Detailed predictions for the  fluxes largely depend on the particle physics model for the DM,
e.g. the mass, the annihilation channels etc.
It is therefore desirable to describe DM annihilations and their products within
a general and model-independent framework and Effective Field Theory (EFT) provides such a tool.

Of course, the EFT is  only applicable whenever there is a separation of scales between the process
to describe (the annihilation of non-relativistic DM at a scale $\sim M_{\rm DM}$) and the underlying
microscopic physics of the interactions (at a scale $\Lambda$). This may not always be the case, as, for example, in
the case of supersymmetry with a compressed spectrum.

If we want to describe the annihilation of two non-relativistic DM particles, whose relative velocity
is $v\sim 10^{-3}$ (in units of $c$) in our Galaxy today, it is convenient to expand the
cross section in powers of $v$
\begin{equation}}   \def\ee   {\end{equation}
v\sigma=a+b\,v^2+\mathcal{O}(v^4)\,,
\ee
where the first term corresponds to annihilation in the state of orbital angular momentum $L=0$ ($s$-wave)
while the second term describes $L=1$ ($p$-wave).
For the annihilation $\textrm{DM\, DM}\to f\bar f$ of a self-conjugate DM particle (real scalar or Majorana fermion)  into SM fermions of mass $m_f$,
helicity arguments lead to
$a\propto (m_f/M_{\rm DM})^2$, and hence a very suppressed s-wave term for
light final state fermions (e.g.~leptons), while the $p$-wave is suppressed by $v^2$.

It is clear that the correct operator expansion must be done in terms of two parameters:
mass dimension of the operator and relative velocity. The effective
lagrangian would be generically given by an infinite series of non-renormalizable operators
\begin{equation}}   \def\ee   {\end{equation}
\mathscr{L}_{\rm eff}=
\sum_{d>4}
{1\over \Lambda^{d-4}}
\left(c^{(d,s)}{\cal O}^{\,(d, s)} + c^{(d,p)}{\cal O}^{\,(d, p)} \right)\,,
\ee
where ${\cal O}^{(d, s \,{\rm or}\, p )}$ indicates that the operator of dimension $d$
describes $s$-wave or $p$-wave annihilations,
and we neglect annihilations in waves higher than $p$.
There is no obvious ordering of the importance of the operators. An important role in this respect
is played by ElectroWeak
(EW) bremsstrahlung, which has received significant attention recently \cite{paper0, bell1, paper1, bell2,
cheung, ibarra1, paper2, iengo, barger, ibarra2,proc} (for earlier studies on the impact of
gauge boson radiation on DM annihilations or cosmic ray physics, see \cite{bergstrom1, bergstrom2, list1}).
Taking into account processes with the inclusion of EW radiation  eludes the helicity suppression and
opens up an $s$-wave contribution to the cross section
\cite{bell1, paper1, ibarra1, paper2, ibarra2, bergstrom1}.

In this paper we will classify the operators (up to dimension 8) according to the $v$-behaviour of the amplitude
connecting two self-conjugate particles -- real scalars or Majorana fermions --  in the initial state  with the final state of two massless fermions and possibly a gauge boson.

The effects of lifting the helicity suppression  by means of EW radiation has so far been studied within
the context of explicit models \cite{bell1, paper1, ibarra1, paper2, ibarra2}.
The results of the present paper provide a model-independent approach to this problem and can
be used to place robust constraints on the new physics responsible for the DM sector.
Although several analyses in the literature place phenomenological constraints on the coefficients of the dimension-six operators  (see e.g. Refs.~\cite{eftconstraints1, eftconstraints2,eftconstraints3, eftconstraints4}),
the important role of higher-dimensional operators is typically underestimated.

The remainder of the paper is organized as follows.
In Section \ref{sec:ops}, we explain our methodology and construct the
effective operators contributing to $s$-wave DM annihilations.
We compute the differential and total annihilation cross sections
for each operator in Section \ref{sec:xsec}, and compare them
to the contribution from the typical lowest-dimensional operator.
We conclude in Section \ref{sec:conclusions}, mentioning phenomenological
applications and possible analyses that could be carried out using these results.
Finally, we collect some useful relations and identities in the Appendices,
for the convenience of the reader.

\section{Effective operators}
\label{sec:ops}

We will carry out our analysis under the following set of assumptions on the DM sector:
\begin{enumerate}
\item  the DM is either
a Majorana fermion or a real scalar field;
\item the DM is neutral under the SM gauge group;
\item  there exists a $Z_2$ symmetry under which the DM is odd and the SM is even;
\item the DM couples only to the fermions in the SM spectrum, which are assumed to be massless.
\end{enumerate}
Of course, any of these assumptions may not hold in reality and if this is the case
our analysis requires modifications.
Assumption 1 specifies the conditions under which the helicity suppression is effective,
while  assumption 2 is there to simplify the discussion, even though it is not strictly necessary
and relaxing this assumption can also lead to interesting effects
(see Refs.~\cite{paper2, paper3}).
Assumption 3 is commonly used in DM phenomenology to ensure the stability of the DM particle.

More clarifications about assumption 4 are in order.
 Considering $m_f\neq 0$
would introduce another mass scale into the problem and would render the operator
classification much less transparent. Our analysis is still valid in the regime $m_f / \Lambda \ll v$,
which may not hold for the third-generation quarks.
The other piece of information in assumption 4 is that the DM particle only couples to the fermion
sector of the SM. This needs not to be true, of course.
Allowing for DM interactions with the other SM particles, namely gauge bosons or the Higgs,
the possibility of having additional operators contributing to  DM annihilations in
$s$-wave opens up.
For instance for Majorana DM, $\chi$, coupling to the Higgs doublet $H$ or a generic  field-strength
$F^{\mu\nu}$ one can have $(\bar\chi\gamma^5\chi)(H^\dag H)$
and
$\left(\bar{\chi}\gamma^5\chi\right)
F_{\mu\nu}F^{\mu\nu}$, which are CP odd, and
$\left(\bar{\chi}\gamma^5\chi\right)
F_{\mu\nu}\tilde F^{\mu\nu}$, which is CP even.
For a list of effective operators connecting DM to vector bosons and/or Higgs bosons see e.g. Refs.\cite{Rajaraman, rizzo}. We will not consider these possibilities here.

We want to classify the operators of dimension $d=6, 7,  8$ according to the $v$-behaviour of the amplitude
connecting two DM particles in the initial state  with two massless SM fermions $f$ and a gauge boson in the final state.
We look for operators which are hermitian, gauge invariant (under $SU(2)_L\otimes U(1)_Y$) and
 giving a non-zero contribution to the amplitude for the annihilation of DM into two fermions
and one gauge boson. We will further classify the operators according to their $CP$ transformation
properties (see Appendix \ref{app:CP} for the relevant transformation properties).
The operators containing $\slashed{D}f$ give no contribution to the process under consideration,
due to the Equation Of Motion (EOM).
In order to ensure manifest gauge invariance,
we first introduce the following notation for the covariant derivative
\begin{eqnarray}
\overrightarrow{D}_{\mu}f&=& \left(\overrightarrow{\partial}_{\mu}+igT^{a}W_{\mu}^{a}+ig^{\prime}Y_fB_{\mu}\right)f,\label{eq:cov1}\\
\bar{f}\overleftarrow{D}_{\mu}&=&
\bar{f}\left(\overleftarrow{\partial}_{\mu}-igT^{a}W_{\mu}^{a}-ig^{\prime}Y_fB_{\mu}\right),\label{eq:cov2}
\end{eqnarray}
where $T^a=\sigma^a/2$, with $\sigma^a$ being the usual Pauli matrices, and the charge $Q_f$ of the fermion $f$
is related to its hypercharge $Y_f$ by $Q_f=T^3+Y_f$.
Furthermore, the field-strength $W_{\mu\nu}$ of the SM gauge fields is related to the covariant derivatives as usual
\begin{eqnarray}
\left[\overrightarrow{D}_{\mu},\overrightarrow{D}_{\nu}\right]f&=&(igT^a)W^a_{\mu\nu}f\,,\label{eq:dualidentity1}\\
\overline{f}\left[\overleftarrow{D}_{\mu},\overleftarrow{D}_{\nu}\right]&=&\overline{f} (ig T^a)W^a_{\mu\nu}\,.\label{eq:dualidentity2}
\end{eqnarray}
In the following, we will separately deal with the cases where the DM is a Majorana fermion or a real scalar.
For simplicity, we will consider only left-handed SM fermions  $f_L$, but the analysis can be
applied to right-handed fermions straightforwardly.

\subsection{Majorana fermion DM}

Let us first suppose the DM particle is a Majorana fermion $\chi$.
It is possible to build several operators containing two DM fields, two SM fermion fields and zero or
one gauge bosons.
They can be built in full generality requiring gauge invariance and hermiticity,
and further classified according to their $CP$ properties.
The Majorana-flip properties and the chiralities of the SM fermions make several structures identically zero.
The only two Majorana fermion bilinears that are non-vanishing in the limit $v\to 0$ are $\bar{\chi}\gamma^5\chi$
and $\bar\chi\gamma^\mu\gamma^5\chi$ (see  Appendix \ref{app:bilinears}).
As we are interested in $s$-waves only, we will limit our analysis to these two bilinears for the Majorana fermions.
For the pseudo-scalar bilinear  $\bar \chi \gamma^5\chi$ no contractions with SM fermions can be built,
as they would vanish either by chirality or by the
EOM of  fermions, so we are left with the axial-vector bilinear.

The lowest-dimensional terms that can be written out of two $\chi$'s and two $f_L$'s are of dimension 6.
At this level, there is only one non-vanishing operator satisfying all criteria:
\begin{equation}
\mathcal{O}_{\rm M}=
\left(\bar{\chi}\gamma^5\gamma^{\mu}\chi\right)
\left[
\bar{f}_L\gamma_{\mu} f_L
\right]\,.
\label{O6}
\end{equation}
The $\mu=0$ component of this operator could {\it a priori} give a $v$-independent
contribution to the scattering amplitude, but it actually vanishes because of the identity
$u_f^\dag({\mathbf p}) v_f(-{\mathbf p})=0$.
Therefore this dimension-6 operator only contributes to the $p$-wave, as  noted in Ref.~\cite{paper1}, due to the helicity suppression, which cannot be removed
by simply radiating a gauge boson from the external final state leg.
In order to look for $s$-wave terms, we need to consider higher-dimensional operators with one EW gauge
boson whose radiation in the annihilation process  lifts the helicity suppression.

At dimenion 7, there are no operators contributing to the $s$-wave cross section.
In fact, all possible structures vanish either due to Majorana-flip properties or the chirality of the SM fermions, or because of the EOM of the $f_L$s.
A priori, the $\mu=0$ component of the operator $(\bar{\chi}\gamma^5{\partial_\mu}\chi)
[\overline{f}_L\gamma^\mu f_L]$ would give a $v$-independent contribution to the
scattering amplitude, but it vanishes again because of the identity
$u_f^\dag({\mathbf p}) v_f(-{\mathbf p})=0$.

At the level of dimension 8, there are several structures that can be built
requiring gauge invariance and hermiticity.
They contain two $\chi$'s (in the bilinear $\bar\chi\gamma^\mu\gamma^5\chi$), two $f_L$'s and two covariant derivatives.
It is possible to reduce the number of independent operators, contributing to the
cross section for the process under consideration, by using EOM and the identities
in Appendix \ref{app:identities}.
In addition, some structures can be related to each other by  terms (like
$\left(\bar{\chi}\gamma^5\gamma^{\mu}\chi\right)\partial^2\left[
\bar{f}_L\gamma_{\mu}f\right]$),
which do not contribute to the $s$-wave annihilation into two fermions and a gauge boson.
Therefore they contribute in exactly the same way to the amplitude for the process we are interested in.
There remain only five independent operators of dimension 8 contributing to the
 $s$-wave annihilation of DM into two SM fermions and a gauge boson,
 listed in Table \ref{tableOM} together with their $CP$ conjugation properties.
We remain agnostic about the presence or absence of $CP$ violation in the Dark Matter
sector, which can possibly induce $CP$ violation in the SM at loop level and therefore be
further constrained.
\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\textbf{Name} & \textbf{Operator} & \textbf{CP}\\
\hline
$\mathcal{O}_{{\rm M}1}$ &
$\left(\overline{\chi}\gamma^5\gamma^{\mu}\chi\right)
\left[\left(\overline{f}_L\overleftarrow{D}_{\rho}\right)\gamma_{\mu}\left(\overrightarrow{D}^{\rho}f_L\right)\right]$ & $+$\\
\hline
$\mathcal{O}_{{\rm M}2}$&
$i \epsilon_{\mu\nu\rho\sigma} \left(\overline{\chi}\gamma^5\gamma^{\mu}\chi\right)
\left[\overline{f}_L\gamma^{\nu}\overrightarrow{D}^{\rho}\left(\overrightarrow{D}^{\sigma}f_L\right)-
\left(\overline{f}_L
\overleftarrow{D}^{\sigma}\right)\overleftarrow{D}^{\rho}\gamma^{\nu}f_L\right]$
& $+$ \\
\hline
$\mathcal{O}_{{\rm M}3}$ &
$i\epsilon_{\mu\nu\rho\sigma}
\left(\overline{\chi}\gamma^5\gamma^{\mu}\chi\right)
\left[\left(\overline{f}_L\overleftarrow{D}^{\nu}\right)\gamma^{\rho}
\left(\overrightarrow{D}^{\sigma}f_L\right)-
\left(\overline{f}_L
\overleftarrow{D}^{\sigma}\right)\gamma^{\rho}\left(\overrightarrow{D}^{\nu}f_L\right)\right]$
& $+$\\
\hline
$\mathcal{O}_{{\rm M}4}$ &
$
i \left(\overline{\chi}\gamma^5\gamma^{\mu}\chi\right) \left[\overline{f}_L\overrightarrow{\slashed{D}}\left(\overrightarrow{D}_{\mu}f_L\right)-\left(\overline{f}_L
\overleftarrow{D}_{\mu}\right)\overleftarrow{\slashed{D}}f_L\right]$
& $-$ \\
\hline
$\mathcal{O}_{{\rm M}5}$ &
$ i \left(\overline{\chi}\gamma^5\gamma^{\mu}\chi\right) \left[\overline{f}_L\gamma_{\mu}\overrightarrow{D}_{\rho}\left(\overrightarrow{D}^{\rho}f_L\right)-\left(\overline{f}_L
\overleftarrow{D}_{\rho}\right)\overleftarrow{D}^{\rho}\gamma_{\mu}f_L\right]$
& $-$\\
\hline
\end{tabular}
\caption{List of dimension-8 operators contributing to the $s$-wave  cross section for the
annihilation of Majorana DM into two fermions and a gauge boson.
\label{tableOM}}
\end{center}
\end{table}

 All other operators have either a larger dimensionality or produce more powers of $v^2$ in
the annihilation cross section.
Notice that we chose to keep a Lorentz-covariant formalism, despite looking at the non-relativistic
limit. This implies that the same operator can lead to both $v$-independent and $v$-dependent
terms in the amplitudes; for example, the operator $\mathcal{O}_{{\rm M}1}$ also gives a contribution
to the $p$-wave cross section, but we will not consider it as it is very suppressed.

To summarize, the dimension-8 operator contributing to the $s$-wave annihilation
cross section of Majorana DM into SM fermions is given by the sum of the operators in Table
\ref{tableOM},
$\mathcal{O}_{{\rm M}}^{(8,s)}=\sum_{i=1}^5 c^{(8,s)}_{i} \mathcal{O}_{{\rm M}i}$;
the leading interactions are therefore described in terms of only a few  operators
\begin{equation}}   \def\ee   {\end{equation}
\mathscr{L}_{\rm eff}=
{1\over \Lambda^2}c^{(6,p)}{\cal O}_{\rm M}+
{1\over \Lambda^4}\sum_{i=1}^5 c^{(8,s)}_{i}  \mathcal{O}_{{\rm M}i}
+\textrm{higher-dim}\,.
\ee
In absence of $CP$-violation in the DM sector, only the first three operators
in the sum need to be considered.

\subsection{Real scalar DM}
\label{subsec:scalarops}

Next, let us consider the case where the DM particle is a real scalar $\phi$.
By angular momentum conservation, two real scalars  cannot annihilate into two massless fermions in the configuration with $L=0$; this process would  only occur through a chirality flip induced by a mass term. For this reason,
in the limit $m_f=0$ we are considering here, no $s$-wave annihilation $\phi\phi \to \bar f_L f_L$ is possible.
One can recover this result in the language of effective operators.
At dimension 5, only a single operator can be constructed which, however, vanishes due to chirality,
$\phi^2 (\overline{f}_L f_L) = 0$.
At dimension 6, we have $\phi^2 \partial_\mu(\overline{f}_L\gamma^\mu f_L) = 0$,
 by the EOM (cf. Eq.~(\ref{eq:ID1bis})).
At dimension 7, there are no possible Lorentz contractions to construct an operator.

Nevertheless,
at the level of dimension 8, several gauge invariant hermitian operators can be built out of two
$\phi$'s,
two $f_L$'s and covariant derivatives. As discussed already for the Majorana case,
it is possible to reduce the number of independent operators, contributing to the
cross section for the process under consideration, by using the EOM and the identities
in Appendix \ref{app:identities}.
We are left with four CP-even operators and three CP-odd operators,  listed in Table \ref{tableOR}.
A $v$-dependent  annihilation $\phi\phi \to \bar f_L f_L$ is mediated by the operator ${\cal O}_{\rm R4}$, while the $s$-wave annihilation of two $\phi$'s can proceed by switching on the emission
of a gauge boson in the final state.

Other operators involving $\phi\partial_\mu\phi$ or a gauge field-strength $F^{\mu\nu}$ can be obtained from the listed ones by integration by parts, using the EOM of the DM particle or the identities
(\ref{eq:dualidentity1})-(\ref{eq:dualidentity2}).
For instance, the operator $\phi^2 \partial_\nu \left[ \bar{f}_L \gamma_\mu f_L \right] F^{\mu\nu}$, considered in Ref.~\cite{barger2}, is expressed in this basis as $1/g ({\cal O}_{\rm R1}-{\cal O}_{\rm R2})$.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
\textbf{Name} & \textbf{Operator} & \textbf{CP}\\
\hline
$\mathcal{O}_{{\rm R1}}$ &
$i \phi^2 \left[\overline{f}_L\overrightarrow{\slashed{D}}\overrightarrow{D}^{\nu}
 (\overrightarrow{D}_{\nu}f_L)-(\overline{f}_L\overleftarrow{D}_{\nu})\overleftarrow{D}^{\nu}\overleftarrow{\slashed{D}}f_L
 \right]
$ & $+$\\
\hline
$\mathcal{O}_{{\rm R2}}$ &
$i\phi^2 \left[\overline{f}_L\overrightarrow{D}^{\nu}\overrightarrow{\slashed{D}}
 (\overrightarrow{D}_{\nu}f_L)-(\overline{f}_L\overleftarrow{D}_{\nu})\overleftarrow{\slashed{D}}\overleftarrow{D}^{\nu}f_L
 \right]
$ & $+$\\
\hline
$\mathcal{O}_{{\rm R3}}$ &
$ \phi^2 \epsilon_{\mu\nu\rho\sigma} \left[
   \overline{f}_L\gamma^{\mu}\overrightarrow{D}^{\nu}\overrightarrow{D}^{\rho}(\overrightarrow{D}^{\sigma}f_L)+
   (\overline{f}_L\overleftarrow{D}^{\sigma})\overleftarrow{D}^{\rho}\overleftarrow{D}^{\nu}\gamma^{\mu}f_L
   \right]
$ & $+$\\
\hline
$\mathcal{O}_{{\rm R4}}$ &
$i (\partial_\mu\phi\partial_\nu\phi)
\left[\bar f_L\gamma^\mu
\overrightarrow{D}^{\nu} f_L-
\bar f_L\overleftarrow{D}^{\nu}\gamma^\mu f_L\right]$
& $+$\\
\hline
$\mathcal{O}_{{\rm R5}}$&
$\phi^2 \left[(\overline{f}_L\overleftarrow{D}^{\nu})\overrightarrow{\slashed{D}}
 (\overrightarrow{D}_{\nu}f_L)+(\overline{f}_L\overleftarrow{D}_{\nu})\overleftarrow{\slashed{D}}(\overrightarrow{D}^{\nu}f_L)
 \right]$
& $-$ \\
\hline
$\mathcal{O}_{{\rm R6}}$ &
$i \phi^2 \epsilon_{\mu\nu\rho\sigma} \left[
   \overline{f}_L\gamma^{\mu}\overrightarrow{D}^{\nu}\overrightarrow{D}^{\rho}(\overrightarrow{D}^{\sigma}f_L)-
   (\overline{f}_L\overleftarrow{D}^{\sigma})\overleftarrow{D}^{\rho}\overleftarrow{D}^{\nu}\gamma^{\mu}f_L
   \right]$
& $-$\\
\hline
$\mathcal{O}_{{\rm R7}}$ &
$
i \phi^2\epsilon_{\mu\nu\rho\sigma}\left[
   (\overline{f}_L\overleftarrow{D}^{\nu})\gamma^{\mu}\overrightarrow{D}^{\rho}(\overrightarrow{D}^{\sigma}f_L)-
   (\overline{f}_L\overleftarrow{D}^{\sigma})\overleftarrow{D}^{\rho}\gamma^{\mu}(\overrightarrow{D}^{\nu}f_L)
   \right]$
& $-$ \\
\hline
\end{tabular}
\caption{List of dimension-8 operators contributing to the $s$-wave  cross section for the
annihilation of real scalar DM into two fermions and a gauge boson.
\label{tableOR}}
\end{center}
\end{table}

To summarize, the dimension-8 operator contributing to the $s$-wave annihilation
cross section of real scalar DM into SM fermions is given by the sum of the operators in Table
\ref{tableOR}:
$\mathcal{O}_{{\rm R}}^{(8,s)}=\sum_{i=1}^7 c^{(8,s)}_{i}  \mathcal{O}_{{\rm R}i}$;
the leading interactions are therefore described in terms of only a few  operators
\begin{equation}}   \def\ee   {\end{equation}
\mathscr{L}_{\rm eff}=
{1\over \Lambda^4}\sum_{i=1}^7 c^{(8,s)}_{i} \mathcal{O}_{{\rm R}i}
+\textrm{higher-dim}\,.
\ee
In absence of $CP$-violation in the DM sector, only the first four operators
in the sum need to be considered.

\section{Annihilation cross sections}
\label{sec:xsec}

In this section we show analytical results for the annihilation cross sections due to the operators
found above.
For simplicity, we restrict ourselves to considering left-handed SM fermions only, but the results can be
easily adapted to account for annihilations into right-handed fermions as well.
We consider  the process
\begin{equation}\label{eq:MainAnnihilation}
{\rm DM}(k_1)\, {\rm DM}(k_2)\to
f_{i,L}(p_{1})\, \bar{f}_{j,L}(p_2)\, V(k),
\end{equation}
where light fermions in the final state - described here by the generic $SU(2)_L$ doublet $F=(f_1,f_2)^{T}$  - can be both leptons and quarks, and where $V=W^{\pm},Z,\gamma$.
Note that diagrams with gauge boson emission from the final state legs have to be included in order to compute a gauge invariant amplitude (see Fig.~\ref{fig:EffectiveAmplitude}).

\begin{figure}[t]
\begin{center}
  \includegraphics[width=12 cm]{EffectiveAmplitude.pdf}
  \caption{\emph{Diagrams for the annihilation process in Eq.~(\ref{eq:MainAnnihilation}).}}\label{fig:EffectiveAmplitude}
  \end{center}
\end{figure}

It is convenient to  introduce the kinematical variables $y, z$ defined by
\begin{eqnarray}
 p_1^0 &=& (1-y)\sqrt{s}/2\, ,\label{eq:p1} \\
 p_2^0 &=& (1-z)\sqrt{s}/2\, ,\label{eq:p2}\\
  k^0 &=& (y+z)\sqrt{s}/2\, , \label{eq:k}
\end{eqnarray}
which are subject to the following phase space constraints
\begin{equation}}   \def\ee   {\end{equation}
\frac{m_{V}^{2}}{s} \leq  y \leq 1\,,\quad
\frac{m_{V}^{2}}{s \, y} \leq  z \leq 1 - y + \frac{m_{V}^{2}}{s} \,,
\label{eq:phasespace}
\ee
where $s=(k_1+k_2)^2= 4 M_{\rm DM}^2 / (1 - v^2/4)$.
The scattering amplitudes will be proportional to a coefficient $A$ containing
the correct gauge couplings according to the different possible final states and their $SU(2)_L\otimes U(1)_Y$ quantum numbers; more explicitly,
\begin{equation}}   \def\ee   {\end{equation}
\begin{tabular}{rclrcl}
$A(f_{1}f_{1}\gamma)$ &=& $- \dfrac{g s_W}{2}(1+y_f)$\,,
&
$A(f_2f_2\gamma) $ &=& $+\dfrac{g s_W}{2}(1-y_f)$\,,\\
$A(f_1f_1Z)$ &=& $-\dfrac{g}{2c_W}[1-(1+y_f)s_W^2]$\,,
&
$A(f_2f_2Z)$ &=&  $+\dfrac{g}{2c_W}[1-(1-y_f)s_W^2]$\,,\\
$A(f_1f_2W^+)$ &=& $-\dfrac{g}{\sqrt{2}}$\,,
&
$A(f_2f_1W^-)$ &=& $-\dfrac{g}{\sqrt{2}}$\,,
\end{tabular}
\ee
where $Y_f=y_f/2$ (e.g. $y_L=-1$, $y_Q=1/3$), $g$ is the gauge coupling
and $s_W, c_W$ are the sin and cos of the
weak angle, respectively.

The double-differential annihilation cross section is given in terms of the
squared amplitude $\overline{|\mathcal{M}|^2}$ (averaged
over the initial spins and summer over the final ones) as
\begin{equation}
v\frac{d^2\sigma}{dy dz}={\overline{|\mathcal{M}|^2}\over 128 \, \pi^3}\,,
\end{equation}
while the total cross section is obtained by integrating over the kinematical
variables on the phase space domain defined by Eq.~(\ref{eq:phasespace}).
Let us now show the results of the computation of these cross sections
 for each of the operators found in the previous section.

\subsection{Majorana fermion DM}

In the zero-velocity approximation, $s=4M_{\rm DM}^2$ and the double-differential cross sections
for the annihilation  process (\ref{eq:MainAnnihilation}) mediated by the dimension-8  operators
$\mathcal{O}_{{\rm M}1}, \ldots, \mathcal{O}_{{\rm M}5}$
in Table \ref{tableOM} are
\begin{eqnarray}
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}1}} &=& |c_{{\rm M}1}^{(8, s)}|^{2}
\frac{A^2\, M_{\rm DM}^6} {2 \pi^{3} \Lambda^{8}}
\left[
1-y-z+\frac{m_V^2}{4M_{\rm DM}^2}
\right]\left[
y^2+z^2-\frac{m_V^2}{2M_{\rm DM}^2}
\right]\,,\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}2}}  &=&
4 \, \frac{|c_{{\rm M2}}^{(8, s)}|^{2}}{|c_{{\rm M1}}^{(8, s)}|^{2}}
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}1}}
 \,,\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}3}}  &=&
4 \, \frac{|c_{{\rm M}3}^{(8, s)}|^{2}}{|c_{{\rm M}1}^{(8, s)}|^{2}}
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}1}}
\,,\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}4}}  &=& |c_{{\rm M}4}^{(8, s)}|^{2}
\frac{2 A^2\,M_{\rm DM}^6}{\pi^{3} \Lambda^{8}} \left[
(1-y-z)(y^2+z^2)+\frac{m_V^2}{4M_{\rm DM}^2}\left[(y+z)^2-2(y+z)+2\right]
\right]
\,,\nonumber\\
&&\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}5}}  &=&
\frac{|c_{{\rm M}5}^{(8, s)}|^{2}}{|c_{{\rm M}4}^{(8, s)}|^{2}}
\left.v\frac{d^2\sigma}{dydz}\right\vert_{{\cal O}_{{\rm M}4}}
\,.
\end{eqnarray}
Notice that in the limit $m_V/M_{\rm DM}\to 0$ the differential cross sections
are the same for all operators, up to an overall numerical factor.

In the indirect searches for DM, the observables measured experimentally
are the fluxes of cosmic rays, which are directly related to the energy spectra of particles
generated by DM annihilations at the production point. By integrating the double-differential cross sections listed above once, one obtains the energy spectra of the SM fermions
and of the gauge bosons at production.
We consider the  distributions of the final fermion energy ($E_f$) and of the final gauge boson energy ($E_V$), defined as
\begin{equation}}   \def\ee   {\end{equation}
\frac{dN}{d\ln x}\equiv \frac{1}{\sigma(\textrm{DM\,DM} \to f \overline f V)}\frac{d\sigma(\textrm{DM\,DM} \to f \overline f V) }{d\ln x}\,,
\label{dNdx}
\ee
where $x\equiv E_{f, V}/M_{\rm DM}$, and shown in Fig.~\ref{fig:spectraM}.

 \begin{figure}[t]
\begin{center}
\includegraphics[width=7.5cm]{dNdx_f_M.pdf}
\hspace{0.5cm}
\includegraphics[width=7.5cm]{dNdx_V_M.pdf}
  \caption{Majorana fermion DM. \emph{The energy distributions $d N/d\ln(E/M_{\rm DM})$ of the final fermion \emph{(left panel)} and of the final gauge boson \emph{(right panel)},  for the different dimenion-8 operators.
We  set $c^{(6,p)}=c^{(8,s)}_{{\rm M}i}=1$ and $M_{\rm DM}=1$ TeV.
 }}
\label{fig:spectraM}
  \end{center}
\end{figure}

The distributions originating from the set of operators ${\cal O}_{\rm M1}, {\cal O}_{\rm M2}, {\cal O}_{\rm M3}$ differ just by an overall factor, which cancels out by normalizing the spectra, and therefore produce the same curves.
The same argument applies to the other set ${\cal O}_{\rm M4}, {\cal O}_{\rm M5}$.
However, the operator ${\cal O}_{\rm M1}$ would instead give a different result with respect to ${\cal O}_{\rm M4}$, but the difference is not visible in the fermion energy distribution (left panel of Fig.~\ref{fig:spectraM})
because of the smallness of $m_V/M_{\rm DM}$.

The dimension-6 operator ${\cal O}_{\rm M}$ mediates the  two-body annihilation
$\chi\chi \to f_L\overline{f}_L$ in $p$-wave.
The processes where a gauge boson is radiated from the final state fermion
are still in $p$-wave,
as discussed in Ref.~\cite{paper1}. Therefore, the energy spectra of fermions and gauge bosons originating from the dimension-6 operator are going to be much less important than those from the dimension-8
operators.
The total cross section for  the two-body process $\chi\chi \to f_L\overline{f}_L$ mediated by ${\cal O}_{\rm M}$ is simply
\begin{equation}}   \def\ee   {\end{equation}
\left.v\sigma(\chi\chi \to f_L\overline{f}_L)\right\vert_{{\rm M}}=  |c^{(6,p)}|^2 \frac{  M_{\rm DM}^{2}}{12 \pi \Lambda^{4}}
v^{2}\,,
\ee
while the total cross sections of the three-body processes (\ref{eq:MainAnnihilation}) mediated by
the dimension-8 operators are obtained by integrating over the full phase space
\begin{eqnarray}}   \def\eea  {\end{eqnarray}
\left.v\sigma\right\vert_{{\rm M}1}&=&|c_{{\rm M}1}^{(8, s)}|^{2} \frac{ A^{2} M_{\rm DM}^{6}}{240 \pi^{3} \Lambda^{8}} \left[ 4  - 15 \frac{m_{V}^{2}}{M_{\rm DM}^{2}} + 5 \frac{m_{V}^{4}}{M_{\rm DM}^{4}} \left( -4 + 6 \ln \frac{2 M_{\rm DM}}{m_{V}}  \right)
+\mathcal{O}\left(\frac{m_V^6}{M_{\rm DM}^6}\right) \right]\,, \\
\left.v\sigma\right\vert_{{\rm M}2}&=&
4 \,\frac{|c_{{\rm M}2}^{(8, s)}|^{2}}{|c_{{\rm M}1}^{(8, s)}|^{2}}
\left.v\sigma\right\vert_{{\rm M}1}\,, \\
\left.v\sigma\right\vert_{{\rm M}3}&=&
4 \, \frac{|c_{{\rm M}3}^{(8, s)}|^{2}}{|c_{{\rm M1}}^{(8, s)}|^{2}}
\left.v\sigma\right\vert_{{\rm M}1} \,, \\
\left.v\sigma\right\vert_{{\rm M}4}&=& |c_{{\rm M}4}^{(8, s)}|^{2} \frac{ A^{2} M_{\rm DM}^{6}}{120 \pi^{3} \Lambda^{8}} \left[8 + 15 \frac{m_{V}^{2}}{M_{\rm DM}^{2}} + \frac{m_{V}^{4}}{M_{\rm DM}^{4}}  \left( 40 - 60 \ln \frac{2 M_{\rm DM}}{m_{V}} \right)+\mathcal{O}\left(\frac{m_V^6}{M_{\rm DM}^6}\right) \right] \,,\\
\left.v\sigma\right\vert_{{\rm M}5}&=&
\frac{|c_{{\rm M}5}^{(8, s)}|^{2}}{|c_{{\rm M}4}^{(8, s)}|^{2}}
\left.v\sigma\right\vert_{{\rm M}4}\,.
\eea
For simplicity, we have only reported here the leading terms in the expansion in powers of
$m_V/M_{\rm DM}$,
but the complete analytical expressions are used in the plots.
The sub-leading terms can  be of the same order as the contributions from higher-dimensional
operators we are neglecting.

The relative importance of the $s$-wave three-body process due to dimension-8 operators with respect
to the $p$-wave two-body annihilation due to the dimenion-6 operator is captured by the
ratios of the total cross sections,  plotted in Fig.~\ref{fig:ratioM}.
Three-body cross sections can be sizeably larger than the two-body ones.
 It is evident that the dimension-8 operators dominate the total cross section,
provided  that the effective operator scale $\Lambda$ is not too large with respect to the DM mass $M_{\rm DM}$.
It is clear that limiting an EFT analysis for DM annihilations to the dimension-six operators
 misses the right result, as the cross section receives important contributions from operators
of dimension higher than six.

One could have expected this result by an order-of-magnitude estimate
 of the two-body and three-body total cross sections originating from the
operators above
\begin{eqnarray}}   \def\eea  {\end{eqnarray}
\left.v\sigma(2\to 2)\right\vert_{{\cal O}_6} &\sim&  |c^{(6,p)}|^2 {v^2}{ M_{\rm DM}^2\over \Lambda^4}\\
\left.v\sigma(2\to 3)\right\vert_{{\cal O}_8} &\sim& |c^{(8,s)}|^2 {\alpha_W\over 4\pi}
{ M_{\rm DM}^6\over \Lambda^8}\,,
\eea
where $\alpha_W=g^2/(4\pi)$,
from which we learn that the three-body cross section due to the dimension-8 operator can be bigger than
 the two-body one due to dimension-6, provided that
\begin{equation}}   \def\ee   {\end{equation}
{|c^{(6,p)}|\over |c^{(8,s)}|}
{\Lambda^2\over M_{\rm DM}^2}\lesssim {1\over v} \sqrt{{\alpha_W\over 4\pi}}\simeq 50
\label{estimate}
\ee
for $v=10^{-3}$.
This estimate is confirmed by the numerical results in Fig.~\ref{fig:ratioM}.

Let us conclude this subsection with a comment on  the toy-model studied in Ref.~\cite{paper1}, consisting of a Majorana spinor $\chi$ with mass $M_{\rm DM}$ and a scalar $S$ with mass $M_{S} > M_{\rm DM}$, in addition to the SM particle content. The added particles are a SM singlet and a $SU(2)_L$-doublet respectively. The Lagrangian is of the form
\begin{equation}
\mathcal{L} = \mathcal{L}_{\rm SM} +  \frac{1}{2} \bar{\chi} ( i \slashed{\partial} - M_{\rm DM} ) \chi +
(D_{\mu} S)^{\dagger}(D_{\mu}S) - M_S^{2} S^{\dagger} S +(y_{\chi} \bar{\chi} (F i \sigma_{2} S) + \textrm{h.c.})\,.
\label{toymodel}
\end{equation}
In the limit $M_S\gg M_{\rm DM}$, the interactions can be accurately described by the effective operator
${\cal O}_{{\rm M}1}$ in Table \ref{tableOM}.
Indeed, in  the limit $v \rightarrow 0$, the amplitudes for the process $\chi\chi\to f_L \bar{f}_L Z$
as due to the toy model  Eq.~(\ref{toymodel}) and the effective operator ${\cal O}_{{\rm M}1}$ match, provided that $4 c^{(8,s)}_{ \rm{M} 1}/\Lambda^4= -  y_\chi^2/M_S^4$.

 \begin{figure}[t]
\begin{center}
 \includegraphics[width=7.5cm]{ratiosigmaVSmDM.pdf}
 \hspace{0.5cm}
 \includegraphics[width=7.5cm]{ratiosigmaVSLambda.pdf}
  \caption{Left panel: \emph{The ratio
  of the total cross sections $\left.v\sigma(\chi\chi\to \ell_L^+\ell_L^-Z)\right\vert_{{\rm Mi}}/
  \left.v\sigma(\chi\chi \to \ell_L^+\ell_L^-)\right\vert_{{\rm M}}$
as a function of $M_{\rm DM}$. We set $c^{(6,p)}=c^{(8,s)}_i=1$, and $\Lambda=5$ TeV.
The curves corresponding to ${\cal O}_{\rm M4}, {\cal O}_{\rm M5}$ (in green) departs from those
of ${\cal O}_{\rm M2}, {\cal O}_{\rm M3}$ at low $M_{\rm DM}$ and are slightly higher.}
Right panel: \emph{The ratio
  of the total cross sections $\left.v\sigma(\chi\chi\to \ell_L^+\ell_L^-Z)\right\vert_{{\rm Mi}}/
  \left.v\sigma(\chi\chi \to \ell_L^+\ell_L^-)\right\vert_{{\rm M}}$
as a function of  $\Lambda/M_{\rm DM}$.
We  set $c^{(6,p)}=c^{(8,s)}_i=1$ and $M_{\rm DM}=1$ TeV. }}
\label{fig:ratioM}
  \end{center}
\end{figure}

\subsection{Real scalar DM}

As discussed in the previous section, there is no two-body annihilation $\phi\phi \to  f_L\overline{f}_L$
in $s$-wave
in the limit $m_f=0$. So the EW bremsstrahlung opens up a three-body annihilation channel
which is otherwise absent. The first non-vanishing contribution to the $s$-wave annihilation cross section
comes from the dimension-8 operators  in Table \ref{tableOR} and
there is no contribution from lower-dimensional operators to compare with.

For the  case of real scalar DM, we computed the double-differential cross sections
 for the annihilation  process (\ref{eq:MainAnnihilation}) mediated by the dimension-8 operators
$\mathcal{O}_{{\rm R1}}, \ldots, \mathcal{O}_{{\rm R7}}$
listed  in Table \ref{tableOR}
\begin{eqnarray}
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R1}} &=&|c_{{\rm R1}}^{(8, s)}|^{2}
\frac{4 A^2\, M_{\rm DM}^6} {\pi^{3} \Lambda^{8}}
\left[
(1-y-z) (y^{2} + z^{2}) + \frac{m_V^2}{4M_{\rm DM}^2}
\left[(y+z)^2-2(y+z)+2\right]
\right]\,,\nonumber\\
&&\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R2}} &=&|c_{{\rm R2}}^{(8, s)}|^{2}
\frac{A^2\, M_{\rm DM}^6} {\pi^{3} \Lambda^{8}}
\left[
(1-y-z) (y^{2} + z^{2}) - \frac{m_V^2}{4M_{\rm DM}^2}
\left[3(y+z)^2-2(y+z+3y z)-2\right]\right.\nonumber\\
&&\qquad\qquad\qquad\qquad
\left. - \frac{m_V^4}{8M_{\rm DM}^4} + \frac{m_V^6}{16M_{\rm DM}^6}
\right]\,,\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R3}} &=&|c_{{\rm R3}}^{(8, s)}|^{2}
\frac{\,A^2\,M_{\rm DM}^8}{\pi^{3} \Lambda^{8}}
\left[\left(1 - y - z +  \frac{m_V^2}{4 M_{\rm DM}^2} \right) \left( y^{2} + z^{2} - \frac{m_V^{2}}{2 M_{\rm DM}^{2}} \left[ 2(y + z) - 1\right] + \frac{m_V^{4}}{4 M_{\rm DM}^{4}} \right)
\right]\,,\nonumber\\
&&\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R4}} &=&
\frac{1}{16}\frac{|c_{{\rm R4}}^{(8, s)}|^{2}}{|c_{{\rm R1}}^{(8, s)}|^{2} }
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R1}}
\,, \\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R5}}  &=& |c_{{\rm R5}}^{(8, s)}|^{2}
\frac{\, A^2\,M_{\rm DM}^8}{ \pi^{3} \Lambda^{8}}
\left[\left(1 - y - z +  \frac{m_V^2}{4 M_{\rm DM}^2} \right) \left( y^{2} + z^{2} -  \frac{m_V^{2}}{2 M_{\rm DM}^{2}} \right)
\right]\,,\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R6}}  &=&
\frac{|c_{{\rm R6}}^{(8, s)}|^{2}}{|c_{{\rm R5}}^{(8, s)}|^{2}}
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R5}}
\,,\\
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R7}}  &=&
\frac{|c_{{\rm R7}}^{(8, s)}|^{2}}{|c_{{\rm R5}}^{(8, s)}|^{2}}
\left.v\frac{d^2\sigma}{dydz}\right\vert_{\textrm{R5}}
\,.
\end{eqnarray}
They are clearly symmetric under the exchange $y\leftrightarrow z$.
The energy distributions of the fermion and gauge boson, defined as in Eq.~(\ref{dNdx}),
are shown in Fig.~\ref{fig:spectraR}.
Some of the distributions (e.g. those due to the operators ${\cal O}_{\rm R1}, {\cal O}_{\rm R4}$)
 differ just by an overall factor, which cancels out by normalizing the spectra, and therefore produce the same curves. In the limit $m_V\ll M_{\rm DM}$ the distributions from all different operators become proportional
 to each other.
Therefore, because of the smallness of $m_V/M_{\rm DM}$,
the differences  in the fermion energy distribution are not visible (left panel of Fig.~\ref{fig:spectraR}).

 \begin{figure}[t]
\begin{center}
\includegraphics[width=7.5cm]{dNdx_f_R.pdf}
\hspace{0.5cm}
\includegraphics[width=7.5cm]{dNdx_V_R.pdf}
  \caption{Real scalar DM. \emph{The energy distributions $dN/d\ln(E/M_{\rm DM})$ of the final fermion \emph{(left panel)} and of the final gauge boson \emph{(right panel)},  for the different dimenion-8 operators.
 We  set $c^{(8,s)}_{{\rm R}i}=1$ and $M_{\rm DM}=1$ TeV.
}}
\label{fig:spectraR}
  \end{center}
\end{figure}
For convenience, we also report the results for the total cross sections, obtained
after integrating over the kinematical variables

\begin{eqnarray}}   \def\eea  {\end{eqnarray}
\left.v\sigma\right\vert_{{\rm R1}}&=&|c_{{\rm R1}}^{(8, s)}|^{2} \frac{ A^{2} M_{\rm DM}^{6}}{60 \pi^{3} \Lambda^{8}} \left[ 8 + 15 \frac{m_{V}^{2}}{M_{\rm DM}^{2}}  + \frac{m_{V}^{4}}{M_{\rm DM}^{4}} \left( 40 - 60 \ln \frac{2 M_{\rm DM}}{m_{V}} \right)
+\mathcal{O}\left(\frac{m_V^6}{M_{\rm DM}^6}\right) \right]\,, \\
\left.v\sigma\right\vert_{{\rm R2}}&=&|c_{{\rm R2}}^{(8, s)}|^{2} \frac{ A^{2} M_{\rm DM}^{6}}{120 \pi^{3} \Lambda^{8}} \left[ 4 + 25 \frac{m_{V}^{2}}{M_{\rm DM}^{2}}  -10 \frac{m_{V}^{4}}{M_{\rm DM}^{4}} \left( 1 + 3 \ln \frac{2 M_{\rm DM}}{m_{V}} \right)
+\mathcal{O}\left(\frac{m_V^6}{M_{\rm DM}^6}\right) \right]\,, \\
\left.v\sigma\right\vert_{{\rm R3}}&=&|c_{{\rm R3}}^{(8, s)}|^{2} \frac{ A^{2} M_{\rm DM}^{6}}{120 \pi^{3} \Lambda^{8}} \left[ 4 - 5 \frac{m_{V}^{2}}{M_{\rm DM}^{2}}  + \frac{m_{V}^{4}}{M_{\rm DM}^{4}} \left( 50 - 30 \ln \frac{2 M_{\rm DM}}{m_{V}} \right)
+\mathcal{O}\left(\frac{m_V^6}{M_{\rm DM}^6}\right) \right]\,, \\
\left.v\sigma\right\vert_{{\rm R4}}&=&{1\over 16}\frac{|c_{{\rm R2}}^{(8, s)}|^{2}}{|c_{{\rm R1}}^{(8, s)}|^{2}}
\left.v\sigma\right\vert_{{\rm R1}}\,,\\
\left.v\sigma\right\vert_{{\rm R5}}&=& |c_{{\rm R5}}^{(8, s)}|^{2} \frac{ A^{2} M_{\rm DM}^{6}}{120 \pi^{3} \Lambda^{8}} \left[ 4  - 15 \frac{m_{V}^{2}}{M_{\rm DM}^{2}} + 5 \frac{m_{V}^{4}}{M_{\rm DM}^{4}} \left( -4 + 6 \ln \frac{2 M_{\rm DM}}{m_{V}}  \right)
+\mathcal{O}\left(\frac{m_V^6}{M_{\rm DM}^6}\right)\right] \\
\left.v\sigma\right\vert_{{\rm R6}}&=&
\frac{|c_{{\rm R6}}^{(8, s)}|^{2}}{|c_{{\rm R5}}^{(8, s)}|^{2}}
\left.v\sigma\right\vert_{{\rm R5}}\,,\\
\left.v\sigma\right\vert_{{\rm R7}}&=&
\frac{|c_{{\rm R7}}^{(8, s)}|^{2}}{|c_{{\rm R5}}^{(8, s)}|^{2}}
\left.v\sigma\right\vert_{{\rm R5}}\,.
\eea
As for the Majorana case,
we have only shown the leading terms in the expansion in $m_V/M_{\rm DM}$.
The sub-leading terms can  be of the same order as the contributions from higher-dimensional
operators we are neglecting.

\section{Conclusions}
\label{sec:conclusions}

We have considered the annihilation of a self-conjugate DM particle
-- either a Majorana fermion or a real scalar -- into light SM fermions.
In these cases, the simple annihilation DM DM $\to f_L\overline{f}_L$
is helicity suppressed.
The radiation of electroweak gauge bosons lifts the suppression and opens up an
 $s$-wave channel, independently of the relative velocity of the annihilating DM particles.
 From the effective operator point of view,
we found that this effect is encoded by dimension-8 operators.

The dimension-8 operators, despite suffering from a high mass dimensionality, enjoy
no suppression by the relative velocity.
We have found all the dimension-8 operators mediating
 $s$-wave annihilations of DM into SM fermions and a gauge boson,
 and for each of them we have computed the differential and total cross sections.
 These operators encode the so-called ``virtual internal bremsstrahlung'' effect.

In the Majorana DM case, the two-body annihilation  DM DM $\to f_L\overline{f}_L$
 proceeds through a dimension-6 operator and is velocity suppressed ($p$-wave).
We have shown explicitly and quantitatively
that the dimension-8 operators actually provide a bigger contribution  to the
annihilation cross section than the dimension-6 operators, for DM masses well above
the weak scale.

For real scalar DM, the $s$-wave  annihilation into massless fermions is forbidden, and becomes
allowed  in the presence of an extra gauge boson in the final state.
We found the dimension-8 operators mediating annihilations into SM fermions and a gauge boson,
which switches on an $s$-wave annihilation process.

In conclusion, our analysis has strengthened the statement that, when dealing with the annihilation of
 heavy self-conjugate DM today,
it would be erroneous to neglect the effects of dimension-8 operators.
The results presented in this paper are of interest also when applied to phenomenological analyses, as they provide a model-independent language into which to translate the experimental data, present exclusions and possible future signals.

\section*{Acknowledgements}
We thank Duccio Pappadopulo and Riccardo Rattazzi for useful discussions.
The work of AT is supported by the
Swiss National Science Foundation under contract 200020-138131
and  by a Marie Curie Early Initial Training Network Fellowship of the European Community's Seventh Framework Programme under contract PITN-GA-2008-237920-UNILHC.
The work of AU is supported by the ERC Advanced Grant n$^{\circ}$ $267985$, ``Electroweak Symmetry Breaking, Flavour and Dark Matter: One Solution for Three Mysteries" (DaMeSyFla).