11institutetext: Department of Physics, University of Cincinnati, Cincinnati,
Ohio 45221,USA22institutetext: Department of Physics and Astronomy, Wayne
State University, Detroit, Michigan 48201, USA33institutetext: Department of
Physics and Astronomy, University of South Carolina, Columbia, South Carolina
29205, USA44institutetext: Jozef Stefan Institute, Jamova 39, Ljubljana,
Slovenia

# Enhanced neutrino polarizability

S. Bansal, 2 G. Paz 2,3 A.A. Petrov 4 M. Tammaro 1 J. Zupan
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We point out that neutrinos can have enhanced couplings to photons, if light
(pseudo)scalar mediators are present, resulting in potentially measurable
neutrino polarizability. We show that the expected suppression from small
neutrino masses can be compensated by the light mediator mass, generating
dimension 7 Rayleigh operators at low scales. We explore the rich
phenomenology of such models, computing in details the constraints on the
viable parameter space, spanned by the couplings of the mediator to neutrinos
and photons. Finally, we build several explicit models that lead to enhanced
neutrino polarizability by modifying the inverse see-saw majoron, i.e., the
pseudo-Nambu-Goldstone boson of the $U(1)_{L}$ global lepton number
responsible for generating small neutrino masses.

##  1 Introduction

Electromagnetic interactions of neutrinos serve as a primary venue for
discovering new physics interactions. It can be viewed as a qualitatively
different pathway to uncover physics beyond the standard model compared to the
observation of neutrino masses two decades ago Super-Kamiokande:1998kpq . For
instance, if neutrinos are of the Majorana type, their masses do point to a
new physical scale, $\Lambda$, since in this case the neutrino masses are
generated through a non-renormalizable dimension-5 Weinberg operator, ${\cal
L}\supset
y^{\prime}_{ij}\big{(}\bar{L}_{i}^{c}H^{c\dagger}H^{\dagger}L_{j}\big{)}/\Lambda$.
However, it is equally possible that neutrinos are of the Dirac type, in which
case the neutrino masses are due to the renormalizable Yukawa interactions,
${\cal L}\supset y_{ij}\big{(}\bar{\nu}_{Ri}H^{\dagger}L_{j}\big{)}$. To be
certain that the neutrino masses imply the existence of a new physics scale,
$\Delta L=2$ neutrinoless double $\beta$ decay needs to be discovered first,
see, e.g., Gonzalez-Garcia:2002bkq ; Gonzalez-Garcia:2007dlo ; DellOro:2016tmg
.

In contrast, if neutrinos are found to couple to photons in the current or
immediately planned experiments, this would unambiguously point to the
existence of a new physical scale. The operators of the lowest dimension,
invariant under $SU(2)_{L}\times U(1)_{Y}$, that couple neutrinos to photons
$F_{\mu\nu}$ are the dipole operators, which for Dirac neutrinos are of
dimension 6,
$\big{(}\bar{\nu}_{Ri}\sigma^{\mu\nu}H^{\dagger}L_{j}\big{)}B_{\mu\nu}/\Lambda^{2}$,
and the dimension 8 Rayleigh operators such as
$\big{(}\bar{\nu}_{Ri}HL_{j}\big{)}B_{\mu\nu}B^{\mu\nu}/\Lambda^{4}$ (similar
operators can be written for the weak isospin fields $W_{\mu\nu}^{a}$ by
direct substitutions of the weak hypercharge fields $B_{\mu\nu}$). After the
Higgs obtains a vev, $H=(0,v)/\sqrt{2}$, these operators lead to neutrino
dipole moments, $\bar{\nu}_{Ri}\sigma_{\mu\nu}\nu_{Lj}F^{\mu\nu}$, and
neutrino polarizability111In the manuscript we use interchangeably neutrino
polarizability and neutrino Rayleigh operators.,
$\bar{\nu}_{Ri}\nu_{Lj}F^{\mu\nu}F_{\mu\nu}$, respectively. The Dirac neutrino
mass term, $m_{\nu}\bar{\nu}_{L}\nu_{R}$, as well as the neutrino dipole
moments and the neutrino polarizability operators, are all chirality flipping.
The new physics that generates at some loop level the neutrino dipole moments
and/or the neutrino polarizability is, therefore, expected to generate at the
same loop level also the contributions to the neutrino masses. Unless there
are large cancellations between tree level and radiatively generated
contributions to the neutrino masses, the dipole moments and polarizability
thus need to be tiny, effectively proportional to the tiny neutrino masses,
$m_{\nu}$, and out of reach of the experiments. In this manuscript, we show
that this is not necessarily the case for Rayleigh operators, for which the
$m_{\nu}$ suppression can be parametrically compensated if the couplings to
photons arise from tree-level exchanges of light new physics.

Similar naive dimensional analysis arguments apply to Majorana neutrinos,
though with several important differences. First, if neutrinos are Majorana,
the same operators, the neutrino mass term, the dipole, and the Rayleigh
operators, require an extra Higgs insertion compared to Dirac neutrinos. That
is, for Majorana neutrinos the mass term is of dimension 5, the dipole
operators are of dimension 7,
$\big{(}\bar{L}_{i}^{c}H\sigma^{\mu\nu}H^{\dagger}L_{j}\big{)}B_{\mu\nu}$,
while Rayleigh operators are of dimension 9,
$\big{(}\bar{L}_{i}^{c}H^{c}{}^{\dagger}H^{\dagger}L_{j}\big{)}B_{\mu\nu}B^{\mu\nu}$
(and similarly for $W^{a}_{\mu\nu}$). More importantly, these operators
violate the lepton number by $\Delta L=2$. This breaking is expected to be
small, explaining why the neutrino masses are small and implying that the
neutrino magnetic moment and neutrino polarizability will be small.

There are, however, exceptions to this general rule. First of all, for
Majorana neutrinos, the tensor and scalar neutrino currents have definite
symmetry under the interchange of the neutrinos (unlike in the case of Dirac
neutrinos). Since
$\bar{\nu}_{iL}^{c}\sigma_{\mu\nu}\nu_{jL}=-\bar{\nu}_{jL}^{c}\sigma_{\mu\nu}\nu_{iL}$
is odd, while $\bar{\nu}_{iL}^{c}\nu_{jL}=\bar{\nu}_{jL}^{c}\nu_{iL}$ is even
under the interchange of the two neutrinos, any new physics that is odd under
the same flavor exchange will only contribute to the neutrino magnetic moments
and not to the neutrino masses Voloshin:1987qy . This has been used in Refs.
Babu:1989wn ; Babu:1990wv ; Babu:2020ivd ; Babu:2021jnu  to build explicit
models of enhanced neutrino magnetic moments.

No such symmetry distinguishes the neutrino mass operator from the Rayleigh
operators since the neutrino currents in both are exactly the same. Neutrino
polarizability is thus inevitably suppressed by the same small $\Delta L=2$
breaking spurion as neutrino masses. That is, neutrino polarizability is
model-independently proportional to tiny neutrino masses. However, it can
still be parametrically enhanced if generated by a tree-level exchange of a
light scalar or pseudo-scalar mediator. A prototypical example is a pseudo-
Nambu-Goldstone boson (pNGB) due to spontaneous breaking of the lepton number
– the majoron, which couples derivatively to the $\Delta L=2$ current,
$i\bar{\nu}_{L}^{c}\nu_{L}\partial_{\mu}\phi/f_{\phi}=-im_{\nu}\nu_{L}^{c}\nu_{L}\phi/f_{\phi}$.
Generically, majoron also couples to photons through a higher dimension
operator, $\phi FF/\Lambda_{\gamma}$. For the minimal majoron, this operator
is additionally suppressed by the majoron mass squared, $m_{\phi}^{2}$, while
this suppression is absent in non-minimal models. At energies below $m_{\phi}$
this then leads to the neutrino polarizability of the form $\nu\nu
FF\times(m_{\nu}/f_{\phi})\times 1/(m_{\phi}^{2}\Lambda_{\gamma})$. The small
majoron mass compensates for the $m_{\nu}$ suppression, leading to
parametrically enhanced neutrino polarizability within reach of astrophysical
and terrestrial experiments. In this manuscript, we perform the first
phenomenological analysis of the existing constraints and possible future
probes of neutrino polarizability over a wide range of mediator masses, from
eV, i.e., comparable to the neutrino masses, up to the GeV scale.

The paper is organized as follows. In Section 2, we introduce the neutrino
dipole, anapole, and polarizability operators within an EFT framework. The
enhanced neutrino polarizability via light mediator exchange is detailed in
Sec. 2.3. In Section 3, we explore the consequences of this interaction for
cosmological observables such as CMB and BBN. In Section 4, we analyze bounds
from anomalous star cooling rates due to the production of light $\phi$
particles. At higher energy scales, the Rayleigh operator can be probed with
neutrino scatterings in terrestrial experiments, including the production of
$\phi$ particles in colliders; these are discussed in Section 5. In Section 6,
we discuss UV complete models that lead to enhanced neutrino polarizability,
focusing on spontaneously broken $U(1)_{L}$. Our conclusions are summarized in
Section 7. Appendix A contains our notation and conventions, while appendix B
contains further details on the calculation of production rates of light
(pseudo)scalars in stellar cores. Appendix C contains further details on
constraints from invisible decays of heavy (pseudo)scalars.

##  2 Neutrino couplings to photons

Neutrino couplings to photons arise from higher dimensional operators. Using
the notation of Ref. Altmannshofer:2018xyo  and restricting the discussion to
low energies, well below the electroweak symmetry breaking scale, the relevant
operators are given by222The dimension six anapole moment operator induces a
contact interaction and can be replaced through the use of the equation of
motion by the four fermion operators, a choice made in the construction of the
complete basis in Ref. Altmannshofer:2018xyo . See Section 2.2 for further
details. (see also Appendix A),

$\begin{split}{\cal L}_{\rm
EFT}&\supset\sum_{i>j}\frac{C_{1,ij}^{(5)}}{\Lambda}\frac{e}{8\pi^{2}}(\bar{\nu}_{i}\sigma^{\mu\nu}P_{L}\nu_{j})F_{\mu\nu}+\frac{1}{2}\sum_{i,j}\biggr{[}\frac{C_{1,ij}^{(7)}}{\Lambda^{3}}\frac{\alpha}{12\pi}(\bar{\nu}_{i}P_{L}\nu_{j})F_{\mu\nu}F^{\mu\nu}\\\
&\qquad+\frac{C_{2,ij}^{(7)}}{\Lambda^{3}}\frac{\alpha}{8\pi}(\bar{\nu}_{i}P_{L}\nu_{j})F_{\mu\nu}\tilde{F}^{\mu\nu}\biggr{]}+{\rm
h.c.}+\cdots,\end{split}$ (1)

with ellipses denoting higher dimension terms. The indices $i,j=e,\mu,\tau$
represent the SM neutrino flavors, while $F_{\mu\nu}$ is the electromagnetic
field strength tensor, with
$\tilde{F}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}$ its
dual. Here, and in the rest of the paper, the neutrinos, $\nu_{i}$, are
assumed to be Majorana fermions. Throughout the manuscript, we also use the
four-component notation with the conventions from Ref. Dreiner:2008tw , so
that $\nu=\nu^{c}$.

The dimension 5 operators in (1) encode the neutrino dipole moments. For
Majorana neutrinos the flavor conserving dipole moments vanish because the
dipole is antisymmetric in flavor indices,
$(\bar{\nu}_{i}\sigma^{\mu\nu}P_{L}\nu_{j})=-(\bar{\nu}_{j}\sigma^{\mu\nu}P_{L}\nu_{i})$.
The dimension-7 Rayleigh operators, on the other hand, are symmetric in flavor
indices, $C_{1,ij}^{(7)}=C_{1,ji}^{(7)}$, and thus mediate also flavor
diagonal transitions. The definitions of the Wilson coefficients in (1)
include the loop factor, anticipating that in many models the operators would
be generated at one loop, while $\Lambda$ is the mass scale associated with
the masses of particles running in the loop (see also the discussion below and
in Section 6).

Below we will also use a short hand notation, where $\Lambda$ is absorbed in
the definitions of the Wilson coefficients that now become dimensionful,

$\hat{C}_{1(2),ij}^{(7)}\equiv C_{1(2),ij}^{(7)}/\Lambda^{3}.$ (2)

Quite often we will also assume that the neutrino polarizability is flavor
diagonal, so that (no summation implied)

$\hat{C}_{1(2),ij}^{(7)}=\hat{C}_{1(2),i}^{(7)}\delta_{ij},$ (3)

and similarly for dimensionless Wilson coefficients, $C_{1(2),ij}^{(7)}$. For
flavor universal case we will denote

$\hat{C}_{1(2),ij}^{(7)}=\hat{C}_{1(2)}^{(7)}\delta_{ij},$ (4)

Finally, we also define

$\hat{C}_{1(2)}^{\rm
Re}=\sum_{ij}2\operatorname{Re}\Big{[}\hat{C}_{1(2),ij}^{(7)}\Big{]}.$ (5)

In the remainder of this section we discuss in more detail the neutrino dipole
moments (Sec. 2.1), neutrino anapole moments (Sec. 2.2), and neutrino
polarizability (Sec. 2.3), including possible enhancements.

###  2.1 Neutrino dipole moments

The neutrino dipole moments are tightly constrained from the searches for
solar neutrino scatterings on electrons by Borexino, which gives at 90%CL
$\mu_{\nu}^{\rm eff}<2.8\cdot 10^{-11}\mu_{B}$ Borexino:2017fbd , while
Xenon1T sees a hint of an excess, $\mu_{\nu}^{\rm eff}\in(1.4,2.9)\cdot
10^{-11}\mu_{B}$ XENON:2020rca  (the $\mu_{\nu}^{\rm eff}$ is a linear
combination of magnetic moments that depends on flavor composition of neutrino
flux on Earth, for details see Refs. Borexino:2017fbd ; XENON:2020rca , and
also Appendix A). Interperting both measurements as bounds and taking
$C_{1,ij}^{(5)}=1$, this translates to $\Lambda\gtrsim 10^{6}$ GeV
Altmannshofer:2018xyo . While the bound on $\Lambda$ is impressive, it is
useful to compare it with the typical sizes of neutrino masses,

${\cal L}\supset-\frac{1}{2}(m_{\nu})_{ij}\bar{\nu}_{i}P_{L}\nu_{j}+{\rm
h.c.}.$ (6)

For concreteness let us assume that the neutrino dipole moments are generated
at one loop, so that parametrically

$(\lambda_{\nu})_{ij}\mu_{B}\equiv\frac{C_{1,ij}^{(5)}}{\Lambda}\frac{e}{4\pi^{2}}\sim
e\frac{y_{i}y_{j}}{16\pi^{2}}\frac{v^{2}}{M^{3}}=2.8\times
10^{-11}\,\mu_{B}\frac{y_{i}y_{j}}{\big{(}M/2.4{\rm~{}TeV}\big{)}^{3}},$ (7)

where $M$ is the typical mass of new physics particles in the loop, $y_{i}$
their couplings to neutrinos, and we included two insertions of the Higgs
electroweak vev, $v=246$ GeV, as required to project out only the neutrino
part of the electroweak leptonic doublet. Sample diagrams for the one loop
radiative corrections are shown in Fig. 1 (top right). Generically, the same
loop, but without attached photon, Fig. 1 (top left), will also contribute to
the neutrino masses

$\displaystyle m_{\nu}$
$\displaystyle\sim\frac{y_{i}y_{j}}{16\pi^{2}}\frac{v^{2}}{M}=0.05{\rm~{}eV}\,\frac{y_{i}y_{j}}{\big{(}M/7.7\cdot
10^{9}{\rm~{}TeV}\big{)}}.$ (8)

In both (7) and (8) we assumed that the $\Delta L=2$ mass insertion, denoted
with red cross in Fig. 1, is of the same size as the typical mass $M$ of the
new particles. Comparing (7) with the experimental bound, $\mu_{\nu}^{\rm
eff}<2.8\cdot 10^{-11}\mu_{B}$ Borexino:2017fbd , shows that the neutrino
magnetic moments can be large enough to be observed in the near future only if
the related radiative corrections to the neutrino masses are suppressed below
the generic expectations given by Eq. (8).

Figure 1: Sample one loop contributions to the neutrino mass matrix (top
left), dipole moments (top right) and polarizability (bottom). The Higgs vev
insertions are denoted with black cross, while the $\Delta L=2$ mass insertion
by a red cross.

Such a suppression of the neutrino masses is possible due to the Voloshin
mechanism Voloshin:1987qy , i.e., exploiting the fact that the operator
$\bar{\nu}_{i}^{c}\sigma^{\mu\nu}P_{L}\nu_{j}$ is antisymmetric, while
$\bar{\nu}_{i}^{c}\nu_{j}$ is symmetric under the exchange of flavor indices.
An explicit realization are models with approximate horizontal $SU(2)_{H}$
symmetry, in which $(\nu_{e},\nu_{\mu})$ form a doublet of $SU(2)_{H}$, while
$\nu_{\tau}$ is a singlet Babu:2020ivd . The $SU(2)_{H}$ allows a nonzero
magnetic dipole term, $\bar{\nu}_{i}\sigma^{\mu\nu}P_{L}\nu_{j}\epsilon^{ij}$,
which is a singlet under the horizontal symmetry. The neutrino mass terms
vanish in the limit of unbroken $SU(2)_{H}$,
$\bar{\nu}_{i}\nu_{j}\epsilon^{ij}=0$, due to the symmetric nature of the mass
term. The neutrino masses are thus proportional to the charged lepton masses
that break the $SU(2)_{H}$, giving rise to small enough neutrino masses
without tuning, while neutrino magnetic moments can be observably large
Babu:1989wn ; Babu:2020ivd ; Babu:1990wv ; Leurer:1989hx . The Voloshin
mechanism can be applied also to the transition dipole moments to sterile
neutrinos Brdar:2020quo .

###  2.2 Neutrino anapole moments

The anapole or toroidal moment of the neutrino is represented by a dimension 6
operator

$\ {\cal L}_{\rm
EFT}=\sum_{i,j}\frac{C_{F,ij}^{(6)}}{\Lambda^{2}}\bar{\nu}_{i}\gamma^{\mu}\gamma_{5}\nu_{j}\partial^{\nu}F_{\mu\nu},$
(9)

which has no classical analogue in the mulitipole expansion. The operator
breaks both charge, $C$, and parity, $P$, but conserves the time-reversal
symmetry $T$. This is immediately apparent in the non-relativistic limit,
where the interaction Hamiltonian is ${\cal
H}_{a}\sim\vec{\sigma}\cdot\vec{J}_{\rm em}$. The anapole moment was first
proposed by Zeldovich in 1958 Zeldovich:1958  and can be viewed as the direct
interaction between neutrino and the electromagnetic current

$\partial^{\nu}F_{\mu\nu}=J_{\mu}^{\rm
em}=\sum_{f}eQ_{f}\bar{f}\gamma_{\mu}f.$ (10)

Here, the sum runs over the SM fermions with charges $Q_{f}$ and mass smaller
then the scale $\mu\lesssim 2$ GeV, at which we define the EFT. The anapole
moment operator does not lead to an emission of a propagating photon, but
rather to a short range interaction described by dimension 6 four fermion
operators

${\cal L}_{\rm
EFT}=\sum_{i,j,f}\frac{eQ_{f}C_{F,ij}^{(6)}}{\Lambda^{2}}\,\bar{\nu}_{i}\gamma^{\mu}\gamma_{5}\nu_{j}\,\bar{f}\gamma_{\mu}f.$
(11)

That is, the anapole operator can be replaced by the sum over four-fermion
operators that are already part of the operator basis in Ref.
Altmannshofer:2018xyo , where we also refer the interested reader for the
discussion of the phenomenology of non-standard neutrino interactions due to
point-like four-fermion interactions.

The anapole moment of the neutrino is related to the neutrino charge radius
Degrassi:1989ip . Defining the effective electromagnetic form factor of the
neutrino by the relation $\langle\nu_{i}|J_{\mu}^{\rm
em}|\nu_{j}\rangle=F_{1,ij}(q^{2})\bar{u}_{i}\gamma_{\mu}P_{L}u_{j}+\ldots$,
where we do not display the $F_{2}$ term, the neutrino charge is $F_{1}(0)=0$,
while its effective mean-square charge radius is

$\langle r^{2}\rangle_{ij}=6\frac{\partial F_{1,ij}(q^{2})}{\partial
q^{2}}\Big{|}_{q^{2}=0}.$ (12)

Evaluating the single photon exchange contribution to the scattering of
charged SM fermions on neutrinos, the $q^{2}$ factor in the
$F_{1}^{\prime}(0)q^{2}$ term cancels the $1/q^{2}$ pole, and results in a
contact contribution of the form (11). The neutrino charge radius is therefore
directly proportional to the neutrino anapole moment

$\langle r^{2}\rangle_{ij}=6\frac{C_{F,ij}^{(6)}}{\Lambda^{2}}.$ (13)

###  2.3 Light scalar mediator model for enhanced neutrino polarizability

In generic new physics models the neutrino polarizability will be highly
suppressed. For instance, if the dimension 7 Rayleigh operators in (1) result
from heavy particles running in a loop, Fig. 1 (bottom), and if we assume that
the neutrino masses are dominated by a similar loop without photons attached,
Fig. 1 (top left), the NDA expectation is

$\frac{C_{1(2),ij}^{(7)}}{\Lambda^{3}}\sim
y_{i}y_{j}\frac{v^{2}}{M^{5}}\sim\frac{m_{\nu}}{M^{4}},$ (14)

where in the last estimate we used the relation (8). The searches for new
charged particles at the LEP and LHC requires
$M\gtrsim{\mathcal{O}}(100~{}{\rm GeV})$. This gives an NDA estimate for the
neutrino polarizability that is orders of magnitudes below the present and
future experimental sensitivities, see Table 1.

Figure 2: The tree-level $\phi$ exchange that leads to neutrino polarizability
once $\phi$ is integrated out.

The crucial assumption in the NDA estimate (14) was that all the relevant new
physics is heavy. If this is not the case, the effective scale $\Lambda$ can
be significantly lower. A simple example is a model with a light scalar
mediator, $\phi$, that couples to neutrinos, and, through dimension 5
operators, also to photons,

${\cal L}_{\rm int}\supset-\frac{\alpha}{12\pi}\frac{c_{\gamma}}{f_{\phi}}\phi
F_{\mu\nu}F^{\mu\nu}-\frac{\alpha}{8\pi}\frac{c_{\gamma}^{\prime}}{f_{\phi}}\phi
F_{\mu\nu}\tilde{F}^{\mu\nu}+\frac{1}{2}c_{\nu}^{ij}\left(\bar{\nu}_{i}P_{L}\nu_{j}\right)\phi+{\rm
h.c.}.$ (15)

Here, $c_{\gamma},c_{\gamma}^{\prime}$ and $c_{\nu}^{ij}$ are dimensionless
couplings, while $f_{\phi}$ is the UV scale (for a pNGB $\phi$ the $f_{\phi}$
is related to the scale of spontaneous symmetry breaking, see Section 6). In
general UV models we expect, $c_{\gamma,\gamma^{\prime}}\sim{\mathcal{O}}(1)$
and $c_{\nu}\propto m_{\nu}$, the latter suppressed either by $f_{\phi}$ or
some other dimensionful scale.

If the scalar mediator is heavy enough to be integrated out, the tree level
exchange of $\phi$ generates the Rayleigh operators in (1), cf. Fig. 2,

$\frac{C_{1,ij}^{(7)}}{\Lambda^{3}}=c_{\nu}^{ij}c_{\gamma}\frac{1}{m_{\phi}^{2}f_{\phi}},\qquad\frac{C_{2,ij}^{(7)}}{\Lambda^{3}}=c_{\nu}^{ij}c_{\gamma}^{\prime}\frac{1}{m_{\phi}^{2}f_{\phi}}.$
(16)

For light mediator the effective scale $\Lambda$ is thus parametrically
smaller then the UV scale $f_{\phi}$. For instance, taking
$c_{\nu}^{ij}=10^{-3}$ and setting the other dimensionless couplings to 1, we
have $\Lambda=\\{1\,{\rm GeV},10\,{\rm MeV},100\,{\rm keV}\\}$ for
$m_{\phi}=\\{1\,{\rm MeV},1\,{\rm keV},1\,{\rm eV}\\}$ and $f_{\phi}=1$ TeV.
Whether or not the EFT description can be used for a particular process
depends on the typical energy and momentum exchange. These range from eV to
GeV for the observables we take into account in the phenomenological analysis,
as listed in Table 1.

Process |  ${\cal C}_{2}^{(7)}/\Lambda^{3}$ (GeV-3) |  EFT thr. (GeV) | Sec.
---|---|---|---
BBN | $-$ | $\sim 10^{-3}$ | 3.2
$\nu$ decay | $1.2\times 10^{11}$ | $\sim 10^{-10}$ | 3.3
$\nu$ self-interaction | $-$ | $\sim 10^{-6}$ | 3.1
HB star | $1.9\times 10^{6}$ | $\sim 10^{-5}$ | 4.1
SN1987a | $-$ | $\sim 10^{-1}$ | 4.2
Borexino | $0.68$ | $\sim 10^{-4}$ | 5.1
Xenon-nT | $1.6$ | $\sim 10^{-4}$ | 5.2
MiniBoone | $4\times 10^{-3}$ | $\sim 1$ | 5.3
$\pi^{0}\to\gamma\gamma\to\bar{\nu}\nu$ | $4.7\times 10^{3}$ | $\sim 0.1$ | 5.4
$B^{0}\to\gamma\gamma\to\bar{\nu}\nu$ | $3.7\times 10^{4}$ | $\sim 5$ | 5.4
$h\to\gamma\gamma\to\bar{\nu}\nu$ | $1.2$ | $\sim 10^{2}$ | 5.4
Table 1: Summary of bounds (2nd column) on the neutrino polarizability Wilson
coefficient ${\cal C}_{2}^{(7)}/\Lambda^{3}$, assuming flavor universality,
Eq. (3). The bounds from processes list in the 1st column were obtained under
the assumption that the EFT framework (1) applies, i.e., that the mediators
are heavier than the process dependent typical energy and momentum exchanges
(3rd column). Further details can be found in sections listed in the fourth
column. For the $h\to\gamma\gamma\to\bar{\nu}\nu$ decay (last row) the bound
is on ${\cal C}_{1}^{(7)}/\Lambda^{3}$.

In the remainder of the paper we discuss different probes of neutrino
polarizability, summarized in Fig. 3 and in Tables 1, 2, 3. We derive bounds
both assuming an EFT, Eq. (1), and assuming the scalar mediator model, Eq.
(15), for four mass benchmarks, $m_{\phi}=1$ eV, 1 keV, 1 MeV, 1 GeV. For
easier comparison with the literature, we take in the numerical analysis the
couplings to neutrinos to be flavor universal,

$c_{\nu}^{ij}=ic_{\nu}\delta_{ij},$ (17)

set $c_{\gamma}=0$, and keep $c_{\gamma}^{\prime}\neq 0$. For such purely
pseudoscalar $\phi$, we adopt the commonly used notation

$g_{\phi\gamma}\equiv\frac{\alpha}{2\pi}\frac{c_{\gamma}^{\prime}}{f_{\phi}}.$
(18)

Most of the experimental probes we consider are not sensitive to photon
polarization and do not distinguish $c_{\gamma}$ from $c_{\gamma}^{\prime}$.
In fact, most of the phenomenology is governed by the partial decay widths for
$\phi$ to photons or neutrinos,

$\displaystyle\Gamma(\phi\to\gamma\gamma)$
$\displaystyle=\left(\frac{\alpha}{8\pi}\right)^{2}\frac{m_{\phi}^{3}}{4\pi
f_{\phi}^{2}}\biggr{[}\frac{4}{9}(c_{\gamma})^{2}+(c_{\gamma}^{\prime})^{2}\biggr{]}\,,$
(19) $\displaystyle\Gamma(\phi\to\nu\nu)$
$\displaystyle=\frac{m_{\phi}}{32\pi}\sum_{ij}\big{|}c_{\nu}^{ij}\big{|}^{2}\,.$
(20)

Numerically, the ratio of $\phi\to\gamma\gamma$ and $\phi\to\nu\nu$ branching
fractions is given by

$\frac{{\cal B}(\phi\to\gamma\gamma)}{{\cal B}(\phi\to\nu\nu)}=6.7\times
10^{-19}\times\biggr{(}\frac{m_{\phi}}{{\rm MeV}}\biggr{)}^{2}\left(\frac{{\rm
TeV}}{f_{\phi}}\right)^{2}\left(\frac{(2c_{\gamma}/3)^{2}+c_{\gamma}^{{}^{\prime}2}}{\big{|}c_{\nu}^{ij}\big{|}^{2}}\right),$
(21)

For comparable values of $c_{\gamma,\gamma^{\prime}}$ and $c_{\nu}$, with
$f_{\phi}$ in the TeV regime therefore $\phi$ predominantly decays to
neutrinos. This is, for instance, a typical situation for the enhanced
neutrino polarizability model in Sec. 6.4, in which $c_{\nu}\sim
m_{\nu}/f_{\phi}^{\prime}$ with $f_{\phi}^{\prime}\ll f_{\phi}$. For the
modified majoron model in Sec. 6.3, on the other hand, $c_{\nu}\sim
m_{\nu}/f_{\phi}$ so that $c_{\nu}\sim 10^{-13}$ for $f_{\phi}\sim 1\,$TeV,
and thus the decays to photons dominate. In the phenomenological analysis in
Sections 3 to 5.4 we treat $c_{\gamma}^{\prime}$ and $c_{\nu}$ as free
parameters (assuming flavor universal neutrino couplings), and set
$c_{\gamma}=0$ (except for bounds from Higgs decays, see below). For the
parameters used in the numerical analysis, Eqs. (17), (18), the ratio of
branching fractions is given by

$\frac{{\cal B}(\phi\to\gamma\gamma)}{{\cal
B}(\phi\to\nu\nu)}=0.17\times\biggr{(}\frac{m_{\phi}}{{\rm
MeV}}\biggr{)}^{2}\left(\frac{g_{\phi\gamma}}{10^{-4}\,\text{GeV}^{-1}}\right)^{2}\left(\frac{10^{-7}}{c_{\nu}}\right)^{2}.$
(22)

Process |  $m_{\phi}=$ eV |  $m_{\phi}=$ keV | Sec.
---|---|---|---
BBN | $c_{\nu}\lesssim 4\times 10^{-5}$ | $c_{\nu}\lesssim 4.4\times 10^{-6}$ | 3.2
$\nu$ self-interaction | $c_{\nu}<2.8\times 10^{-7}$ | $c_{\nu}<2.8\times 10^{-4}$ | 3.1
HB star | $g_{\phi\gamma}\in[3.5\times 10^{-3},10^{-11}]$ | $g_{\phi\gamma}\in[3.3\times 10^{-3},10^{-11}]$ | 4.1
SN1987a | $g_{\phi\gamma}\in[10^{-2},5\times 10^{-6}]$ | $g_{\phi\gamma}\in[10^{-2},5\times 10^{-6}]$ | 4.2
$c_{\nu}\in[10^{-3},1]$ | $c_{\nu}\in[10^{-6},10^{-2}]$
Borexino | $c_{\nu}g_{\phi\gamma}<5.3\times 10^{-8}$ | $c_{\nu}g_{\phi\gamma}<5.3\times 10^{-8}$ | 5.1
Xenon-nT | $c_{\nu}g_{\phi\gamma}<8.5\times 10^{-9}$ | $c_{\nu}g_{\phi\gamma}<8.5\times 10^{-9}$ | 5.2
MiniBoone | $c_{\nu}g_{\phi\gamma}<4\times 10^{-6}$ | $c_{\nu}g_{\phi\gamma}<4\times 10^{-6}$ | 5.3
$M/\tau$ rare dec. | $c_{\nu}<4\times 10^{-3}$ | $c_{\nu}<4\times 10^{-3}$ | 5.4
$0\nu 2\beta$ | $c_{\nu}<8\times 10^{-6}$ | $c_{\nu}<8\times 10^{-6}$ | 5.4
Beam dump | $-$ | $g_{\phi\gamma}<10^{-2}$ | 5.4
$e^{+}e^{-}\to 3\gamma$ | $-$ | $-$ | 5.4
$\pi^{0}\to\gamma\gamma\to\bar{\nu}\nu$ | $c_{\nu}g_{\phi\gamma}<2\times 10^{-2}$ | $c_{\nu}g_{\phi\gamma}<2\times 10^{-2}$ | 5.4
$B^{0}\to\gamma\gamma\to\bar{\nu}\nu$ | $c_{\nu}g_{\phi\gamma}<180$ | $c_{\nu}g_{\phi\gamma}<180$ | 5.4
BaBar | $g_{\phi\gamma}<1.5\times 10^{-4}$ | $g_{\phi\gamma}<1.5\times 10^{-4}$ | 5.4
$h\to\gamma\gamma\to\bar{\nu}\nu$ | $c_{\nu}g_{\phi\gamma}<2.4$ | $c_{\nu}g_{\phi\gamma}<2.4$ | 5.4
Table 2: The bounds on flavor universal pseudoscalar mediator couplings to neutrinos, $c_{\nu}$, Eq. (17), and photons, $g_{\phi\gamma}$ (in GeV-1), Eq. (18), or the product of the two, for various processes (1st column), for two mass benchmarks, $m_{\phi}=1\,$eV, $1\,$keV (2nd and 3rd columns). Further details are given in sections listed in the last column. The BBN bounds also require $\Gamma_{\phi\to\gamma\gamma}\ll\Gamma_{\phi\to\nu\nu}$, while the quoted SN bounds refer to the regions where one of the coupling dominates, either $c_{\nu}\gg g_{\phi\gamma}$ or $g_{\phi\gamma}\gg c_{\nu}$ (for intermediate region see the main text). When the mass benchmark is in the EFT regime for the corresponding process, we use Eq. (16) for the matching.  Process |  $m_{\phi}=$ MeV |  $m_{\phi}=$ GeV | Sec.
---|---|---|---
BBN | $c_{\nu}\lesssim 5\times 10^{-9}$ | $-$ | 3.2
$\nu$ self-interaction | $c_{\nu}<2.8\times 10^{-1}$ | $c_{\nu}<2.8\times 10^{2}$ | 3.1
HB star | $-$ | $-$ | 4.1
SN1987a | $g_{\phi\gamma}\in[10^{-2},5\times 10^{-6}]$ | $-$ | 4.2
$c_{\nu}\in[10^{-9},10^{-5}]$ | $-$
Borexino | $c_{\nu}g_{\phi\gamma}<7.9\times 10^{-10}$ | $c_{\nu}g_{\phi\gamma}<7.9\times 10^{-4}$ | 5.1
Xenon-nT | $c_{\nu}g_{\phi\gamma}<3.5\times 10^{-10}$ | $c_{\nu}g_{\phi\gamma}<3.5\times 10^{-4}$ | 5.2
MiniBoone | $c_{\nu}g_{\phi\gamma}<4\times 10^{-6}$ | $c_{\nu}g_{\phi\gamma}<1.2\times 10^{-5}$ | 5.3
$M/\tau$ rare dec. | $c_{\nu}<4\times 10^{-3}$ | $c_{\nu}<0.3$ | 5.4
$0\nu 2\beta$ | $c_{\nu}<2\times 10^{-5}$ | $-$ | 5.4
Beam dump | $g_{\phi\gamma}<10^{-5}$ | $-$ | 5.4
$e^{+}e^{-}\to 3\gamma$ | $-$ | $g_{\phi\gamma}<10^{-2}$ | 5.4
$\pi^{0}\to\gamma\gamma\to\bar{\nu}\nu$ | $c_{\nu}g_{\phi\gamma}<2\times 10^{-2}$ | $c_{\nu}g_{\phi\gamma}<0.9$ | 5.4
$B^{0}\to\gamma\gamma\to\bar{\nu}\nu$ | $c_{\nu}g_{\phi\gamma}<180$ | $c_{\nu}g_{\phi\gamma}<180$ | 5.4
BaBar | $g_{\phi\gamma}<1.5\times 10^{-4}$ | $g_{\phi\gamma}<1.5\times 10^{-4}$ | 5.4
$h\to\gamma\gamma\to\bar{\nu}\nu$ | $c_{\nu}g_{\phi\gamma}<2.4$ | $c_{\nu}g_{\phi\gamma}<2.4$ | 5.4
Table 3: Same as Table 2, but for mass benchmarks, $m_{\phi}=1\,$MeV,
$1\,$GeV.

Figure 3:  Summary of the bounds on the pseudoscalar mediator model for
enhanced neutrino polarizability, Eq. (15), for four mass benchmarks,
$m_{\phi}=1\,$eV, $1\,$keV, $1\,$MeV, $1\,$GeV, as indicated. The green, blue
and red regions indicate exclusion from cosmological (Sec. 3), stellar cooling
(Sec. 4), and terrestrial (Sec. 5) constraints, respectively. The transition
region between neutrino dominated (above) to photon dominated (below) bounds
is shown as hatched, with the upper (lower) boundary corresponding to ${\cal
B}(\phi\to\gamma\gamma)=1\%$ (${\cal B}(\phi\to\nu\nu)=1\%)$.

##  3 Cosmological constraints

It is well known that the precision cosmological data impose some of the
strongest constraints on the light mediator models, such as ALPs and majoron
models Cadamuro:2011fd ; Millea:2015qra ; Depta:2020wmr ; Blinov:2019gcj .
These constraints come from a variety of cosmological measurements including
those from the measurements of cosmic microwave background (CMB) and
abundances of heavier nuclei. In this section, we apply such constraints to
the case where the light mediator can couple to both photons and neutrinos,
Eq. (15).

The four most relevant processes that determine the evolution history of
$\phi$ are pair annihilation of neutrinos ($\nu\nu\to\phi\phi$), neutrino
coalescence ($\nu\nu\to\phi$), Primakoff conversion ($\gamma\to\phi$) and
photon coalescence ($\gamma\gamma\to\phi$). The pair annihilation processes
are dominant at early times when the number densities are large, whereas
coalescence processes are dominant at temperature $T\sim m_{\phi}$.

Figure 4: Left: Constraints on $g_{\phi\gamma}$ from Planck bounds on $N_{\rm
eff}$ (blue), setting $c_{\nu}=0$. Right: Constraints on $c_{\nu}$ from BBN
(red), $\nu$SI (green) and $\phi\leftrightarrow\nu\nu$ (brownish green) from
Planck data, setting $g_{\phi\gamma}=0$. Dotted lines correspond to constant
values of $c_{\nu}$ (left) and $g_{\phi\gamma}$ (right) for which
$\Gamma(\phi\to\gamma\gamma)=\Gamma(\phi\to\nu\nu)$.

Figure 5: Constraints from cosmology on pseudoscalar mediator couplings to
neutrinos, $c_{\nu}$, and photons, $g_{\phi\gamma}$, for three mass
benchmarks, $m_{\phi}=1$ eV, keV, MeV, as denoted (for the considered
parameter ranges there are no cosmological bounds on $m_{\phi}=1$ GeV). The
red, blue and green regions are excluded from BBN (from $\nu-\phi$ coupling),
from Plank measurement of $N_{\rm eff}$, and by CMB bounds, respectively. The
brown bands, labeled $\nu\nu\to\nu\nu$, are excluded by neutrino self-
scattering bounds, derived in the EFT regime, $m_{\phi}\gg 100$ eV, see the
text for details. The dotted lines show iso-countours of $\phi$ lifetimes,
while the hatched transition region is defined as in Fig. 3.

###  3.1 Constraints from Planck

The decay of $\phi$ to photons and neutrinos can change the number of
relativistic degrees of freedom in the Universe, $N_{\rm eff}$. If the decay
happens before photon-neutrino decoupling, $T\sim 1~{}\text{MeV}$, the $\phi$
energy and entropy are quickly distributed between the different species
present in the thermal bath; the standard cosmology scenario is thus
preserved. However, if $\phi$ decays after decoupling, there is no energy and
entropy exchange between neutrinos and photons. As a consequence, the value of
$N_{\rm eff}$ during recombination can be different from the $\Lambda$CDM
prediction, $N_{\rm eff}^{\rm SM}=3.044$ Gariazzo:2019gyi . The effect can be
two-fold: the $\phi\to\nu\nu$ decays will increase $N_{\rm eff}$, whereas the
$\phi\to\gamma\gamma$ decays will decrease it. This happens because, as $\phi$
decays to photons, the photon temperature increases, while the effective
neutrino energy density is diluted.

The Planck measurement of CMB observables limits the allowed range of
relativistic degrees of freedom to $N_{\rm eff}=2.99^{+0.34}_{-0.33}$ at 95%
C.L. Planck:2018vyg . If $\phi$ decays only to neutrinos, one finds $N_{\rm
eff}=3.57$ using conservation of entropy Huang:2017egl , while in the opposite
regime, $\phi$ decaying only to photons, one finds $N_{\rm eff}=2.4$
Cadamuro:2011fd . However, the latter case presents additional complications.
Firstly, the presence of extra photons can also shift the time of matter-
radiation equality. Secondly, if $\phi$ decays close to or after
recombination, $z_{\rm rec}\sim 1089$, the standard assumptions of cosmology
no longer hold. Lastly, additional complications may occur when photons keep
$\phi$ in equilibrium, but $\phi$ predominantly decays to neutrinos. In these
cases, the Planck bound on $N_{\rm eff}$ is not reliable and cannot be used.
In this work we do not attempt to compute bounds in parts of the parameter
space which can lead to these effects, and only employ Planck’s $N_{\rm eff}$
bounds when $\phi$ predominantly decays to neutrinos after decoupling from the
thermal bath.

Another effect that can change $N_{\rm eff}$ is the re-thermalization of
$\phi$. A light $\phi$, $m_{\phi}\lesssim\rm{MeV}$, can thermalize again with
photons (neutrinos) due to the inverse decay process, $\gamma\gamma\to\phi$
($\nu\nu\to\phi$). If this rethermalization occurs at $T\gtrsim m_{\phi}$, the
abundance of $\phi$ will increase to its thermal equilibrium value. However,
as the Universe cools down, $\phi$ becomes non-relativistic at $T\sim
m_{\phi}/3$ and decays out of equilibrium to photons (neutrinos). As shown in
Cadamuro:2011fd , this effect is subdominant for ALPs that only couple to
photons, we expect the same for $\phi$ as well.

In Fig. 4 (left), we show the constraints on $\phi$ coupling to photons
(setting $c_{\nu}=0$) as a function of $\phi$ mass that were obtained for ALPs
using Planck data in Ref. Depta:2020wmr  (see also Ref. Millea:2015qra ). Note
that there are no constraints for $m_{\phi}<1\,\text{keV}$: for these masses
$\phi$ decays to photons close to time of recombination. For masses larger
than a few MeV, $\phi$ decays before photon-neutrino decoupling and therefore
only very small couplings are constrained. The dotted lines in Fig. 4 (left)
show constant values of $c_{\nu}$ for which ${\cal
B}(\phi\to\gamma\gamma)={\cal B}(\phi\to\nu\nu)$. The couplings to neutrinos
need to be below these values for the $N_{\rm eff}$ bounds to apply unchanged.

We discuss next the CMB bounds on $\phi$ couplings to neutrinos. The
$\nu-\phi$ coupling suppresses the neutrino anisotropic stress energy tensor,
which leads to distortions in the CMB power spectrum. In the $\Lambda$CDM
model, neutrinos are free-streaming particles and this description agrees very
well with the Planck measurements. The new coupling with $\phi$ induces
neutrino self-scattering ($\nu$SS), $\nu\nu\to\nu\nu$. Self-scattering
neutrinos behave like a fluid, rather than free-streaming radiation, and thus
leave their imprints in the CMB power spectrum. The processes that define the
latter spectrum happen at the typical temperature $T\sim 100$ eV, which
defines two different $m_{\phi}$ regimes.

For $m_{\phi}\gg 100$ eV, the flavor universal $\nu$SS mediated by $\phi$ is
effectively described by a dimension 6 scalar operator

$\mathcal{L}_{\nu{\rm
SS}}\supset\frac{1}{4}\frac{c_{\nu}^{2}}{m_{\phi}^{2}}(\bar{\nu}_{i}P_{L}\nu_{i})(\bar{\nu}_{j}P_{L}\nu_{j})\equiv
G_{\nu}(\bar{\nu}\nu)(\bar{\nu}\nu)\,,$ (23)

where we already assumed that the process is flavor-universal. Fits to
cosmological data find that $\log_{10}(G_{\nu}\times\rm{MeV}^{2})<-3.57$ at
95% C.L. or $\log_{10}(G_{\nu}\times\rm{MeV}^{2})=-1.711\pm 0.099$ at 68% C.L.
Lancaster:2017ksf . Stronger neutrino self-interactions are ruled out by
terrestrial experiments: from meson decays, $\tau$ decays, and double beta
decays Blinov:2019gcj . We show these bounds in the right panel of Figure 4.

In the opposite regime, $m_{\phi}\ll 100$ eV, $\phi$ is a relativistic degree
of freedom at the CMB formation temperature and one needs to include
$\nu-\phi$ interactions in the Boltzmann equation. Ref. Escudero:2019gvw
computed such constraints on light $\phi$ from the Planck 2018 data, which we
show as a green band in Fig. 5. The green horizontally hatched region
indicates the parts of the parameter space in which $\phi$ and $\gamma$
decouple around the time of recombination and the standard assumptions of
cosmology may not hold.

Finally, we comment on the effect of the $\nu\nu\gamma\gamma$ interaction due
to Rayleigh operators, i.e., for $m_{\phi}$ induced interactions, but in the
EFT regime, cf. Eq. (16). The Rayleigh operators can keep neutrinos and
photons in thermal equilibrium and modify the CMB power spectrum measured by
Planck. As already mentioned, the relevant scale for CMB is $T\sim 100$ eV.
Thus we need to estimate if thermal equilibrium can be achieved at $T\lesssim
100$ eV. The scattering rate for the process can be approximated as

$\Gamma_{\nu\nu\gamma\gamma}\sim\left(\frac{\alpha}{2\pi}\frac{{\cal
C}_{2}^{(7)}}{\Lambda^{3}}\right)^{2}T^{7}.$ (24)

At the decoupling temperature $T_{\rm fo}$, the scattering rate should satisfy
the condition

$\frac{\Gamma_{\nu\nu\gamma\gamma}}{H}\lesssim
1\quad\implies\left(\frac{\alpha}{2\pi}\frac{{\cal
C}_{2}^{(7)}}{\Lambda^{3}}\right)^{2}T_{\rm fo}^{5}M_{\rm Pl}\lesssim 1\,,$
(25)

where $H$ is the Hubble rate at decoupling. Setting $T_{\rm fo}=100\rm{~{}eV}$
leads to the following bound on the neutrino polarizability operator

$\Lambda^{3}/{\cal C}_{2}^{(7)}\gtrsim(0.2\,\text{MeV})^{3},$ (26)

or in terms of the $\phi$ mediator model,

$c_{\nu}^{2}\,(g_{\phi\gamma}\times\text{GeV})^{2}\left(\frac{\rm{keV}}{m_{\phi}}\right)^{4}\lesssim
10^{-7}\,.$ (27)

The bound in (26) is model independent as long as the mediators generating the
Rayleigh operator are heavier than about 100 eV. However, this EFT bound is
also relatively weak. For instance, if the neutrino polarizability is induced
by the pseudoscalar mediator, Eq. (15), the other cosmological bounds for
$m_{\phi}>100$ eV are more stringent than Eq. (27), cf. Figs. 3, 4. The
neutrino-photon interaction in this case therefore freezes-out much before
recombination and does not lead to any new constraint.

###  3.2 Constraints from BBN

In order to estimate the impact of $\phi$ on the abundances of primordial
elements produced during the BBN, we consider three parameter regimes. Below,
we compute the BBN bounds for the case when decays to neutrinos dominate, by
considering the extreme case of no coupling to photons, $g_{\phi\gamma}=0$. In
the opposite regime, $c_{\nu}=0$, the bounds from Planck data dominate and the
effects on BBN can be neglected Depta:2020wmr . The constraints in the
intermediate regime, shown as the diagonally hatched regions in Figs. 3 and 5,
are more involved to estimate and go beyond the scope of this paper. The upper
(lower) boundary of the intermediate regime region are defined by requiring
that the neutrino (photon) channel accounts for 99% of the total width.

The BBN bounds on $\phi-\nu$ coupling, in the limit of no couplings to
photons, is shown in red color in Fig. 5, indicated as “BBN$\nu$”, and is
obtained as follows. The abundance of $\phi$ during BBN can increase $N_{\rm
eff}$, which in turn modifies the expansion rate of the Universe and thus the
abundance of heavy elements. The two processes that can keep $\phi$ in thermal
equilibrium are the neutrino pair annihilation ($\nu\nu\to\phi\phi$) and the
neutrino coalescence ($\nu\nu\to\phi$). For very light $\phi$, $m_{\phi}\ll
1~{}{\rm MeV}$, the pair annihilation process dominates and keeps $\phi$ in
thermal equilibrium during BBN, whereas for $m_{\phi}\lesssim 1~{}{\rm MeV}$,
the inverse decay dominates. The BBN bounds on $\phi-\nu$ coupling can be
written as Escudero:2019gvw

$c_{\nu}<\frac{1}{\left(5\times
10^{-9}\frac{\mathrm{MeV}}{m_{\phi}}\right)^{-1}+\left(4\times
10^{-5}\right)^{-1}}\,,$ (28)

where the two parts in the denominator come from the two aforementioned
processes. The above bound is shown as excluded red regions in Fig. 5.

Both contributions can be obtained by requiring that $\phi$ is not in thermal
equilibrium at the photon-neutrino decoupling temperature, $T_{D}\sim
1~{}\text{MeV}$, as this would otherwise result in $\Delta N_{\rm eff}\sim
0.5$ at the time of BBN. We then require

$\displaystyle\frac{\Gamma_{\nu\nu\to\phi\phi}(T_{D})}{H(T_{D})}\lesssim
1\implies c_{\nu}^{4}\lesssim\frac{T_{D}}{M_{pl}}\sim
10^{-5},\rm{~{}~{}~{}for~{}m_{\phi}\ll 1{\rm MeV}}\,,$ (29)

while for $m_{\phi}\lesssim 1~{}{\rm MeV}$,

$\frac{\Gamma_{\nu\nu\to\phi}(T_{D})}{H(T_{D})}\sim\frac{c_{\nu}^{2}m_{\phi}^{2}M_{pl}}{T_{D}^{3}}\lesssim
1\implies c_{\nu}m_{\phi}\lesssim 10^{-10}~{}{\rm MeV}.$ (30)

In the opposite limit ($c_{\nu}=0$), where $\phi$ behaves like an ALP coupling
to photons, Ref. Depta:2020wmr  found that the BBN constraints on ALPs are
weaker than those from the Planck data. Therefore, in the parameter regime
where decays to photons dominate, we indicate in Fig. 5 only the Planck
constraints.

###  3.3 Neutrino decay

The Rayleigh operators can induce the decay of a neutrino into a lighter mass
eigenstate, along with two photons, $\nu_{i}\to\nu_{j}\gamma\gamma$, where
$m_{\nu_{i}}>m_{\nu_{j}}$. The sum of the neutrinos masses is bounded from CMB
Planck data to be $\sum_{i}{m_{\nu_{i}}}<\rm{0.12~{}eV}$ (95% CL, Planck
TT,TE,EE+lowE +lensing+BAO Planck:2018vyg ). Therefore, the photon energy
spectrum will follow the typical 3-body decay distribution, with a maximum
energy $\lesssim 0.1~{}\rm{eV}$. Depending on the time of their injection,
these photons may leave their imprints in the CMB power spectrum measured by
Planck Planck:2018vyg  or the CMB blackbody spectrum measured by COBE/FIRAS
Fixsen:2009ug ; Fixsen_2002 .

We assume for simplicity that the final state neutrino is massless,
$m_{\nu_{j}}=0$, and that the mediator is heavy, $m_{\phi}\gg m_{\nu_{i}}$,
which is true for all the mass benchmarks considered. For our purposes it
suffices to estimate the decay width using naive dimensional analysis,

$\Gamma_{\nu_{i}\to\nu_{j}\gamma\gamma}\sim\frac{1}{8\pi}\left(\frac{g_{\phi\gamma}c_{\nu}^{ba}}{4\pi}\right)^{2}\frac{m_{\nu}^{7}}{m_{\phi}^{4}}\,,$
(31)

which corresponds to a lifetime

$\tau_{\nu_{i}}\simeq 8\times
10^{17}\left(\frac{1}{c_{\nu}\left(g_{\phi\gamma}\cdot{\rm
GeV}\right)}\right)^{2}\left(\frac{m_{\phi}}{{\rm
keV}}\right)^{4}\left(\frac{0.1~{}{\rm eV}}{m_{\nu}}\right)^{7}~{}{\rm
years}\,.$ (32)

If the neutrino lifetime becomes comparable or smaller than the age of the
Universe, $t_{0}\simeq 1.4\times 10^{10}$ years, the emitted photon would
affect the observed CMB spectrum. We then require

$c_{\nu}^{ij}g_{\phi\gamma}\leq 8\times 10^{3}\times\left(\frac{m_{\phi}}{{\rm
keV}}\right)^{2}\left(\frac{0.1~{}{\rm eV}}{m_{\nu}}\right)^{7/2}~{}{\rm
GeV}^{-1}\,.$ (33)

Numerically, taking $m_{\nu}=0.1$ eV, and $m_{\phi}=1$ eV, this gives
$c_{\nu}^{ij}g_{\phi\gamma}\lesssim 8\times 10^{-3}~{}{\rm GeV}^{-1}$, while
for $m_{\phi}=1$ GeV, the constraint is $c_{\nu}^{ij}g_{\phi\gamma}\lesssim
8\times 10^{15}~{}{\rm GeV}^{-1}$.

##  4 Stellar cooling constraints

If stellar dynamics is able to produce light new physics states that
efficiently escape from its core, it can lead to excessively large stellar
cooling rates. Requiring that the additional cooling does not exceed the
standard model one, typically leads to very stringent bounds on light new
physics sectors. In this section we evaluate the stellar cooling bounds for
the light scalar $\phi$ that couples to photons and neutrinos, Eq. (15), with
the results summarized in Fig. 6. The stellar cooling rates are controlled by
the $\phi$ production rates, as well as its decay length and/or mean free
path.

In the analysis we distinguish two cases. In the first category are the
cooling rates for Horizontal Branch stars (HB), Red Giants (RG) and White
Dwarves (WD), for which the core temperature is low, $T\sim$few keV. In the
second category are the Supernova (SN) cooling constraints, for which the core
of the proto-neutron star is much denser and hotter, $T_{\rm SN}\sim 30$ MeV.
The $\phi$ production mechanisms are the Primakoff conversion,
$\gamma\to\phi$, the photon coalescence, $\gamma\gamma\to\phi$, and for SN
also the neutrino coalescence, $\nu\nu\to\phi$. The rates for these processes
are given in Appendix B.

Figure 6: Stellar cooling bounds on pseudoscalar couplings to neutrinos,
$c_{\nu}$ and photons, $g_{\phi\gamma}$, for three mass benchmarks,
$m_{\phi}=1$ eV, keV, MeV, while there are no bounds for $m_{\phi}=1$ GeV. The
blue and purple regions show the bounds from HB stars and SuperNovae,
respectively.

The inverse reactions constitute the main decay channels: the decays to two
neutrinos, $\phi\to\nu\nu$, Eq. (19), the decays to two photons,
$\phi\to\gamma\gamma$, Eq. (19), and the inverse Primakoff process,
$\phi\to\gamma$, Eq. (106). The $\phi$ decays to photons, as well as the
Primakoff process, reduce the stellar cooling rates, since they reduce the
number of $\phi$ particles that escape the interior of the star. The same is
true for $\phi\to\nu\nu$ decays, whenever these occur inside the SN, since the
neutrinos get trapped inside the proto-neutron star. The HB, RG and WD cooling
rates, on the other hand, are not affected by the $\phi\to\nu\nu$ decays,
since the neutrinos escape from these types of stars in the same way as the
undecayed $\phi$ would have.

For processes that involve photons we need to take into account finite
temperature effects. To first approximation, the effect of a photon
propagating in plasma instead of in a free space can be taken into account by
replacing the dispersion relation for a massless photon with the one for a
collective excitation – plasmon, i.e., a massive photon with the mass equal to
the plasma frequency

$\omega_{P}^{2}=\frac{4\pi\alpha n_{e}}{E_{F}}\,,$ (34)

where $n_{e}$ is the electron number density in the core and $E_{F}$ the Fermi
energy of the electron gas. For rough numerical estimates we can use the non-
relativistic Fermi gas expression for the plasma frequency $\omega_{P}\sim
28.7\,{\rm eV}\big{[}{\rho}/({{\rm g/cm^{3}}})\big{]}^{1/2}$.

The $\phi$ emissivity, i.e., the energy emitted in the production of the final
state $\phi$ per unit volume and time, is given by

$Q_{\phi}=g_{I}\int\frac{{d}^{3}p}{(2\pi)^{3}}\Gamma_{I\to\phi}E_{\phi}f(E_{\phi})=\int_{m_{\phi}}^{\infty}dE_{\phi}E_{\phi}\frac{{d}^{2}N_{\phi}}{{d}E_{\phi}{d}t}.$
(35)

Here, $p$ is the $\phi$ momentum, $\Gamma_{I\to\phi}$ the $\phi$ production
rate from initial state $I$, with $g_{I}$ number of degrees of freedom (the
number of polarization states), and $f(E_{\phi})$ the thermal distribution of
particles $\phi$ in the stellar core. The number of emitted $\phi$ particles
per unit of energy and time, ${{d}^{2}N_{\phi}}/{{d}E_{\phi}{d}t}$ is a sum of
different $\phi$ production mechanisms: Primakoff conversion, photon
coalescence, and neutrino coalescence (for SN), see Appendix B.

###  4.1 Horizontal branch stars

For cold cores, i.e., for WD and HB stars, we can safely approximate the
stellar core with a homogeneous sphere of radius $R_{c}$. The total luminosity
is then given by

$L_{\phi}=V_{c}\left[g_{I}\int\frac{{d}^{3}p}{(2\pi)^{3}}\Gamma_{I\to\phi}E_{\phi}f(E_{\phi})S(R_{c})\right]\,,$
(36)

where $V_{c}=4\pi R_{c}^{3}/3$ is the volume of the stellar core, and we used
Eq. (35) for the emissivity. The factor $S(R_{c})$ takes into account the
suppression of the $\phi$ luminosity due to decays of $\phi$ into photons.
Namely, if $\phi$ decays via $\phi\to\gamma\gamma$ while still inside the
stellar core, then its energy is reabsorbed by the plasma. If $\phi$ decays
instead via the $\phi\to\nu\nu$ channel, its energy is still carried away by
the two neutrinos, thus contributing to the exotic cooling of the star. The
suppression factor $S(R_{c})$ is then given by

$S(R_{c})=1-\exp\left(-\Gamma_{\nu}{R_{c}}\right)+\exp\left(-\Gamma_{\rm
tot}{R_{c}}\right)\,,$ (37)

where

$\Gamma_{\nu}=\frac{1}{\beta_{\phi}}\frac{\Gamma_{\phi\to\nu\nu}}{\gamma_{\phi}},\qquad\Gamma_{\rm
tot}=\frac{1}{\beta_{\phi}}\Big{(}\frac{\Gamma_{\phi\to\gamma\gamma}+\Gamma_{\phi\to\nu\nu}}{\gamma_{\phi}}+{\Gamma_{\phi\to\gamma}}\Big{)}\,,$
(38)

with $\gamma_{\phi}=E_{\phi}/m_{\phi}$ the Lorentz factor, and $\beta_{\phi}$
the velocity of $\phi$. In writing (37) we approximated the suppression factor
for $\phi$ originating from any given point inside the star to be the same as
when traversing distance $R_{c}$, i.e., the typical linear dimension. Beside
total luminosity, $L_{\phi}$, the quantity often considered in the literature
is the luminosity per unit mass, $\epsilon_{\phi}=L_{\phi}/M$, where $M$ is
the mass of the stellar core.

For cold cores the main energy-loss process, within the standard model, is the
neutrino emission via plasmon decay into two photons Haft:1993jt . The
measurements of RG and HB cooling rates give a typical value for the observed
luminosity per unit mass of $|\epsilon_{\nu}|\simeq 100~{}{\rm erg}~{}{\rm
g}^{-1}~{}{\rm s}^{-1}$, with a $\sim 10\%$ uncertainty. Requiring that the
extra cooling due to emission of $\phi$ is smaller than the experimental error
on the measurement gives PhysRevD.37.549

$|\epsilon_{\phi}|\leq 10~{}{\rm erg}~{}{\rm g}^{-1}~{}{\rm s}^{-1}\,.$ (39)

The excluded regions in the $c_{\nu}-g_{\phi\gamma}$ plane are shown in Fig. 6
as blue bands for two mass benchmarks, $m_{\phi}=1\,$eV, $1\,$keV, while for
$m_{\phi}=1\,$MeV, $1\,$GeV benchmarks $\phi$ is too heavy to be produced in a
cold stellar core. We only show the bounds obtained from HB stars, due to
their higher density and thus higher luminosity than RG stars. The benchmark
values are $T_{\rm HB}\sim 10^{8}{\rm K}\sim 8.6$ keV, $\rho_{\rm HB}\sim
10^{4}{\rm g/cm^{3}}$ and $R_{\rm HB}\sim 0.03R_{\odot}$ for temperature,
density and helium-burning core radius respectively, where
$R_{\odot}=6.96\times 10^{5}$ km is the radius of the Sun. The associated
plasma frequency is $\omega_{\rm HB}\sim 3$ keV.

The lower boundary of excluded $g_{\phi\gamma}$ values in Fig. 6,
$g_{\phi\gamma}\sim 10^{-11}~{}{\rm GeV}^{-1}$, indicates the value at which
$\phi$ starts to be produced efficiently in the star, thus exceeding the bound
on exotic cooling, Eq. (39). The upper boundary of $g_{\phi\gamma}$ exclusion,
on the other hand, indicates the onset of parameter region for which $\phi$ is
trapped: the coupling to photon is strong enough that $\phi$ always decays
inside the stellar core. In the parameter space scanned in this work, the
trapping regime is never reached; the decays to neutrinos compensate the
exponential suppression from photon decay, Eq. (37), and thus the production
of $\phi$ still contributes to stellar cooling. For very small values of
couplings to neutrinos, below the values shown in the figure, there are
however regions where neutrino coupling is not able to overcome the trapping.

In the EFT regime, where $\phi$ is too heavy to be produced on-shell, the
cooling mechanism is due to a Primakoff conversion
$\gamma^{*}\gamma_{L}\to\nu\nu$, where $\gamma^{*}$ is the plasmon and the
transition is induced by the longitudinal plasmon of the external electric
field of the plasma, $\gamma_{L}$. In the limit where the latter is taken as
static, the transition can be approximated as the two body decay
Raffelt:1996wa , and the rate scales roughly as
$\Gamma_{\gamma^{*}\gamma_{L}\to\nu\nu}\sim\omega_{HB}^{7}/\Lambda^{6}$.
Demanding that the cooling is smaller than the error on the measured rates
gives the lower bound on the effective suppression scale of the Rayleigh
operator

$\Lambda^{3}/{\cal C}_{2}^{(7)}\gtrsim\big{(}8\,\text{MeV}\big{)}^{3}.$ (40)

###  4.2 Supernova cooling

There are several important differences between cooling rates deduced from HB
stars and the SN. Firstly, at proto-neutron star densities and temperatures
the neutrinos produced inside the core are efficiently trapped, leading to a
thermal population of neutrinos with a chemical potential $\mu_{\nu}\sim
200\,$MeV. The decays of a propagating $\phi$ into neutrinos, therefore no
longer lead to enhanced cooling rates. Secondly, the total luminosity is quite
sensitive to the exact radial profile of the SN core after the start of the
explosion. The SM temperature and density profiles as a function of the
distance from the center, $r$, can be obtained via numerical simulations, and
depend both on the initial conditions of the progenitor star and the explosion
mechanism. In the numerics we use the profile from Ref. Fischer:2018kdt  at
the benchmark time $t=1$ s after the start of the explosion, leading to
constraints that are in general agreement with the results of Ref.
Fiorillo:2022cdq .

The production of $\phi$ is $r$ dependent, since the production rates depend
on $T(r),\omega_{P}(r),\ldots$ After the production, $\phi$ propagates inside
the core and contributes to the SN cooling, if it escapes the neutrino-sphere
of radius $R_{\nu}\sim 23\,$km, i.e., the region where neutrino production
rate is higher than the absorption. In the opposite case, the energy taken by
the $\phi$ is re-deposited into neutrinos. The probability that $\phi$ reaches
a distance $R_{\nu}$ is controlled by the optical depth Lucente_2020 ,

$\tau_{\phi}\left(r,E_{\phi},R_{\nu}\right)=\int_{r}^{R_{\nu}}\Gamma_{\rm
tot}(\tilde{r})\frac{{\rm d}\tilde{r}}{d_{\phi}(\tilde{r})}\,,$ (41)

with $\Gamma_{\rm tot}$ given in (38) now dependent on radial distance. The
total luminosity is then given by

$L_{\phi}=4\pi\int_{0}^{R_{\nu}}{\rm d}r~{}r^{2}\left(g_{I}\int\frac{{\rm
d}^{3}p}{(2\pi)^{3}}\Gamma_{I\to\phi}E_{\phi}f(E_{\phi})~{}\exp\left[-\tau_{\phi}\left(r,E_{\phi},R_{\nu}\right)\right]\right)\,.$
(42)

The bounds on neutrino and photon couplings can be obtained by imposing that
$L_{\phi}$ does not exceed the measured neutrino luminosity,

$L_{\phi}\lesssim L_{SN}^{\nu}\sim 3\times 10^{52}~{}{\rm erg}~{}{\rm
s}^{-1}\,,$ (43)

i.e., which is the usual rule of thumb prescription that translates the
absence of large cooling effects in the observed neutrino flux from SN1987a to
a bound on the production of light particles.

The resulting excluded region is shown as a purple band in Fig. 6. In contrast
to HB stars, the SN core is hot enough to produce new particles with masses up
to $m_{\phi}\sim 100$ MeV. In the parameter range we are interested in, the
bound is affected both by $\phi$ couplings to photon and neutrinos. The
excluded region in Fig. 6 is a horizontal band when $c_{\nu}$ dominates, and
is a vertical band when $g_{\phi\gamma}$ is more important. The change from
one regime to the other can be roughly understood through Eq. (21).

For $\phi$ heavier than $m_{\phi}\gtrsim 100$ MeV the $\phi$ particle is not
produced on shell in the proto-neutron star, and thus does not contribute to
cooling. The situation is different from HB stars, where the Rayleigh operator
created by integrating out a heavy $\phi$ can still enhance the cooling rates
via the production of neutrinos through the Primakoff transition,
$\gamma^{*}\gamma_{L}\to\nu\nu$. The neutrinos then escape and lead to
enhanced stellar cooling rates. For SN, the neutrinos are instead trapped
inside the dense SN core. Increased coupling between photons and neutrinos,
due to a new off-shell degree of freedom, therefore has no visible effect.

Finally, we comment on the possibility that the neutrinos produced in the
SN1987a core would interact with CMB photons and modify the observed CMB
spectrum. The scattering length is given by
$\lambda_{\nu\gamma}\sim(n_{\gamma}\sigma_{\nu\gamma})^{-1}$, where
$\sigma_{\nu\gamma}$ is the cross section for $\nu\gamma\to\nu\gamma$
scattering, mediated by the $s-$ and $u-$channel tree level $\phi$ exchange.
Taking the neutrinos to have fixed energy $E_{\nu}\sim 30\,$MeV, and the CMB
photon the typical energy $E_{\gamma}\sim T_{\rm CMB}\sim 2\times
10^{-13}\,$GeV, with the number density $n_{\gamma}\sim 2.2\times 10^{8}\,{\rm
m}^{-3}$, we find

$\lambda_{\nu\gamma}\sim 1.2\times 10^{53}\left(\frac{m_{\phi}}{{\rm
GeV}}\right)^{4}\Big{(}\frac{1}{c_{\nu}\left(g_{\phi\gamma}\times{\rm
GeV}\right)}\Big{)}^{2}~{}{\rm m}\,.$ (44)

The condition that the scattering length $\lambda_{\nu\gamma}$ is less than
the distance of SN1987a from Earth, $d_{\rm SN}\sim 2\times 10^{21}\,$m, is
achieved for

$c_{\nu}\left(g_{\phi\gamma}\times{\rm GeV}\right)\gtrsim 2.5\times
10^{-3}\left(\frac{m_{\phi}}{{\rm eV}}\right)^{2}\,.$ (45)

This can exclude part of the parameter space we are intersted in for the mass
benchmark $m_{\phi}=1$ eV. The bound is shown as a red region in Fig. 6. For
$m_{\phi}>{\mathcal{O}}(\text{eV})$ only relatively large values of
$c_{\nu},g_{\phi\gamma}$ are covered, which are already well excluded by
cosmology and laboratory searches, and thus do not appear in the plots.

##  5 Bounds from terrestrial experiments

Next we discuss the bounds on neutrino polarizability from terrestrial
detectors. In neutrino and dark matter experiments, the incoming neutrinos can
scatter on electrons or nuclei in the detector. The Rayleigh operators induce
at 1-loop the $\nu X\to\nu X$ scattering, where $X=e,N$ is either an electron
or a nucleon, and at tree level $\nu X\to\nu X\gamma$, i.e., neutrino
scattering with an emission of an extra photon. The scatterings on nucleons,
such as the coherent neutrino nucleus scattering COHERENT:2017ipa ;
COHERENT:2020iec ; CONUS:2020skt ; CONUS:2021dwh ; nuGeN:2022bmg , leads to
less stringent bounds than scattering on electrons Altmannshofer:2018xyo . The
resulting bounds from Borexino, Xenon-nT and MiniBoone are given in Sections
5.1, 5.2, and 5.3, respectively. In Section 5.4 we discuss collider
constraints. Summary of the terrestrial constraints on the pseudoscalar
coupling to neutrinos and photons is given in Fig. 7, for the four mass
benchmarks, $m_{\phi}=1$ eV, 1 keV, 1 MeV, 1 GeV.

Figure 7:  Constraints on neutrino and photon couplings of pseudoscalar
mediator from terrestrial experiments: from Xenon-nT (purple), Borexino (red),
MiniBoone (blue), BaBar monophoton (orange), Belle II $e^{+}e^{-}\to 3\gamma$
(light red) and beam dump (light red) searches, for the four mass benchmarks
$m_{\phi}=1\,$eV, $1\,$keV, $1\,$MeV, $1\,$GeV. The dotted lines are iso-
contours of $\phi$ lifetimes.

###  5.1 Bounds from Borexino

Borexino measured scattering of solar neutrinos on electrons Borexino:2017fbd
, where the detector response does not distinguish between $\nu e\to\nu e$ and
$\nu e\to\nu e\gamma$ scattering events, and thus the two need to be added
when comparing with the measured rates. We calculate the tree-level $\nu
e\to\nu e\gamma$ scattering numerically by implementing the Rayleigh operator
in MadGraph Alwall:2011uj . For the much smaller 1-loop induced $\nu e\to\nu
e$ scattering we use the NDA estimate,
$\sigma_{i}\simeq\big{(}\big{|}\hat{C}_{1,i}^{(7)}\big{|}^{2}+3\big{|}\hat{C}_{2,i}^{(7)}\big{|}^{2}/2\big{)}\alpha^{4}s^{2}/(48\pi^{4})$,
where $s=m_{e}^{2}+2m_{e}E_{\nu}$ is the center of mass energy of the
scattering process, and assume flavor diagonal couplings, cf. Eq. (3).

The event rate per day per 100 tons of detector is given by

$\begin{split}R_{i}(\hat{C}_{1}^{(7)},\hat{C}_{2}^{(7)})=TN_{e}\int_{E_{\nu,{\rm
min}}}^{E_{\nu,{\rm
max}}}dE_{\nu}\,\phi_{i}(E_{\nu})&\Big{[}P_{i}^{e}\sigma_{\nu_{e}e}(E_{\nu},\hat{C}_{1,e}^{(7)},\hat{C}_{2,e}^{(7)})\\\
&+P_{i}^{\mu}\sigma_{\nu_{\mu}e}(E_{\nu},\hat{C}_{1,\mu}^{(7)},\hat{C}_{2,\mu}^{(7)})\Big{]}\,,\end{split}$
(46)

with $T=1\,{\rm day}=8.64\cdot 10^{4}\,{\rm s}$ the exposure time and
$N_{e}=3.307\cdot 10^{31}$ the number of target electrons in 100 tons of
detector mass, while $\ell=e,\mu$ are the incoming neutrino flavors. The label
$i$ in (46) denotes the main components of the solar neutrino flux on Earth,
$\phi_{i}(E_{\nu})$, Ref. Vitagliano_2020 ; due to proton-proton fusion
($i=pp$), Berillium 7 electron capture ($i={}^{7}{\rm Be}$), and proton
electron capture ($i=pep$). The $\nu_{e}$ from $pp$ have a continuous energy
spectrum with the maximal energy $E_{\nu,{\rm max}}=0.423$ MeV, while 7Be and
$pep$ neutrinos are monochromatic, with energies $E_{{}^{7}{\rm Be}}=0.863$
MeV and $E_{pep}=1.445$ MeV, respectively. The minimal incoming neutrino
energy that can still produce the threshold $\sim 50$ keV recoil in Borexino
is $E_{\nu,{\rm min}}=0.139$ MeV. The $\nu_{e}$ neutrinos produced in the Sun
undergo flavor oscillations while propagating to Earth. The $\nu_{e}$ survival
probabilities are $P_{i}^{e}=\\{0.554\,,0.536\,,0.529\,\\}$ for
$i=\\{pp,^{7}{\rm Be},pep\\}$, once matter effects are taken into account
Khan_2020 , while $P_{i}^{\mu}=(1-P_{i}^{e})/2$, assuming maximal
$\theta_{23}$ for simplicity.

The bounds on $\hat{C}_{j,\ell}^{(7)}$ are obtained using the following chi-
squared function

$\chi^{2}(\vec{\alpha},\hat{C}_{j,\ell})=\sum_{i}\frac{\Big{[}R_{{\rm
meas},i}-R_{i}(\hat{C}_{1,\ell}^{(7)},\hat{C}_{2,\ell}^{(7)})(1+\alpha_{i})\Big{]}^{2}}{\sigma_{i}^{2}}+\left(\frac{\alpha_{i}}{\sigma_{\alpha_{i}}}\right)^{2},$
(47)

where the sum is over the three types of solar neutrino fluxes. The measured
event rates in Borexino phase-I and their statistical uncertainties are
Khan_2020  $R_{{\rm meas},i}\pm\sigma_{i}=\\{134\pm 10,48.3\pm 1.1,2.43\pm
0.36\\}$, $i=\\{pp,^{7}{\rm Be},pep\\}$, to be compared with the SM rates
$R_{i}(0)=\\{131.4,48.1,2.8\\}$, where the theoretical errors on the
predictions are accounted for by marginalizing over the parameters
$\alpha_{i}$, with $\sigma_{\alpha_{i}}=\\{1.1\%,5.8\%,1.5\%\\}$. The
resulting $1\sigma$ allowed ranges on Rayleigh operators are

$\left|\frac{2}{3}\hat{C}_{1,\ell}^{(7)}+\hat{C}_{2,\ell}^{(7)}\right|\leq\left\\{0.71,~{}2.55,~{}0.68\right\\}~{}\text{GeV}^{-3}\,,\quad\text{for}\quad\ell=\\{e,\mu,\text{univ.}\\},$
(48)

assuming photon couplings to either only $\nu_{e}$ or $\nu_{\mu}$ or both
(with universal couplings). In Table 1 we list the result for the universal
couplings, assuming only the CP-odd Rayleigh operator is nonzero.

For light mediators, with mass much lower than a typical momentum exchange in
Borexino, $m_{\phi}^{2}\ll|q^{2}|$, where $|q|\sim 100$ keV, the EFT framework
no longer applies, and we include the full $\phi$ propagator in the MadGraph
calculation of the cross sections. Comparison with the measurements then gives
for light $\phi$

$c_{\nu}g_{\phi\gamma}\leq\left\\{0.56,~{}1.91,~{}0.53\right\\}\times
10^{-7}~{}\text{GeV}^{-1}\,,\quad\text{for}\quad\ell=\\{e,\mu,\text{univ.}\\}.$
(49)

In Table 2 and 3 we quote only the bound for the flavor universal case. The
corresponding bounds for the four benchmark masses are shown as excluded red
regions in Fig. 7.

Figure 8: Neutrino flux at the MiniBoone detector for each of the neutrino
flavor components, as denoted in the legend.

###  5.2 Bounds from dark matter detectors

Dark matter direct detection experiments are sensitive to enhanced $\nu
e\to\nu e$ scattering rates from solar neutrinos, similar to the Borexino
bounds discussed in the previous section, but with a lower recoil energy
threshold of $\sim 2\,$keV, which translates to a lower required minimal
energy of the incoming neutrinos, $E_{\nu,{\rm min}}\simeq 16$ keV, in Eq.
(46). The strongest constraints come from the recent measurement of electron
recoil events by XENONnT XENON:2022mpc  with exposure of 1.16 tonne-years,
with no excess observed over the background rate $R_{\rm bg}=(16.1\pm
1.3,\text{stat})$ events/(t$\cdot$y$\cdot$keV) in the $(1,30)\,$keV recoil
energy search window (see also CDEX:2019exx ; PandaX-II:2020udv ;
EDELWEISS:2018tde ; LUX:2017glr ). The resulting constraints on Rayleigh
operators are

$\left|\frac{2}{3}\hat{C}_{1,\ell}^{(7)}+\hat{C}_{2,\ell}^{(7)}\right|\leq\left\\{1.6,~{}9.6,~{}1.6\right\\}~{}\text{GeV}^{-3}\,,$
(50)

for $m_{\phi}\gg q^{2}$, and

$c_{\nu}g_{\phi\gamma}\leq\left\\{2.4,~{}9.9,~{}2.6\right\\}\times
10^{-8}~{}\text{GeV}^{-1}\,,$ (51)

for $m_{\phi}\ll q^{2}$. In Tables 1, 2, and 3 we list only the constraint for
the pseudoscalar case with flavor universal couplings, with the corresponding
excluded parameter regions for the four mass benchmarks shown as blue regions
in Fig. 7.

In deriving the above bounds we included only $\nu e\to\nu e$ scattering as
the signal. Dark matter detectors have in principle the possibility to probe
also the subleading $\nu e\to\nu e\gamma$ process, by searching for an extra
photon. It would be interesting to explore, if this signature can give
enhanced sensitivity to neutrino polarizability.

###  5.3 Bounds from MiniBoone

MiniBoone is an electron neutrino appearance experiment in which $\nu_{e}$ are
detected through quasi-elastic charged current interaction, with a typical
momentum exchange $q^{2}\sim-2~{}{\rm GeV}^{2}$. The initial neutrino flux is
composed at $99.4\%$ by $\nu_{\mu}+\bar{\nu}_{\mu}$ and peaks at
${\mathcal{O}}(500~{}{\rm MeV})$, with two modes of operation, the neutrino
and antineutrino modes. For details on neutrino fluxes we use Ref.
MiniBooNE:2008hfu ; MiniBooNE:2021bgc . The MiniBoone detector is filled with
pure mineral oil, CH2, which acts as both a target and a scintillator.

The signal of neutrino polarizability interactions is the Rayleigh operator
induced $\nu A\to\nu X_{A}+\gamma$ scattering, where $A$ is the initial
nucleus, while $X_{A}$ denotes the final states from either elastic or
inelastic scattering. Experimentally, the signature is similar to the
radiative up-scattering Bolton:2021pey ; Schwetz:2020xra ; Coloma:2017ppo ;
Atkinson:2021rnp ; Acero:2022wqg  and thus the same type of analyses would
also be sensitive to neutrino scattering through polarizability operators. In
the MiniBoone detector, however, this signature is indistinguishable from the
SM quasi-elastic charged current scattering of electron neutrino, $\nu_{e}A\to
eX_{A}$. The only difference is that the Rayleigh induced process leads to a
softer deposited energy spectrum due to the final state neutrino that escapes
the detector. For scattering on carbon we assume that it is dominated by
quasi-elastic scattering, i.e., by neutrino scattering on a single nucleon
bound inside the carbon nucleus, which then gets kicked out of the nucleus. In
the calculation of the total scattering rates we also include the scattering
on hydrogen, $\nu p\to\nu p+\gamma$, which constitutes a subdominant component
of the signal.

Figure 9: The number of measured MiniBoone events (blue dots and error bars)
in each photon/electron energy bin Aguilar-Arevalo:2020nvw , compared to the
expected SM signal+background (black line). Red line denotes the predicted
signal from a Rayleigh operator setting the NP scale to $\Lambda=6.3$ GeV (and
dimensionless coupling to unity), multiplying it by a factor 10 for clarity.

The photon spectrum is given by

$\frac{{d}N}{{d}E_{\gamma}}=N_{p}\int~{}{d}E_{\nu}\epsilon(E_{\gamma})\phi_{\rm
tot}(E_{\nu})\frac{{d}^{2}\sigma}{{d}E_{\gamma}{d}E_{\nu}}\,,$ (52)

where $N_{p}=2.8\times 10^{32}$ is the total number of target protons in the
detector, $\epsilon(E_{\gamma})$ is the acceptance as function of the photon
energy and $\phi_{\rm tot}(E_{\nu})$ is the total neutrino flux at the
detector. The integration is performed over the initial neutrino energies
interval, $E_{\nu}\in(0.05,7)$ GeV, with $\sim 98\%$ of the flux below
$E_{\nu}\leq 2$ GeV. The differential cross section
${d}^{2}\sigma/{d}E_{\gamma}{d}E_{\nu}$ for a neutrino scattering on a single
proton in the nucleus is calculated with MadGraph. We use this simple
approximation of scattering on single protons to evaluate our bounds; a more
refined calculation would involve modelling the nuclear responses of $C$ and
$H$ in the detector.

The bounds on neutrino polarizability is obtained using the chi-squared
function

$\chi^{2}(X)=\sum_{i}\frac{\left(N_{i}^{\rm MB}-(N_{i,bg}^{\rm MB}+N_{i}^{\rm
NP})\right)^{2}}{\sigma_{i}^{2}}\,,$ (53)

where the sum is over photon energy bins, with $N_{i}^{\rm MB}$ the number of
measured events in bin $i$, $N_{i,bg}^{\rm MB}$ the expected number of SM
background events, and $\sigma_{i}$ the experimental uncertainty Aguilar-
Arevalo:2018gpe ; Aguilar-Arevalo:2020nvw . The maximal allowed number of NP
events, $N_{i}^{\rm NP}$, places the bounds on the Rayleigh operator Wilson
coefficients.

$\left|\tfrac{2}{3}\hat{C}_{1}^{(7)}+\hat{C}_{2}^{(7)}\right|\lesssim 4\times
10^{-3}~{}{\rm GeV}^{-3}.$ (54)

In Fig. 9 we show the photon energy spectrum induced by the Rayleigh
interaction with the saturated bound and enlarged by a factor of 10 (red),
compared to the expected SM background (blue). The bound in (54) is dominated
by the large photon energy ”tail”, since the photon spectrum from the Rayleigh
interaction is broad and peaks at $\sim 500$ MeV.

For light mediator we keep the full $\phi$ propagator in the calculation of
the scattering cross sections, which then leads to the bounds

$c_{\nu}g_{\phi\gamma}\lesssim 1.2\times 10^{-5}\,,~{}~{}m_{\phi}=1~{}{\rm
GeV}\,;\qquad c_{\nu}g_{\phi\gamma}\lesssim 4\times
10^{-6}\,,~{}~{}m_{\phi}=1~{}{\rm MeV,keV,eV}\,.$ (55)

For the 1 GeV benchmark the $\phi$ mass is comparable to the typical momentum
exchange, while to the other three benchmarks the $\phi$ mass can be neglected
in the propagator.

While an enhanced $\nu A\to\nu X_{A}+\gamma$ scattering rate is an intriguing
possibility in view of the longstanding MiniBoone anomaly Aguilar-
Arevalo:2018gpe ; Aguilar-Arevalo:2020nvw , we note in passing that for a
massless neutrino in the final state the photon spectrum does not match the
observed low energy anomaly, see Fig. 9. For $\nu A\to NX_{A}+\gamma$
scattering, on the other hand, where $N$ is a heavier sterile neutrino, the
final state photon would be softer and could potentially match the MiniBoone
measurements (for sample of other new physics explanations of the MiniBoone
anomaly, see, e.g., Refs. Batell:2022xau ; Foppiani:2022qsi ;
Arguelles:2018mtc ; Kamp:2022bpt ; Vergani:2021tgc ; Dentler:2019dhz ;
Fischer:2019fbw ). We leave the full investigation of such a possibility for
future work.

###  5.4 Collider constraints

There are a number of constraints on neutrino polarizability from measurements
of higher energy processes, mostly from producing $\phi$ mediator on shell.

##### Rare meson and tau decays.

Couplings of $\phi$ to neutrinos generate the three body $M\to\ell\nu\phi$
decays of mesons $M=K,D,B,\ldots,$ via the emission of $\phi$ from the
neutrino leg. The decay is kinematically allowed for
$m_{\phi}<m_{M}-m_{\ell}$. For
$m_{\phi}=1\,\text{eV},1\,\text{keV},1\,\text{MeV}$ benchmarks the most
stringent constraint of this type comes from bounds on $K^{+}\to
e^{+}\nu_{e}\phi$ decays, giving $c_{\nu}\lesssim 4\times 10^{-3}$, while for
$m_{\phi}=1\,\text{GeV}$ the most stringent bound is from
$\tau\to\ell\nu\bar{\nu}\phi$ decays, leading to $c_{\nu}\lesssim 0.3$, see
Ref. Blinov:2019gcj  and references within. These constraints are shown as
dark gray shaded excluded regions in Fig. 7.

##### Neutrinoless double $\beta$ decay.

The neutrinoless double $\beta$ decay ($0\nu 2\beta$) experiments can be used
to also search for $0\nu 2\beta\phi$ transitions, where $\phi$ is emitted from
one of the neutrino lines. Present experimental bounds translate to a
constraint $c_{\nu}\lesssim 10^{-5}$ Blum:2018ljv  for $m_{\phi}\lesssim 2$
MeV, which we show as light gray shaded excluded region in Fig. 7.

##### Beam dump experiments.

Light pseudoscalars coupling to photons can be produced in electron and proton
beam dump experiments via Primakoff process, and are then searched for via
their decays to two photons. In our case $\phi$ has an additional invisible
decay channel $\phi\to\nu\nu$, which dilutes the $\gamma\gamma$ signal, if the
decays to neutrinos dominate. Rescaling the bounds on ALP couplings to
photons, $g_{\phi\gamma}(m_{\phi})$, from Ref. Dolan:2017osp  with the
diphotonic branching ratio ${\text{Br}}(\phi\to\gamma\gamma)$, Eq. (21), gives
the bound

$g_{\phi\gamma}(m_{\phi})<{g_{\phi\gamma}^{\rm
b.d.}(m_{\phi})}/{\sqrt{{\text{Br}}(\phi\to\gamma\gamma)}}\,.$ (56)

where $g_{\phi\gamma}^{\rm b.d.}(m_{\phi})$ is the bound quoted in
Dolan:2017osp . The constraint in (56) is shown in Fig. 7 as the dark red
excluded regions. The exception to rescaling rule in (56) is the newer NA64
analysis NA64:2020qwq  that included both recoils due to invisible $\phi$
escaping the detector as well as the $\phi\to\gamma\gamma$ events. Assuming
$\phi$ only has photon couplings the NA64 analysis is less sensitive than the
other beam dump experiments for the four $\phi$ mass benchmarks we consider.
For this reason and because the NA64 result is difficult to recast for the
more general case of an arbitrary invisible branching ratio, we do not include
it in Fig. 7.

##### Peripheral heavy ion collisions.

The production of ALP with couplings to photons is $Z^{4}$ coherently enhanced
in ultra-peripheral ion collisions Knapen:2016moh . The ALP was then searched
for in the two photon decay channel ATLAS:2019azn . This leads to a
significant bound only for heavy ALPs, with mass $m_{\phi}\gtrsim 10$ GeV,
thus we do not show it in Fig. 7.

##### Search for $e^{+}e^{-}\to 3\gamma$.

Belle II collaboration performed a search for ALPs decaying to two photons,
$e^{+}e^{-}\to\gamma(\phi\to\gamma\gamma)$ PhysRevLett.125.161806 , and set
stringent bounds on $g_{\phi\gamma}$ in the mass range $0.2\lesssim
m_{\phi}\lesssim 9.5$ GeV, assuming $Br(\phi\to\gamma\gamma)=100\%$. At the
mass benchmark $m_{\phi}=1$ GeV, the coupling to photons is constrained to be
$g_{\phi\gamma}\lesssim 10^{-3}~{}{\rm GeV}^{-1}$ for values $c_{\nu}$ small
enough that $Br(\phi\to\nu\nu)\ll Br(\phi\to\gamma\gamma)$, while for larger
$c_{\nu}$ we rescale the Belle II bound as in (56). The excluded region is
shown with dark red in Fig. 7. Note that a similar bound on the $m_{\phi}=1$
GeV mass benchmark follows from searches for anomalous $2\gamma$ and $3\gamma$
signal at LEP Jaeckel:2015jla .

##### Invisible decays of spin-0 particles.

The $\nu\nu\gamma\gamma$ interaction would induce $S\to\gamma\gamma\to\nu\nu$
decays, i.e., the $S\to\gamma\gamma$ transition leads at one loop to
$S\to\nu\nu$ decays, where for the initial spin-0 particle we consider
$S=\pi^{0},B^{0}$ and the Higgs boson, $h$, and assume that the EFT limit for
$\nu\nu\gamma\gamma$ interaction applies. The SM rates to $S\to\nu\nu$ are
negligible Bhattacharya:2018msv . If just one combination of Rayleigh Wilson
coefficients, $\hat{C}_{1}^{\rm Re}$ or $\hat{C}_{2}^{\rm Re}$ in Eq. (5),
contributes, then (cf. Appendix C)

$\big{|}\hat{C}_{a}^{\rm Re}\big{|}\leq
k_{a}\,\frac{128\pi^{2}}{\alpha}\frac{1}{m_{S}^{3}}\sqrt{\frac{\mbox{Br($S\to$
inv.)}}{\mbox{Br}(S\to\gamma\gamma)}}\,,$ (57)

where $k_{1}=3/2$ and $k_{2}=1$. This is the case for $\pi^{0}$ and $h$,
leading to

$\pi^{0}:\quad\big{|}\hat{C}_{2}^{\rm Re}\big{|}\leq 4.7\cdot 10^{3}\mbox{
GeV}^{\,-3},\qquad\qquad h:\quad\big{|}\hat{C}_{1}^{\rm Re}\big{|}\leq
1.2\mbox{ GeV}^{\,-3},$ (58)

from experimental bounds $\text{Br}(\pi^{0}\to\text{inv})\leq 4.4\cdot
10^{-9}$, $\text{Br}(h\to\text{inv})\leq 0.19$ Workman:2022ynf , along with
$\text{Br}(\pi^{0}\to\gamma\gamma)\approx 0.99$ and the SM prediction ${\rm
Br}_{\rm SM}(h\to\gamma\gamma)=2.3\cdot 10^{-3}$
LHCHiggsCrossSectionWorkingGroup:2013rie , using the fact that the Higgs
properties are consistent with the SM. From $\text{Br}(B^{0}\to\text{inv})\leq
2.4\cdot 10^{-5}$Workman:2022ynf  we obtain, on the other hand,

$B^{0}:\qquad\big{|}\tfrac{2}{3}\hat{C}_{1}^{\rm Re}+\hat{C}_{2}^{\rm
Re}\big{|}\leq 3.7\cdot 10^{4}\mbox{ GeV}^{\,-3},$ (59)

assuming $B^{0}\to\gamma\gamma$ is as predicted in the SM, see Appendix C for
details. The constraint in (59) may therefore change, if new physics affects
$B^{0}\to\gamma\gamma$ decays. Note that at present the above bounds are quite
weak, and the use of EFT may be questioned. For light mediator $\phi$, with
$m_{S}\neq m_{\phi}$, the above results still apply, but with replacement
$\hat{C}_{1(2)}^{\rm Re}\to
2\sum_{ij}\operatorname{Re}c_{\nu}^{ij}c_{\gamma}^{(^{\prime})}[f_{\phi}(m_{S}^{2}-m_{\phi}^{2})]^{-1}$
(the case $m_{S}\simeq m_{\phi}$ is more involved, and we do not attempt it
here).

##### Monophoton searches.

Neutrino polarizability leads to a monophoton signature in $e^{+}e^{-}$
collisions. This is either due to
$e^{+}e^{-}\to\gamma^{*}\to\gamma\,\nu\bar{\nu}$ scattering, generated by
Rayleigh operators in the EFT limit, or by an on-shell production of the light
$\phi$ mediator, $e^{+}e^{-}\to\gamma^{*}\to\gamma\,\phi$, where $\phi$ then
decays to two neutrinos or escapes the detector.

We recast the BaBar monophoton search BaBar:2017tiz  for the case of light
mediator $\phi$. The results in BaBar:2017tiz  were interpreted in terms of
the bounds on dark photon $A^{\prime}$ mixing parameter, $\varepsilon$. The
cross section for $e^{+}e^{-}\to A^{\prime}\gamma$ is in the
$m_{e}\ll\sqrt{s}$ limit given by Boehm:2003hm ; Borodatchenkova:2005ct

$\frac{d\sigma_{A^{\prime}\gamma}}{d\cos\theta}=\frac{2\pi\alpha^{2}\varepsilon^{2}}{s^{2}(s-m_{A^{\prime}}^{2})}\Big{[}\frac{s^{2}+m_{A^{\prime}}^{4}}{\sin^{2}\theta}-\frac{(s-m_{A^{\prime}}^{2})^{2}}{2}\Big{]},$
(60)

where $\theta$ is the emerging angle of the photon, restricted to
$|\cos\theta|<0.6$ in the BaBar analysis. For each benchmark value of
$m_{\phi}$ we compare the above cross section, integrated over
$\cos\theta\in[-0.6,0.6]$ and taking $m_{A^{\prime}}=m_{\phi}$, with the
$e^{+}e^{-}\to\gamma\phi$ cross section Dolan:2017osp ,

$\frac{d\sigma_{\phi\gamma}}{d\cos\theta}=\frac{\alpha^{3}}{512\pi^{2}}\left(\left|\frac{2}{3}\frac{c_{\gamma}}{f_{\phi}}\right|^{2}+\left|\frac{c^{\prime}_{\gamma}}{f_{\phi}}\right|^{2}\right)\left(1-\frac{m_{\phi}^{2}}{s}\right)^{3}\left(3+\cos
2\theta\right)\,,$ (61)

integrated over the same range $\cos\theta\in[-0.6,0.6]$. Using the bound
$\varepsilon<9.5\times 10^{-4}$ from BaBar:2017tiz , valid for all four
benchmark masses, then gives

$\left(\left|\frac{2}{3}\frac{c_{\gamma}}{f_{\phi}}\right|^{2}+\left|\frac{c^{\prime}_{\gamma}}{f_{\phi}}\right|^{2}\right){\cal
B}(\phi\to\nu\nu)\leq 0.018\mbox{ GeV}^{\,-2}\,,$ (62)

with branching ratio to neutrinos given by Eqs. (19), (20). The corresponding
excluded region is denoted with orange in Fig. 7.

##  6 UV models of enhanced neutrino polarizability

Next, we discuss several UV models that lead to enhanced contributions to
neutrino Rayleigh operators. In Section 6.1 we first review the minimal
singlet majoron model. This does not predict large neutrino polarizability,
but can be used as a useful benchmark. The other models we consider below, the
majoron as a QCD axion and the majoron from non-minimal inverse see-saw
models, discussed in Sections 6.2-6.4, have enhanced neutrino polarizabitily
relative to the minimal majoron model. All of the models rely on spontaneously
broken global lepton number, $U(1)_{L}$, and the associated Goldstone boson,
the majoron Chikashige:1980ui ; Gelmini:1980re , whose tree level exchanges
lead to enhanced neutrino polarizability, in the same way as for the
simplified model in Section 2.3.

###  6.1 Minimal singlet Majoron

The minimal singlet majoron model assumes that the SM neutrinos, $\nu_{i}$,
are Majorana fermions, and that their masses are suppressed as the result of
the type I see-saw, with spontaneously broken lepton number Chikashige:1980ui
. The SM is supplemented by three right-handed neutrinos $N_{R,i}$ and a
singlet scalar, $S$, that carries lepton number $L=-2$.333We use a short-
handed notation $N_{R,i}=P_{R}N_{i}$, where $N_{i}$ is the four component
Majorana fermion field, where as in the rest of the paper we use the
notational conventions from Ref. Dreiner:2008tw . The terms in the Lagrangian
relevant for the neutrino masses are thus given by

${\cal L}=-\bar{L}yN_{R}H-\frac{1}{2}\bar{N}_{R}^{c}\lambda N_{R}S+{\rm
h.c.},$ (63)

where $y,\lambda$ are $3\times 3$ complex matrices. The lepton number is
spontaneusly broken once $S$ obtains a vev,
$S=(f_{\phi}+\sigma+i\phi)/\sqrt{2}$. The radial mode $\sigma$ is assumed to
be heavy with mass $m_{\sigma}\sim{\mathcal{O}}(f_{\phi})$ and not relevant
for our discussion. The majoron, $\phi$, is the pNGB of the spontaneously
broken lepton number. Its mass term, ${\cal L}\supset-m_{\phi}^{2}\phi^{2}/2$,
represents a (small) explicit breaking of the shift symmetry, where $m_{\phi}$
is taken as a free parameter Gu:2010ys ; Frigerio:2011in  (it could arise from
Planck scale physics since gravity is expected to break global symmetries
Akhmedov:1992hi ; Rothstein:1992rh ; Alonso:2017avz ).

The SM neutrino masses induced by the interactions (63) are parametrically
given by $m_{\nu_{i}}\sim y^{2}v^{2}/\lambda f_{\phi}$ so that for $\lambda
f_{\phi}\gg v$ (the see-saw limit) the couplings $y$ can be large, which is
one of the main motivations for contemplating the see-saw models. The
couplings of majoron to the SM fermions are given by Heeck:2019guh

${\cal
L}_{\phi}\supset\frac{i\phi}{2f_{\phi}}m_{f_{i}}c_{f_{i}}\bar{f}_{i}\gamma_{5}f_{i}+\cdots.$
(64)

The majoron–SM-fermion couplings are thus suppressed by the SM fermion masses,
$m_{f_{i}}$, while the dimensionless coefficients are, in the see-saw limit,

$c_{\nu_{i}}=1,\qquad
c_{d_{i}}=-c_{u_{i}}=\frac{1}{16\pi^{2}}\operatorname{Tr}(yy^{\dagger}),\qquad
c_{\ell_{i}}=\frac{1}{16\pi^{2}}\big{[}\operatorname{Tr}(yy^{\dagger})-(yy^{\dagger})_{ii}\big{]},$
(65)

for the couplings to the SM neutrinos, down and up quarks, and the charged
leptons, respectively. The ellipses in (64) denote the flavor off-diagonal
terms for charged leptons that we do not display, but can be found in
Heeck:2019guh . Even though the majoron couplings to the SM neutrinos are
generated at tree level, while the couplings to quarks and charged leptons are
generated at one loop, cf. Eq. (64), the couplings to quarks and charged
leptons can still be larger, if the Yukawa couplings $y$ are sizable. That is,
the couplings to quarks and charged leptons are parametrically enhanced for
$f_{\phi}\gg v$ by $\sim\lambda f_{\phi}m_{f_{i}}/(4\pi v)^{2}$ compared to
the $\phi$ couplings to neutrinos.

Couplings of majoron to photons and gluons are generated at two loops. In the
limit of a light majoron, $m_{\phi}\ll m_{\ell_{i},u_{i},d_{i}}$, they match
onto the dimension 7 operators
$\big{(}\partial^{2}\phi\big{)}F_{\mu\nu}\tilde{F}^{\mu\nu}$ and
$\big{(}\partial^{2}\phi\big{)}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}$ and are thus
suppressed both by the two loop factors $\alpha_{{\rm em}}/(16\pi^{2})^{2}$
and $\alpha_{s}/(16\pi^{2})^{2}$, respectively, as well as by the majoron
mass, $m_{\phi}^{2}/m_{f_{i}}^{2}$, for a contribution from a SM charged
fermion $f_{i}$ running in the loop. For a heavier majoron the latter
suppression is lifted and in the corresponding transition amplitude
$m_{\phi}^{2}/m_{f_{i}}^{2}$ is replaced by an ${\mathcal{O}}(1)$ factor.

###  6.2 Majoron as a QCD axion

If the spontaneously broken $U(1)_{L}$ is anomalous under QCD, the majoron
will act as the QCD axion and solve the strong CP problem, with $U(1)_{L}$
identified as the Peccei-Quinn symmetry Mohapatra:1982tc ; Langacker:1986rj ;
Shin:1987xc ; Ma:2001ac ; Ma:2017vdv . Recent concrete realization of this
idea can be found in Ma:2017vdv , where the SM was supplemented by a set of
color octet fermions $\Psi_{R}^{A}$ (with lepton number charge $L=1$,
electroweak singlets), color octet scalars, $\Phi^{A}$ ($L=0$, doublets of
$SU(2)_{L}$), and an electroweak and color singlet scalar $S$ ($L=2)$. The
lepton number is spontaneously broken once $S$ obtains a vev, $\langle
S\rangle=f_{\phi}/\sqrt{2}$, giving Majorana mass contribution,
$M_{\Psi_{k}}=y_{\Psi}^{k}\langle S\rangle$, to the color octet fermion mass
matrix as a result of the term ${\cal
L}\supset-\frac{1}{2}y_{\Psi}^{i}S\bar{\Psi}_{iR}^{c}\Psi_{iR}+{\rm h.c.}$, in
the Lagrangian.

At one loop the breaking of $U(1)_{L}$ generates the Majorana mass matrix for
the SM neutrinos. In the limit of heavy color octet scalars, with almost
degenerate mass $M_{0}\gg M_{\Psi_{k}}\gg v$, the radiatively generated
neutrino mass matrix is given by Ma:2017vdv ,

$\big{(}m_{\nu}\big{)}_{ij}=\frac{1}{4\pi^{2}}\sum_{k}h_{\Psi}^{ik}h_{\Psi}^{jk}y_{\Psi}^{k}\langle
S\rangle\frac{\Delta M_{0}^{2}}{M_{0}^{2}},$ (66)

where the summation is over color octet fermion, and $\Delta
M_{0}^{2}=2\lambda_{5}v^{2}$ is the mass splitting between the CP-even and CP-
odd components of the color octet scalar due to the term ${\cal
L}\supset-\frac{1}{2}\lambda_{5}(H^{\dagger}\Phi^{A})^{2}+{\rm h.c.}$ in the
scalar potential. Note that the neutrino masses are proportional to the
$\Delta L=2$ breaking vev, $m_{\nu}\propto\langle S\rangle$, to two insertions
of electroweak vev, $m_{\nu}\propto v^{2}$, as well as to the Yukawa
interactions between the SM lepton doublets and the new color octet fields,
${\cal L}\supset
h_{\Psi}^{ij}{}^{*}\bar{L}_{i}\Psi_{jR}^{A}\Phi^{A\dagger}+{\rm h.c.}$. Taking
$f_{\phi}\sim{\mathcal{O}}({\rm TeV})$ and all couplings to be
${\mathcal{O}}(10^{-2})$ the observed neutrino masses are obtained for
$M_{0}\sim{\mathcal{O}}(10\,{\rm TeV})$.

The majoron $\phi$ is part of the singlet,
$S(x)=\frac{1}{\sqrt{2}}\big{(}f_{\phi}+\sigma(x)\big{)}\exp\big{(}i\phi(x)/f_{\phi}\big{)}$,
and has interactions to neutrinos suppressed by the neutrino masses,

${\cal
L}_{\phi}\supset-\frac{i}{2}(m_{\nu})_{ij}\big{(}\bar{\nu}_{i}P_{L}\nu_{j}\big{)}\frac{\phi}{f_{\phi}}+{\rm
h.c.}.$ (67)

In the notation of Eq. (15) the coupling of $\phi$ to neutrinos is thus given
by

$c_{\nu}^{ij}=-i\frac{(m_{\nu})_{ij}}{f_{\phi}}.$ (68)

The triangle anomaly induces couplings of the majoron to gluons,

${\cal
L}_{\phi}\supset-\frac{3n_{\Psi}\alpha_{s}}{8\pi}\frac{\phi}{f_{\phi}}\tilde{G}_{\mu\nu}^{A}G^{A\mu\nu},$
(69)

where $n_{\Psi}$ is the number of color octet fermions. The couplings to
gluons are not suppressed by the majoron mass, since they are generated from
the QCD anomaly, unlike the minimal singlet majoron case, Section 6.1. At low
energies, below QCD confinement, the above interaction induces interactions of
$\phi$ with nucleons and pions, and other hadronic states. It also leads to
the majoron mass, in the same way as for the standard QCD axion,

$m_{\phi}\simeq 6\,{\rm keV}\times\biggr{(}\frac{1\,{\rm
TeV}}{f_{\phi}/(3n_{\Psi})}\biggr{)}.$ (70)

as well as to the couplings to photons, Eq. (15), with

$c_{\gamma}^{\prime}\simeq 2.0\times 3n_{\Psi}.$ (71)

At energies below $m_{\phi}$ the majoron can be integrated out, giving rise to
the neutrino polarizability with

$\frac{C_{2,ij}^{(7)}}{\Lambda^{3}}=i\frac{(m_{\nu})_{ij}c_{\gamma}^{\prime}}{f_{\phi}^{2}m_{\phi}^{2}}\simeq
i\biggr{(}\frac{1}{81\,{\rm
GeV}}\biggr{)}^{3}\times\frac{1}{n_{\Psi}}\times\frac{(m_{\nu})_{ij}}{0.1\,{\rm
eV}}.$ (72)

Note that the dependence on $f_{\phi}$ drops out, due to the relation between
$m_{\phi}$ and $f_{\phi}$, with $m_{\phi}f_{\phi}\sim m_{\pi}f_{\pi}$ a
constant that is fixed entirely by the QCD dynamics apart from the
$n_{\Psi}\sim{\mathcal{O}}(1)$ factors that depends on the UV physics. The PQ
breaking scale $f_{\phi}$ thus does not determine directly the effective
suppression scale of the Rayleigh operator, $\Lambda$, but rather just
determines the range of the validity of the EFT, via the requirement $E,q\gg
m_{\phi}\propto 1/f_{\phi}$. The effective scale $\Lambda$ is given by a
combination of QCD and neutrino mass scales,
$\Lambda\sim\big{(}m_{\pi}^{2}f_{\pi}^{2}/m_{\nu}\big{)}^{1/3}$, which
accidentally turns out to be close to the weak scale.

While the majoron that is the QCD axion is an example of the model that leads
to enhanced neutrino polarizability, with the effective scale $\Lambda$ much
smaller than the UV scales $f_{\phi}$ and $M_{0}$, it does not map
straightforwardly onto the phenomenological analyses forming the bulk of the
present paper. A QCD axion couples to gluons, Eq. (69), which was not taken
into account in our analysis.

In the final two examples, the non-minimal majoron models in Sections 6.3 and
6.4 the parametric suppressions encountered in the minimal see-saw majoron and
the majoron as the QCD axion model are lifted.

###  6.3 Majoron from inverse see-saw with extra triplet fermions

In this model both the neutrino mass generation sector is enlarged as is the
sector that leads to couplings of majoron to the photons.

####  6.3.1 The inverse see-saw sector

The sector relevant for the generation of the SM neutrino masses contains
three generations of left-handed and right-handed Weyl fermions, $N_{R,i}$ and
$N_{L,i}$, $i=1,2,3$, singlets under the SM gauge group, and carry a global
lepton number $L=+1$. The interaction Lagrangian is given by (suppressing
generation indices)

$-{\cal
L}_{Y}=y_{\nu}\bar{L}\tilde{H}^{\dagger}N_{R}+\bar{N}_{L}M_{N}N_{R}+\frac{\lambda_{R}}{2}\bar{N}_{R}^{c}SN_{R}+\frac{\lambda_{L}}{2}\bar{N}_{L}SN_{L}^{c}+{\rm
h.c.},$ (73)

where $L_{i}$ are the SM lepton doublets, and $H$ is the SM Higgs doublet. The
Yukawa coupling $y_{\nu}$ and neutrino Dirac mass matrix $M_{N}$ are general
$3\times 3$ complex matrices, while $\lambda_{L,R}$ are symmetric $3\times 3$
complex matrices.

The lepton number is spontaneusly broken once $S$ obtains a vev,

$S=(f_{\phi}+\sigma+i\phi)/\sqrt{2}.$ (74)

This then gives the following neutrino mass matrix in the basis
$\\{\nu_{L}^{c},N_{R},N_{L}^{c}\\}$, see, e.g. Dev:2012sg ,

${\cal M}=\begin{pmatrix}0&M_{D}&0\\\ M_{D}^{\top}&\mu_{R}&M_{N}^{\top}\\\
0&M_{N}&\mu_{L}\end{pmatrix},$ (75)

where we shortened the notation to $M_{D}=y_{\nu}v/{\sqrt{2}}$,
$\mu_{L,R}=\lambda_{L,R}f_{\phi}/\sqrt{2}$. Without loss of generality we can
work in the basis, where $M_{N}$ is diagonal. We will assume the hierarchy
$\mu_{R,L}\ll M_{D}\ll M_{N}$, to be understood as the hierarchy among all the
eigenvalues of the corresponding $3\times 3$ complex matrices.

We start the analysis with the one-generation case in which case the three
parameters $\mu_{L},M_{D},M_{N}$ can be made real via phase redefinitions of
$\nu_{L}^{c},N_{R},N_{L}^{c}$, while $\mu_{R}$ is a complex parameter. The
lightest neutrino mass eigenstate is then, up to corrections of higher order
in $\mu_{L,R}/M_{N,D}$, given by

$\nu=c_{\alpha}\nu_{L}-s_{\alpha}N_{L},\qquad
t_{\alpha}\simeq\frac{M_{D}}{M_{N}},$ (76)

where we abbreviated $c_{\alpha}\equiv\cos\alpha$,
$s_{\alpha}\equiv\sin\alpha$, $t_{\alpha}\equiv\tan\alpha$. In order for $\nu$
to be predominantly composed of the neutral component of the electroweak
doublet, with only a small admixture of the sterile neutrino, we require
$M_{D}\ll M_{N}$. The corresponding mass is

$m_{\nu}\simeq\mu_{L}s_{\alpha}^{2}.$ (77)

The two heavy neutrinos are mass degenerate, with masses equal to
$\sqrt{M_{N}^{2}+M_{D}^{2}}$, up to $\mu_{L,R}$ suppressed corrections. Note
that the SM neutrino masses are proportional to the lepton number breaking
parameter $\mu_{L}=\lambda_{L}f_{\phi}/\sqrt{2}$, and vanish in the limit
$f_{\phi}\to 0$, as expected. The neutrino masses can now be small either due
to a small value of $f_{\phi}$, the smallness of the mixing angle $s_{\alpha}$
(i.e., $M_{D}/M_{N}\ll 1$), or a combination of the two.

These results extend trivially to the case of three generations. To linear
order in $M_{D}/M_{N}$ the mass eigenstates, the left-handed fields
$\nu,N_{1},N_{2}$, are expressed in terms of the initial states as

$\displaystyle\nu_{L}^{c}$
$\displaystyle=\nu^{c}+\big{(}M_{N}^{-1}\big{)}^{T}M_{D}\frac{1}{\sqrt{2}}\big{(}N_{1}^{c}+N_{2}^{c}\big{)},$
(78) $\displaystyle N_{L}^{c}$
$\displaystyle=-M_{N}^{-1}M_{D}^{T}\nu^{c}+\frac{1}{\sqrt{2}}\big{(}N_{1}^{c}+N_{2}^{c}\big{)},$
(79) $\displaystyle N_{R}$
$\displaystyle=\frac{1}{\sqrt{2}}\big{(}N_{1}^{c}-N_{2}^{c}\big{)}.$ (80)

The light neutrino mass terms are then given by

${\cal
L}\supset\frac{1}{2}\bar{\nu}M_{D}\big{(}M_{N}^{-1}\big{)}^{T}\mu_{L}M_{N}^{-1}M_{D}^{T}P_{R}\nu^{c}+{\rm
h.c.}$ (81)

and are proportional to $\mu_{L}$ lepton number violating parameter, while the
dependence on $\mu_{R}$ only enters at higher orders.

The interactions with the majoron can be obtained by replacing $\mu_{L}\to
i\mu_{L}\phi/f_{\phi}$ in the mass terms, which then gives for the majoron
couplings to neutrinos,

${\cal
L}_{\phi}\supset-\frac{i}{2}(m_{\nu})_{ij}\big{(}\bar{\nu}_{i}P_{L}\nu_{j}\big{)}\frac{\phi}{f_{\phi}}+{\rm
h.c.}.$ (82)

In the notation of Eq. (15) the coupling of $\phi$ to neutrinos is thus given
by

$c_{\nu}^{ij}=-i\frac{(m_{\nu})_{ij}}{f_{\phi}}.$ (83)

On the face of it, this is the same result as for the minimal majoron, cf.
Eqs. (67), (68). However, there is a major difference between the two, namely
that in the inverse see-saw model the smallness of neutrino masses, $m_{\nu}$,
can be due to the smallness of $f_{\phi}$. In principle, $f_{\phi}$ can be as
small as $m_{\nu}$ and thus $c_{\nu}$ as large as
$c_{\nu}\sim{\mathcal{O}}(1)$. Numerically, the bounds on self-interactions of
neutrinos limit the value of $c_{\nu}$ to be well below 1, cf. Section 3.

####  6.3.2 Couplings to photons via heavy electroweak triplets

We assume that the field content of the theory contains a set of heavy
$SU(2)_{L}$ triplet fermions $\Psi_{R}^{a}$, charged under lepton number,
$L=-1$. They obtain their masses through interactions with the scalar $S$,

${\cal
L}_{\Psi_{R}}=-\frac{1}{2}y_{\Psi}S\overline{\big{(}\Psi_{R}^{a}\big{)}^{c}}\Psi_{R}^{a}+{\rm
h.c.},$ (84)

after $S$ obtains a vev, Eq. (74), giving $M_{\Psi}=y_{\Psi}\langle
S\rangle=y_{\Psi}f_{\phi}/\sqrt{2}$. This generates a coupling of PNGB $\phi$
with the $SU(2)_{L}$ gauge fields through anomaly,

${\cal
L}_{\phi}\supset-\frac{9\alpha}{64\pi}n_{\Psi}\frac{\phi}{f_{\phi}}\tilde{W}_{\mu\nu}^{a}W^{a\mu\nu},$
(85)

where $n_{\Psi}$ is the number of $\Psi_{R}^{a}$ generations, and
$W^{a}_{\mu\nu}$ is the $SU(2)_{L}$ field strength. This gives couplings of
$\phi$ to $W^{\pm}$, $Z^{0}$ and photons, where for the latter

${\cal
L}_{\phi}\supset-\frac{9\alpha}{64\pi}n_{\Psi}s_{w}^{2}\frac{\phi}{f_{\phi}}\tilde{F}_{\mu\nu}F^{\mu\nu},$
(86)

where $s_{w}=\sin\theta_{w}$, with $\theta_{w}$ the weak mixing angle.

Since $\Psi_{R}^{a}$ carry electroweak charges they could be produced in
$e^{+}e^{-}$ collisions at LEP or in $pp$ collisions at the LHC. The bounds on
their mass depends on the $\Psi_{R}^{a}$ decay channels, and is thus model
dependent. In general, we expect the bound to be in the range of a few 100
GeV. Taking this as the typical lower bound on $f_{\phi}$, this would then
translate to the following typical size of the Rayleigh operator,

$\begin{split}\frac{C_{2,ij}^{(7)}}{\Lambda^{3}}&=i\frac{(m_{\nu})_{ij}c_{\gamma}^{\prime}}{f_{\phi}^{2}m_{\phi}^{2}}=i\frac{9}{8}\frac{(m_{\nu})_{ij}n_{\Psi}s_{w}^{2}}{f_{\phi}^{2}m_{\phi}^{2}}\\\
&\simeq
i\biggr{(}\frac{1}{7.3\,\text{GeV}}\biggr{)}^{3}\biggr{(}\frac{100\,\text{GeV}}{f_{\phi}}\biggr{)}^{2}\times\biggr{(}\frac{1\,\text{keV}}{m_{\phi}}\biggr{)}^{2}\times
n_{\Psi}\times\frac{(m_{\nu})_{ij}}{0.1\,{\rm eV}}.\end{split}$ (87)

That is, for majoron mass $m_{\phi}\sim{\mathcal{O}}(\text{keV})$ the
effective suppression scale of the Rayleigh operator is only
$\Lambda\sim\mathcal{O}(10\,\text{GeV})$ for
$m_{\phi}\sim\mathcal{O}(1\,\text{keV})$, and is parametrically smaller than
the $U(1)_{L}$ breaking scale, $f_{\phi}$. For instance, for
$m_{\phi}\sim\mathcal{O}(1\,\text{eV})$ the effective suppression scale would
be $\Lambda\sim\mathcal{O}(0.1\,\text{GeV})$. The majoron mass is an explicit
$U(1)_{L}$ breaking term and is treated as a free parameter. On general
grounds one expects $m_{\phi}\ll f_{\phi}$ so that the explicit breaking is
smaller than the spontaneous symmetry breaking, and thus $U(1)_{L}$ is a good
approximate symmetry.

###  6.4 Enhanced neutrino polarizability from $U(1)_{L}\times U(1)^{\prime}$

In non-minimal versions of the above model one can obtain even larger
parametric enhancements of the neutrino polarizability. Let us consider an
example of a model with two global $U(1)$ factors, $U(1)_{L}\times
U(1)^{\prime}$, a simple modification of the inverse see-saw model in Section
6.3, but with two different scalars in the inverse see-saw sector and the
electroweak triplet sector. That is, the model contains two scalar SM gauge
singlets, $S$ and $S^{\prime}$, where the scalar $S$ carries a global charge
$L=2$ under $U(1)_{L}$, while $S^{\prime}$ carries a charge of $+2$ under
$U(1)^{\prime}$. Otherwise the field content is the same as in Section 6.3.
The SM is extended by three generations of left-handed and right-handed
sterile neutrinos, $N_{R,i}$ and $N_{L,i}$, $i=1,2,3$, singlets under the SM
gauge group, with global lepton number $L=+1$, and a set of $n_{\Psi}$ weak
triplets that carry a charge $-1$ under $U(1)^{\prime}$.

The $U(1)_{L}$ and $U(1)^{\prime}$ are broken by $S$ and $S^{\prime}$ once
these obtain a vev, $S=(f_{\phi}+\sigma+i\phi)/\sqrt{2}$ and
$S^{\prime}=(f_{\phi}^{\prime}+\sigma^{\prime}+i\phi^{\prime})/\sqrt{2}$. In
general, the two vevs can be very different, which can be phenomenologically
beneficial. In particular, it is possible to have $f_{\phi}\ll
f_{\phi}^{\prime}$, which would explain the smallness of neutrino masses.444We
set aside the question of a hierarchy problem in the scalar potential. In
general one would need to assume that some terms, such as
$S^{\dagger}SS^{\prime\dagger}S^{\prime}$ or $S^{\dagger}SH^{\dagger}H$, are
suppressed by small couplings. Such a hierarchy of vevs would then also lead
to an enhancement of Rayleigh operators.

More explicitly, the Lagrangian of the model is the same as in Section 6.3,
except that in subsection 6.3.2 we should replace $S\to S^{\prime}$,
$\phi\to\phi^{\prime}$, $f_{\phi}\to f_{\phi}^{\prime}$. The interactions of
$\phi$ and $\phi^{\prime}$ with the SM particles are thus given by

${\cal L}_{\rm
int}\supset-\frac{i}{2}(m_{\nu})_{ij}\big{(}\bar{\nu}_{i}P_{L}\nu_{j}\big{)}\frac{\phi}{f_{\phi}}+{\rm
h.c.},$ (88)

and

${\cal L}_{\rm
int}\supset-\frac{9\alpha}{64\pi}n_{\Psi}\frac{\phi^{\prime}}{f_{\phi}^{\prime}}\tilde{W}_{\mu\nu}^{a}W^{a\mu\nu}\to-\frac{9\alpha}{64\pi}n_{\Psi}s_{w}^{2}\frac{\phi^{\prime}}{f_{\phi}^{\prime}}\tilde{F}_{\mu\nu}F^{\mu\nu}.$
(89)

As in Section 6.3, we assume that the dominant explicit breaking of the global
symmetry $U(1)_{L}\times U(1)^{\prime}$ is given by the masses of the two
PNGBs. Completely generally, the mass term is given

${\cal L}_{m}=-\frac{1}{2}m_{\phi}^{2}\phi^{2}-\frac{1}{2}m_{\phi}^{\prime
2}\phi^{\prime 2}-m_{\phi\phi^{\prime}}^{2}\phi\phi^{\prime},$ (90)

leading to two mass eigenstates, $m_{\phi_{1},\phi_{2}}$, with
$\phi_{1}=c_{\theta}\phi+s_{\theta}\phi^{\prime}$,
$\phi_{2}=-s_{\theta}\phi+c_{\theta}\phi^{\prime}$, where
$c_{\theta}=\cos\theta$, $s_{\theta}=\sin_{\theta}$, and $\theta$ the mixing
angle. In principle, all the explicit symmetry breaking terms in (90) can be
of comparable size, and thus the mixing angle large.

It is instructive to calculate the $\nu\nu\to\gamma\gamma$ amplitude due to
tree level exchanges of $\phi_{1,2}$,

${\cal
M}_{\nu\nu\to\gamma\gamma}=i\frac{9}{8}\frac{(m_{\nu})_{ij}}{f_{\phi}f_{\phi}^{\prime}}n_{\Psi}s_{w}^{2}c_{\theta}s_{\theta}\biggr{(}\frac{1}{q^{2}-m_{\phi_{1}}^{2}}-\frac{1}{q^{2}-m_{\phi_{2}}^{2}}\biggr{)},$
(91)

where $q^{\mu}$ is the sum of the initial neutrino momenta. In the center of
mass of the $\nu\nu$ collision it is given by $q^{\mu}=(2E_{\nu},\vec{0})$,
where $E_{\nu}$ is the neutrino energy. For high energy collisions, $q^{2}\gg
m_{\phi_{1,2}}^{2}$, the scattering amplitude vanishes as ${\cal M}\propto
m_{\phi_{1,2}}^{2}/q^{4}$, while for low energy processes, $q^{2}\ll
m_{\phi_{1,2}}^{2}$, it matches onto the neutrino polarizability operator,
with

$\frac{C_{2,ij}^{(7)}}{\Lambda^{3}}=i\frac{9}{8}\frac{(m_{\nu})_{ij}}{f_{\phi}f_{\phi}^{\prime}}n_{\Psi}s_{w}^{2}c_{\theta}s_{\theta}\biggr{(}\frac{1}{m_{\phi_{2}}^{2}}-\frac{1}{m_{\phi_{1}}^{2}}\biggr{)}.$
(92)

Numerically, since $f_{\phi}$ and $f_{\phi}^{\prime}$ can take very different
values, the Rayleigh operator can be suppressed by a light effective scale.
For instance, taking $m_{\phi_{1}}\ll m_{\phi_{2}}$ for simplicity,

$\frac{\big{|}C_{2,ij}^{(7)}\big{|}}{\Lambda^{3}}\simeq\biggr{(}\frac{1}{16\,\text{MeV}}\biggr{)}^{3}\biggr{(}\frac{100\,\text{GeV}}{f_{\phi}^{\prime}}\biggr{)}\biggr{(}\frac{1\,\text{keV}}{f_{\phi}}\biggr{)}\biggr{(}\frac{1\,\text{keV}}{m_{\phi_{1}}}\biggr{)}^{2}\times
n_{\Psi}c_{\theta}s_{\theta}\times\biggr{(}\frac{(m_{\nu})_{ij}}{0.1\,{\rm
eV}}\biggr{)}.$ (93)

Even lower effective scales than shown in the above numerical example, can
thus be obtained if $f_{\phi}$ and $m_{\phi_{1}}$ are smaller. For instance,
if they are comparable with the neutrino masses, $f_{\phi}\sim
m_{\phi_{1}}\sim{\mathcal{O}}(\text{eV})$, the effective scale would be
$\Lambda\sim{\mathcal{O}}(\text{keV})$.

Assuming $m_{\phi_{1}}\ll m_{\phi_{2}}$ the above model matches onto the
simplified model singlet mediator models for enhanced neutrino polarizability,
Section 2.3, with $\phi_{1}$ playing the role of the light mediator,
$\phi$,555Note that $\phi$ from now on denotes the single mediator from
Section 2.3 and not the field from the beginning of this subsection, Eqs.
(88)-(90). and

$c_{\nu}^{ij}=i\frac{(m_{\nu})_{ij}}{f_{\phi}}c_{\theta},\qquad
g_{\phi\gamma}=\frac{9\alpha}{16\pi}\frac{n_{\Psi}s_{w}^{2}}{f_{\phi}^{\prime}}s_{\theta},$
(94)

the nonzero coefficients in the Lagrangian (15) (we use the notation in Eq.
(18)). The contributions from heavier $\phi_{2}$ state are suppressed. Note
that $\phi=\phi_{1}$ has both flavor diagonal and off-diagonal couplings to
neutrinos, and thus only approximately matches onto the constraints shown in
Fig. 3 in which flavor universal neutrino couplings were assumed, cf. Eq.
(17). Nevertheless, the constraints on $c_{\nu}$ shown in Fig. 3 should
approximate well the constraints on $c_{\nu}^{ij}$ from Eq. (94). Numerically,

$\big{|}c_{\nu}^{ij}|=10^{-4}c_{\theta}\biggr{(}\frac{(m_{\nu})_{ij}}{0.1\,\text{eV}}\biggr{)}\biggr{(}\frac{1\,\text{keV}}{f_{\phi}}\biggr{)},\qquad
g_{\phi\gamma}=3\times
10^{-9}\text{GeV}^{-1}\biggr{(}\frac{100\,\text{GeV}}{f_{\phi}^{\prime}}\biggr{)}\biggr{(}\frac{s_{\theta}n_{\Psi}}{10^{-3}}\biggr{)},$
(95)

and thus for $m_{\phi}=1\,\text{keV},1\,\text{MeV}$ mass benchmarks the
$U(1)_{L}\times U(1)^{\prime}$ model can cover the whole experimentally still
available parameter space in Fig. 3, even when imposing $f_{\phi}>m_{\phi}$.
For the $m_{\phi}=1\,\text{GeV}$ mass benchmark the main constraint on
couplings to photons would be from searches for on shell production of
electroweak triplets, limiting $f_{\phi}^{\prime}$ to be above several 100
GeV, and thus $g_{\phi\gamma}\lesssim 10^{-6}\,\text{GeV}^{-1}$ even for large
mixing angles, $\sin\theta\sim{\mathcal{O}}(1)$. For couplings to the
neutrinos the requirement $f_{\phi}>m_{\phi}$ leads to
$\big{|}c_{\nu}^{ij}|\lesssim 10^{-7}$.

##  7 Conclusions

In this paper, we examined the theory and phenomenology of New Physics sources
of neutrino polarizability, that is, the electromagnetic interaction of
neutrinos with two photons. We present the latter in an EFT framework, where
these interactions are described by the Rayleigh operators; see Eq. (1). While
naively, one would expect these chirality-flipping operators to be suppressed
by the neutrino mass, we show that such a suppression can be compensated if
the interaction is mediated by a light scalar or pseudo-scalar particle.

Such models can have a wide variety of phenomenological consequences,
depending on the new particle’s mass. We fix four mass benchmarks,
$m_{\phi}=1$ eV, keV, MeV, GeV, and explore the constraints in the parameter
space defined by the coupling to photons, $g_{\phi\gamma}$, and neutrinos,
$c_{\nu}$. To provide results easily comparable with existing literature, we
limit ourselves to the case of a pseudoscalar mediator, with flavor-universal
couplings to neutrinos. The main results are summarized in Tables 1, 2 and 3,
and in Fig. 3.

The first three mass benchmarks, $m_{\phi}=1$ eV, keV, MeV, are largely
excluded by cosmological and astrophysical observables. This is to be expected
since CMB spectrum measurements strongly constrain the number of relativistic
degrees of freedom in the Universe, $N_{{\rm eff}}$, up to scales of order
$T\sim 100$ eV. Similarly, $N_{{\rm eff}}$ affects the abundance of primordial
elements produced during BBN, which takes place when the Universe temperature
is $T\sim 1-2$ MeV. Finally, exotic emission of neutrino and light new
particles can be excluded by measurements of Horizontal Branch star cooling
rates and neutrino fluxes from SN1987a. The former has a typical temperature
of $T_{HB}\sim 8$ keV, while the latter can reach $T_{SN}\sim 30$ MeV in the
inner core, largely setting the $\phi$ mass reach of the corresponding bounds.

We performed comprehensive analysis for coupling values in rather large
ranges, $c_{\nu}\in[10^{-12},1]$ and
$g_{\phi\gamma}\in[10^{-12}\,\text{GeV}^{-1},1\,\text{GeV}^{-1}]$. In this
parameter space, the eV benchmark is completely excluded. The heavier
benchmarks, keV and MeV, allow for very small coupling values, mainly due to
the disappearance of CMB bounds. The former benchmark is not excluded for
$c_{\nu}\lesssim 10^{-6}$ and $g_{\phi\gamma}\cdot{\rm GeV}\lesssim 10^{-11}$,
while the latter is not excluded for $c_{\nu}\lesssim 10^{-9}$ and
$g_{\phi\gamma}\cdot{\rm GeV}\lesssim 10^{-9}$.

When the mediator mass is heavier than the typical scales of cosmological and
astrophysical processes, we expect these bounds to disappear or become
negligible. This is evident for the last mass benchmark, $m_{\phi}=1$ GeV,
where only terrestrial experiments are able to probe parts of the parameter
space. Rare $\tau$ lepton decays bound $c_{\nu}\lesssim 0.3$, while monophoton
search in $e^{+}e^{-}$ collision leads to $g_{\phi\gamma}\cdot{\rm
GeV}\lesssim 10^{-4}$.

We also discussed UV complete models that lead to enhanced neutrino
polarizability, all of which are based on the appearance of a pseudo-Nambu-
Goldstone boson associated with the spontaneous breaking of the global lepton
number symmetry, $U(1)_{L}$. Such a pNGB, the majoron, automatically couples
to neutrinos since it participates in neutrino mass generation. In cases where
majoron has enhanced couplings to photons, the tree-level exchanges of the
majoron result in a parametrically enhanced neutrino polarizability, large
enough to saturate the present experimental bounds. Such enhanced couplings to
photons are, for instance, generated if majoron couples to a separate sector
of heavy charged states. A concrete example of a model in which all such
parametric enhancements are present is the non-minimal inverse see-saw model
discussed in detail in Section 6.4, which has two global symmetries, the
lepton number $U(1)_{L}$ and the anomalous $U^{\prime}(1)$. The neutrino
polarizability is then generated with very low effective scale suppression,
despite being suppressed by the neutrino masses.

Acknowledgements: We thank R. Harnik, J. Kopp, P. Machado and E. Vitagliano
for useful discussions. MT acknowledges the financial support from the
Slovenian Research Agency (research core funding No. P1-0035). JZ acknowledges
support in part by the DOE grant DE-SC0011784. GP and AAP were supported in
part by the DOE grant DE-SC0007983.

##  Appendix A Notations and conventions

Throughout the manuscript we use the four-component notation following the
conventions of Ref. Dreiner:2008tw . For Majorana neutrinos we thus have,

$\nu=\nu^{c}=\begin{pmatrix}\xi_{\alpha}\\\
\xi^{\dagger\dot{\alpha}}\end{pmatrix},\qquad\bar{\nu}=\bar{\nu}^{c}=\big{(}\xi_{\alpha},\xi^{\dagger\dot{\alpha}}\big{)},$
(96)

where $\xi_{\alpha}$ is a two-component Weyl spinor, so that, for instance,

$\bar{\nu}_{i}P_{L}\nu_{j}=\xi_{i}\xi_{j}.$ (97)

The normalization we use for the dimension 5 dipole operators and the
dimension 7 operators in the EFT Lagrangian, Eq. (1), is straightforwardly
related to other notations commonly used in the literature. The neutrino
dipole moments are conventionally defined as

${\cal L}_{\rm
eff}\supset\sum_{i>j}\frac{1}{2}(\lambda_{\nu})_{ij}(\bar{\nu}_{i}\sigma^{\mu\nu}P_{L}\nu_{j})F_{\mu\nu}+{\rm
h.c.},$ (98)

where the antisymmetric $3\times 3$ matrix, $\lambda=\mu-id$, decomposes into
the magnetic ($\mu$) and electric ($d$) dipole moments, see, e.g., Ref.
Miranda:2019wdy . In terms of the Wilson coefficients in (1) we have

$(\lambda_{\nu})_{ij}\mu_{B}\equiv\frac{C_{1,ij}^{(5)}}{\Lambda}\frac{e}{4\pi^{2}},$
(99)

where $\mu_{B}$ is the Bohr magneton.

For polarizabilities, we can follow the conventions used for nucleon
polarizabilities, see, e.g., Ref. Hagelstein:2015egb , and define for non-
relativistic neutrinos

${\cal L}_{\rm
NR}=2\pi\left(\alpha_{E1,i}\vec{E}^{2}+\beta_{M1,i}\vec{B}^{2}\right)\otimes
1_{\nu_{i}}.$ (100)

Here, $\alpha_{E1}$ is the electric, and $\beta_{M1}$ the magnetic scalar
polarizability, with $E_{i}=F_{0i}$ the electric, and
$B_{i}=\epsilon_{ijk}F_{jk}/2$ the magnetic field, and $1_{\nu_{i}}$ the
neutrino number operator for neutrinos of flavor $i$ (the non-relativistic
version of the $\bar{\nu}_{i}\nu_{i}$ operator). At dimension 7 in the EFT
expansion in $1/\Lambda$, Eq. (1), we have

$\alpha_{E1,i}=\beta_{M1,i}=\frac{\alpha}{24\pi^{2}\Lambda^{3}}C_{1,ii}^{(7)}.$
(101)

The relation (101) is broken by $m_{\nu}^{2}/\Lambda^{2}$ suppressed
contributions from higher-order operators, for instance from
$(\partial_{\mu}\bar{\nu}_{i}P_{L}\partial_{\tau}\nu_{j})F^{\mu\nu}F^{\tau}_{\nu}$.
These can become important only if the effective scale $\Lambda$ is low, i.e.,
if there are light mediators with mass $m_{\phi}$ comparable to the neutrino
mass that get integrated out in the construction of (101). In our numerical
examples, however, we always have $m_{\phi}\gg m_{\nu}$.

##  Appendix B Further details on stellar cooling rate calculations

In this appendix we collect the $\phi$ production rates relevant for the
stellar cooling bounds discussed in Section 4.

###  B.1 Primakoff conversion

The rate for the Primakoff conversion of a photon $\gamma$ (more precisely,
the transverse plasmon, $\gamma_{T}$) to a pseudoscalar $\phi$ in the field of
the nucleus in a plasma, is given by Lucente_2020 ; Carenza:2020zil ,

$\begin{split}\Gamma_{\gamma\to\phi}&=\left(\frac{c_{\gamma}^{\prime}\alpha}{8\pi
f_{\phi}}\right)^{2}\frac{T\kappa^{2}}{2\pi}\frac{p}{E_{\phi}}\Big{[}\frac{\left((k+p)^{2}+\kappa^{2}\right)\left((k-p)^{2}+\kappa^{2}\right)}{4kp\kappa^{2}}\ln\left(\frac{(k+p)^{2}+\kappa^{2}}{(k-p)^{2}+\kappa^{2}}\right)-\\\
&-\left(\frac{(k^{2}-p^{2})^{2}}{4kp\kappa^{2}}\right)\ln\left(\frac{(k+p)^{2}}{(k-p)^{2}}\right)-1\Big{]}\,,\end{split}$
(102)

where $k$ and $p$ are the incoming photon and outgoing $\phi$ momenta,
respectively, while $\kappa$ is the Debye screening length,

$\kappa^{2}=\frac{4\pi\alpha}{T}\biggr{(}n_{e}^{\rm
eff}+\sum_{j}Z_{j}^{2}n_{j}^{\rm eff}\biggr{)}\,.$ (103)

Here $T$ is the temperature of the star at the radius where the Primakoff
conversion occurs, while $n_{e}^{\rm eff}$ and $n_{j}^{\rm eff}$ are the
effective number densities of electrons and ions, the latter with charge
$Z_{j}e$.

The core of a HB star is a non-relativistic, non-degenerate gas of electrons
and helium ions, thus $n_{X}^{\rm eff}=n_{X}$ for $X=e,\text{He}$. The
electrons forms an ideal Fermi gas, so that the number density is given by
$n_{e}=p_{F}^{3}/(3\pi^{2})$, where the Fermi momentum is $p_{F}=88$ keV
Raffelt:1996wa . The number density of Helium ions is $n_{\rm
He}=\rho/(m_{u}A_{He})$, where $\rho_{HB}\sim 10^{4}{\rm g/cm^{3}}$ is the HB
star density, $m_{u}=0.932$ GeV the atomic mass unit, and $A_{\rm He}=4$ the
atomic number of Helium. Combining the two terms in Eq. (103), we get
$\kappa_{\rm HB}\simeq 27$ keV.

The electrons in the SN core form instead a highly degenerate relativistic gas
and, as such, do not contribute to the screening. The plasma is composed of a
degenerate gas of protons; the degeneracy reduces the effective number of
proton targets in the Primakoff process, thus, more care is needed in
calculating $n_{p}^{\rm eff}$. The number density of degenerate non-
relativistic protons at a distance $r$ from the SN core center is given by

$n_{p}(r)=2\int\frac{{\rm
d}^{3}p}{(2\pi)^{3}}\exp\left[\left(\frac{p^{2}}{2m_{p}^{\rm
eff}(r)T(r)}-\eta_{p}(r)\right)+1\right]^{-1}\,,$ (104)

where $m_{p}^{\rm eff}$ is the effective proton mass and $\eta_{p}$ is the
proton degeneracy parameter Payez:2014xsa : protons are (non) degenerate for
$\eta_{p}>1~{}(\eta_{p}<1)$. Neglecting the proton recoil, the effect of
degeneracy on the number of targets can be calculated as Payez:2014xsa  (see
also Eq. (D.26) in Raffelt:1996wa )

$\frac{n_{p}^{\rm eff}}{n_{p}}=\frac{2}{n_{p}}\int\frac{{\rm
d}^{3}p}{(2\pi)^{3}}\hat{f}_{p}\left(1-\hat{f}_{p}\right)\,,$ (105)

where $\hat{f}_{p}$ is the Fermi-Dirac distribution for protons.

The numerical inputs are functions of the SN core profile, i.e., how the
temperature and density change with the radial distance $r$ from the center.
For numerical results in Sec. 4.2 we use the SN core profile from
Fischer:2018kdt  at the benchmark time $t=1\,$s after the start of the
explosion. For rough numerical estimates we can take $r\sim 10\,$km, $T\sim
30\,$MeV, $\eta_{p}\sim 1$, $m_{p}^{\rm eff}\sim 800\,$MeV, which gives
$n_{p}^{\rm eff}\sim 0.6~{}n_{p}$ and $\kappa_{\rm
SN}\sim{\mathcal{O}}(40\,\text{MeV})$.

Inside dense stellar cores the dispersion relation for the transverse plasmon
of energy $\omega$ is well approximated by introducing an effective photon
thermal mass, $\omega_{P}$, Eq. (34), so that
$k=\big{(}\omega^{2}-\omega_{P}^{2}\big{)}^{1/2}$. The rate for the inverse
Primakoff conversion, $\phi\to\gamma$, is then given by

$\Gamma_{\phi\to\gamma}=\frac{2\beta_{\gamma}}{\beta_{\phi}}\Gamma_{\gamma\to\phi}\,,$
(106)

where $\beta_{i}$ is the velocity of particle $i$, that is,
$\beta_{\gamma}=(1-\omega_{P}^{2}/\omega^{2})^{1/2}$ for plasmon of energy
$\omega$ and $\beta_{\phi}=(1-m_{\phi}^{2}/E_{\phi}^{2})^{1/2}$ for $\phi$
with energy $E_{\phi}$. The factor of two in (106) is due to the two possible
polarizations of a transverse plasmon.

###  B.2 Photon and neutrino coalescence

The rates for photon coalescence, $\gamma\gamma\to\phi$, and for neutrino
coalescence, $\nu\nu\to\phi$, can be calculated by first considering a generic
$\phi$ production process $a_{1}+\dots+a_{n}\to\phi$, where
$\\{a_{1},\dots,a_{n}\\}$ is a set of initial states, e.g., the two photons in
photon coalescence. The phase space production rate for the (pseudo)scalars
$\phi$ can be extracted from the corresponding Boltzmann equation for the
phase space density $\hat{f}_{\phi}$. In the limit $\hat{f}_{\phi}\ll 1$ this
gives, see, e.g., Ref. Carenza:2020zil ,

$\frac{\partial\hat{f}_{\phi}}{\partial
t}=\frac{1}{2E_{\phi}}\int\prod_{i\to\\{a\\}}\frac{{d}^{3}k_{i}}{(2\pi)^{3}2E_{i}}(2\pi)^{4}\delta^{(4)}\left(p-k_{i}\right)|\overline{\cal
M}|^{2}\prod_{i\to\\{a\\}}\hat{f}_{i}(E_{i})\,,$ (107)

where $p$ is the $\phi$ four-momentum and $k_{i}$ the momenta of the initial
states, $\overline{\cal M}$ is the spin averaged matrix element for the
process, while the integration is performed over phase space of the initial
states. In writing (107) we assumed that, once produced, the $\phi$ escapes
the SN, and thus we can take $\hat{f}_{\phi}\simeq 0$ in the possible
collision terms on the right and ignore them. The number of $\phi$ produced is
then given by

${d}N_{\phi}=\hat{f}_{\phi}\frac{{d}^{3}p}{(2\pi)^{3}}\quad\Rightarrow\quad\frac{{d}^{2}N_{\phi}}{{d}t\,{d}E_{\phi}}=\frac{\partial\hat{f}_{\phi}}{\partial
t}|\vec{p}\,|E_{\phi}\frac{{d}\Omega}{(2\pi)^{3}}\,.$ (108)

##### Photon coalescence.

In the production of on-shell $\phi$ via annihilation of two photons,
$\gamma(k_{1})+\gamma(k_{2})\to\phi(p)$, the two photon momenta need to
satisfy $(k_{1}+k_{2})^{2}=m_{\phi}^{2}$. The emissivity for the case when
$\phi$ is a pseudoscalar, i.e., an ALP, is well known in the literature, see
for example Ref. Lucente_2020 . If $\phi$ is a scalar the amplitude squared
changes to

$|\overline{\cal
M}|^{2}=\left(\frac{\alpha}{3\sqrt{2}\pi}\frac{c_{\gamma}}{f_{\phi}}\right)^{2}m_{\phi}^{4}\biggr{[}\biggr{(}1-\frac{4\omega_{P}^{2}}{m_{\phi}^{2}}\biggr{)}^{2}+\frac{6\omega_{P}^{4}}{m_{\phi}^{4}}\biggr{]}\,.$
(109)

The last term in the square brackets does not appear in the pseudoscalar case;
the size of it is, however, relevant only in proximity of the kinematical
threshold, $m_{\phi}=2\omega_{P}$, and becomes quickly negligible for heavier
$m_{\phi}$ masses. Performing the integration over the initial momenta in
(107) gives

$\frac{\partial\hat{f}_{\phi}}{\partial
t}=\left(\frac{\alpha}{3\sqrt{2}\pi}\frac{c_{\gamma}}{f_{\phi}}\right)^{2}\frac{m_{\phi}^{4}}{32\pi
E_{\phi}}\sqrt{1-\frac{4\omega_{P}^{2}}{m_{\phi}^{2}}}\biggr{[}\biggr{(}1-\frac{4\omega_{P}^{2}}{m_{\phi}^{2}}\biggr{)}^{2}+\frac{6\omega_{P}^{4}}{m_{\phi}^{4}}\biggr{]}\exp\left(-\frac{E_{\phi}}{T}\right)\,,$
(110)

where we assumed that the initial photons follow the Maxwell-Boltzmann
distribution, where, from conservation of energy, the sum of the two photon
energies satisfies $E_{1}+E_{2}=E_{\phi}$. The emissivity due to the photon
coalescence is then given by

$\begin{split}Q_{\phi,\gamma\gamma}&=\left(\frac{\alpha}{3\sqrt{2}\pi}\frac{c_{\gamma}}{f_{\phi}}\right)^{2}\frac{1}{32\pi^{3}}\int_{m_{\phi}}^{\infty}{\rm
d}E_{\phi}E_{\phi}|\vec{p}|m_{\phi}^{4}\\\
&\times\sqrt{1-\frac{4\omega_{P}^{2}}{m_{\phi}^{2}}}\biggr{[}\biggr{(}1-\frac{4\omega_{P}^{2}}{m_{\phi}^{2}}\biggr{)}^{2}+\frac{6\omega_{P}^{4}}{m_{\phi}^{4}}\biggr{]}\exp\left(-\frac{E_{\phi}}{T}\right)\,,\end{split}$
(111)

where $|\vec{p}|=\sqrt{E_{\phi}^{2}-m_{\phi}^{2}}$. The result for
pseudoscalar $\phi$ is obtained by neglecting the second term in the square
bracket and by replacing $c_{\gamma}\to 3/2~{}c_{\gamma}^{\prime}$.

##### Neutrino coalescence.

The amplitude squared for the $\phi$ production via scattering of two
neutrinos, $\nu(k_{1})+\nu(k_{2})\to\phi(p)$, is given by

$|\bar{\cal M}|^{2}=\frac{c_{\nu}^{2}}{2}\left(k_{1}\cdot
k_{2}\right)=\frac{c_{\nu}^{2}m_{\phi}^{2}}{4}\,,$ (112)

where we used the momentum conservation and assumed massless neutrinos. The
neutrinos follow the Fermi-Dirac thermal distribution,
$\hat{f}_{\nu}(E)=(\exp[(E-\mu)/T]+1)$, where $\mu$ is the neutrino chemical
potential.

The Boltzmann equation can be written as

$\begin{split}\frac{\partial\hat{f}_{\phi}}{\partial
t}&=\frac{1}{2E_{\phi}}\int\frac{{\rm d}^{3}k_{1}}{2E_{1}(2\pi)^{3}}\frac{{\rm
d}^{3}k_{2}}{2E_{2}(2\pi)^{3}}(2\pi)^{4}\delta^{(4)}\left(p-k_{1}-k_{2}\right)\frac{c_{\nu}^{2}m_{\phi}^{2}}{4}\hat{f}_{\nu}(E_{1})\hat{f}_{\nu}(E_{2})\\\
&=\frac{c_{\nu}^{2}m_{\phi}^{2}}{16E_{\phi}(2\pi)^{2}}\int\frac{{\rm
d}^{3}k_{1}}{E_{1}}\delta\left((p-k_{1})^{2}\right)\hat{f}_{\nu}(E_{1})\hat{f}_{\nu}(E_{2})\,,\end{split}$
(113)

where in the second line we integrated over ${\rm d}^{4}k_{2}$ using the on-
shell condition. The delta funtion can be written as

$\delta\left((p-k_{1})^{2}\right)=\delta\left(m_{\phi}^{2}-2p\cdot
k_{1}\right)=\frac{1}{2|\vec{p}|E_{1}}\delta\biggr{(}\cos\theta-\frac{(m_{\phi}^{2}-2E_{\phi}E_{1})}{2|\vec{p}|E_{1}}\biggr{)}\,,$
(114)

where $\theta$ is the angle between $\vec{p}$ and $\vec{k}_{1}$. Requiring
that $-1<\cos\theta<1$ gives the integration limits

$E_{1,{\rm min}}=\frac{E_{\phi}-\sqrt{E_{\phi}^{2}-m_{\phi}^{2}}}{2}\,,\qquad
E_{1,{\rm max}}=\frac{E_{\phi}+\sqrt{E_{\phi}^{2}-m_{\phi}^{2}}}{2}\,.$ (115)

Changing to spherical coordinates for the integration over ${\rm d}^{3}k_{1}$
then gives,

$\frac{\partial\hat{f}_{\phi}}{\partial
t}=\frac{c_{\nu}^{2}m_{\phi}^{2}}{64\pi E_{\phi}}\int_{E_{1,{\rm
min}}}^{E_{1,{\rm max}}}\frac{{\rm
d}E_{1}}{|\vec{p}|}\hat{f}_{\nu}(E_{1})\hat{f}_{\nu}(E_{\phi}-E_{1})\,.$ (116)

Using the above expression in $\phi$ production rate, Eq. (108), and then in
the general expression for emissivity, Eq. (35), gives the emissivity due to
the neutrino coalescence,

$Q_{\phi,\nu\nu}=\frac{c_{\nu}^{2}m_{\phi}^{2}}{16(2\pi)^{3}}\int_{m_{\phi}}^{\infty}{\rm
d}E_{\phi}\int_{E_{1,min}}^{E_{1,max}}{\rm
d}E_{1}E_{\phi}\hat{f}_{\nu}(E_{1})\hat{f}_{\nu}(E_{\phi}-E_{1})\,.$ (117)

##  Appendix C Rare decays of heavy (pseudo)scalars

In this appendix we provide further details on the derivation of the
constraints on the $\nu\nu\gamma\gamma$ effective interactions via
$S\to\gamma\gamma\to\nu\nu$ decays that were given in Section 5.4. Throughout,
we assume that we can use EFT to describe the neutrino-photon interactions,
Eq. (1).

To construct the bound, one needs first to compute the discontinuity of the
$S\to\nu\bar{\nu}$ amplitude. While many possible intermediate states can
contribute to $S\to\nu\bar{\nu}$, only on-shell intermediate states generate
the absorptive part. The total amplitude can then be obtained from a
dispersion relation Donoghue:1996fv . Selecting one intermediate state,
$\gamma\gamma$, one obtains a bound on the parameters of the
$\nu\nu\gamma\gamma$ interaction from the bounds on the invisible widths of
the $\pi$ and $B$-mesons and the Higgs bosons (see also Bhattacharya:2018msv
).

To obtain the bound, one needs to parameterize the $S\to\gamma\gamma$ decay
amplitude, which can be written as Bosch:2002bv

${\cal A}(S\to\gamma\gamma)\equiv{\cal
A}_{+}\langle\gamma(q_{1},\epsilon_{1})\gamma(q_{2},\epsilon_{2})|F^{\mu\nu}F_{\mu\nu}|0\rangle+{\cal
A}_{-}\langle\gamma(q_{1},\epsilon_{1})\gamma(q_{2},\epsilon_{2})|\tilde{F}^{\mu\nu}F_{\mu\nu}|0\rangle\,,$
(118)

where $q_{1},q_{2}$ are the photon momenta and $\epsilon_{1},\epsilon_{2}$
their polarizations. The matrix elements of the operators
$F^{\mu\nu}F_{\mu\nu}$ and $\tilde{F}^{\mu\nu}F_{\mu\nu}$ in Eq. (118) are
given by

$\begin{split}\langle\gamma(q_{1},\epsilon_{1})\gamma(q_{2},\epsilon_{2})|F^{\mu\nu}F_{\mu\nu}|0\rangle&=-4\big{(}q_{1}\cdot
q_{2}g^{\alpha\beta}-q_{1}^{\beta}q_{2}^{\alpha}\big{)}\epsilon_{1\alpha}\epsilon_{2\beta}\,,\\\
\langle\gamma(q_{1},\epsilon_{1})\gamma(q_{2},\epsilon_{2})|\tilde{F}^{\mu\nu}F_{\mu\nu}|0\rangle&=4\epsilon^{\alpha\beta\rho\sigma}\epsilon_{1\alpha}\epsilon_{2\beta}q_{1\rho}q_{2\sigma}\,.\end{split}$
(119)

Note that the CP-even (CP-odd) matrix element is symmetric (anti-symmetric) in
$\epsilon_{1}\leftrightarrow\epsilon_{2}$ interchange.

The discontinuity of the $S\to\gamma\gamma\to\bar{\nu}\nu$ amplitude reads

$\mbox{Disc }{\cal
A}(S\to\gamma\gamma\to\bar{\nu}\nu)=i\frac{m_{S}^{4}}{2\pi}\left[{\cal
A}_{+}\left({\cal A}^{R}_{+}\right)^{*}+{\cal A}_{-}\left({\cal
A}^{R}_{-}\right)^{*}\right],$ (120)

where the contributions from the Rayleigh operators are encoded in (cf. Eq.
(5))

${\cal A}^{R}_{+}=\frac{1}{2}\left(\frac{\alpha}{12\pi}\right)\hat{C}_{1}^{\rm
Re},\qquad{\cal
A}^{R}_{-}=\frac{1}{2}\left(\frac{\alpha}{8\pi}\right)\hat{C}_{2}^{\rm Re},$
(121)

Notice that while ${\cal A}^{R}_{\pm}$ are real since they arise from an
effective Lagrangian, ${\cal A}^{S}_{\pm}$ are not necessarily real, as they
might receive contributions from other on-shell intermediate states. In fact,
please note that the discontinuity computed in Eq. (120) assumes that the on-
shell transitions $S\to\gamma\gamma$ and $\gamma\gamma\to\nu\bar{\nu}$ are
dominated by local interactions, i.e. there are no discontinuities generated
by on-shell contributions in $A_{\pm}^{(R)}$. This implies that one cannot
obtain a meaningful bound from the invisible decays of the $D^{0}$, as
$D^{0}\to\gamma\gamma$ is dominated by the non-local contributions
Burdman:2001tf .

Since $|\mbox{Im }{\cal A}(S\to\gamma\gamma\to\bar{\nu}\nu)|\leq|{\cal
A}(S\to\gamma\gamma\to\bar{\nu}\nu)|$, the decay rate calculated using
$|\mbox{Im }{\cal A}|$ is smaller or equal to the decay rate calculated using
$|{\cal A}|$ which in turn is smaller than $|{\cal A}(S\to{\rm invisible})|$.
The rate calculated using $|\mbox{Disc }{\cal A}|$ is

$\begin{split}\Gamma_{{\rm Im}}&=\frac{|\mbox{Im }{\cal A}|^{2}}{16\pi
m_{S}}=\frac{m_{S}^{9}}{256\pi^{3}}\frac{1}{4}\left(\frac{\alpha}{8\pi}\right)^{2}\left|\frac{2}{3}{\cal
A}_{+}\hat{C}_{1}^{\text{Re}*}+{\cal
A}_{-}\hat{C}_{2}^{\text{Re}*}\right|^{2}\leq\Gamma(S\to{\rm
inv}).\end{split}$ (122)

Note that summing over the neutrino spins gives a factor of $m_{S}^{2}$ in the
rate. If only $\hat{C}_{1}^{\text{Re}}$ or $\hat{C}_{2}^{\text{Re}}$ are
nonzero, the above relation can be rewritten in the form of the bound in (57),
that was used to obtain bounds on Rayleigh operators from the bounds on
invisible $\pi^{0}$ and Higgs decays.

For the $B^{0}$ meson decays both ${\cal A}_{+}$ and ${\cal A}_{-}$ are
nonzero. At leading order in the $1/m_{b}$ expansion Bosch:2002bv  they are
given by

${\cal A}_{+}={\cal
A}_{-}=\frac{G_{F}}{\sqrt{2}}\frac{\alpha}{3\pi}\frac{f_{B}}{2}V_{td}^{*}V_{tb}C_{7\gamma}\frac{m_{B^{0}}}{\lambda_{B}}\,,$
(123)

where $f_{B}=0.190$ GeV FlavourLatticeAveragingGroup:2019iem  is the $B$-meson
decay constant, $\lambda_{B}=0.35$ GeV Bosch:2002bv  is the inverse of the
first inverse moment of the $B$-meson light-cone distribution amplitude, while
$C_{7\gamma}=-0.38$ Benzke:2010tq  is the Wilson coefficient in the effective
weak Hamiltonian. Higher order corrections to the expansion are $\sim 10\%$,
and we neglect them in the following. Using (122) together with the total
decay width, $\Gamma_{B^{0}}=4.3\times 10^{-13}$ GeV and
$\text{Br}(B^{0}\to\text{ inv})\leq 2.4\cdot 10^{-5}$ gives (59).

## References

  * (1) Super-Kamiokande collaboration, Y. Fukuda et al., _Evidence for oscillation of atmospheric neutrinos_ , _Phys. Rev. Lett._ 81 (1998) 1562–1567, [hep-ex/9807003].
  * (2) M. C. Gonzalez-Garcia and Y. Nir, _Neutrino Masses and Mixing: Evidence and Implications_ , _Rev. Mod. Phys._ 75 (2003) 345–402, [hep-ph/0202058].
  * (3) M. C. Gonzalez-Garcia and M. Maltoni, _Phenomenology with Massive Neutrinos_ , _Phys. Rept._ 460 (2008) 1–129, [0704.1800].
  * (4) S. Dell’Oro, S. Marcocci, M. Viel and F. Vissani, _Neutrinoless double beta decay: 2015 review_ , _Adv. High Energy Phys._ 2016 (2016) 2162659, [1601.07512].
  * (5) M. B. Voloshin, _On Compatibility of Small Mass with Large Magnetic Moment of Neutrino_ , _Sov. J. Nucl. Phys._ 48 (1988) 512.
  * (6) K. S. Babu and R. N. Mohapatra, _Model for Large Transition Magnetic Moment of the $\nu_{e}$_, _Phys. Rev. Lett._ 63 (1989) 228.
  * (7) K. S. Babu and R. N. Mohapatra, _Large transition magnetic moment of the neutrino from horizontal symmetry_ , _Phys. Rev. D_ 42 (1990) 3778–3793.