# Fermion pair radiation from accelerating classical systems

Margarita Gavrilova<EMAIL_ADDRESS>Department of Physics, LEPP, Cornell
University, Ithaca, NY 14853, USA  Mitrajyoti Ghosh<EMAIL_ADDRESS>Department of Physics, LEPP, Cornell University, Ithaca, NY 14853, USA  Yuval
Grossman<EMAIL_ADDRESS>Department of Physics, LEPP, Cornell University,
Ithaca, NY 14853, USA  Walter Tangarife<EMAIL_ADDRESS>Department of
Physics, Loyola University Chicago, Chicago, IL 60660, USA  Tien-Hsueh Tsai
<EMAIL_ADDRESS>Department of Physics, National Tsing Hua
University, Hsinchu 300, Taiwan

###### Abstract

Accelerating classical systems that couple to a fermion-antifermion pair at
the microscopic level can radiate pairs of fermions and lose energy in the
process. In this work, we derive the generalization of the Larmor formula for
fermion pair radiation. We focus on the case of a point-like classical source
in an elliptical orbit that emits fermions through vector and scalar
mediators. Ultra-light fermion emission from such systems becomes relevant
when the mass of the mediator is larger than the frequency of the periodic
motion. This enables us to probe regions of the parameter space that are
inaccessible in on-shell bosonic radiation. We apply our results to pulsar
binaries with mediators that couple to muons and neutrinos. Using current data
on binary period decays, we extract bounds on the parameters of such models.

###### Contents

  1. I Introduction
  2. II Fermion pair radiation by a point-like object
    1. II.1 General formalism
    2. II.2 Power loss formulae
  3. III Discussion of the power-loss formula
    1. III.1 General features of the power-loss formula
    2. III.2 Asymptotic behavior for the case of circular orbits
    3. III.3 Fermion-pair radiation in the SM
  4. IV Fermion pair radiation by pulsar binaries
    1. IV.1 Pulsar binaries as a classical source
    2. IV.2 Neutrino pair radiation by pulsar binaries in the SM
    3. IV.3 New physics constraints from the neutrino pair radiation by pulsar binaries
  5. V Conclusion
  6. A Derivation of the power loss formula
    1. A.1 The case of a vector boson mediator
    2. A.2 The case of the scalar mediator

##  I Introduction

Radiation by a classical system is a well-known phenomenon. Probably the most
familiar example is the radiation of electromagnetic waves by an accelerating
point-like particle. The power loss, in this case, is calculated using the
famous Larmor formula [1, 2], which, in natural units, is given by

$P_{\text{loss}}=\frac{1}{6\pi}q^{2}\bm{a}^{2},$ (1)

where $q$ is the electric charge of the particle and $\bm{a}$ is its
acceleration. The Larmor formula in Eq. (1) has been also generalized to other
types of radiation by accelerating classical sources, such as radiation of
massive vector and scalar bosons [3, 4, 5, 6, 7, 8, 9, 4].

Generalizations of the Larmor formula to exotic types of radiation are
motivated, among other things, by their applications to new physics searches.
The basic idea is that if a new physics radiation accompanies an accelerating
astrophysical object, the power loss effect can be enhanced thanks to the
large number density of an object, even if the coupling between the new
physics and the Standard Model (SM) is very small. This expected enhancement
can be used to obtain constraints on various new physics scenarios using
astrophysical observations. One example is the radiation of an ultra-light
gauged $L_{\mu}-L_{\tau}$ vector boson [10, 11, 12, 13] by pulsar binaries.
The measurement of the orbital period decay, when compared to the prediction
due to the gravitational wave (GW) radiation, was used to constrain the mass
of the $L_{\mu}-L_{\tau}$ gauge boson and its couplings to the SM [6, 5, 14,
15, 16, 17, 18].

In this paper, we extend the previous work and derive the generalization of
the Larmor formula to the case of fermion-antifermion pair radiation by
classical systems. The interest in this scenario is twofold. First, it is
interesting theoretically since it is one more example of a case where a
fermion pair behaves like a boson (other cases are Cooper pairs in
superconductors and the mediation of forces between objects via 2-fermion
forces [19, 20, 21, 22]). Thus we can study the coherent radiation of
fermions. The key point is that single-fermion emission changes the source and
thus can not be treated classically. Fermion-pair emission, however, can take
place without changing any quantum degrees of freedom of the emitting system
(such as spin). Thus, fermion-pair emission (or emission of any even number of
fermions) can be treated classically.

The second aspect is phenomenological. In particular, we consider radiation by
astrophysical systems. In the SM, as we show below, the effect of the fermion
pair radiation is negligible. In beyond the SM (BSM) theories, however, such
processes can be enhanced, enabling us to probe various new physics scenarios
using astrophysical observations. In particular, fermion-pair radiation can
become significant in models with a new light mediator (a vector or scalar
boson) that couples to some light fermionic degrees of freedom. These
fermionic degrees of freedom can be the well-known neutrinos or some new BSM
fermions. The effects of this radiation can become relevant when the mediator
is too heavy to be produced on-shell, but the fermions are much lighter and
can be radiated out. Since fermion pairs can be produced via off-shell
mediators, the fermion pair radiation can be used to probe broader regions of
the parameter space of such models.

As a particular application of our result for the fermion-pair radiation, we
consider two models: (_i_) a model with a gauged $L_{\mu}-L_{\tau}$ symmetry
and (_ii_) a model with a muonophilic scalar that couples to the muon and the
muon neutrino. We study the implications of these scenarios for the power loss
by pulsar binaries and compare our results to the cases of on-shell vector
boson radiation [6, 5, 3] and on-shell scalar radiation [3]. A stark
difference is that the emission of neutrino pairs in a particular harmonic
mode of the periodic system is not kinematically forbidden when the mediator
mass becomes larger than the frequency of that particular mode. In the case of
on-shell bosonic radiation, radiation from a harmonic mode is cut off once the
boson mass exceeds the frequency of that particular mode due to energy
conservation. We use the available period decay data for pulsar binaries to
demonstrate how neutrino pair radiation, mediated by BSM bosons, can be used
to probe a broader parameter space than the on-shell boson emission. We,
however, do not perform a comprehensive study of other bounds on the models we
consider.

This paper is organized as follows: In Sec. II, we discuss the general
machinery required for calculating fermion-pair radiation from a classical
system. In Sec. III, we discuss the main features of the power-loss formula.
In Sec. IV, we perform the computation for the particular case where the
classical system is a binary system. We then use available data to place
constraints on the parameters of a few models. We conclude in Sec. V. The
detailed calculations are shown in the appendix.

##  II Fermion pair radiation by a point-like object

In this section, we outline the calculation of the power of fermion-pair
radiation that accompanies a non-relativistic point-like object. We formulate
a general approach to the derivation of the power loss formula with a focus on
the case of elliptical orbits. The fermion pair radiation is realized in our
analysis via the coupling of the classical object to a massive boson: a
vector, or a scalar, which is unstable and decays into a fermion pair. We
consider the emission of Dirac fermions and generalize our result to the case
of Weyl fermions when we discuss the application of our result to the SM in
Section III.3. While a point-like object is a purely theoretical entity, it is
worthwhile to perform this calculation since the approximation of a radiating
extended object as a point is valid in the limit of long-wavelength radiation.

###  II.1 General formalism

We describe a point-like object as a classical source using classical current,
$J_{\text{cl}}^{\mu}(x)$ and classical density, $\rho_{\text{cl}}(x)$, which
are given by

$J_{\text{cl}}^{\mu}(x)=Q\delta^{3}({\bm{x}}-{\bm{x}}(t))u^{\mu},$ (2)
$\rho_{\text{cl}}(x)=N\delta^{3}({\bm{x}}-{\bm{x}}(t)).$ (3)

Here, $Q$ is the total charge of the object under the symmetry of interest,
$N$ is the number of the relevant microscopic constituents, ${\bm{x}}(t)$ is
its position as a function of time, $t$, and $u^{\mu}$ is its four-velocity.

Assuming motion in the $x-y$ plane, in the non-relativistic limit, the four-
velocity of the object is given by

$u^{\mu}=\left(1,\dot{x},\dot{y},0\right).$ (4)

We focus on the case of the elliptical motion in the $x-y$ plane, which can be
parametrically described by

$x=a(\cos\xi-e),\qquad y=a\sqrt{1-e^{2}}\sin\xi,\qquad\Omega\,t=\xi-e\sin\xi,$
(5)

where $e$ is the eccentricity, $a$ is the semi-major axis of the ellipse, and
$\Omega$ is the fundamental frequency of revolution. One full revolution
around the ellipse corresponds to changing the parameter $\xi$ from $0$ to
$2\pi$.

The power loss due to the fermion-pair radiation is calculated using

$P_{\text{loss}}=\int(\omega_{1}+\omega_{2})\,\mathrm{d}\Gamma,$ (6)

where $\omega_{1}$ and $\omega_{2}$ are the energies of the emitted fermion
and anti-fermion, respectively, and $\mathrm{d}\Gamma$ is the differential
rate of the fermion-pair emission. The rate depends on the type of mediator,
i.e., a scalar or a vector, and the specific form of the classical current or
density.

In general, the acceleration is not constant. In the case of periodic orbits,
the motion can be decomposed into harmonic modes with frequencies
$\Omega_{n}=n\Omega$, where $\Omega$ is the fundamental frequency of
revolution. The total emission rate can then be written as a sum of emission
rates at different harmonics $n$,

$\mathrm{d}\Gamma=\sum_{n}\mathrm{d}\Gamma_{n}\,.$ (7)

The sum goes over all kinematically allowed harmonics $n>2m_{\psi}/\Omega$,
where $m_{\psi}$ is the mass of the emitted fermions. The emission rate at
harmonic $n$ is found using

$\mathrm{d}\Gamma_{n}=\sum_{s_{1},s_{2}}|{\cal
M}_{n}(s_{1},s_{2})|^{2}(2\pi)\delta(\Omega_{n}-\omega_{1}-\omega_{2})\frac{\mathrm{d}^{3}{{\bm{k}}_{1}}}{(2\pi)^{3}\omega_{1}}\frac{\mathrm{d}^{3}{{\bm{k}}_{2}}}{(2\pi)^{3}\omega_{2}}.$
(8)

Here, $k_{1}=\left(\omega_{1},\bm{k}_{1}\right)$ and
$k_{2}=\left(\omega_{2},\bm{k}_{2}\right)$ are the four-momenta of the fermion
and anti-fermion respectively, and $s_{1}$($s_{2}$) is the spin of the fermion
(anti-fermion). The microscopic physics enters via ${\cal
M}_{n}\left(s_{1},s_{2}\right)$, which is the matrix element of the fermion-
pair emission at harmonic $n$. At leading order, this matrix element is
obtained from the diagram in Fig. 1. In the diagram, $\otimes$ denotes the
classical source, which is given by the classical current,
$J_{\text{cl}}^{\mu}(x)$, in the case of vector mediator and by the density,
$\rho_{\text{cl}}(x)$, in the case of the scalar mediated radiation.

$\psi,k_{1}$$\overline{\psi},k_{2}$ Figure 1: Feynman diagram for a fermion
pair emission by a classical current.

The total power loss via fermion-pair radiation is simply a sum of power
losses over all harmonics

$P_{\rm loss}=\sum_{n}P_{n},\qquad
P_{n}=\int(\omega_{1}+\omega_{2})\,\mathrm{d}\Gamma_{n}.$ (9)

Here, $P_{n}$ is the power loss of the $n$th harmonic.

In what follows, we consider two types of mediators: a massive gauge boson and
a massive scalar. We only consider $s$-channel exchange and remark on
$t$-channel exchange at the end of this subsection.

First, we consider a vector mediator, $A_{\mu}$, that corresponds to a broken
$\operatorname{U}(1)^{\prime}$ and has mass $m_{A}$. This gauge boson couples
to a classical current $J^{\mu}_{\text{cl}}(x)$, which has charge $Q$ under
$\operatorname{U}(1)^{\prime}$. The gauge boson $A_{\mu}$ is unstable and
decays into a fermion pair. The terms in the effective Lagrangian, relevant
for the fermion-pair radiation via $A_{\mu}$, are

${\cal L}_{\text{eff}}\supset
gA_{\mu}J_{\text{cl}}^{\mu}+gq_{\psi}\bar{\psi}\gamma^{\mu}A_{\mu}\psi\,,$
(10)

where $q_{\psi}$ is the $\operatorname{U}(1)^{\prime}$ charge of the fermion
$\psi$, $g$ is a dimensionless coupling constant, and $J_{\text{cl}}^{\mu}(x)$
is the classical current defined in Eq. (2). Both the vector boson and the
fermions are assumed to be massive with masses $m_{A}$ and $m_{\psi}$,
respectively. The leading order matrix element for the emission, at the $n-$th
harmonic, is given by

${\cal
M}_{n}(s_{1},s_{2})=g^{2}q_{\psi}\,\bar{u}(k_{1},s_{1})\gamma^{\mu}v(k_{2},s_{2})\,\frac{i(-\eta_{\mu\nu}+(k_{1}+k_{2})_{\mu}(k_{1}+k_{2})_{\nu}/m_{A}^{2})}{(k_{1}+k_{2})^{2}-m_{A}^{2}+im_{A}\Gamma_{A}}\,J_{\text{cl}}^{\nu}(\Omega_{n})\,,$
(11)

where $J^{\nu}_{\text{cl}}(\Omega_{n})$ is the Fourier transform of
$J_{\text{cl}}^{\nu}(x)$, given by

$J^{\nu}_{\text{cl}}(\Omega_{n})={\Omega\over
2\pi}\int_{0}^{2\pi/\Omega}\mathrm{d}t\int\mathrm{d}^{3}\bm{x}\ e^{i(n\Omega
t-\bm{p}\cdot\bm{x})}J^{\nu}_{\text{cl}}(x)$ (12)

with $\bm{p}=\bm{k}_{1}+\bm{k}_{2}$, $\Gamma_{A}$ is the decay width of the
gauge boson, and $2\pi/\Omega$ is the period. We assume that the decay into a
$\bar{\psi}\psi$ pair is the only decay channel for the gauge boson $A_{\mu}$,
and that the fermion mass $m_{\psi}$ is negligible compared to the gauge boson
mass $m_{A}$. Under these assumptions, the decay width of $A_{\mu}$ is given
by

$\Gamma_{A}={g^{2}q_{\psi}^{2}m_{A}\over 12\pi}.$ (13)

The other case we consider is that of a scalar mediator, $\phi$, for which the
relevant terms in the Lagrangian are

${\cal L}\supset g\phi\rho_{\text{cl}}+g^{\prime}\phi\bar{\psi}\psi,$ (14)

where $g$ is the dimensionless coupling between the scalar $\phi$ and the
classical source, $g^{\prime}$ is the Yukawa coupling of the fermion $\psi$ to
the scalar $\phi$, and $\rho_{\text{cl}}(x)$ is the number density of relevant
particles in the classical source. Both the scalar and the fermions are
assumed to be massive with masses $m_{\phi}$ and $m_{\psi}$, respectively. The
matrix element in this case is given by

${\cal
M}_{n}(s_{1},s_{2})=gg^{\prime}\bar{u}(k_{1},s_{1})v(k_{2},s_{2})\frac{i\rho_{\text{cl}}(\Omega_{n})}{(k_{1}+k_{2})^{2}-m_{\phi}^{2}+im_{\phi}\Gamma_{\phi}},$
(15)

where $\rho_{\text{cl}}(\Omega_{n})$ is the Fourier transform of
$\rho_{\text{cl}}(x)$,

$\rho_{\text{cl}}(\Omega_{n})={\Omega\over
2\pi}\int_{0}^{2\pi/\Omega}\mathrm{d}t\int\mathrm{d}^{3}\bm{x}\ e^{i(n\Omega
t-\bm{p}\cdot\bm{x})}\rho_{\text{cl}}(x),$ (16)

and the decay width of the scalar is $\Gamma_{\phi}$. As in the case of the
vector mediator, we assume that the fermionic decay mode is the only available
mode, and the fermion mass $m_{\psi}$ can be neglected compared to the mass of
a scalar $m_{\phi}$. Thus we have

$\Gamma_{\phi}=\frac{g^{\prime 2}m_{\phi}}{8\pi}.$ (17)

So far, we have only considered the $s$-channel contribution to the fermion
pair radiation. Fermion pair radiation via $t-$channel process mediated by a
vector or scalar is also a possibility. Such contributions, however, are
highly suppressed for $m_{S}\gg\Omega,m_{M}$, where $m_{S}$ is the mass of the
particles in the source that couple to the fermion pairs $\bar{\psi}\psi$ at
the microscopic level, and $m_{M}$ is the mediator mass. Since the emitted
fermions have energy of the order of $\Omega$, the fundamental frequency of
the system, the $t$-channel contribution to the momentum entering the
propagator is of the order of $m_{S}-\Omega$. Thus the $t$-channel propagator
is schematically given by

$\Pi\sim{1\over(m_{S}-\Omega)^{2}-m_{M}^{2}}.$ (18)

In the case where $m_{S}$ is much larger than both $\Omega$ and $m_{M}$, the
propagator is dominated by the mass of the source particles, and the process
is heavily suppressed. In this paper, we assume that the mass hierarchy
$m_{S}\gg\Omega,m_{M}$ and neglect the $t-$channel contributions to the
fermion pair radiation everywhere.

###  II.2 Power loss formulae

Using Eqs. (8)–(15), we can calculate the power loss via fermion-pair
radiation from a point-like object moving in an elliptical orbit. The detailed
derivations are shown in Appendix A, and here we only quote the final result.
The power loss due to fermion-pair emission in harmonic $n>2m_{\psi}/\Omega$,
for the cases of the vector and scalar mediator, can be written as

$\displaystyle P^{A}_{n}$
$\displaystyle=\frac{g^{4}q_{\psi}^{2}Q^{2}}{12\pi^{3}}\,a^{2}\Omega^{4}\,B_{n}^{A}(n_{A},n_{\psi},n_{\Gamma}),$
(19) $\displaystyle P^{\phi}_{n}$ $\displaystyle=\frac{g^{2}g^{\prime
2}N^{2}}{12\pi^{3}}\,a^{2}\Omega^{4}\,B_{n}^{\phi}(n_{\phi},n_{\psi},n_{\Gamma}).$
(20)

The functions $B_{n}^{M}(n_{A},n_{\psi},n_{\Gamma})$, where $M=A,\phi$, are
given by

$B_{n}^{M}(n_{M},n_{\psi},n_{\Gamma})\,\equiv\,\left(J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}(ne)^{2}\right)\int_{n_{\psi}}^{n-n_{\psi}}\mathrm{d}x\
F^{M}(x,n,n_{M},n_{\psi},n_{\Gamma}).$ (21)

Here

$n_{M}\equiv m_{M}/\Omega,\qquad n_{\psi}\equiv m_{\psi}/\Omega,\qquad
n_{\Gamma}\equiv\Gamma_{M}/\Omega,$ (22)

and $J_{n}(ne)$ is a Bessel function of order $n$ with argument $ne$. The
integration variable in Eq. (21) is defined by $x\equiv\omega_{1}/\Omega$,
where $\omega_{1}$ is the energy of one of the final-state fermion. In what
follows, for brevity, we use the notation

$F^{M}(x)\equiv F^{M}(x,n,n_{M},n_{\psi},n_{\Gamma}),\qquad B_{n}^{M}\equiv
B_{n}^{M}(n_{M},n_{\psi},n_{\Gamma}).$ (23)

The functions $F^{M}(x)$ have the general form

$\displaystyle F^{M}(x)=F^{M}_{0}(x)$ $\displaystyle+$
$\displaystyle{F^{M}_{1}(x)\over
n_{M}n_{\Gamma}}\left[\tan^{-1}\left({a(x)+b(x)\over
n_{M}n_{\Gamma}}\right)-\tan^{-1}\left({a(x)-b(x)\over
n_{M}n_{\Gamma}}\right)\right]$ (24) $\displaystyle+$ $\displaystyle
F^{M}_{2}(x)\tanh^{-1}\left({2a(x)b(x)\over
a(x)^{2}+b(x)^{2}+n_{M}^{2}n_{\Gamma}^{2}}\right),$

with $a(x)$ and $b(x)$ being universal for both gauge boson and scalar
mediators,

$\displaystyle a(x)$ $\displaystyle=$ $\displaystyle
2n_{\psi}^{2}-n_{M}^{2}+2nx-2x^{2}\,,$ $\displaystyle b(x)$ $\displaystyle=$
$\displaystyle 2\sqrt{x^{2}-n_{\psi}^{2}}\sqrt{(n-x)^{2}-n_{\psi}^{2}}\,.$
(25)

The functions $F_{0}^{M}(x)$, $F_{1}^{M}(x)$, and $F^{M}_{2}(x)$ are different
for the two cases. For a gauge boson mediator, we obtain

$\displaystyle F^{A}_{0}(x)$ $\displaystyle=$ $\displaystyle b(x)/2n\,,$
$\displaystyle F^{A}_{1}(x)$ $\displaystyle=$ $\displaystyle{1\over
4n}\left(n_{A}^{4}+4n^{2}n_{\psi}^{2}-n_{A}^{2}n_{\Gamma}^{2}+2n_{A}^{2}n^{2}-4nxn_{A}^{2}+4x^{2}n_{A}^{2}\right)\,,$
$\displaystyle F^{A}_{2}(x)$ $\displaystyle=$ $\displaystyle{1\over
2n}\left(n_{A}^{2}+n^{2}-2nx+2x^{2}\right)\,,$ (26)

while for a scalar mediator,

$\displaystyle F^{\phi}_{0}(x)$ $\displaystyle=$ $\displaystyle-b(x)/2n\,,$
$\displaystyle F^{\phi}_{1}(x)$ $\displaystyle=$ $\displaystyle{1\over
4n}\left(n_{\phi}^{2}n_{\Gamma}^{2}+(n^{2}-n_{\phi}^{2})(n_{\phi}^{2}-4n_{\psi}^{2})\right)\,,$
$\displaystyle F^{\phi}_{2}(x)$ $\displaystyle=$ $\displaystyle{1\over
4n}\left(n^{2}+4n_{\psi}^{2}-2n_{\phi}^{2}\right).$ (27)

Eqs. (19)–(II.2) are the main results of our work. Analytical integration of
$F^{A}(x)$ and $F^{\phi}(x)$ is challenging, but it still can be performed in
certain limits. In Sec. III.2, we consider two limiting cases: the case of
$n_{M}\ll 1$, which reproduces the Larmor formula, and $n_{M}\gg 1$, which is
relevant for the fermion pair radiation in the SM. In general, however,
calculating the power loss requires numerical analysis. We perform such an
analysis in Sec. IV when we discuss a particular phenomenological application
of our result.

##  III Discussion of the power-loss formula

The power loss due to fermion-pair emission by a classical source on an
elliptical orbit is given by Eqs. (19)-(II.2). Below we discuss the main
features and the asymptotic behavior of this result.

###  III.1 General features of the power-loss formula

We start with the general features that hold for both the vector and scalar
cases.

$\bullet$ The radiation rate is proportional to the charge-squared; that is,
the functions $P^{A}_{n}$ and $P^{\phi}_{n}$ depend on $Q^{2}$ and $N^{2}$,
respectively. This is a manifestation of the fact that the fermion-pair
radiation that we are considering is coherent.

$\bullet$ The form of $F^{M}(x)$, with $M=A,\phi$, in Eq. (24) is somewhat
general. We show in Appendix A that the overall form of $F^{M}(x)$, at the
tree level, is the same for any renormalizable theory that couples fermions to
a classical source moving in an elliptical orbit. Note that the functions
$a(x)$ and $b(x)$ defined in Eq. (II.2) are purely kinematic and thus have the
same form for any theory of fermion pair emission, while the form of
$F_{0}^{M}(x),F_{1}^{M}(x)$, and $F_{2}^{M}(x)$ vary with the theory
considered. For instance, considering non-renormalizable interactions would
lead to a different momentum dependence of the matrix element that could, in
principle, change the form of $F^{M}(x)$.

$\bullet$ The power loss for both vector and scalar mediators behaves
qualitatively the same way despite the different functional forms of
$F_{i}^{A}(x)$ vs. $F_{i}^{\phi}(x)$, with $i=0,1,2$. This is not surprising
since there is nothing fundamentally different between the matrix elements for
the vector and scalar cases.

$\bullet$ Energy conservation implies that the functions $F^{M}(x)$ are
invariant under $x\to(n-x)$ exchange. The reason is that the total energy
radiated in fermion pairs in the $n$-th harmonic is $n\Omega$. The
transformation $x\to(n-x)$ exchanges the energies of the emitted fermion and
anti-fermion, and the emission rate is the same regardless of the order in
which the integrals are carried out. This invariance results from the fact
that the fermion-antifermion emission from a classical system is essentially a
2-body decay. Note that this has nothing to do with the details of the
considered model.

$\bullet$ For $n_{A}<n$, the power loss has a very weak dependence on $n_{A}$.
This is true for the particular models that we chose here but is not expected
to be true in general. For an example when this is not the case, see the
discussion of Proca fields in Ref. [6], where dependence on $n_{A}$ appears
due to the absence of gauge symmetry.

$\bullet$ There is an interplay of three energy scales: The mass of the
mediator, $m_{M}$, the mass of the fermion, $m_{\psi}$, and the frequency of
the harmonics, $n\Omega$. The fermions cannot be produced when
$2m_{\psi}>n\Omega$. In the opposite limit, when $2m_{\psi}<n\Omega$, the
production rate depends strongly on the mediator mass. For
$m_{M}<2m_{\psi}<n\Omega$, fermion production is strongly suppressed since the
on-shell boson is kinematically forbidden from decaying into fermions. (Note
that strictly speaking, our result cannot be straightforwardly applied in this
case as everywhere we assume $\Gamma_{M}>0$.) For $2m_{\psi}<m_{M}<n\Omega$,
the fermions are produced via decay of the on-shell mediator. Thus the power
loss in the fermion-pair radiation is equal to that of the on-shell boson
radiation. The region of the parameter space where $m_{M}>n\Omega>2m_{\psi}$
is of the most interest to us, as in this region the fermions are
kinematically allowed, the mediator is off-shell, and therefore the fermion
pair emission is most significant.

$\bullet$ As an example that illustrates the qualitative features of the power
loss, consider Fig. 2. It shows $B^{A}_{n}$, defined in Eq. (21), as a
function of $n_{A}$ for massless fermions for the first four harmonics. The
most striking feature of the plots is a sharp drop at $n_{A}\sim n$. This
behavior follows from the fact that at $n_{A}\sim n$, the radiation regime
switches from the radiation dominated by on-shell boson production
($n_{A}<n$), which is proportional to $g^{2}$ to the off-shell production
($n_{A}>n$) proportional to $g^{4}$. The power loss in the regime dominated by
fermion-pair radiation is thus suppressed by $g^{2}$ compared to the power
loss in the regime dominated by the on-shell boson radiation. The power loss
in the case of the scalar mediator exhibits the same behavior.

Figure 2: $B^{A}_{n}~{}{\rm vs}~{}n_{A}$ for fixed eccentricity, $e=10^{-3}$,
coupling constant $g=10^{-15}$, and massless final state fermions,
$m_{\psi}=0$. See Eqs. (21)-(26) for the definition of $B_{n}^{A}$.

$\bullet$ Comparing our results to the cases of vector [6, 5, 3] and scalar
radiation [3], we note that from kinematic considerations alone, boson
radiation drops to zero as soon as $n_{M}=n$. This is not what we observe for
the fermion-pair emission. In the case of fermion-pair radiation, off-shell
boson production is possible, even though there is an extra suppression by
$g^{2}$ for a vector and ${g^{\prime}}^{2}$ for a scalar compared to on-shell
boson radiation. As a result, the regime $n_{M}>n$ opens up new regions of the
parameter space for each harmonic $n$ and is of particular phenomenological
interest to us.

$\bullet$ Next, we remark on the dependence of the power loss on the
eccentricity in the case of orbits close to circular. For that, we note that
the eccentricity only enters the power loss through the Bessel function
prefactor of $B^{M}_{n}$ in Eq. (21), which we denote as $K(n,e)$,

$K(n,e)=J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}(ne)^{2}\,.$ (28)

We recall that $J_{n}(z)$ and $J_{n}^{\prime}(z)$ behave asymptotically, in
the limit $z\ll 1$, as

$J_{n}(z)\approx\frac{1}{n!}\left(\frac{z}{2}\right)^{n},\qquad
J_{n}^{\prime}(z)\approx\frac{n}{n!}\frac{1}{2}\left(\frac{z}{2}\right)^{n-1}\approx\frac{n}{z}J_{n}(z),\qquad
z\ll 1.$ (29)

Using Eq. (29), we find for the eccentricity dependent prefactor $K(n,e)$, in
the limit $ne\ll 1$, that

$\displaystyle K(n,e)$
$\displaystyle=J^{\prime}_{n}(ne)^{2}+\frac{n^{2}-(ne)^{2}}{(ne)^{2}}J_{n}(z)^{2}\approx
J^{\prime}_{n}(ne)^{2}+\frac{n^{2}}{(ne)^{2}}J_{n}(z)^{2}$
$\displaystyle=2\frac{n^{2}}{z^{2}}J_{n}(ne)^{2}=\frac{(ne)^{2n-2}}{2^{2n-1}}\frac{n^{2}}{\left(n!\right)^{2}}=\frac{(ne)^{2n-2}}{2^{2n-1}\left((n-1)!\right)^{2}}.$
(30)

Thus we learn that in the limit $ne\ll 1$, prefactor $K(n,e)$ scales with the
eccentricity as

$K(n,e)\propto\left(ne\right)^{2n-2}.$ (31)

This shows that for small eccentricities (and thus orbits close to circular
ones), the contributions from higher harmonics die away very fast as $n$
increases. For $n=1$ and $e\ll 1$, we have $K(1,e)\approx 1/2$. For each
subsequent harmonic power drops by a factor of order $e^{2}$, until the
factorial in the denominator of $K(n,e)$ (see Eq. (III.1)) starts to dominate.
Then the contributions from the higher harmonics start to decay away even
faster. Fig. 2 illustrates the behavior of the power loss for the first four
harmonics in the case of small eccentricity $e=10^{-3}$.

$\bullet$ The case of highly eccentric orbits $e\sim 1$ is significantly more
involved. First, the contributions from different modes do not follow the
simple hierarchy of the low eccentricity case. The contributions from higher
modes can be of the same order or even larger than the first mode depending on
the values of other parameters. See the left panel of Fig. 3 to compare the
$n$-dependence of $B_{n}^{A}$ for different eccentricity values. Second, as
Fig. 3 demonstrates, the hierarchy of modes in the on-shell dominated part of
the parameter space does not carry into the off-shell dominated region.
Consider the green line corresponding to a highly eccentric orbit with
$e=0.6$. For $n_{A}=10^{-1}$ (left panel), the maximum contribution to the
power loss comes from the mode with $n=2$ and the first 5 modes contribute at
about the same order. The situation is drastically different for $n_{A}=50$
(right panel). The maximum contribution to the power loss comes from the $n=8$
mode. We learn that for $e\sim 1$, generally speaking, the power loss per mode
first increases as we increase $n$ and then starts decreasing after reaching a
certain value of $n$. Where this maximum occurs depends on other parameters.

Figure 3: _Left:_ $B_{n}^{A}$ as a function of $n$ in the regime where the
radiation is dominated by on-shell boson production. Different colors
correspond to different values of eccentricity. The values of $n_{\psi}$,
$n_{A}$ and $g$ are fixed. _Right:_ $B_{n}^{A}$ as a function of $n$ for a
highly eccentric orbit with $e=0.6$ in the regime where the radiation is
dominated by off-shell boson production.

###  III.2 Asymptotic behavior for the case of circular orbits

We now move to the discussion of the asymptotic behaviour of the power loss in
two limiting cases $m_{M}\ll\Omega$ and $m_{M}\gg\Omega$, where $m_{M}$ is the
mass of the mediator, $M=A,\,\phi$. In this subsection, for simplicity we
consider the straightforward case of circular orbits ($e=0$) and massless
fermions ($m_{\psi}=0$). For the eccentricity dependent part of the power
loss, $K(n,e)$, we have

$\lim_{e\,\to\,0}K(n,e)=\lim_{e\,\to\,0}\left(J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}(ne)^{2}\right)={1\over
2}\delta_{n,1}.$ (32)

Thus the only mode that contributes to the power loss in the circular orbit
limit is the mode with $n=1$.

First, let us consider the regime of light mediators, $m_{M}\ll\Omega$, or
equivalently $n_{M}\ll 1$. In this limit, $F^{M}(x)$ defined in Eq. (24) is
dominated by the second term. We thus neglect the first and the third terms of
$F^{M}(x)$ and take the second term’s limit $n_{M}\to 0$. After that, the
integral in (21) can be performed analytically, yielding the following
asymptotic expressions for the power radiated via vector and scalar,
respectively:

$P^{A}(m_{A}\ll\Omega)\approx\frac{g^{2}}{6\pi}Q^{2}a^{2}\Omega^{4},$ (33)
$P^{\phi}(m_{\phi}\ll\Omega)\approx\frac{g^{2}}{12\pi}N^{2}a^{2}\Omega^{4}.$
(34)

The asymptotic behavior that we find for $P^{A}$ and $P^{\phi}$ reproduces the
known results for the on-shell vector [6, 5, 3] and scalar [3] radiation. This
is expected as, in the regime $m_{M}\ll\Omega$, the fermion pair radiation is
dominated by on-shell boson production. Additionally, Eq. (33) also reproduces
the Larmor formula for the power of the electromagnetic wave radiation given
in Eq. (1). To see this, recall that the acceleration on a circular orbit is
equal to $a\Omega^{2}$, where $a$ is the radius of the orbit and $\Omega$ is
the frequency of revolution.

Next, we study the regime when on-shell boson production is kinematically
forbidden, and the fermion pair radiation takes place through the off-shell
mediator. This is the limit of heavy mediators, $m_{M}\gg\Omega$, or
equivalently $n_{M}\gg 1$. As in the case of the light mediators, we take the
$n_{M}\to\infty$ limit of $F^{M}(x)$ and find that the resulting expression
can be integrated analytically. Upon performing the integration, we find that
the vector and scalar-mediated radiation behave as

$P^{A}(m_{A}\gg\Omega)\approx\frac{g^{4}q_{\psi}^{2}Q^{2}}{210\pi^{3}}\frac{a^{2}\Omega^{8}}{m_{A}^{4}}=\frac{1}{35\pi^{2}}\frac{g^{2}q_{\psi}^{2}\Omega^{4}}{m_{A}^{4}}\times
P^{A}(m_{A}\ll\Omega),$ (35)
$P^{\phi}(m_{\phi}\gg\Omega)\approx\frac{g^{2}g^{\prime
2}N^{2}}{840\pi^{3}}\frac{a^{2}\Omega^{8}}{m_{\phi}^{4}}=\frac{1}{70\pi^{2}}\frac{g^{\prime
2}\Omega^{4}}{m_{\phi}^{4}}\times P^{\phi}(m_{\phi}\ll\Omega).$ (36)

We learn that in the limit of heavy mediators, the fermion pair radiation is
suppressed compared to on-shell boson radiation by the following factors:

  1. 1.

A factor of $g^{2}q_{\psi}^{2}$ or $g^{\prime 2}$, which, at the amplitude
level, comes from the coupling of the mediator to the fermion pair.

  2. 2.

A factor of $\Omega^{4}/m_{\phi}^{4}$, which comes from the propagator of the
mediator.

  3. 3.

A phase space factor of $1/35\pi^{2}$ or $1/70\pi^{2}$, which arises from the
fact that there are more particles in the final state in the case of the off-
shell pair production than in the case of the on-shell boson production.

Note that Eqs. (35) and (36) can be interpreted as integrating out the heavy
mediator, resulting in an effective 4-Fermi interaction with a coefficient
proportional to $g^{2}/m_{A}^{2}$ or $gg^{\prime}/m_{\phi}^{2}$. Thus, it is
also valid for $t$-channel and $u$-channel interactions.

Last, we compare the results of the vector to that of the scalar mediators.
Consider $m_{A}=m_{\phi}$, $Q^{2}=N^{2}$ and $g^{\prime}=gq_{\psi}$. In this
case, the power radiated via the vector mediator is greater than the power
radiated via the scalar mediator in both radiation regimes. In particular, we
have

$\frac{P^{A}(m_{A}\ll\Omega)}{P^{\phi}(m_{\phi}\ll\Omega)}\approx
2,\qquad\frac{P^{A}(m_{A}\gg\Omega)}{P^{\phi}(m_{\phi}\gg\Omega)}\approx 4.$
(37)

These factors are related to the different number of degrees of freedom
between the vector and scalar cases. There are two polarization states for an
on-shell massless vector, while the scalar has only one. For the deeply off-
shell mediator, the correspondence is not so clear, but it seems to us that it
is related to the fact that off shell gauge boson, $A_{\mu}$, has four degrees
of freedom

###  III.3 Fermion-pair radiation in the SM

The expression in Eq. (35) can be used to estimate the power loss due to
fermion pair radiation by classical sources within the SM. In this subsection,
we consider neutrino pair radiation mediated by $Z$-boson. The contribution
due to $W$-boson mediated pair emission is qualitatively the same as the
$Z$-boson contribution and is expected to be of the same order. The main
difference between the two contributions is due to the fact that $W$-boson
mediated radiation is only relevant for leptons in the source while $Z$-boson
contribution is present for all types of fermions.

Consider a source made of $N_{\Psi}$ fermions of type $\Psi$ with the total
weak charge $Q=N_{\Psi}q_{\Psi}$. To apply Eq. (35) to the neutrino pair
radiation in the SM, we need to recall that Eq. (35) was derived under the
assumption of vectorial couplings, while the SM is a chiral theory. The
relevant parts of the SM Lagrangian are different from the Lagrangian in Eq.
(10); in particular, in the SM we have

$\mathcal{L}_{\text{SM}}\supset-i\frac{g}{2\cos\theta_{W}}\left(\bar{\Psi}\gamma^{\mu}(c_{V}^{\Psi}-c_{A}^{\psi})\Psi+\bar{\nu}\gamma^{\mu}(c_{V}^{\nu}-c_{A}^{\nu})\nu\right)Z_{\mu}.$
(38)

Thus Eq. (35) yields the following expression for the $Z$-boson mediated power
loss due to the neutrino pair radiation in the SM

$P^{Z}(m_{Z}\gg\Omega)\approx\frac{1}{210\pi^{3}}\frac{g^{4}q_{\nu}^{2}q_{\Psi}^{2}N_{\Psi}^{2}}{16\cos^{4}\theta_{W}}\frac{a^{2}\Omega^{8}}{m_{Z}^{4}},$
(39)

where we perform the replacement $g\rightarrow g/(2\cos\theta_{W})$ in Eq.
(35) and define

$q_{\psi}^{2}=q_{\nu}^{2}=(c_{V}^{\nu})^{2}+(c_{A}^{\nu})^{2},\qquad
q_{\Psi}=c_{V}^{\Psi},\qquad m_{A}=m_{Z}\,.$ (40)

Note that, for the source, only vectorial coupling $c_{V}^{\Psi}$ enters the
power loss. This is because we consider coherent radiation.

The expression in Eq. (39) can be rewritten as

$P^{Z}(m_{Z}\gg\Omega)\approx
G_{\text{eff}}^{2}q_{\Psi}^{2}q_{\nu}^{2}N_{\Psi}^{2}\frac{a^{2}\Omega^{8}}{210\pi^{3}}\,,$
(41)

where $G_{\text{eff}}=\sqrt{2}G_{F}$ and $G_{F}$ is the Fermi constant. When
the power loss is written in the form of Eq. (41), it becomes clear that it is
the same as what one would obtain by performing the calculation for the
effective Fermi theory with the effective Lagrangian given by

$\mathcal{L}_{\text{eff}}^{Z}\supset
G_{\text{eff}}[\bar{\Psi}\gamma^{\mu}(c_{V}^{\Psi}-c_{A}^{\Psi}\gamma^{5})\Psi][\bar{\nu}\gamma_{\mu}(c_{V}^{\nu}-c_{A}^{\nu}\gamma^{5})\nu].$
(42)

This, of course, is not surprising as we consider radiation at the energy
$\Omega$, which is much less than the electroweak scale, $\Omega\ll m_{Z}$. In
fact, the result in Eq. (41) applies to any effective 4-Fermi interaction.
While we derive our results for $s$-channel exchange, in the limit where the
mediator is much heavier than the orbit frequency, we do not need to
distinguish between $s$-channel and $t$-channel. Thus, Eqs. (39) and (41) can
also be used for $t$-channel $W$-exchange in the SM.

Finally, we discuss the situation when there are several different types of
fermions in the source. In this case, we need to first add all the amplitudes
that correspond to the radiation by different fermions $\Psi$ (for leptons, we
add both $Z$-boson and $W$-boson contributions). Then, we square the sum of
the relevant amplitudes to obtain the total emission rate.

We end this subsection with the following remark. The power loss due to
neutrino pair radiation in the SM was estimated in Ref. [4] to be
$P^{Z}_{SM}\sim G_{F}^{2}\Omega^{6}$. Using the explicit calculation, however,
we find that $P^{Z}_{SM}\sim G_{F}^{2}a^{2}\Omega^{8}$. That is, there is an
extra factor of $a^{2}\Omega^{2}$ compared to the estimation of Ref. [4]. In
fact, our result includes the semi-major axis $a$ as an additional energy
scale of the system.

##  IV Fermion pair radiation by pulsar binaries

We now move to discuss the phenomenological applications of our results to
astrophysical systems. We focus on the neutrino-pair emission from pulsar
binaries [23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. A pulsar binary is a
binary system of a pulsar and companion. This choice is motivated by the
availability of extensive period decay data for such systems. In particular,
we apply our results to two binaries: Hulse-Taylor binary PSR B1913+16 [34,
35, 36] (a system of a pulsar and a neutron star) and PSR J1738+0333 [29, 37]
(a system of a pulsar and a white dwarf). The parameters characterizing the
two systems are summarized in Table 1.

In what follows, we first discuss the applicability of our results of Section
II.2 to pulsar binaries in general. Then we estimate the contribution to the
power loss due to neutrino pair emission in the SM and show that it is
negligible compared to the gravitational wave radiation. We then consider
neutrino pair radiation in two BSM scenarios via ultralight vector and scalar
mediators and apply our results to the pulsar binaries with the parameters in
Table 1.

Binary system |  PSR B1913+16 [36] |  PSR J1738+0333 [29]
---|---|---
Eccentricity $e$ | 0.6171340(4) | $3.4(11)\times 10^{-7}$
Pulsar mass $m_{1}$ ($M_{\odot}$) | $1.438(1)$ | $1.46(6)$
Companion mass $m_{2}$ ($M_{\odot}$) | $1.390(1)$ | $0.181(8)$
Binary period $T_{b}$ ($\text{GeV}^{-1}$) | $4.240\times 10^{28}$ | $4.657\times 10^{28}$
Intrinsic period decay $\dot{T}_{b}$ | $-2.398(4)\times 10^{-12}$ | $-2.59(32)\times 10^{-14}$
Predicted period decay due to GW $\dot{T}_{GW}$ | $-2.40263(5)\times 10^{-12}$ | $-2.77(19)\times 10^{-14}$
Ratio of period decays $\mathcal{R}=\dot{T}_{b}/\dot{T}_{GW}$ | $0.9983(16)$ | $0.94(13)$
Orbital frequency $\Omega=2\pi/T_{b}$ (GeV) | $1.482\times 10^{-28}$ | $1.349\times 10^{-28}$
Semi-major axis $a$ ($\text{GeV}^{-1}$) | $9.878\times 10^{24}$ | $8.77\times 10^{24}$
Table 1: The relevant parameters for the PSR B1913+16 and PSR J1738+0333
binary systems. Figures in parenthesis are the $1\sigma$ uncertainties in the
last quoted digit, where all the uncertainties are symmetrized. $M_{\odot}$ is
the mass of the sun. The relative experimental error of the binary period
$T_{b}$ is $\sim 10^{-12}$ for PSR B1913+16, and $\sim 10^{-11}$ for PSR
J1738+0333. The double line separates binary parameters quoted in Ref. [36,
29] and the ones we derive. Values of the semi-major axis $a$ are calculated
using Eq. (47).

###  IV.1 Pulsar binaries as a classical source

The results for the fermion pair radiation, summarized in Eqs. (19)-(II.2),
were derived for the case of classical current describing non-relativistic
point-like object following an elliptical orbit. To justify the application of
our results to pulsar binaries, we note the following:

  1. 1.

A pulsar binary can be treated as a classical source. The typical size of a
pulsar binary can be estimated as the size of the semi-major axis which varies
between $10^{6}$ and $10^{8}$ km, that is, $a\sim 10^{24}-10^{26}\text{
GeV}^{-1}$. The wavelength of the radiation is determined by the fundamental
frequency of the orbit, and for a typical pulsar binary with periods in the
range of $10^{-1}-10^{3}$ days, the wavelength is $\lambda\sim
10^{28}-10^{32}\text{ GeV}^{-1}$. Thus, $\lambda\gg a$ and we conclude that
pulsar binaries can be treated as classical radiation sources.

  2. 2.

Stars of the pulsar binary can be treated as point-like objects. Typical sizes
of stars in a binary vary from $r\sim 10\,{\text{km}}\sim 10^{19}\text{
GeV}^{-1}$, for neutron stars, and $r\sim 10^{3}\,{\text{km}}\sim
10^{21}\text{ GeV}^{-1}$, for white dwarfs. Thus $r\ll a,\,\lambda$ and both
pulsar and its companion can be treated as point-like objects. Moreover,
$r\ll\lambda$ implies the coherence of the radiation.

  3. 3.

The motion of the pulsar and its companion in the binary system is non-
relativistic. We can roughly estimate the orbital velocity of the stars in a
binary as $v\sim a\Omega$, which for characteristic values quoted above
implies $v\lesssim 10^{-2}$.

  4. 4.

For a wide range of pulsar binary systems, the observed power loss is such
that it has no significant effect on the eccentricity of the orbit. Thus we
can treat the orbit as elliptical over the time of observation. For example,
the Hulse-Taylor binary has $e\sim 1$, with
$T_{b}(\mathrm{d}e/\mathrm{d}t)\lesssim 10^{-11}$, where $T_{b}$ is the binary
period and $\mathrm{d}e/\mathrm{d}t$ is the time derivative of the
eccentricity [36].

Now that we have established that the results of Section II.2 can be applied
to pulsar binaries, we proceed in two steps. First, we modify our expressions
for the classical current and number density in Eqs. (2) and (3) to the case
of two point-like objects on an elliptical orbit. Second, we perform the
standard reduction of the two-body problem to a one-body problem.

We write the classical current and number density as

$J_{\text{cl}}^{\mu}(x)=\sum_{b=1,2}Q_{b}\,\delta^{3}({\bm{x}}-{\bm{x}}_{b}(t))u_{b}^{\mu},$
(43)

and

$\rho_{\text{cl}}(x)=\sum_{b=1,2}N_{b}\,\delta^{3}({\bm{x}}-{\bm{x}}_{b}(t)\,,$
(44)

respectively. Here, $b=1,2$ is the index that labels the stars of the binary
system, ${\bm{x}}_{b}(t)$ is the position of the $b$-th star at time $t$, and
$u_{b}^{\mu}$ is its four-velocity.

Next, we move to the binary system’s Center-of-Mass (CoM) frame. For that, we
define $\bm{R}$, the coordinate of center of mass, and $\bm{r}$, the distance
between the two stars,

$\displaystyle\bm{R}=\frac{m_{1}}{m_{1}+m_{2}}\bm{x}_{1}+\frac{m_{2}}{m_{1}+m_{2}}\bm{x}_{2},\qquad\bm{r}=\bm{x}_{1}-\bm{x}_{2}\,,$
(45)

where $m_{1}$ and $m_{2}$ are the masses of the two stars.

As we are not concerned with the translational motion of the system as a
whole, which is described by $\bm{R}$, we can solely focus on $\bm{r}$. This
is the standard two-body to one-body problem reduction for central force
motion. The non-relativistic classical trajectory of the stars in the CoM
frame can thus be described by the vector $\bm{r}=(x,y,0)$ and is given by
elliptical orbits as in Eq. (4):

$x=a(\cos\xi-e),\ \ y=a\sqrt{1-e^{2}}\sin\xi,\ \ \ \Omega t=\xi-e\sin\xi,$
(46)

where $e$ is the eccentricity, $a$ is the semi-major axis of the elliptical
orbit, and the fundamental frequency of revolution is given by

$\Omega=\sqrt{\frac{G_{N}(m_{1}+m_{2})}{a^{3}}}\,.$ (47)

The results of Eqs. (19)-(II.2) generalize to the case of binary systems via
the following replacements that follow from the 2-body to 1-body reduction
procedure:

$Q^{2}\to M^{2}\left({Q_{1}\over m_{1}}-{Q_{2}\over m_{2}}\right)^{2},\qquad
N^{2}\to M^{2}\left({N_{1}\over m_{1}}-{N_{2}\over m_{2}}\right)^{2},$ (48)

where

$M={m_{1}m_{2}\over m_{1}+m_{2}}$ (49)

is the reduced mass of the binary system. As a result we obtain the following
expressions for the power loss in $n$-th harmonic for a vector and scalar
mediators respectively:

$\displaystyle P^{A}_{n}$ $\displaystyle={g^{4}q_{\psi}^{2}\over
12\pi^{3}}M^{2}\left({Q_{1}\over m_{1}}-{Q_{2}\over
m_{2}}\right)^{2}a^{2}\Omega^{4}\,B_{n}^{A}(n_{A},n_{\psi},n_{\Gamma}),$ (50)
$\displaystyle P^{\phi}_{n}$ $\displaystyle=\frac{g^{2}g^{\prime
2}}{12\pi^{3}}M^{2}\left({N_{1}\over m_{1}}-{N_{2}\over
m_{2}}\right)^{2}a^{2}\Omega^{4}\,B_{n}^{\phi}(n_{\phi},n_{\psi},n_{\Gamma}),$
(51)

where the functions $B_{n}^{A}$ and $B_{n}^{\phi}$ are defined in Eqs.
(21)-(II.2).

###  IV.2 Neutrino pair radiation by pulsar binaries in the SM

In the SM, for the pulsar binary, the power loss via electroweak mediators is
discussed in Sec. III.3. Here, we simply generalize it to the case of 2-body
motion using Eq. (48). We obtain the following expression for the power loss
in neutrino pair radiation via $Z$-exchange in the SM

$P_{\text{SM}}\approx\frac{G_{F}^{2}\left({c_{V}^{\nu}}^{2}+{c_{A}^{\nu}}^{2}\right)}{105\pi^{3}\cos^{2}\theta_{W}}M^{2}a^{2}\Omega^{8}\left(\frac{1}{m_{1}}\sum_{i=n,p,e,\dots}c_{V}^{i}N_{1i}Q_{1i}-\frac{1}{m_{2}}\sum_{i=n,p,e,\dots}c_{V}^{i}N_{2i}Q_{2i}\right)^{2}$
(52)

where the sum goes over all microscopic constituents of binary stars, such as
neutrons ($n$), protons ($p$), electrons ($e$), etc. To perform a numerical
estimate, we consider a pulsar binary with a neutron star companion and assume
that all of the neutron star mass is in the form of neutrons. We consider a
typical pulsar-neutron star binary with

$m_{1,2}\sim M_{\odot}\sim 10^{57}\text{GeV},\qquad a\sim
10^{25}\;\text{GeV}^{-1},\qquad\Omega\sim 10^{-28}\ \text{GeV},$ (53)

and non-zero dipole moment

$M^{2}\left(\frac{Q_{1}}{m_{1}}-\frac{Q_{2}}{m_{2}}\right)^{2}\sim
Q_{1,2}^{2}\sim 10^{114},$ (54)

where $Q_{b}=N_{b}(n)-N_{b}(\bar{n})\approx N_{b}(n)\approx
M_{\odot}/m_{n}\approx 10^{57}$, with $b=1,2$, are the neutron charges of the
neutron stars, $N_{b}(n)$ and $N_{b}(\bar{n})$ are the numbers of neutrons and
anti-neutrons respectively, $m_{n}$ is the neutron mass. Using
$c_{V}^{\nu}=c_{A}^{\nu}=1/2$, $c_{V}^{n}=-1/2$, and the measured values of
$m_{n}$, $G_{F}$, and $\theta_{W}$, we find the following numerical estimate
for the radiated power

$P_{\text{SM}}\sim 10^{-56}\text{eV}^{2}.$ (55)

To see if the above result is significant, we compare it to the power loss in
the form of gravitational wave (GW) radiation. Using the quadrupole formula
for the GW radiation [38] for the case of circular orbit ($e=0$) we have

$P_{\text{GW}}=\frac{32}{5}G_{N}M^{2}a^{4}\Omega^{6}\sim
10^{8}\,\text{GeV}^{2}$ (56)

where $G_{N}$ is Newton’s gravitational constant. The rough estimates in Eqs.
(55) and (56) show that, in the SM, the fermion-pair radiation by
astrophysical objects is completely negligible compared to the gravitational
wave radiation.

We close the subsection with one remark. Within the SM, neutron stars also
emit synchrotron radiation of fermion-antifermion pairs in their self-produced
magnetic fields, as shown in Ref. [39]. This phenomenon is different from the
one we consider here. Synchrotron radiation is an incoherent effect. Thus, the
power loss, in this case, scales as $N$, the number of neutrons in the star.
In the case we are considering, the radiation is coherent and comes from the
star’s acceleration as a whole. Then, the net power that is radiated is
proportional to $N^{2}$.

###  IV.3 New physics constraints from the neutrino pair radiation by pulsar
binaries

Since extra radiation in the SM is negligible, any observed deviation from the
gravitational wave radiation would be strong evidence for the physics beyond
the SM. In particular, fermion-pair radiation can be enhanced in BSM models
with light vector or scalar mediators, with $m_{A,\phi}\ll m_{Z}$. To explain
why such light bosonic states have evaded detection so far, we must require
that they have small couplings, thus evading all the available constraints.
The smallness of couplings, however, still can be compensated in cases where
the object has a large charge under the new symmetries. This can be the case
for astrophysical objects. Thus, such objects are our prime focus in the rest
of this work.

In particular, in this subsection, we demonstrate how our results can be used
to derive new physics bounds from the neutrino pair radiation by pulsar
binaries. As we mentioned above, we use two distinct pulsar binary systems,
the Hulse-Taylor binary PSR B1913+16 and PSR J1738+0333. The relevant
properties of the two systems are summarized in Table. 1. The Hulse-Taylor
binary is a pulsar binary with a neutron star companion, it is highly
eccentric, and the mass ratio of the two stars is close to 1. The PSR
J1738+0333, on the other hand, is a pulsar-white dwarf binary with an almost
circular orbit and a high pulsar-to-companion mass ratio. For both systems,
the data on the orbital period decay is shown in Table 1. Both binaries lie
within $1\sigma$ of the general relativity prediction.

In our analysis, we exploit the fact that typical neutron stars contain a very
large number of muons, $N(\mu)\sim 10^{55}$ [40, 41, 42, 43]. Thus, the
effects of muonophilic new physics can be significantly enhanced. The presence
of the large muon number in neutron stars is attributed to the fact that when
the electron chemical potential, $\mu_{e}$, is larger than the muon mass
$\mu_{e}>m_{\mu}$, it becomes energetically favorable for relativistic
electrons at the Fermi surface to decay into muons via
$e^{-}\to\mu^{-}+\bar{\nu}_{\mu}+\nu_{e}$. Moreover, the muonic beta-decay
$n\to p+\mu^{-}+\bar{\nu}_{\mu}$ and inverse beta-decay $p+\mu^{-}\to
n+\nu_{\mu}$ reactions become energetically favorable, while the muon decay
$\mu^{-}\rightarrow e^{-}+\bar{\nu}_{e}+\nu_{\mu}$ is forbidden by Fermi
statistics.

Being motivated by the neutron star muonic content, we consider neutrino pair
emission by pulsar binaries via the following two types of BSM mediators:

  * •

$\operatorname{U}(1)_{L_{\mu}-L_{\tau}}$ massive gauge boson with

$\mathcal{L}\supset
gA_{\alpha}\left(\bar{\mu}\gamma^{\alpha}\mu-\bar{\tau}\gamma^{\alpha}\tau+\bar{\nu}_{\mu}\gamma^{\alpha}\nu_{\mu}-\bar{\nu}_{\tau}\gamma^{\alpha}\nu_{\tau}\right),$
(57)

  * •

Massive muonophilic scalar with

$\mathcal{L}\supset
g\phi\bar{\mu}\mu+g^{\prime}\phi\bar{\nu}_{\mu}\nu_{\mu}\,.$ (58)

It is known that at least two of the SM neutrinos are massive, while the third
neutrino can be very light or massless. This means that only one neutrino mass
eigenstate can be radiated in the two scenarios we consider here. A realistic
treatment of neutrino emission would include insertions of the corresponding
PMNS matrix elements [44], resulting in an additional factor of order one.
Since we already neglecting an ${\cal O}(1)$ factor coming from the estimate
of the muon number density in the neutron stars, we also ignore any PMNS
factors in the rest of this section.

Note also that in a theory with general couplings to the left and right-handed
neutrinos, i.e.,
$gA_{\alpha}\bar{\nu}\gamma^{\alpha}(c_{V}-c_{A}\gamma^{5})\nu$, the results
for the power loss are qualitatively similar. Moreover, in the case of
massless neutrinos, the power loss for the case of the general coupling is the
same as the power loss for the case of purely vectorial coupling up to
$g^{2}\rightarrow g^{2}(c_{A}^{2}+c_{V}^{2})$ replacement. This is why in what
follows, for simplicity, we consider the case of the vectorial coupling only.

These two BSM models imply the possibility for the neutrino pair radiation at
rates enhanced compared to the SM. Our results from Eqs. (50) and (51) thus
can be used to set bounds on the coupling constants and masses of the new
bosons.

The presence of the muonophilic new physics, however, not only alters the
radiation patterns of pulsar binaries, but it also has important implications
for the neutron star’s equation of state. In particular, the presence of a
repulsive (vector) or attractive (scalar) interaction between muons could
affect the muon number, which depends on the coupling $g$ to the new physics.
In the following, we write the muon number as $N(\mu,g)$ to keep the
dependence on $g$ explicit.

The number of muons becomes $g$-dependent as the interactions change the muon
chemical potential. The muon interaction due to the $L_{\mu}-L_{\tau}$ vector
boson is repulsive, and thus the chemical potential is increased compared to
its SM value by $\varepsilon\sim g^{2}N(\mu,g)/R$, where $R$ is the radius of
the neutron star the boson mass is neglected. When the coupling $g$ is small,
such that $\varepsilon\ll m_{\mu}$, the effect of the new interaction is
insignificant, and the number of muons is approximately given by its value in
the limit of no interaction $N(\mu,g=0)$. When the interaction is strong, such
that $\varepsilon\gg m_{\mu}$, it becomes energetically less favorable to have
muons inside the neutron star and thus $N(\mu,g)<N(\mu,g=0)$.

Similar reasoning applies to the case of the scalar mediator. The only
difference is the sign of the interaction. In the scalar case, the interaction
between muons is attractive. Thus the muon chemical potential is decreased by
$\varepsilon$. This leads to the increase of the muon number for larger
couplings $N(\mu,g)>N(\mu,g=0)$. In both cases, the change from the regime
when $N(\mu,g)\approx N(\mu,g=0)$ to the situation when the interaction starts
to affect the muon number happens for couplings such that $\varepsilon\sim
m_{\mu}$, or numerically $g\sim 10^{-18}$ for a typical neutron star [5].

However, in what follows, we ignore the effect of the new physics on the muon
number. Everywhere in our analysis, we use the muon number in the limit of no
new physics interaction, that is we set $N(\mu)=N(\mu,g=0)\sim 10^{55}$ [40,
41, 42, 43]. In principle, $g$-independence of muon number can be achieved in
models with both vector and scalar mediators with fine-tuned coupling
constants such that the repulsive and attractive interactions cancel each
other.

To apply Eqs. (50) and (51), we define $N_{b}(\mu)$ and $N_{b}(\bar{\mu})$ as
the number of muons and antimuons respectively in neutron star labeled by
$b=1,2$. Then, as there are almost no tau leptons in neutron stars,
$Q_{b}=N_{b}(\mu)-N_{b}(\bar{\mu})$ is the total charge of the neutron star
under the $L_{\mu}-L_{\tau}$ gauge symmetry, and
$N_{b}=N_{b}(\mu)+N_{b}(\bar{\mu})$ is the total number of muons and anti-
muons in the star. Additionally, since $N_{b}(\bar{\mu})\approx 0$, we have
$Q_{b}\approx N_{b}$.

The energy lost through radiation in a binary star system can be directly
probed by measuring the decay of the orbital period. Assuming that the
attractive gravitational force between the two stars is such that their orbits
stay Keplerian, the decay rate of the period of revolution $T_{b}$ is related
directly to the energy lost via radiation [6]:

$\dot{T}_{b}=-6\pi
a^{5/2}G_{N}^{-3/2}(m_{1}m_{2})^{-1}(m_{1}+m_{2})^{-1/2}\times
P_{\text{loss}},$ (59)

where $\dot{T}_{b}$ is the time derivative of the binary period, $G_{N}$ is
the gravitational constant, $m_{1}$ and $m_{2}$ are the masses of the stars in
the binary system, $a$ is the semi-major axis of the elliptical orbit, and
$P_{\text{loss}}$ is the total power radiated. The decay of the period per
unit of time is dimensionless and is measured experimentally.

GW emission is the dominant source of power loss in a binary star system.
Assuming that the GW emission and neutrino pair emission are the only sources
of energy loss, we have

$P_{\text{loss}}=P_{\text{GW}}+P_{\bar{\nu}\nu},$ (60)

where $P_{\bar{\nu}\nu}$ is the power loss due to the neutrino pair radiation
and $P_{GW}$ is the power loss due to GW emission, which, to the leading
order, is given by the GW quadrupole radiation formula [38],

$P_{\text{loss}}^{GW}={32\over
5}G\Omega^{6}M^{2}a^{4}(1-e^{2})^{-7/2}\left(1+{73\over 24}e^{2}+{37\over
96}e^{4}\right),$ (61)

where $M$ is the reduced mass of the system, as defined in Eq. (49). The
binary period decay $\dot{T}_{b}$ thus can be written as a sum of two
contributions,

$\dot{T}_{b}=\dot{T}_{\text{GW}}+\dot{T}_{\bar{\nu}\nu}\,.$ (62)

We next introduce the period decay ratio $\mathcal{R}$ as the ratio of the
measured period decay to the theoretical prediction of the period decay due to
GW radiation,

$\mathcal{R}=\frac{\dot{T}_{b}}{\dot{T}_{\text{GW}}}=1+\frac{\dot{T}_{\bar{\nu}\nu}}{\dot{T}_{\text{GW}}}.$
(63)

We use the measured value of $\mathcal{R}$ to set $2\sigma$ limits on the
masses and couplings of the BSM mediators of neutrino pair radiation as

$\frac{\dot{T}_{\bar{\nu}\nu}}{\dot{T}_{\text{GW}}}\leq\left(\mathcal{R}-1\right)+2\sigma\,.$
(64)

Figure 4: Left: Constraints on $g\,{\text{vs}}\,m_{A}$ from the highly
eccentric PSR B1913+16 (Hulse-Taylor) Bounds from the neutrino pair radiation
(solid) and vector boson radiation (dashed) are shown such that the region
above the curves is excluded by the measurements of the period decay. The
system parameters are taken from Table 1. Right: Constraints on $g\,{\rm
vs}\,m_{\phi}$ from PSR J1738+033. The dashed gray line corresponds to the
bound set by the emission of the scalar boson only, while the solid lines show
the bounds from including a coupling $g^{\prime}$ to the neutrinos.

The resulting constraints on the parameter space $(g,m_{A})$ and
$(g,m_{\phi})$ that we derive from the period decay data for the Hulse-Taylor
binary and PSR J1738+033 are shown in Fig. 4. When deriving the constraints,
we use $Q_{b}=N_{b}=10^{55}$ with $b=1,2$ and $q_{\nu}=1$. For the gauge boson
mediator (left panel), we calculate the period decay due to neutrino pair
emission, $\dot{T}_{\bar{\nu}\nu}$, using Eqs. (50) and (59). As we take all
three neutrinos to be massless, and as $L_{\mu}-L_{\tau}$ boson couples to two
neutrino types, there is an extra factor of $2$ in Eq. (50). Similarly, for
the case of the scalar mediator (right panel), we use Eqs. (51) and (59). As
there is no symmetry that requires equality of $g$ and $g^{\prime}$ in the
case of the scalar mediator, we present our results for the scalar case in the
$(g,m_{\phi})$ plane for four different values of $g^{\prime}$ that vary from
$10^{-7}$ to $10^{-1}$.

First, let us discuss the left panel of Fig. 4, which shows constraints on the
mass and coupling of the gauge boson. For the PSR J1738+0333 (red line), whose
orbit is very close to circular, the effect of neutrino pair radiation becomes
significant for the mediator masses greater than the second harmonic
frequency, $m_{A}>2\Omega$. For the highly eccentric Hulse-Taylor binary, off-
shell radiation dominates for $m_{A}>85\Omega$. In the region $m_{A}>2\Omega$
($m_{A}>85\Omega$) for PSR J1738+0333 (Hulse-Taylor binary), the boundary of
the excluded region is approximately quadratic in the mediator mass. This is
in stark contrast with the case of the on-shell boson emission discussed in
Ref. [6, 5, 3], where the boundary of the excluded region jumps in steps at
$m_{A}=n\Omega$, with $n$ being an integer. For comparison, the dashed lines
in Fig. 4 show the bounds due to the on-shell boson radiation.

Finally, we comment on the right panel of Fig. 4, which shows the constraints
on the mass $m_{\phi}$ and coupling $g$ for different values of $g^{\prime}$
in the case of the scalar mediated radiation. We only demonstrate the
constraints for PSR J1738+0333; the results for the Hulse-Taylor binary are
qualitatively the same. Depending on the value of $g^{\prime}$ the off-shell
scalar radiation starts to dominate for $m_{\phi}>\Omega$ ($g^{\prime}\gtrsim
10^{-4}$) or $m_{\phi}>2\Omega$ ($g^{\prime}\lesssim 10^{-8}$). As one can see
from the plot, $g^{\prime}=10^{-1}$ provides the strongest bound.

We conclude this section by noting that we do not perform a detailed analysis
of the bounds on muonophilic light states. We only remark that very strong
bounds on light states are derived from fifth force searches. Most of these
bounds do not apply in our case as these experiments are done using materials
made out of protons, neutrons, and electrons.

##  V Conclusion

It is well known that fermion pairs can behave as bosons in several
circumstances. In this work, we show that fermion pairs can also constitute
classical radiation just like bosonic states do. We use this understanding to
derive the generalization of the Larmor formula for the case of the fermion
pair emission.

Being motivated by the potential of applying fermion pair radiation to
astrophysical objects, we consider the case of classical sources following
elliptical orbits. The most interesting regime of fermion pair radiation is
when the mediator is off-shell, which takes place when the mass of the
mediator is much smaller than the frequency of the periodic motion of the
source. In this regime, the fermion pair emission takes over from on-shell
boson production. This opens up a window into a broader region of parameter
space for various models that allow for the fermion pair radiation by
classical sources.

Subsequently, we apply our results to neutrino-antineutrino emission by two
pulsar binary systems PSR B1913+16 and PSR J1738+0333. Neutrino pair emission
by binary systems is highly suppressed in the SM compared to GW radiation, but
can be significantly enhanced in various BSM scenarios. In particular, we
consider two possibilities: light muonophilic vector and scalar mediators that
couple to the SM neutrinos. Using period decay data for the two binary
systems, we derive bounds on the parameters of the two models. While we did
not perform a comprehensive study of the relevance of these bounds, the key
point is that they provide a demonstration of the fact that fermion pair
radiation can be used to enhance BSM probes using astrophysical data.

There are several future directions to go from here. Here are a few that we
find particularly interesting:

  * •

A thorough and detailed study of the bounds that we find on specific models is
called for. This, however, is complicated by the large uncertainties that come
from the estimates on the neutron star constituents. In particular, new
physics interactions alter the equation of state of a neutron star and,
currently, there is no precise quantitative understanding of how this affects
its content.

  * •

It also would be interesting to see if we can find more systems to which our
results can be applied. In particular, exotic astrophysical systems and exotic
types of new physics models.

  * •

In this work, we only consider fermion pair radiation; however, the results
can be modified to also include bosonic pair radiation. All that needs to be
done is to calculate the relevant matrix elements. It is expected to result in
a different kinematic dependence.

We conclude with the main message of our paper: If nature includes new light
states, fermion pair radiation can be one more tool in our toolbox to probe
them.

## Acknowledgements

We are grateful to Kfir Blum, Jeff Dror, Toby Opferkuch, Nadav Outmezguine,
Ira Rothstein, Ryosuke Sato, and Kohsaku Tobioka for useful discussions. The
work of YG is supported in part by the NSF grant PHY1316222. The research of
WT is supported by NSF Grant No. PHY-2013052.

##  Appendix A Derivation of the power loss formula

We present below an explicit derivation of the power loss formula for the
fermion pair radiation by a point-like classical object on an elliptical
orbit. We perform the calculation separately for the case of vector and scalar
mediators. In our calculation, we follow closely the analysis in Ref. [6].

###  A.1 The case of a vector boson mediator

The power loss is a sum over different harmonics, as given by Eqs. (8) and
(9). The matrix element, at leading order, for a vector boson mediator, is
given by Eq. (11). It includes the Fourier Transform of the classical current
$J_{\text{cl}}^{\mu}(x)$ defined in Eq. (2). We rewrite it here for
convenience:

${\cal
M}_{n}(s_{1},s_{2})=g^{2}Q_{\psi}\,\bar{u}(k_{1},s_{1})\gamma^{\mu}v(k_{2},s_{2})\,\frac{i(-\eta_{\mu\nu}+(k_{1}+k_{2})_{\mu}(k_{1}+k_{2})_{\nu}/m_{A}^{2})}{(k_{1}+k_{2})^{2}-m_{A}^{2}+im_{A}\Gamma_{A}}\,J_{\text{cl}}^{\nu}(\Omega_{n})\,,$
(65)

where $\eta_{\mu\nu}$ is the Minkowski metric tensor. Note that the
contribution from the $(k_{1}+k_{2})_{\mu}(k_{1}+k_{2})_{\nu}$ term vanishes
by means of the Dirac equation since the fermions are on-shell, that is,

$\displaystyle\bar{u}(\not{k}_{1}+\not{k}_{2})v=\bar{u}(m_{\psi}-m_{\psi})v=0.$
(66)

Squaring the amplitudes corresponding to different harmonics and summing over
spins, we find

$\displaystyle|\overline{{\cal M}_{n}}|^{2}$ $\displaystyle=$
$\displaystyle\sum_{s_{1},s_{2}}|{\cal
M}_{n}|^{2}=\frac{g^{4}Q_{\psi}^{2}}{\left((k_{1}+k_{2})^{2}-m_{A}^{2}\right)^{2}+m_{A}^{2}\Gamma_{A}^{2}}J_{\text{cl}}^{\mu}(\Omega_{n})J_{\text{cl}}^{*\nu}(\Omega_{n})\operatorname{Tr}\left[(\not{k}_{1}+m_{\nu})\gamma_{\mu}(\not{k}_{2}-m_{\nu})\gamma_{\nu}\right]$
(67) $\displaystyle=$
$\displaystyle\frac{4g^{4}Q_{\psi}^{2}}{\left((k_{1}+k_{2})^{2}-m_{A}^{2}\right)^{2}+m_{A}^{2}\Gamma_{A}^{2}}J_{\text{cl}}^{\mu}(\Omega_{n})J_{\text{cl}}^{*\nu}(\Omega_{n})\left(k_{1\mu}k_{2\nu}+k_{1\nu}k_{2\mu}-\frac{1}{2}(k_{1}+k_{2})^{2}\eta_{\mu\nu}\right).$

Finally, we are ready to write the expression for the rate of energy loss due
to $\psi\bar{\psi}$ emission at harmonic $n$ by the classical source as

$\displaystyle P_{n}=\left(\frac{dE}{dt}\right)_{n}$
$\displaystyle=\int\Omega_{n}\,d\Gamma_{n}$
$\displaystyle=\Omega_{n}\int\frac{\mathrm{d}^{3}{{\bm{k}}_{1}}}{(2\pi)^{3}(2\omega_{1})}\frac{\mathrm{d}^{3}{\bm{k}}_{2}}{(2\pi)^{3}(2\omega_{2})}(2\pi)\delta(\Omega_{n}-\omega_{1}-\omega_{2})|\overline{\mathcal{M}_{n}}|^{2}$
$\displaystyle=\Omega_{n}\int\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}\frac{|{\bm{k}}_{1}|\mathrm{d}{\omega_{1}}}{2(2\pi)^{3}}\frac{|{\bm{k}}_{2}|\mathrm{d}{\omega_{2}}}{2(2\pi)^{3}}(2\pi)\delta(\Omega_{n}-\omega_{1}-\omega_{2})|\overline{\mathcal{M}_{n}}|^{2}\,,$
(68)

where $|{\bm{k}}_{1,2}|=\sqrt{\omega_{1,2}^{2}-m_{\psi}^{2}}$, we used
$\Omega_{n}=\omega_{1}+\omega_{2}$ for the total energy carried away by the
fermion pair, $d\Phi_{1,2}$ are the differential elements of solid angles in
the fermion’s direction of flight, and
$\left|\overline{\mathcal{M}_{n}}\right|^{2}$ is given in Eq. (67). The total
power radiated is found by summing over all kinematically allowed harmonics:

$P=\sum_{n}P_{n}.$ (69)

To calculate the power radiated in fermion pairs by a point-like source in an
elliptical orbit, we need to evaluate the integrals in Eq. (A.1), after
substituting in the explicit form of $J_{\text{cl}}^{\mu}(\Omega_{n})$ in Eq.
(67). Using Eqs. (2) and (5), we find the Fourier Transform
$J_{\text{cl}}^{\mu}(\Omega_{n})$ as:

$J_{\text{cl}}^{i}{(\Omega_{n}})=a\Omega Qj^{i}_{n},\qquad
J_{\text{cl}}^{0}(\Omega_{n})=a\Omega
Q\left(\frac{{\bm{j}}_{n}\cdot{\bm{p}}}{n\Omega}\right),$ (70)

where the 3-vector ${\bm{j}}_{n}$ is defined as

${\bm{j}}_{n}=\left(-iJ_{n}^{\prime}(ne),\frac{\sqrt{1-e^{2}}}{e}J_{n}(ne),0\right),$
(71)

with $J_{n}(z)$ denoting a Bessel function, and
${\bm{p}}={\bm{k}}_{1}+{\bm{k}}_{2}$.

The terms in the numerator of $|\overline{\mathcal{M}}|^{2}$ in Eq. (67), are
then given by

$\left(J_{\text{cl}}^{\mu}(\Omega_{n})k_{1\mu}\right)\left(J_{\text{cl}}^{\nu*}(\Omega_{n})k_{2\nu}\right)=a^{2}\Omega^{2}Q^{2}j_{n}^{i}j_{n}^{j*}\left[\frac{\omega_{1}\omega_{2}}{\left(n\Omega\right)^{2}}p^{i}p^{j}-\frac{\omega_{1}}{n\Omega}p^{i}k_{2}^{j}-\frac{\omega_{2}}{n\Omega}k_{1}^{i}p^{j}+k_{1}^{i}k_{2}^{j}\right],$
(72)

and

$|J^{\mu}_{\text{cl}}(\Omega_{n})|^{2}=|J_{\text{cl}}^{0}(\Omega_{n})|^{2}-|\bm{J}_{\text{cl}}(\Omega_{n})|^{2}=a^{2}\Omega^{2}Q^{2}j_{n}^{i}j_{n}^{j*}\left[\frac{p^{i}p^{j}}{(\Omega
n)^{2}}-\delta^{ij}\right],$ (73)

where we used $\Omega_{n}=n\Omega$. Note that all quantities above are
3-vectors with Latin indices $i=1,2,3$, and a sum over $i$ and $j$ is
implicit. The expression for
$\left(J_{\text{cl}}^{\mu}(\Omega_{n})k_{2\mu}\right)\left(J_{\text{cl}}^{\nu*}(\Omega_{n})k_{1\nu}\right)$
is obtained from Eq. (72) via complex conjugation.

Next we note that the denominator of $|\overline{\mathcal{M}_{n}}|^{2}$, see
Eq. (67), depends only on $m_{A}$, $\Gamma_{A}$, $\omega_{1,2}$, the
magnitudes $|{\bm{k}}_{1,2}|$ and the relative angle between the two momenta
${\bm{k}}_{1}$, and ${\bm{k}}_{2}$ that we denote as $\gamma$. Because of
this, it is convenient to perform the change of coordinates in the integral in
Eq. (A.1) from the integration over the solid angles
$\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}$ to the integration over
$\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}^{r}$ where the solid angle of the second
neutrino is measured relative to the direction of ${\bm{k}}_{1}$, hence the
super index $r$. (Equivalently, one can also choose to integrate over
$\mathrm{d}\Phi_{1}^{r}\mathrm{d}\Phi_{2}$.) The Jacobian of this coordinate
change is unity since the transformation is simply a coordinate rotation, and
thus

$\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}=\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}^{r}.$
(74)

Defining

$\mathrm{d}\Phi_{b}=\sin\theta_{b}\mathrm{d}\theta_{b}\mathrm{d}\phi_{b},\qquad\mathrm{d}\Phi_{2}^{r}=\sin\gamma\mathrm{d}\gamma\mathrm{d}\delta,\qquad
b=1,2\,,$ (75)

we find the following relations between the two sets of integration variables

$\displaystyle\cos\gamma$ $\displaystyle=$
$\displaystyle\cos\theta_{1}\cos\theta_{2}+\sin\theta_{1}\sin\theta_{2}\cos\left(\phi_{2}-\phi_{1}\right),$
$\displaystyle\sin\delta$ $\displaystyle=$
$\displaystyle\frac{\sin\theta_{2}\sin\left(\phi_{2}-\phi_{1}\right)}{\sin\gamma}\,.$
(76)

Since, out of all the angular variables, the denominator only depends on the
relative angle $\gamma$, the integrals over $\theta_{1}$, $\phi_{1}$ and
$\delta$ can be taken easily using the following relations

$\displaystyle\int\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}k_{a}^{i}k_{a}^{j}$
$\displaystyle=$
$\displaystyle\int\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}^{r}k_{a}^{i}k_{a}^{j}=\delta^{ij}\frac{8\pi^{2}}{3}{\bm{k}}_{a}^{2}\int\sin\gamma\mathrm{d}\gamma,$
$\displaystyle\int\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}k_{1}^{i}k_{2}^{j}$
$\displaystyle=$
$\displaystyle\int\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{1}^{r}k_{1}^{i}k_{2}^{j}=\delta^{ij}\frac{8\pi^{2}}{3}({\bm{k}}_{1}\cdot{\bm{k}}_{2})\int\sin\gamma\mathrm{d}\gamma,$
$\displaystyle\int\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}$ $\displaystyle=$
$\displaystyle\int\mathrm{d}\Phi_{1}\mathrm{d}\Phi_{2}^{r}=8\pi^{2}\int\sin\gamma\mathrm{d}\gamma\,.$
(77)

Using this and the results of Eqs. (72) and (73), we perform the integration
over $\theta_{1}$, $\phi_{1}$ and $\delta$ in Eq. (A.1), and find the
following expression for the power radiated in harmonic $n$,

$\displaystyle P_{n}$
$\displaystyle=\frac{g^{4}\left(n\Omega\right)}{12\pi^{3}}a^{2}\Omega^{2}Q_{\psi}^{2}Q^{2}\left|\bm{j}_{n}\right|^{2}\int\frac{\delta(n\Omega-\omega_{1}-\omega_{2})}{\left((k_{1}+k_{2})^{2}-m_{A}^{2}\right)^{2}+m_{A}^{2}\Gamma_{A}^{2}}\times$
$\displaystyle\left[-\frac{1}{2}\left(k_{1}+k_{2}\right)^{2}\left[\left({\bm{k}}_{1}+{\bm{k}}_{2}\right)^{2}/\left(\left(n\Omega\right)^{2}-3\right)\right]+2\frac{\omega_{1}\omega_{2}}{\left(n\Omega\right)^{2}}\left({\bm{k}}_{1}+{\bm{k}}_{2}\right)^{2}\right.$
$\displaystyle\left.-2\frac{\omega_{1}}{n\Omega}\left({\bm{k}}_{2}^{2}+{\bm{k}}_{1}\cdot{\bm{k}}_{2}\right)-2{\omega_{2}\over
n\Omega}\left({\bm{k}}_{1}^{2}+{\bm{k}}_{1}\cdot{\bm{k}}_{2}\right)+2{\bm{k}}_{1}\cdot{\bm{k}}_{2}\right]\times$
$\displaystyle\omega_{1}\omega_{2}\left(1-\frac{m_{\psi}^{2}}{\omega_{1}^{2}}\right)^{1/2}\left(1-\frac{m_{\psi}^{2}}{\omega_{2}^{2}}\right)^{1/2}\sin\gamma\,\mathrm{d}\gamma\mathrm{d}\omega_{1}\mathrm{d}\omega_{2}\,,$
(78)

where the only integrals left are the integrals over $\gamma$, $\omega_{1}$
and $\omega_{2}$.

Next, we introduce the following dimensionless variables and parameters

$x_{1}=\frac{\omega_{1}}{\Omega},\quad x_{2}=\frac{\omega_{2}}{\Omega},\quad
n_{\psi}=\frac{m_{\psi}}{\Omega},\quad n_{A}=\frac{m_{A}}{\Omega},\quad
n_{\Gamma}=\frac{\Gamma_{A}}{\Omega}.$ (79)

Performing the change of variables in Eq. (A.1) from $(\omega_{1},\omega_{2})$
to $(x_{1},x_{2})$, we rewrite the expression for the power radiated in
harmonic $n$ as follows:

$P_{n}=\frac{g^{4}}{12\pi^{3}}a^{2}\Omega^{4}Q_{\psi}^{2}|{\bm{j}}_{n}|^{2}\int\sin\gamma\,\mathrm{d}\gamma\
\mathrm{d}x_{1}\,\mathrm{d}x_{2}\,\delta(n-x_{1}-x_{2})\,\mathcal{F}(\cos\gamma,x_{1},x_{2})\,.$
(80)

Upon taking the integral over $x_{2}$ and performing the replacement
$x_{1}\rightarrow x$, we obtain

$P_{n}=\frac{g^{4}}{12\pi^{3}}a^{2}\Omega^{4}Q_{\psi}^{2}Q^{2}|{\bm{j}}_{n}|^{2}\int_{n_{\psi}}^{n-n_{\psi}}\mathrm{d}x\int_{-1}^{1}\mathrm{d}(\cos\gamma)\,\mathcal{F}(\cos\gamma,x)\,,$
(81)

where function $\mathcal{F}(\cos\gamma,x)$ is given by

$\mathcal{F}(\cos\gamma,x)=\frac{b(x)}{2n}\frac{\frac{1}{2}b^{2}(x)\cos^{2}\gamma+b(x)c(x)\cos\gamma+d(x)}{\left(a(x)-b(x)\cos\gamma\right)^{2}+g^{2}}\,,$
(82)

with

$\displaystyle a(x)$ $\displaystyle=$ $\displaystyle
2n_{\psi}^{2}+2x(n-x)-n_{A}^{2}\,,$ $\displaystyle b(x)$ $\displaystyle=$
$\displaystyle 2\sqrt{x^{2}-n_{\psi}^{2}}\sqrt{(n-x)^{2}-n_{\psi}^{2}}\,,$
$\displaystyle c(x)$ $\displaystyle=$
$\displaystyle-\left(n^{2}+2n_{\psi}^{2}\right),$ $\displaystyle d(x)$
$\displaystyle=$ $\displaystyle
2(x(n^{3}-2n^{2}x+2nx^{2}-x^{3})+2n^{2}n_{\psi}^{2}+n_{\psi}^{4}),$
$\displaystyle g^{2}$ $\displaystyle=$ $\displaystyle
n_{A}^{2}n_{\Gamma}^{2}\,.$ (83)

The variable $x$ here is the ratio of the energy of one of the fermions to the
fundamental oscillation frequency. It can be at least $n_{\psi}$ or at most
$n-n_{\psi}$, hence the limits on the integral. Also note that $F$ also
depends on the parameters of the problem namely $n_{A}$, $n_{\psi}$,
$n_{\Gamma}$ defined in Eq. (79), but we do not write them explicitly for
brevity. Lastly, note that the $\gamma$-dependence of the numerator of
function ${\cal F}$ is through a term quadratic in $\cos\gamma$ and a term
linear in $\cos\gamma$. This behavior is attributed to the theory that we pick
– renormalizable theories such as in the case considered here would only
contribute at most two powers of momentum in the matrix element, leading to a
$\cos\gamma$ dependence that is at most quadratic. However non-renormalizable
theories have more momenta in the matrix element, and will give us a different
$\cos\gamma$ dependence in the ${\cal F}$.

Now, we define

$F^{A}(x)\equiv
F^{A}(n,x,n_{\psi},n_{A},n_{\Gamma})=\int_{-1}^{1}\mathrm{d}\left(\cos\gamma\right)\mathcal{F}\left(\cos\gamma,x,n\right),$
(84)

where the superscript $A$ denotes the vector boson mediator.

The integral over $\cos\gamma$ can be taken analytically. Then, we find that
the function $F^{A}(x)$, has the form:

$\displaystyle F^{A}(x)=F^{A}_{0}(x)$ $\displaystyle+$
$\displaystyle{F^{A}_{1}(x)\over
n_{M}n_{\Gamma}}\left[\tan^{-1}\left({a(x)+b(x)\over
n_{M}n_{\Gamma}}\right)-\tan^{-1}\left({a(x)-b(x)\over
n_{M}n_{\Gamma}}\right)\right]$ (85) $\displaystyle+$ $\displaystyle
F^{A}_{2}(x)\tanh^{-1}\left({2a(x)b(x)\over
a(x)^{2}+b(x)^{2}+n_{M}^{2}n_{\Gamma}^{2}}\right),$

with:

$\displaystyle F^{A}_{0}(x)$ $\displaystyle=$ $\displaystyle b(x)/2n\,,$
$\displaystyle F^{A}_{1}(x)$ $\displaystyle=$ $\displaystyle{1\over
4n}\left(n_{A}^{4}+4n^{2}n_{\psi}^{2}-n_{A}^{2}n_{\Gamma}^{2}+2n_{A}^{2}n^{2}-4nxn_{A}^{2}+4x^{2}n_{A}^{2}\right)\,,$
$\displaystyle F^{A}_{2}(x)$ $\displaystyle=$ $\displaystyle{1\over
2n}\left(n_{A}^{2}+n^{2}-2nx+2x^{2}\right)\,.$ (86)

Consequently, the power loss formula of each mode with $n>2n_{\psi}$ becomes

$\displaystyle
P_{n}=\frac{2g^{4}Q_{\psi}^{2}Q^{2}}{3(2\pi)^{3}}a^{2}\Omega^{4}\left(J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}(ne)^{2}\right)\int_{n_{\psi}}^{n-n_{\psi}}\mathrm{d}xF^{A}(x),$
(87)

which gives us Eq. (19) for the case $M=A$, where we define for mediator $M$

$B_{n}^{M}(n_{M},n_{\nu},n_{\Gamma})\,\equiv\,\left(J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}(ne)^{2}\right)\int_{n_{\psi}}^{n-n_{\psi}}\mathrm{d}x\
F^{M}(x,n,n_{M},n_{\nu},n_{\Gamma}),$ (88)

where $J_{n}(z)$ is a Bessel function of order $n$ in the variable $z$.

###  A.2 The case of the scalar mediator

The derivation for the power loss in the scalar mediator is similar to the
vector case, but the matrix element is different, as shown in Eq. (15). This
matrix element contains the number density $\rho_{\text{cl}}(x)$ of source
particles, instead of a current. As such, the difference in the calculation in
this case comes from the calculation of the squared matrix element, which in
this case, is given by:

$\displaystyle\sum_{s_{1},s_{2}}|{\cal M}_{n}(s_{1},s_{2})|^{2}$
$\displaystyle=$ $\displaystyle\frac{g^{2}g^{\prime
2}}{((k_{1}+k_{2})^{2}-m_{\phi}^{2})^{2}+m_{\phi}^{2}\Gamma_{\phi}^{2}}\operatorname{Tr}((\not{k_{1}}+m_{\psi})(\not{k_{2}}-m_{\psi}))|\rho_{\text{cl}}(\Omega_{n})|^{2}$
(89) $\displaystyle=$ $\displaystyle\frac{4g^{2}g^{\prime
2}}{((k_{1}+k_{2})^{2}-m_{\phi}^{2})^{2}+m_{\phi}^{2}\Gamma_{\phi}^{2}}(k_{1}\cdot
k_{2}-m_{\psi}^{2})|\rho_{\text{cl}}(\Omega_{n})|^{2}\,].$

The power loss is again given by Eq. (A.1).

Using Eqs. (3) and (5), we find the Fourier Transform
$\rho_{\text{cl}}^{\mu}(\Omega_{n})$ as:

$\rho_{\text{cl}}^{0}(\Omega_{n})=a\Omega
N\left(\frac{{\bm{j}}_{n}\cdot{\bm{p}}}{n\Omega}\right),$ (90)

where, like in the vector case, we define the 3-vector $j_{n}^{i}$ as follows:

${\bm{j}}_{n}=\left(-iJ_{n}^{\prime}(ne),\frac{\sqrt{1-e^{2}}}{e}J_{n}(ne),0\right),$
(91)

with $J_{n}(z)$ denoting a Bessel function, amd
${\bm{p}}={\bm{k}}_{1}+{\bm{k}}_{2}$.

After performing all the steps analogous to Eqns. (A.1)–(84) in the previous
section, i.e, after performing the angular integration, we get:

$P_{n}=\frac{g^{2}g^{\prime
2}}{12\pi^{3}}a^{2}\Omega^{4}N^{2}|{\bm{j}}_{n}|^{2}\int_{n_{\psi}}^{n-n_{\psi}}\mathrm{d}x\int_{-1}^{1}\mathrm{d}\cos\gamma\,\mathcal{F}(\cos\gamma,x)\,,$
(92)

where function $\mathcal{F}(\cos\gamma,x)$ is given by

$\mathcal{F}(\cos\gamma,x)=-\frac{b(x)}{2n}\frac{\frac{1}{2}b^{2}(x)\cos^{2}\gamma+b(x)c(x)\cos\gamma+d(x)}{\left(a(x)-b(x)\cos\gamma\right)^{2}+g^{2}}\,,$
(93)

with

$\displaystyle a(x)$ $\displaystyle=$ $\displaystyle
2n_{\psi}^{2}+2x(n-x)-n_{\phi}^{2}\,,$ $\displaystyle b(x)$ $\displaystyle=$
$\displaystyle 2\sqrt{x^{2}-n_{\psi}^{2}}\sqrt{(n-x)^{2}-n_{\psi}^{2}}\,,$
$\displaystyle c(x)$ $\displaystyle=$ $\displaystyle{(n-2x)^{2}\over 2},$
$\displaystyle d(x)$ $\displaystyle=$
$\displaystyle(n_{\psi}^{2}-nx+x^{2})(n^{2}-2n_{\psi}^{2}-2nx+2x^{2}),$
$\displaystyle g^{2}$ $\displaystyle=$ $\displaystyle
n_{\phi}^{2}n_{\Gamma}^{2}\,.$ (94)

Like before, we define

$F^{\phi}(x)\equiv
F^{\phi}(n,x,n_{\psi},n_{\phi},n_{\Gamma})=\int_{-1}^{1}\mathrm{d}\left(\cos\gamma\right)\mathcal{F}\left(\cos\gamma,x,n\right),$
(95)

where the superscript $\phi$ denotes the scalar mediator.

The integral over $\cos\gamma$ can be taken analytically to find a form for
$F^{\phi}$:

$\displaystyle F^{\phi}(x)=F^{\phi}_{0}(x)$ $\displaystyle+$
$\displaystyle{F^{\phi}_{1}(x)\over
n_{M}n_{\Gamma}}\left[\tan^{-1}\left({a(x)+b(x)\over
n_{M}n_{\Gamma}}\right)-\tan^{-1}\left({a(x)-b(x)\over
n_{M}n_{\Gamma}}\right)\right]$ (96) $\displaystyle+$ $\displaystyle
F^{\phi}_{2}(x)\tanh^{-1}\left({2a(x)b(x)\over
a(x)^{2}+b(x)^{2}+n_{M}^{2}n_{\Gamma}^{2}}\right),$

with:

$\displaystyle F^{\phi}_{0}(x)$ $\displaystyle=$ $\displaystyle-b(x)/2n\,,$
$\displaystyle F^{\phi}_{1}(x)$ $\displaystyle=$ $\displaystyle{1\over
4n}\left(n_{\phi}^{2}n_{\Gamma}^{2}+(n^{2}-n_{\phi}^{2})(n_{\phi}^{2}-4n_{\nu}^{2})\right)\,,$
$\displaystyle F^{\phi}_{2}(x)$ $\displaystyle=$ $\displaystyle{1\over
4n}\left(n^{2}+4n_{\nu}^{2}-2n_{\phi}^{2}\right).$ (97)

Consequently, the power loss formula of each mode with $n>2n_{\psi}$ becomes

$\displaystyle P_{n}=\frac{2g^{2}g^{\prime
2}}{3(2\pi)^{3}}a^{2}\Omega^{4}N^{2}\left(J^{\prime}_{n}(ne)^{2}+\frac{1-e^{2}}{e^{2}}J_{n}(ne)^{2}\right)\int_{n_{\psi}}^{n-n_{\psi}}\mathrm{d}xF^{\phi}(x),$
(98)

which gives us Eq. (20) for the case $M=\phi$

$P^{\phi}_{n}=\frac{g^{2}g^{\prime
2}}{12\pi^{3}}a^{2}\Omega^{4}\left({N_{1}\over m_{1}}-{N_{2}\over
m_{2}}\right)^{2}B_{n}^{\phi}(n_{A},n_{\nu},n_{\Gamma}).$ (99)

We find that the form of the function $F^{M}$ is general for the two types of
mediators, the difference lying in the explicit forms of the functions
$F_{0}^{M},F_{1}^{M}$ and $F_{2}^{M}$. This is due to the fact that the
$\cos\gamma$ dependence of the function ${\cal F}$ is the same in both cases,
as in both cases, the theory considered is a renormalizable one. As we
explained in the previous sub-section, this general form of $F^{M}$ is not
what we will have when we consider non-renormalizable theories that give us
higher powers of momenta in the numerator of ${\cal F}$.

## References

  * b. Larmor [1897] J. S. b. Larmor, Philosophical Magazine Series 1 44, 503 (1897).
  * Jackson [1998] J. D. Jackson,  _Classical Electrodynamics_ (Wiley, 1998).
  * Krause  _et al._ [1994] D. Krause, H. T. Kloor, and E. Fischbach, Phys. Rev. D 49, 6892 (1994).
  * Mohanty and Kumar Panda [1996] S. Mohanty and P. Kumar Panda, Phys. Rev. D 53, 5723 (1996), arXiv:hep-ph/9403205 .
  * Dror  _et al._ [2020] J. A. Dror, R. Laha, and T. Opferkuch, Phys. Rev. D 102, 023005 (2020), arXiv:1909.12845 [hep-ph] .
  * Kumar Poddar  _et al._ [2019] T. Kumar Poddar, S. Mohanty, and S. Jana, Phys. Rev. D 100, 123023 (2019), arXiv:1908.09732 [hep-ph] .
  * Huang  _et al._ [2019] J. Huang, M. C. Johnson, L. Sagunski, M. Sakellariadou, and J. Zhang, Phys. Rev. D 99, 063013 (2019), arXiv:1807.02133 [hep-ph] .
  * Kumar Poddar  _et al._ [2020] T. Kumar Poddar, S. Mohanty, and S. Jana, Phys. Rev. D 101, 083007 (2020), arXiv:1906.00666 [hep-ph] .
  * Hook and Huang [2018] A. Hook and J. Huang, JHEP 06, 036 (2018), arXiv:1708.08464 [hep-ph] .
  * Foot [1991] R. Foot, Mod. Phys. Lett. A 6, 527 (1991).
  * He  _et al._ [1991] X.-G. He, G. C. Joshi, H. Lew, and R. R. Volkas, Phys. Rev. D 44, 2118 (1991).
  * Foot  _et al._ [1994] R. Foot, X. G. He, H. Lew, and R. R. Volkas, Phys. Rev. D 50, 4571 (1994), arXiv:hep-ph/9401250 .
  * Heeck and Rodejohann [2011] J. Heeck and W. Rodejohann, Phys. Rev. D 84, 075007 (2011), arXiv:1107.5238 [hep-ph] .
  * Kopp  _et al._ [2018] J. Kopp, R. Laha, T. Opferkuch, and W. Shepherd, JHEP 11, 096 (2018), arXiv:1807.02527 [hep-ph] .
  * Alexander  _et al._ [2018] S. Alexander, E. McDonough, R. Sims, and N. Yunes, Class. Quant. Grav. 35, 235012 (2018), arXiv:1808.05286 [gr-qc] .
  * Choi and Jung [2019] H. G. Choi and S. Jung, Phys. Rev. D 99, 015013 (2019), arXiv:1810.01421 [hep-ph] .
  * Fabbrichesi and Urbano [2020] M. Fabbrichesi and A. Urbano, JCAP 06, 007 (2020), arXiv:1902.07914 [hep-ph] .
  * Seymour and Yagi [2020] B. C. Seymour and K. Yagi, Class. Quant. Grav. 37, 145008 (2020), arXiv:1908.03353 [gr-qc] .
  * Feinberg and Sucher [1968] G. Feinberg and J. Sucher, Phys. Rev. 166, 1638 (1968).
  * Hsu and Sikivie [1994] S. D. H. Hsu and P. Sikivie, Phys. Rev. D 49, 4951 (1994), arXiv:hep-ph/9211301 .
  * Ghosh  _et al._ [2020] M. Ghosh, Y. Grossman, and W. Tangarife, Phys. Rev. D 101, 116006 (2020), arXiv:1912.09444 [hep-ph] .
  * Ghosh  _et al._ [2022] M. Ghosh, Y. Grossman, W. Tangarife, X.-J. Xu, and B. Yu, (2022), arXiv:2209.07082 [hep-ph] .
  * van Straten  _et al._ [2001] W. van Straten, M. Bailes, M. C. Britton, S. R. Kulkarni, S. B. Anderson, R. N. Manchester, and J. Sarkissian, Nature 412, 158 (2001), arXiv:astro-ph/0108254 .
  * Kramer  _et al._ [2006] M. Kramer _et al._ , Science 314, 97 (2006), arXiv:astro-ph/0609417 .
  * Stairs  _et al._ [2002] I. H. Stairs, S. E. Thorsett, J. H. Taylor, and A. Wolszczan, Astrophys. J. 581, 501 (2002), arXiv:astro-ph/0208357 .
  * Shannon  _et al._ [2014] R. M. Shannon, S. Johnston, and R. N. Manchester, Mon. Not. Roy. Astron. Soc. 437, 3255 (2014), arXiv:1311.0588 [astro-ph.SR] .
  * Antoniadis  _et al._ [2013] J. Antoniadis _et al._ , Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Bhat  _et al._ [2008] N. D. R. Bhat, M. Bailes, and J. P. W. Verbiest, Phys. Rev. D 77, 124017 (2008), arXiv:0804.0956 [astro-ph] .
  * Freire  _et al._ [2012] P. C. C. Freire, N. Wex, G. Esposito-Farese, J. P. W. Verbiest, M. Bailes, B. A. Jacoby, M. Kramer, I. H. Stairs, J. Antoniadis, and G. H. Janssen, Mon. Not. Roy. Astron. Soc. 423, 3328 (2012), arXiv:1205.1450 [astro-ph.GA] .
  * Ferdman  _et al._ [2014] R. D. Ferdman _et al._ , Mon. Not. Roy. Astron. Soc. 443, 2183 (2014), arXiv:1406.5507 [astro-ph.SR] .
  * van Leeuwen  _et al._ [2015] J. van Leeuwen _et al._ , Astrophys. J. 798, 118 (2015), arXiv:1411.1518 [astro-ph.SR] .
  * Jacoby  _et al._ [2006] B. A. Jacoby, P. B. Cameron, F. A. Jenet, S. B. Anderson, R. N. Murty, and S. R. Kulkarni, Astrophys. J. Lett. 644, L113 (2006), arXiv:astro-ph/0605375 .
  * Davoudiasl and Denton [2019] H. Davoudiasl and P. B. Denton, Phys. Rev. Lett. 123, 021102 (2019), arXiv:1904.09242 [astro-ph.CO] .
  * Hulse and Taylor [1975] R. A. Hulse and J. H. Taylor, Astrophys. J. Lett. 195, L51 (1975).
  * Taylor and Weisberg [1982] J. H. Taylor and J. M. Weisberg, Astrophys. J. 253, 908 (1982).
  * Weisberg and Huang [2016] J. M. Weisberg and Y. Huang, The Astrophysical Journal 829, 55 (2016).
  * Kilic  _et al._ [2015] M. Kilic, J. Hermes, A. Gianninas, and W. R. Brown, Monthly Notices of the Royal Astronomical Society: Letters 446, L26 (2015).
  * Peters and Mathews [1963] P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963).
  * Yakovlev  _et al._ [2001] D. G. Yakovlev, A. D. Kaminker, O. Y. Gnedin, and P. Haensel, Phys. Rept. 354, 1 (2001), arXiv:astro-ph/0012122 .
  * Garani and Heeck [2019] R. Garani and J. Heeck, Phys. Rev. D 100, 035039 (2019), arXiv:1906.10145 [hep-ph] .
  * Pearson  _et al._ [2018] J. M. Pearson, N. Chamel, A. Y. Potekhin, A. F. Fantina, C. Ducoin, A. K. Dutta, and S. Goriely, Mon. Not. Roy. Astron. Soc. 481, 2994 (2018), [Erratum: Mon.Not.Roy.Astron.Soc. 486, 768 (2019)], arXiv:1903.04981 [astro-ph.HE] .
  * Harry and Hinderer [2018] I. Harry and T. Hinderer, Class. Quant. Grav. 35, 145010 (2018), arXiv:1801.09972 [gr-qc] .
  * Zhang and Chen [2001] F. Zhang and L. Chen, Chin. Phys. Lett. 18, 142 (2001), arXiv:nucl-th/0011017 .
  * Maki  _et al._ [1962] Z. Maki, M. Nakagawa, and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).