aainstitutetext: Department of Physics, Oklahoma State University, Stillwater,
OK 74078, USAbbinstitutetext: Department of Physics and McDonnell Center for
the Space Sciences, Washington University,
St. Louis, MO 63130, USA

# Scalar Non-Standard Interactions of Neutrinos in Earth, Sun, Supernovae and
Early Universe

K.S. Babu b,a Garv Chauhan b and P. S. Bhupal Dev<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Non-standard interactions (NSI) of neutrinos with matter mediated by a scalar
field would induce medium-dependent neutrino masses which can modify
oscillation probabilities. Generating observable effects requires an ultra-
light scalar mediator. We derive a general expression for the scalar NSI using
techniques of quantum field theory at finite density and temperature and
discuss various limiting cases applicable to the neutrino propagation in
different media, such as the Earth, Sun, supernovae and early universe. We
also analyze various terrestrial and space-based experimental constraints, as
well as astrophysical and cosmological constraints on these NSI parameters,
applicable to either Dirac or Majorana neutrinos. By combining all these
constraints, we show that observable scalar NSI effects, although precluded in
terrestrial experiments, are still possible in future solar and supernovae
neutrino data, and in cosmological observations such as cosmic microwave
background and big bang nucleosynthesis data.

††preprint: OSU-HEP-19-11

##  1 Introduction

The discovery of neutrino oscillations implies that at least two of the three
neutrinos must have small but non-zero masses Tanabashi:2018oca . The global
neutrino oscillation program is now entering a new era, where the known mixing
angles and mass-squared differences are being measured with an ever-increasing
accuracy. Next-generation of long-baseline oscillation experiments like DUNE
are poised to resolve the sub-dominant effects in oscillation data sensitive
to the currently unknown oscillation parameters, namely the Dirac CP phase,
sign of the atmospheric mass-squared difference and the octant of the
atmospheric mixing angle. These analyses are usually performed within the
$3\times 3$ neutrino mixing scheme under the assumption that neutrinos
interact with matter only through the weak interactions mediated by Standard
Model (SM) $W$ and $Z$ bosons. On the other hand, the origin of neutrino mass
clearly requires some new physics beyond the SM, which often comes with
additional non-standard interactions (NSI) of neutrinos with matter fermions
(i.e. electrons and/or nucleons). Allowing for these NSI in neutrino
production, propagation and/or detection can in principle change the whole
picture and crucially affect the interpretation of the experimental data in
terms of the relevant $3\times 3$ oscillation parameters. It is, therefore, of
paramount importance to understand all possible kinds of NSI effects, and to
see how large these effects could be, while being consistent with other
theoretical and experimental constraints. The study of NSI also opens up the
possibility of using neutrino oscillations to probe the origin of neutrino
mass.

Following the SM interactions of neutrinos with matter via either charged-
current (CC) or neutral-current (NC), which after Fierz transformation can be
written in the form
$(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\alpha})(\bar{f}\gamma_{\mu}Pf)$
(with $f,f^{\prime}\in\\{e,u,d\\}$ the matter fermions and
$P\in\\{P_{L},P_{R}\\}$ the chirality projection operators), NSI induced by
either a vector or charged-scalar mediator can be parametrized in terms of
vector and axial-vector currents Wolfenstein:1977ue :

$\displaystyle{\cal L}_{\rm eff}^{\rm V,NC}\ =\
-2\sqrt{2}G_{F}\sum_{f,P,\alpha,\beta}\varepsilon_{\alpha\beta}^{f,P}(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\nu_{\beta})(\bar{f}\gamma_{\mu}Pf)\,,$
(1) $\displaystyle{\cal L}_{\rm eff}^{\rm V,CC}\ =\
-2\sqrt{2}G_{F}\sum_{f,P,\alpha,\beta}\varepsilon_{\alpha\beta}^{f,P}(\bar{\nu}_{\alpha}\gamma^{\mu}P_{L}\ell_{\beta})(\bar{f}\gamma_{\mu}Pf^{\prime})\,,$
(2)

where $G_{F}$ is Fermi’s constant and the $\varepsilon$ terms quantify the
size of the new interactions. The vector component of NSI given by Eq. (1) or
(2) affect neutrino oscillations during propagation in matter by providing a
new flavor-dependent matter potential. The size of vector NSI is governed by
the parameter $\varepsilon\sim g_{X}^{2}m_{W}^{2}/(g^{2}m_{X}^{2})$, where
$g_{X}$ and $m_{X}$ are respectively the coupling and mass of the mediator
$X$, and $g$ is the $SU(2)_{L}$ gauge coupling. There are two possibilities to
realize experimentally observable vector NSI with $\varepsilon\gtrsim
10^{-2}$: (i) heavy mediator case with $m_{X}\sim{\cal O}(100)$ GeV and
$g_{X}\sim{\cal O}(1)$; and (ii) light mediator case with $m_{X}\ll m_{W}$ and
$g_{X}\ll 1$ such that $g_{X}^{2}/m_{X}^{2}\sim G_{F}$, while evading the low-
energy experimental constraints. For concrete ultraviolet-complete model
realizations, see e.g. Ref. Babu:2019mfe  for the heavy mediator case and
Refs. Farzan:2015doa ; Farzan:2015hkd ; Babu:2017olk  for the light mediator
case. For a recent review on different aspects of vector NSI, see Ref.
Dev:2019anc . For the current global status of the constraints on
$\varepsilon$, see Ref. Coloma:2019mbs .

On the other hand, NSI induced by a neutral scalar mediator is no longer
composed of vector current as in Eq. (1) or (2), but a scalar interaction for
Dirac neutrinos given by Ge:2018uhz

$\displaystyle{\cal L}_{\rm eff}^{\rm S}\ =\
\frac{y_{f}y_{\alpha\beta}}{m_{\phi}^{2}}(\bar{\nu}_{\alpha}\nu_{\beta})(\bar{f}f)\,,$
(3)

where $y_{f}$ and $y_{\alpha\beta}$ are respectively the Yukawa couplings of
the matter fermion and neutrinos to the scalar mediator $\phi$. This cannot be
converted to vector currents, and hence, does not contribute to the matter
potential.111The same is true for tensor NSI of the form
$(\bar{\nu}_{\alpha}\sigma^{\mu\nu}\nu_{\beta})(\bar{f}\sigma_{\mu\nu}f)$.
Instead, it appears as a medium-dependent correction to the neutrino mass
term, with the correction factor $\Delta m_{\nu,\alpha\beta}$ being inversely
proportional to the square of the mediator mass. As we will explicitly show
below, large enough scalar NSI effect is possible only for a sufficiently
light scalar mediator,222Eq. (3) is equally applicable for both light and
heavy mediator, since we are dealing with coherent forward scattering of
neutrinos with $q^{2}\to 0$. since we need $G_{\rm eff}\equiv
y_{f}y_{\alpha\beta}/m_{\phi}^{2}\sim 10^{10}G_{F}$ to have any observable
effect for neutrino propagating in Earth with $\Delta m_{\nu}\sim{\cal
O}(0.1m_{\nu})$. Nevertheless, this could potentially lead to significantly
different phenomenological consequences in reactor, solar, atmospheric and
accelerator neutrino oscillations, as well as for supernovae and early-
universe neutrino interactions.

In this paper, we derive a general formula for evaluating the scalar NSI of
the neutrinos which is applicable to different environments. We perform a
systematic study of the scalar NSI in presence of a light scalar mediator
$\phi$. We consider both Dirac and Majorana neutrino possibilities. The main
objective of our paper is to provide a general field-theoretic derivation of
the scalar NSI effect at finite temperature and density, which can be applied
to different environments, such as Earth, Sun, supernovae and early universe.
Then we go on to derive various constraints on the couplings in Eq. (3) as a
function of the mediator mass $m_{\phi}$ from fifth force experiments, solar
and supernova neutrino data, stellar cooling constraints from red giants (RG)
and horizontal branch (HB) stars, and big bang nucleosynthesis (BBN). We have
considered scalar interactions with electrons and nucleons separately to show
the differences in the constraints. We find that the fifth force experiments
constrain masses of $\phi$ below $0.1$ eV and couplings up to $10^{-24}$.
RG/HB stars constrain couplings up to $10^{-12}$ for nucleons and up to
$10^{-16}$ for electrons coupling to $\phi$. Bounds from BBN constrain
couplings up to $10^{-9}$ for the light scalar mediators. After taking into
account all these constraints, we conclude that any prospects of observing
scalar NSI in Earth matter have been ruled out, while these effects are still
measurable with future solar neutrino data, supernova neutrino bursts or in
cosmological observations of extra relativistic degrees of freedom.

The rest of the paper is organized as follows: In Sec. 2, we present a general
field-theoretic derivation of scalar NSI and discuss various limiting cases
that are applicable to Earth, Sun, supernovae and early universe. In Sec. 3,
we discuss the long-range force effects of a light scalar. In Sec. 4, we
summarize the current experimental constraints on the Yukawa couplings
relevant for scalar NSI as a function of the mediator mass. In Sec. 5, we
discuss the thermal mass of the mediator. In Sec. 6, we derive a quantum-
mechanical bound on the effective in-medium mediator mass. In Sec. 7, we
present our main results and discussions. In Sec. 8, we present a UV-complete
model for scalar NSI. Our conclusions are given in Sec. 9. In Appendix A, we
give the detailed derivation of various limiting cases for the scalar NSI
discussed in Sec. 2. In Appendix B, we provide details of the calculation of
the neutrino self-energy in neutrino background. In Appendix C, we present the
calculation for thermal mass of the scalar mediator.

##  2 Field theoretic origin of scalar NSI

In this section we derive expressions for medium-dependent neutrino mass and
energy when the neutrinos have scalar NSI with matter fermions in the
propagating medium. The results derived here are equally applicable for Dirac
and Majorana neutrinos. As we will see later, for observable scalar NSI it
will be required that the scalar field is very light, which we assume here.
Consider the interaction of fermions $f$ and Dirac neutrinos $\nu$ with a
light scalar $\phi$, with the relevant interaction terms given by the
Lagrangian:

$\mathcal{L}\ \supset\
-y_{\alpha\beta}\overline{\nu}_{\alpha}\phi\nu_{\beta}-y_{f}\bar{f}\phi
f-m_{\alpha\beta}\overline{\nu}_{\alpha}\nu_{\beta}-\frac{m_{\phi}^{2}}{2}\phi^{2}\,.$
(4)

In the case of Majorana neutrinos, the relevant Lagrangian has the form:

$\mathcal{L}\ \supset\
-\frac{y_{\alpha\beta}}{2}\overline{\nu^{c}_{\alpha}}\phi\nu_{\beta}-y_{f}\bar{f}\phi
f-\frac{m_{\alpha\beta}}{2}\overline{\nu^{c}_{\alpha}}\nu_{\beta}-\frac{m_{\phi}^{2}}{2}\phi^{2}~{}.$
(5)

We shall focus primarily on the Dirac neutrinos, but essentially all of our
results will apply for Majorana neutrinos as well, provided that the
normalization of couplings is as in Eq. (5). We shall comment on differences
when they arise between the two cases.

A neutrino with four-momentum $p^{\mu}$ propagating through matter obeys the
Dirac equation given by:

$\left[\not{p}-\Sigma(p)\right]\psi\ =\ 0~{}.$ (6)

In a general medium, the self energy $\Sigma$ of the neutrino gets modified.
We apply real time formalism of field theory at finite temperature and density
in our derivations, which is manifestly Lorentz covariant Das:1997gg . With
pure scalar interactions of the type given in Eqs. (4) and (5), the neutrino
self-energy takes the general form

$\Sigma(p)\ =\ m-(a\not{p}+b\not{u}+c[\not{p},\not{u}]),$ (7)

where $m$ is the neutrino mass inside the medium, $u^{\mu}$ is the four-
velocity of the medium and $a,\,b,\,c$ are functions of only two Lorentz
scalars, viz., $p^{2}$ and $p.u$. In a Lorentz covariant description of field
theory at finite temperature and density one introduces a medium four-velocity
vector $u^{\mu}$ as in Eq. (7) obeying $u^{2}=1$. In real time formalism of
thermal field theory, the finite temperature and density correction to self-
energy of a fermion can be calculated with the help of finite temperature
Green’s function for a free Dirac field Das:1997gg  (for applications to
neutrino propagation in matter see Refs. Notzold:1987ik ; Pal:1989xs ;
Langacker:1992xk ):

$S_{f}(p)\ =\
(\not{p}+m_{f})\left[\frac{1}{p^{2}-m_{f}^{2}+i\epsilon}+i\Gamma_{f}(p)\right]\\\
$ (8)

where

$\Gamma_{f}(p)\ =\
2\pi\delta(p^{2}-m_{f}^{2})[n_{f}(p)\Theta(p_{0})+n_{\bar{f}}(p)\Theta(-p_{0})]~{}.$
(9)

Here $\Theta$ is the Heaviside step function and $n_{f}$ ($n_{\bar{f}}$) is
the Fermi-Dirac distribution function for the fermion (anti-fermion)
occupation number of the medium given by

$n_{f}(p)\ =\ \frac{1}{e^{(|p.u|-\mu)/T}+1}\,,\qquad n_{\bar{f}}(p)\ =\
\frac{1}{e^{(|p.u|+\mu)/T}+1}~{},$ (10)

where $\mu$ is the chemical potential and $T$ is the temperature. Integrating
the occupation number over all possible momentum states yields the total
number density of the fermions (or anti-fermions) in the medium:

$N_{f(\bar{f})}\ =\ g_{f}\int\frac{d^{3}p}{(2\pi)^{3}}n_{f(\bar{f})}(p)\,.$
(11)

Here $g_{f}$ denotes the number of internal degrees of freedom and is equal to
two for electrons, nucleons and neutrinos for the two different spin states.

(a)

(b)

Figure 1: Neutrino self-energy diagrams: (a) Tadpole with background of $f$
and $\nu$, and (b) Self-energy in a neutrino background.

###  2.1 Neutrino self-energy from tadpole diagram

The one-loop thermal self energy corrections for the neutrinos arising from
Eq. (4) or Eq. (5) are shown in Fig. 1. We first compute the one-loop neutrino
thermal mass correction induced by the tadpole diagram in Fig. 1(a). The
Lorentz-invariant form of $\Sigma$ as given in Eq. (7) can be conveniently
evaluated by going to the rest frame of the medium, where the amplitude takes
a simple form:

$\displaystyle-i\Sigma_{\alpha\beta}$ $\displaystyle\ =\
iy_{\alpha\beta}\>\frac{i}{q^{2}-m_{\phi}^{2}}\int\frac{d^{4}k}{(2\pi)^{4}}{\rm
Tr}\left[iy_{f}\>iS_{f}(k)\right]\,.$ (12)

In Eq. (12), we can set $q^{2}=0$ for the momentum transfer because we are
only interested in the coherent forward scattering of neutrinos in matter for
NSI effect. Only retaining the finite temperature and density part of the
self-energy, we obtain

$\displaystyle\Sigma_{\alpha\beta}$ $\displaystyle\ =\
\frac{y_{\alpha\beta}y_{f}m_{f}}{\pi^{2}m_{\phi}^{2}}\int_{0}^{\infty}dk_{0}\int_{0}^{\infty}dk^{2}\>k\>\delta(k^{2}-k_{0}^{2}+m_{f}^{2})\left[n_{f}(k_{0})+n_{\bar{f}}(-k_{0})\right]\,.$
(13)

Integrating over $k^{2}$ using the delta function yields the final result:

$\displaystyle\Sigma_{\alpha\beta}\ =\
\frac{y_{\alpha\beta}y_{f}m_{f}}{\pi^{2}m_{\phi}^{2}}\int_{m_{f}}^{\infty}dk_{0}\>\sqrt{k_{0}^{2}-m_{f}^{2}}\left[n_{f}(k_{0})+n_{\bar{f}}(k_{0})\right]\
\equiv\ \Delta m_{\nu,{\alpha\beta}}~{}.$ (14)

While deriving Eq. (14), it has been assumed that the background medium
contains both fermions and anti-fermions. Thus, Eq. (14) is the complete
expression for scalar NSI of neutrinos at any finite temperature and density
in a background without neutrinos. We have provided details of evaluating the
integral of Eq. (14) in various useful limits in Appendix A.

Note that the scalar NSI of Eq. (14) appears as a medium-dependent mass of the
neutrino. The relevant integral can be evaluated analytically in the high
temperature as well as low temperature regimes. We find:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}\ =\ $
$\displaystyle\frac{y_{f}y_{\alpha\beta}}{m_{\phi}^{2}}\left(N_{f}+N_{\bar{f}}\right)$
($\mu,T\ll m_{f}$) (15) $\displaystyle\Delta m_{\nu,{\alpha\beta}}\ =\ $
$\displaystyle\frac{y_{f}y_{\alpha\beta}}{m_{\phi}^{2}}\frac{m_{f}}{2}\left(\frac{3}{\pi}\right)^{\frac{2}{3}}\left(N_{f}^{2/3}+N_{\bar{f}}^{2/3}\right)$
($\mu>m_{f}\gg T$) (16) $\displaystyle\Delta m_{\nu,{\alpha\beta}}\ =\ $
$\displaystyle\frac{y_{f}y_{\alpha\beta}m_{f}}{3\>m_{\phi}^{2}}\left(\frac{\pi^{2}}{12\>\zeta(3)}\right)^{\frac{2}{3}}\left(N_{f}^{2/3}+N_{\bar{f}}^{2/3}\right)$
($\mu<m_{f}\ll T$) . (17)

If the medium does not contain either fermions or anti-fermions of a certain
type, the corresponding number density should be set to zero in the final
result. If the background has more than one type of fermion, the various
contributions should be added. Eq. (15) for $\mu,T\ll m_{f}$ is the non-
relativistic limit for the scalar NSI expression and matches the result stated
in Ref. Ge:2018uhz . It is most useful in case of Earth and Sun. The limiting
case Eq. (16) is useful for relativistic medium backgrounds as with electrons
in supernovae. For effect of scalar NSI in the early universe, Eq. (17) is the
most relevant. Detailed application of these results is carried out in Sec. 7.

###  2.2 Neutrino self-energy in a neutrino background

There is another important diagram that might contribute to effect of neutrino
propagation in a medium, as shown in Fig. 1(b). This diagram contributes to
neutrino self-energy only in media with a neutrino or an anti-neutrino
background. This situation is realized in early universe cosmology. Here we
evaluate the contribution of Fig. 1(b) in such a background. Again using the
real-time formalism of thermal field theory, we can write this contribution
for a Dirac neutrino as:

$\Sigma_{\alpha\beta}^{\nu}\ =\
-y_{\beta\gamma}y_{\gamma\alpha}\int\frac{d^{4}k}{(2\pi)^{4}}\left(\not{k}+\frac{\not{p}}{2}+m_{\nu}\right)\left[\frac{\Gamma_{\phi}\left(k-\frac{p}{2}\right)}{\left(k+\frac{p}{2}\right)^{2}-m_{\nu}^{2}}+\frac{\Gamma_{\nu}\left(k+\frac{p}{2}\right)}{\left(k-\frac{p}{2}\right)^{2}-m_{\phi}^{2}}\right]\,,$
(18)

where we have used the finite temperature Green’s function for a free bosonic
field given by:

$S_{b}(p)\ =\
\left[\frac{1}{p^{2}-m_{b}^{2}+i\epsilon}-i\Gamma_{b}(p)\right]\\\ $ (19)

where

$\Gamma_{b}(p)\ =\ 2\pi\delta(p^{2}-m_{b}^{2})n_{b}(p)\Theta(p_{0})\,,$ (20)

with the Bose-Einstein distribution function given by

$n_{b}(p)\ =\ \frac{1}{e^{(|p.u|)/T}-1}~{},$ (21)

noting that the chemical potential of the real scalar field $\phi$ is zero. We
have carried out the evaluation of the self energy integral of Eq. (18) in
Appendix B; here we summarize our main results. The contribution of Eq. (18)
can be written as

$\Sigma_{\alpha\beta}^{\nu}\ =\
-\frac{y_{\alpha\gamma}y_{\gamma\beta}}{8\pi^{2}|\textbf{p}|}J\,,$ (22)

with $J$ identified as the integral of Eq. (18), except for an overall factor,
and can be decomposed as

$J\ =\ a\not{p}+b\not{u}+c+d[\not{p},\not{u}]~{}.$ (23)

By taking traces of the integral in Eq. (18) multiplied by
$(1,\,\not{p},\,\not{u},\,\not{p}\not{u})$ we can solve for the Lorentz
scalars $(a,\,b,\,c,\,d)$. Defining

$J_{p}\ =\ \text{Tr}(J\not{p})\,,\quad J_{u}\ =\
\text{Tr}(J\not{u})\,,\quad\text{and}\quad J_{m}\ =\ \text{Tr}(J)\,,$ (24)

we find

$a\ =\ \frac{J_{u}\>(p.u)-J_{p}}{4[(p.u)^{2}-p^{2}]}\,,\quad b\ =\
\frac{J_{p}\>(p.u)-J_{u}p^{2}}{4[(p.u)^{2}-p^{2}]}\,,\quad c\ =\
\frac{J_{m}}{4}\,,\quad\text{ and}\quad d\ =\ 0~{}.$ (25)

It is clear that the coefficient $c$ contributes to the neutrino mass in the
medium. But this effect is negligible in our case, because there is no
$1/m_{\phi}^{2}$ enhancement.

There is also a matter potential that is caused by the neutrino self-
interactions. To arrive at it we examine the pole in the neutrino propagator:

$\displaystyle iS_{\nu}^{-1}(p)\ =\ i(\not{p}-\Sigma^{\nu})\ =\
i[\not{p}(1-A)-B\not{u}]\,,$ (26)

where

$A_{\alpha\beta}\ =\ -\frac{y_{\beta\gamma}y_{\gamma\alpha}}{8\pi^{2}|{\bf
p}|}a\,,\qquad B_{\alpha\beta}\ =\
-\frac{y_{\beta\gamma}y_{\gamma\alpha}}{8\pi^{2}|{\bf p}|}b~{}.$ (27)

$S_{\nu}$ is then obtained in terms of $A$ and $B$ as

$\displaystyle iS_{\nu}(p)\ =\
i\frac{[(1-A)\not{p}-B\not{u}]}{\\{(1-A)p-Bu\\}^{2}}~{}.$ (28)

We define energy and momentum of the neutrino in the rest frame of the medium
as Pal:1989xs

${\cal E}\ =\ p.u,\qquad{\cal P}\ =\ \sqrt{{\cal E}^{2}-p^{2}}~{}.$ (29)

The pole in the neutrino propagator of Eq. (28) occurs at energy values given
by

${\cal E}\ =\ \frac{B}{1-A}\pm{\cal P}~{}.$ (30)

The energy shift for neutrinos is thus $B/(1-A)$, while the shift in
antineutrino energy is $-B/(1-A)$.

For significant regions of the Yukawa couplings $y_{\alpha\beta}$ and $y_{f}$,
the scalar $\phi$ does not get thermalized. In this case, there is no $\phi$
background and the term proportional to $\Gamma_{\phi}(k-p/2)$ should be set
to zero. We present our results in this case here. The contribution from
$\Sigma^{\nu}_{\alpha\beta}$ can then be written as:

$\Sigma_{\alpha\beta}^{\nu}\ =\
-y_{\alpha\gamma}y_{\gamma\beta}\int_{\frac{-p_{0}}{2}}^{\infty}dk_{0}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{(\not{k}+\frac{\not{p}}{2}+m_{\nu})}{\left(k-\frac{p}{2}\right)^{2}-m_{\phi}^{2}}\delta\left[{\left(k+\frac{p}{2}\right)^{2}-m_{\nu}^{2}}\right]n_{\nu}\left(k_{0}+\frac{p_{0}}{2}\right)~{}.$
(31)

We defer the details of evaluating this integral to Appendix B. Here we
present the results in the high temperature limit, assuming that the chemical
potential is vanishing. This condition is generally true in the early universe
when neutrinos propagate in a background of neutrinos. Furthermore, we set the
neutrino mass to be zero, which is a consistent approximation as the medium-
induced $m_{\nu}$ is proportional to the original $m_{\nu}$. In the absence of
neutrino mass, we can set $p_{0}^{2}-|{\bf p}|^{2}=0$. Under these conditions,
our results are as follows:

$\displaystyle\frac{a}{|{\bf p}|}$ $\displaystyle\ =\ $
$\displaystyle-\frac{\pi^{2}T^{2}}{24}\left[12\zeta^{\prime}(-1)-{\rm
ln}\left(\frac{16\pi|{\bf p}|T}{m_{\phi}^{2}}\right)\right]$
$\displaystyle\frac{b}{|{\bf p}|}$ $\displaystyle\ =\ $
$\displaystyle\frac{\pi^{2}T^{2}}{12|{\bf p}|}~{}.$ (32)

Here $\zeta^{\prime}(-1)=-0.165421$ is the derivative of Riemann zeta function
evaluated at argument equal to $-1$. Using these results along with Eq. (30),
we arrive at the energy shift experienced by the neutrino in a background of
neutrinos:

$\Delta{\cal E}_{+,\alpha\beta}\ =\
-\frac{y_{\beta\gamma}y_{\gamma\alpha}T^{2}}{96|{\bf
p}|\left(1-\frac{y_{\beta\gamma}y_{\gamma\alpha}T^{2}}{192|{\bf
p}|^{2}}\left\\{12\zeta^{\prime}(-1)-{\rm ln}\left(\frac{16\pi|{\bf
p}|T}{m_{\phi}^{2}}\right)\right\\}\right)}~{}.$ (33)

While we do not use these results explicitly in our analyses, these are part
of the neutrino scalar NSI which may find use in early universe cosmology
where there is a thermal background of neutrinos.

##  3 Long-range force effects

A light scalar coupling to fermions can lead to long-range forces. This
applies to charged fermions as well as neutrinos propagating through a medium.
Even when the neutrino propagates outside of the medium, such long-range
forces can affect its propagation. Thus, calculating the neutrino energy using
point interactions with a very light mediator does not provide a complete
picture. In this section, we sketch a heuristic derivation to account for
these long-range force effects. Long range effects in non-relativistic media
has been studied in Ref. Wise:2018rnb ; Smirnov_2019 . Here, we have extended
the analysis for all background media, i.e. both non-relativistic and
relativistic cases. This will be especially useful in relativistic media such
as in supernovae and in early universe.

We use the Euler-Lagrange equations for the Lagrangian in Eq. (4) to obtain
equations of motion for $\nu$ and $\phi$:

$\displaystyle(i\not{\partial}-m_{\alpha\beta}-y_{\alpha\beta}\phi)\nu_{\beta}$
$\displaystyle\ =\ 0$ (34) $\displaystyle(\partial^{2}+m_{\phi}^{2})\>\phi-
y_{\alpha\beta}\overline{\nu}_{\alpha}\nu_{\beta}-y_{f}\bar{f}f$
$\displaystyle\ =\ 0~{}.$ (35)

As can be seen from Eq. (34), the interaction vertex
$y_{\alpha\beta}\overline{\nu}_{\alpha}\phi\nu_{\beta}$ leads to an extra
contribution to neutrino mass:

$\Delta m_{\nu,{\alpha\beta}}\ =\
y_{\alpha\beta}\>\langle\phi\rangle_{\text{medium}}~{}.$ (36)

To calculate the mass correction for a neutrino propagating in a medium, we
will need to calculate the expectation value of the operators at finite
temperature and density, appearing in Eqs. (34) and (35).

For a medium in thermal equilibrium with fermion number density $N_{f}$ and
anti-fermion number density $N_{\bar{f}}$ can be represented as a Fock state
$|\Psi\rangle$. This state contains information about particle and anti-
particle distribution in different momentum states. Since the system is
assumed to be in thermal equilibrium, the fermion and anti-fermion density in
each momentum state does not change in time. Thus, we can set $t=0$ and the
state $|\Psi\rangle$ is normalized, i.e, $\langle\Psi|\Psi\rangle=1$. The
field operators for the fermion and anti-fermion fields can be written as
Peskin:1995ev :

$\displaystyle f(x)$ $\displaystyle\ =\
\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{\sqrt{2E_{\textbf{p}}}}\;\;\sum_{s}\left[a^{s}_{\textbf{p}}\>u^{s}(p)\>e^{-ip.x}+b^{s\dagger}_{\textbf{p}}\>v^{s}(p)\>e^{ip.x}\right]\,,$
(37) $\displaystyle\bar{f}(x)$
$\displaystyle=\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{\sqrt{2E_{\textbf{p}}}}\;\;\sum_{s}\left[b^{s}_{\textbf{p}}\>\>\overline{v}^{s}(p)e^{-ip.x}+a^{s\dagger}_{\textbf{p}}\>\overline{u}^{s}(p)\>e^{ip.x}\right]~{}.$
(38)

We need to calculate the expectation value of the operator $\bar{f}f$. While
trying to interpret these quantities classically, we first need to normal
order the product of the quantum fields:

$\displaystyle\langle:\bar{f}f:\rangle$ $\displaystyle\ =\
\langle\Psi|:\bar{f}f:|\Psi\rangle\ =\
\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\frac{1}{\sqrt{2E_{\textbf{p}_{1}}}}\frac{1}{\sqrt{2E_{\textbf{p}_{2}}}}\;$
$\displaystyle\times\sum_{s,s^{\prime}}\left[\langle
a^{s\dagger}_{\textbf{p}_{1}}a^{s^{\prime}}_{\textbf{p}_{2}}\rangle\>\overline{u}^{s}(p_{1})u^{s^{\prime}}(p_{2})\>e^{-i(\textbf{p}_{1}-\textbf{p}_{2})\cdot\textbf{x}}+\langle
a^{s\dagger}_{\textbf{p}_{1}}b^{s^{\prime}\dagger}_{\textbf{p}_{2}}\rangle\>\overline{u}^{s}(p_{1}){v}^{s}(p_{2})\>e^{-i(\textbf{p}_{1}+\textbf{p}_{2})\cdot\textbf{x}}\right.$
$\displaystyle\left.+\langle
b^{s}_{\textbf{p}_{1}}a^{s^{\prime}}_{\textbf{p}_{2}}\rangle\>\overline{v}^{s}(p_{1})u^{s^{\prime}}(p_{2})\>e^{i(\textbf{p}_{1}+\textbf{p}_{2})\cdot\textbf{x}}+\langle
b^{s\dagger}_{\textbf{p}_{1}}b^{s^{\prime}}_{\textbf{p}_{2}}\rangle\>\overline{v}^{s}(p_{1})v^{s^{\prime}}(p_{2})\>e^{i(\textbf{p}_{1}-\textbf{p}_{2})\cdot\textbf{x}}\right]\,,$
(39)

where we have used $\langle A\rangle=\langle\Psi|A|\Psi\rangle$ for brevity
and the symbol : : signifies normal ordering of the product. In Eq. (3), terms
like $a^{\dagger}b^{\dagger}$ and $a\>b$ vanish, since they cannot be
contracted because they act on different subspaces. It is well known from
quantum field theory at zero temperature that $a^{\dagger}a$ and
$b^{\dagger}b$ are the number operators for fermions and anti-fermions
respectively. This can be generalized to finite temperature and density using
the Fermi-Dirac distribution:

$\displaystyle\langle\Psi|a^{s\dagger}_{\textbf{p}_{1}}a^{s^{\prime}}_{\textbf{p}_{2}}|\Psi\rangle\
$ $\displaystyle=\
n_{f}(\textbf{p}_{1})\>\delta(\textbf{p}_{1}-\textbf{p}_{2})\delta_{s,s^{\prime}}\,,$
(40)
$\displaystyle\langle\Psi|b^{s\dagger}_{\textbf{p}_{1}}b^{s^{\prime}}_{\textbf{p}_{2}}|\Psi\rangle\
$ $\displaystyle=\
n_{\bar{f}}(\textbf{p}_{1})\>\delta(\textbf{p}_{1}-\textbf{p}_{2})\delta_{s,s^{\prime}}\,.$
(41)

Eq. (40) can be understood by integrating it over all momentum states which
yields the total number density $N_{f}$:

$\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\int\frac{d^{3}p_{2}}{(2\pi)^{3}}\langle\Psi|a^{s\dagger}_{\textbf{p}_{1}}a^{s^{\prime}}_{\textbf{p}_{2}}|\Psi\rangle\
=\ N_{f}~{}.$ (42)

Using the normalization of states $\overline{u}^{s}(p)u^{s}(p)=2m_{f}$, we
obtain:

$\langle\bar{f}f\rangle\ \equiv\ \langle\Psi|\bar{f}f|\Psi\rangle\ =\
g_{f}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{m_{f}\>}{E_{\textbf{p}}}\left[n_{f}(\textbf{p})+n_{\bar{f}}(\textbf{p})\right]\,.$
(43)

Converting Eq. (43) into an energy integral, we have:

$\displaystyle\langle\bar{f}f\rangle\ =\
\frac{g_{f}m_{f}}{2\pi}\int_{m_{f}}^{\infty}dk_{0}\>\sqrt{k_{0}^{2}-m_{f}^{2}}\left[n_{f}(k_{0})+n_{\bar{f}}(k_{0})\right]\,.$
(44)

Notice that the integral form of Eq. (44) matches Eq. (14) except for the pre-
factors. This implies that generalizing the limiting cases for Eq. (44) is
straightforward.

Now to calculate $\Delta m_{\nu,{\alpha\beta}}$ in Eq. (36), we need to solve
Eq. (35) for $\phi$. Considering $y_{f}\bar{f}f$ as a source term and
neglecting the second term assuming low neutrino number density, we can write
the solution as:

$\langle\phi\rangle(\textbf{x})\ =\ -y_{f}\int
d^{3}\textbf{x}^{\prime}\frac{\langle\bar{f}f\rangle(\textbf{x}^{\prime})}{4\pi|\textbf{x}-\textbf{x}^{\prime}|}e^{-m_{\phi}(|\textbf{x}-\textbf{x}^{\prime}|)}~{}.$
(45)

Under assumptions of spherical symmetry of the medium, integrating over the
angular variables yields the solution of the form:

$\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{f}y_{\alpha\beta}}{m_{\phi}\>r}\left(e^{-m_{\phi}r}\int_{0}^{r}x\>\langle\bar{f}f\rangle\>\sinh{(m_{\phi}\>x)}\>dx+\sinh{(m_{\phi}\>r)}\int_{r}^{\infty}x\>\langle\bar{f}f\rangle\>e^{-m_{\phi}\>x}\>dx\right)~{}.$
(46)

We will work out Eq. (46) in the relativistic limit for two different density
profile distributions, as discussed below. While we do not use these analytic
results in our numerical analysis, these special cases can give insight. We
use actual density profiles of the Sun and supernovae in our numerical
calculations, integrating the relevant integrals exactly.

###  3.1 Constant density distribution

For a relativistic medium like electron background in supernova, the quantity
$\langle\bar{f}f\rangle$ takes the form:

$\langle\bar{f}f\rangle_{\text{SN}}\ =\
\frac{m_{f}}{2}\left(\frac{3N_{f}}{\pi}\right)^{\frac{2}{3}}\,.$ (47)

Given a constant density distribution such that

$N_{f}(r)\ =\ N_{f}(0)\>\Theta(R-r)\,,$ (48)

where $R$ is the radius of the constant-density spherical body. Plugging the
$\langle\bar{f}f\rangle$ in Eq. (46) yields a general form for scalar NSI in
relativistic media with $\mu>m_{f}\gg T$ :

$\displaystyle\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}}{m_{\phi}\>r}\frac{m_{f}}{2}\left(\frac{3}{\pi}\right)^{\frac{2}{3}}$
$\displaystyle\left(e^{-m_{\phi}r}\int_{0}^{r}x\>N_{f}^{2/3}\>\sinh{(m_{\phi}\>x)}\>dx\right.$
$\displaystyle\left.+\sinh{(m_{\phi}\>r)}\int_{r}^{\infty}x\>N_{f}^{2/3}\>e^{-m_{\phi}\>x}\>dx\right)~{}.$
(49)

For number density profile in consideration, the above equation yields:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}m_{f}}{2m_{\phi}\>r}\left(\frac{3N_{f}(0)}{\pi}\right)^{\frac{2}{3}}\times$
$\displaystyle F_{<}$ ($r\leq R$) , (50) $\displaystyle\Delta
m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}m_{f}}{2m_{\phi}\>r}\left(\frac{3N_{f}(0)}{\pi}\right)^{\frac{2}{3}}\times$
$\displaystyle F_{>}$ ($r>R$) , (51)

where

$\displaystyle F_{<}\ $ $\displaystyle=\
1-\frac{m_{\phi}R+1}{m_{\phi}\>r}e^{-m_{\phi}\>R}\sinh{(m_{\phi}\>r)}\,,$ (52)
$\displaystyle F_{>}\ $ $\displaystyle=\
\frac{e^{-m_{\phi}\>r}}{m_{\phi}\>r}[m_{\phi}\>R\cosh{(m_{\phi}\>R)}-\sinh{(m_{\phi}\>R)}]~{}.$
(53)

Note that the pre-factor in Eq. (3.1) matches the scalar NSI contribution
calculated in Eq. (16) assuming point contact interaction.

For the non-relativistic case our formalism gives the same result derived in
Ref. Smirnov_2019  and given below:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}\>N_{f}(0)}{m_{\phi}^{2}}\times$ $\displaystyle
F_{<}$ ($r\leq R$) (54) $\displaystyle\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}\>N_{f}(0)}{m_{\phi}^{2}}\times$ $\displaystyle
F_{>}$ ($r>R$) (55)

where the functions ($F_{<},\>F_{>}$) are identical to the ones in Eqs. (52)
and (53).

###  3.2 Exponential density distribution

Given a relativistic medium ($\mu>m_{f}\gg T$) with the following number
density profile:

$N_{f}(r)\ =\ N_{f}(0)\>e^{-\lambda\>r}\>\Theta(R-r)$ (56)

where $R$ is the radius of the spherical body in consideration, Eq. (3.1)
yields:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}}{2m_{\phi}\>r}\left(\frac{3N_{f}(0)}{\pi}\right)^{\frac{2}{3}}\times$
$\displaystyle G_{<}$ ($r\leq R$) , (57) $\displaystyle\Delta
m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}}{2m_{\phi}\>r}\left(\frac{3N_{f}(0)}{\pi}\right)^{\frac{2}{3}}\times$
$\displaystyle G_{>}$ ($r>R$) , (58)

where

$\displaystyle G_{<}\ $ $\displaystyle=\ \frac{2\lambda
m_{\phi}}{3}\left(\frac{e^{m_{\phi}r}\left(\frac{3m_{\phi}^{2}r}{2\lambda}-\frac{2\lambda\>r}{3}-2\right)+2e^{\frac{2\lambda\>r}{3}}}{(m_{\phi}^{2}-\frac{4\lambda^{2}}{9})^{2}}\right)e^{-r(\frac{2\lambda}{3}+m_{\phi})}$
$\displaystyle\qquad-\left(\frac{\sinh(m_{\phi}r)(m_{\phi}R+\frac{2\lambda\>R}{3}+1)}{(m_{\phi}+\frac{2\lambda}{3})^{2}}\right)e^{-R(\frac{2\lambda}{3}+m_{\phi})}\,,$
(59) $\displaystyle G_{>}\ $ $\displaystyle=\
\sinh(m_{\phi}R)\left(\frac{m_{\phi}^{2}(\frac{2\lambda\>R}{3}-1)-\frac{4\lambda^{2}}{9}(\frac{2\lambda\>R}{3}+1)}{\left(m_{\phi}^{2}-\frac{4\lambda^{2}}{9}\right)^{2}}\right)e^{-(m_{\phi}r+\frac{2\lambda\>R}{3})}+\frac{4\lambda
m_{\phi}}{3\left(m_{\phi}^{2}-\frac{4\lambda^{2}}{9}\right)^{2}}e^{-m_{\phi}r}$
$\displaystyle\qquad+\cosh(m_{\phi}R)\left(\frac{m_{\phi}^{3}R-\frac{4\lambda^{2}\>Rm_{\phi}}{9}-\frac{4\lambda\>m_{\phi}}{3}}{\left(m_{\phi}^{2}-\frac{4\lambda^{2}}{9}\right)^{2}}\right)e^{-(m_{\phi}r+\frac{2\lambda\>R}{3})}~{}.$
(60)

Similar analyses can be done for other relativistic cases such as for early
universe cosmology ($\mu<m_{f}<T$) albeit with a different pre-factor.

For an exponential density distribution with a cut-off in the non-relativistic
case we obtain:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}\>N_{f}(0)}{m_{\phi}\>r}\times$ $\displaystyle
K_{<}$ ($r\leq R$) , (61) $\displaystyle\Delta m_{\nu,{\alpha\beta}}(r)\ =\
\frac{y_{\alpha\beta}\>y_{f}\>N_{f}(0)}{m_{\phi}\>r}\times$ $\displaystyle
K_{>}$ ($r>R$) , (62)

where we can obtain the functions $K_{<}$ and $K_{>}$ by replacing
$\lambda\rightarrow\frac{3\lambda}{2}$ in $G_{<}$ and $G_{>}$ respectively,
i.e., $K(\lambda)_{>(<)}=G(3\lambda/2)_{>(<)}$. This expression is in full
agreement with the result of Ref. Smirnov_2019 .

##  4 Experimental constraints on couplings

In this section we explore two specific scenarios:

  1. (i)

scalar $\phi$ coupling only to electrons and neutrinos.

  2. (ii)

scalar $\phi$ coupling only to nucleons and neutrinos.

Here neutrinos can be either Dirac or Majorana in nature. In this section, we
discuss experimental constraints on the couplings and mass of $\phi$ in the
aforementioned scenarios.

In accordance to Eq. (4), the scalar coupling to electron is denoted by
$y_{e}$. On the other hand, the scalar coupling to quark cannot be probed
directly but only measurable through their effect with scalar-nucleon
interaction. Thus, we present the experimental constraints on scalar-nucleon
coupling labeled as $y_{N}$. The conversion from quark level couplings $y_{q}$
to $y_{N}$ is discussed later in Sec. 4.1.4.

###  4.1 Constraints on $y_{e}$ and $y_{N}$

####  4.1.1 Anomalous electron magnetic moment

A scalar coupling with the electrons will contribute to the electron anomalous
magnetic moment $(g-2)_{e}$ given by Liu:2018xkx :

$\Delta a_{e}\ =\
\frac{1}{8\pi^{2}}\int_{0}^{1}dx\frac{(1-x)^{2}(1+x)y_{e}^{2}}{(1-x)^{2}+x(m_{\phi}/m_{e})^{2}}~{}.$
(63)

There is currently a $2.4\sigma$ discrepancy between the experimentally-
inferred value and SM prediction for $\Delta a_{e}=a_{e}^{\rm exp}-a_{e}^{\rm
SM}=(-88\pm 36)\times 10^{-14}$ Tanabashi:2018oca . A light scalar can
potentially make this discrepancy worse, as it gives a positive contribution,
and thus provides a useful limit on scalar NSI parameters. Using the $3\sigma$
value for the $\Delta a_{e}$, the allowed region in the $y_{e}-m_{\phi}$ plane
is obtained. This constraint is shown in Figs. 4 and 6, labeled as
$(g-2)_{e}$. This constraint yields an almost constant upper bound of
$y_{e}<3.4\times 10^{-6}$ for light scalar mediators.

####  4.1.2 Fifth force experiments

These experiments measure the presence of fifth forces as deviation from the
Newtonian gravitational potential between a given source mass and a test mass,
which is parametrized as follows:

$V(r)\ =\ -\frac{Gm_{1}m_{2}}{r}\left(1+\alpha e^{-r/\lambda}\right)~{}.$ (64)

Given an interaction vertex of the form $y_{f}\bar{f}\phi f$ as in Eqs. (4)
and (5), consider the scattering of two distinguishable fermions in the non-
relativistic limit. The corresponding Yukawa potential for the interaction is
given by (see Sec. 4.7 of Ref. Peskin:1995ev ):

$V(r)\ =\ -\frac{y_{f}^{2}}{4\pi r}\>e^{-m_{\phi}r}\,,$ (65)

where $r$ is the distance between the scattering particles.

For experiments detailed in Refs. Irvine85 ; Colorado02 ; EotWash07 ;
Stanford08 ; HUST12 ; HUST+16  in the range $\lambda=10^{-6}\text{ to }10^{2}$
m, the constraints provided on $\alpha$ in Eq. (64) are not directly
applicable to $y_{f}$ in Eq. (65). In our case, assuming a particle (e.g.
lepton, quark) couples to the scalar mediator with strength $q$ and each
interacting body contains $N$ number of these particles, the potential between
two extended bodies can be written as:

$V_{\phi}(r)\ =\ -\frac{N_{1}q_{1}\>N_{2}q_{2}}{4\pi r}\>e^{-m_{\phi}r}~{}.$
(66)

We identify the inverse of the length scale $\lambda$ as mass of the scalar
particle $\phi$. Thus, we have:

$\alpha\ =\ \frac{N_{1}q_{1}\>N_{2}q_{2}}{4\pi Gm_{1}m_{2}}\ =\
\frac{q_{1}\>q_{2}}{4\pi GA_{1}A_{2}u^{2}}\ =\ \frac{1}{4\pi
Gu^{2}}\frac{q_{1}}{A_{1}}\frac{q_{2}}{A_{2}}\,,$ (67)

where we have used the relation $m=NAu$ ($A$= mass number, $u$ = 1 atomic mass
unit) and $G$ is the gravitational constant. For bounds on $y_{e}$, the
coupling strength will be proportional to the lepton number ($L$) or atomic
number ($Z$), i.e., $q=Zy_{e}$, leading to

$\alpha\ =\ \frac{y_{e}^{2}}{4\pi
Gu^{2}}\frac{Z_{1}}{A_{1}}\frac{Z_{2}}{A_{2}}~{}.$ (68)

Values for charge to mass number ratio for test and source masses can be
obtained from the experimental setups as given in Table 1. These are shown in
Figs. 4 and 6 by the labels I to VII. Similar results follow for coupling to
the nucleons $y_{N}$ by replacing the atomic number by mass number. This
implies that constraints on $y_{N}$ will be independent of the material used
in the experiment, as shown in Figs. 4 and 6.

Label | References | Source Mass Composition | Test Mass Composition | $\frac{Z_{1}}{A_{1}}\frac{Z_{2}}{A_{2}}$
---|---|---|---|---
I | Stanford Stanford08  | Gold, Silicon | Gold | 0.1804
II | Colorado Colorado02  | Tungsten | Tungsten | 0.1621
III | Eot-Wash’07 EotWash07  | Molybdenum, Tantalum | Molybdenum | 0.1839
IV | HUST’12 HUST12  | Tungsten | Tungsten | 0.1621
V | HUST+ ’16 HUST+16  | Tungsten | Tungsten | 0.1621
VI | Irvine A Irvine85  | Copper | Copper | 0.2159
VII | Irvine B Irvine85  | Stainless Steel | Copper | 0.2116
Table 1: The compositions of source and test masses used in the experiment
listed and the corresponding values of ratio
$\frac{Z_{1}}{A_{1}}\frac{Z_{2}}{A_{2}}$.

Additional constraints on $y_{e}$ and $y_{N}$ can be directly obtained from
Ref. ADELBERGER2009102  which used experiments in the range
$\lambda=10^{-1}\text{ to }10^{13}$ m and the corresponding limits on

$\widetilde{\alpha}\ =\ \frac{y_{e(N)}^{2}}{4\pi Gu^{2}}\,.$ (69)

This constraint is labeled as “Torsional Balances” in Fig. 4, 4, 6, and 6. It
can be seen from these figures that fifth-force experiments constrain both
couplings $y_{e}$ and $y_{N}$ with an upper bound in the range
$10^{-25}-10^{-15}$ for $m_{\phi}<0.1$ eV.

####  4.1.3 Constraints from Stellar and Supernova Cooling

$\phi-e$ coupling: The production of the light scalar $\phi$ in stellar bodies
can lead to a new channel for energy loss leading to rapid cooling. This can
help severely constrain the interaction of scalar with electrons. The dominant
production of this scalar is via its resonant mixing with the longitudinal
component of the photon in the plasma Hardy:2016kme . The extra energy loss
processes in red giants (RG) can delay their onset of helium ignition and can
change the helium-burning lifetime of the horizontal branch (HB) stars, in
disagreement with the stellar models that match observations. For bounds from
supernova, the energy loss from production of scalar is required to be less
than that of SN1987A neutrino burst. The energy loss rate from resonant
production of scalar with plasmon is given by Hardy:2016kme ; Knapen:2017xzo

$Q_{\rm res}\ \simeq\
\frac{\omega_{L}}{4\pi}\left(\frac{\omega_{L}}{m_{\phi}}\Pi^{\phi
L}\right)\frac{1}{e^{\frac{\omega_{L}}{T}}-1}\,,$ (70)

where $\omega_{L}$ is the resonant frequency and $\Pi^{\phi L}$ is the mixing
of the scalar with the longitudinal component of the photon in the medium,
given by

$\Pi^{\phi L}\ \simeq\ \frac{y_{e}em_{e}^{\rm
eff}m_{\phi}}{\pi^{2}k}\int_{0}^{\infty}dp\>v^{2}\left[n_{e}(E_{p})+n_{\overline{e}}(E_{p})\right]\left[\frac{\omega_{L}}{v_{k}}\log\left(\frac{\omega_{L}+vk}{\omega_{L}-vk}\right)-\frac{2m_{\phi}^{2}}{\omega_{L}^{2}-k^{2}v^{2}}\right]\,,$
(71)

where $v=p/E_{p}$ is the electron velocity, $m_{e}^{\rm eff}$ is the effective
thermal mass of the electron and $k=\sqrt{\omega_{L}^{2}-m_{\phi}^{2}}$ is the
3-momentum of the scalar mediator $\phi$, where $E_{\phi}=\omega_{L}$ due to
the resonant production of scalar. Ref. Knapen:2017xzo  considers the resonant
production process as dominant over the Compton scattering or electron-ion
interactions.

For large values of the coupling, the scalar can get trapped inside the
star/supernova. This capture would help alleviate the stringent upper bound on
the coupling $y_{e}$. To derive the trapping limit, the detailed balance of
production and absorption rates is used, i.e.,

$\displaystyle\Gamma_{\rm prod}(E_{\phi})\ =\
e^{-\frac{E_{\phi}}{T}}\Gamma_{\rm abs}(E_{\phi})\,.$ (72)

Since we are only interested in ultra-light mediators with $m_{\phi}<1\text{
MeV}$, the absorption through the decay channel $\phi\rightarrow e^{+}e^{-}$
is absent for our purposes. Thus, the absorption rate from the resonant mixing
yields a mean free path length $\lambda$ given by:

$\lambda\ =\ \frac{1}{\Gamma_{\rm abs}(E_{\phi})}\ \sim\
\frac{E_{\phi}^{4}}{Q_{\rm res}}~{}.$ (73)

By requiring the mean free path to be shorter than $R=10\text{ km}$, which is
the typical size of a supernova core, we derive a bound on the coupling
$y_{e}$, as shown in Figs. 4 and 6, labeled “SN1987A”.

In case of SN1987A, constraints on $y_{e}$ range from $10^{-9}$ to $10^{-7}$
for scalar mediators lighter than electron. Even stronger constraints are
obtained from HB/RG stars with an upper bound of $y_{e}\sim 10^{-15}$ for
light scalar mediators.

$\phi-N$ coupling : The constraints are similar to the $\phi-e$ coupling case.
In HB and RG stars with typical temperatures of $10$ keV, the main constraints
for scalar coupling to nucleon in the literature are derived using Compton
scattering, $\gamma+\text{He}\rightarrow\text{He}+\phi$, as the dominant
process. It is required that the new energy loss per unit mass should be less
than $\epsilon<10$ erg/g/s Raffelt:1996wa . As shown in Ref. Hardy:2016kme ,
resonant production through $\phi$ mixing with photon can increase the energy
loss for low scalar masses and therefore are highly constrained.

The constraints from supernova comes from scalar production through
bremsstrahlung process $N+N\rightarrow N+N+\phi$ Ishizuka:1989ts . Bounds on
coupling can be obtained by requiring the energy loss to be less than the
energy contained in the neutrino burst, i.e., $\epsilon<10^{19}$ erg/g/s
Raffelt:1996wa . Similarly, the trapping regime of the scalar being reabsorbed
can be derived using the detailed balance between the absorption and
production rates. Requiring the mean free path
$\lambda\propto\epsilon\rho/T^{4}$ to be smaller than 10 km yields the
constraint on $y_{N}$ Knapen:2017xzo , as shown in Figs. 4 and 6.

In case of SN1987A, constraints on $y_{N}$ range from $10^{-10}$ to $10^{-7}$
for scalar mediators lighter than electron. Similar to $y_{e}$, stronger
constraints are obtained from HB/RG stars with an upper bound of $y_{N}\sim
10^{-12}$ for light scalar mediators.

####  4.1.4 Meson decays

A light scalar coupling to nucleons can be produced in meson decays. The only
process of interest in this case is charged Kaon decay to charged pion and the
scalar: $K^{+}\rightarrow\pi^{+}\phi$. This production cross section is highly
constrained from the measurement of branching ratios from charged Kaon decay :
$\text{Br }(K^{+}\rightarrow\pi^{+}\overline{\nu}\nu)<1.7\times 10^{-10}$
Tanabashi:2018oca .

Using the low-energy effective Lagrangian formalism presented in Ref.
Batell:2018fqo , the branching ratio for the process in consideration is given
by

${\rm BR}(K^{+}\rightarrow\pi^{+}\>\phi)\ =\
\frac{(3y_{u}G_{F}f_{\pi}f_{K}B)^{2}}{32\pi
m_{K^{+}}\Gamma_{K^{+}}}|V_{ud}V_{us}|^{2}\>\lambda^{1/2}\left(1,\frac{m_{\phi}^{2}}{m_{K^{+}}^{2}},\frac{m_{\pi^{+}}^{2}}{m_{K^{+}}^{2}}\right)\,,$
(74)

where $B=\frac{m_{\pi}^{2}}{m_{u}+m_{d}}$ and
$\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2bc-2ac$. Matching the nucleon level
interaction to the effective Lagrangian:

$\mathcal{L}\ \supset\ y_{N}\overline{N}N\phi\,,$ (75)

where $N=p,\,n$, the nucleon coupling $y_{N}$ can be written in terms of
fundamental quark -level couplings $y_{u}(y_{d})$ as :

$y_{N}\ =\ \sum_{q}y_{q}g^{q}_{S}\,,$ (76)

where $g^{q}_{S}$ is the nucleon scalar charge. We assume that the scalar
couples equally to the up and down quark i.e. $y_{u}=y_{d}$. The effective
nucleon couplings to scalar is then given by

$\displaystyle y_{N}$ $\displaystyle\ =\
y_{u}\left(g^{u}_{S}+g^{d}_{S}\right)\simeq 9.47\>y_{u}\,,$ (77)

where we have used $g^{u}_{S}=5.20$ and $g^{d}_{S}=4.27$ Alexandrou:2017qyt .
This constraint is labeled as “$K^{+}\rightarrow\pi^{+}\phi$” in Figs. 4 and
6. It yields an almost constant upper bound of $y_{N}\sim 2.3\times 10^{-5}$
for light scalar mediators.

####  4.1.5 Big Bang Nucleosynthesis

$\phi-e$ coupling: In early universe, the scalar mediator $\phi$ can be in
thermal equilibrium with the SM particles through annihilation
($e^{+}e^{-}\rightarrow\gamma\phi$) and Compton scattering
($e^{-}\gamma\rightarrow e^{-}\phi$). In the limit $s\gg
m_{\phi}^{2},m_{e}^{2}$, the cross sections for these processes are
Knapen:2017xzo

$\displaystyle\sigma_{e\gamma\rightarrow e\phi}\ $ $\displaystyle\approx\
\frac{\alpha_{e}y_{e}^{2}}{s}\left[\log\left(\frac{s}{m_{e}^{2}+m_{\phi}^{2}}\right)+\frac{5}{2}\right]\,,$
(78) $\displaystyle\sigma_{ee\rightarrow\gamma\phi}\ $ $\displaystyle\approx\
\frac{2\alpha_{e}y_{e}^{2}}{s}\log\left(\frac{s}{4m_{e}^{2}}\right)~{},$ (79)

where $\alpha_{e}\equiv e^{2}/4\pi$ is the fine-structure constant. The
thermally averaged cross section for these two processes are given below:

$\displaystyle\langle\sigma_{e\gamma\rightarrow e\phi}\>v\rangle\ $
$\displaystyle=\
\frac{1}{16m_{e}^{2}T^{3}K_{2}(m_{e}/T)}\int_{m_{e}^{2}}^{\infty}ds\>\sigma(s-m_{e}^{2})\sqrt{s}K_{1}\left(\frac{\sqrt{s}}{T}\right)\,,$
(80) $\displaystyle\langle\sigma_{ee\rightarrow\gamma\phi}\>v\rangle\ $
$\displaystyle=\
\frac{1}{8m_{e}^{4}T(K_{2}(m_{e}/T))^{2}}\int_{4m_{e}^{2}}^{\infty}ds\>\sigma(s-4m_{e}^{2})\sqrt{s}K_{1}\left(\frac{\sqrt{s}}{T}\right)~{}.$
(81)

If $\phi$ enters equilibrium with electrons before $T\sim 1\text{ MeV}$, it
can decrease the deuterium abundance which is in conflict with observations
Knapen:2017xzo . In our case, the mediator thermalizes if the thermally
averaged cross section exceeds the Hubble expansion rate $H(T)\sim
1.66\sqrt{g^{*}}T^{2}/M_{\text{Pl}}$ (where $g_{*}$ is the number of
relativistic degrees of freedom and $M_{\rm Pl}$ is the Planck mass) at
$T=1\text{ MeV}$. This yields an upper bound of $y_{e}=5\times 10^{-10}$ for
ultra-light scalar mediators, independent of $m_{\phi}$.

$\phi-N$ coupling : In this case, we require that the scalar $\phi$
thermalizes around QCD phase transition temperature. This will help dilute the
relativistic degrees of freedom ($N_{\rm eff}$) until the nucleosynthesis
phase is reached. Otherwise, the scalar $\phi$ will be in equilibrium with SM
and will have a significant contribution to relativistic degrees of freedom
($\Delta N_{\rm eff}=4/7$) at the time of BBN, in tension with the current
measurements from Planck Aghanim:2018eyx . Thus, we require that the
interaction rate should be lower than the Hubble rate at $T=200$ MeV. We can
estimate the rate of $\phi$ production from the processes like
$u\bar{u}\rightarrow\phi$ (and $d\bar{d}\rightarrow\phi$) as
$\Gamma_{\phi}\sim y_{u}^{2}T$. This should be compared with Hubble rate
$H(T)\sim 1.66\sqrt{g^{*}}T^{2}/M_{\text{Pl}}$. This condition leads to a
stringent constraint on $y_{u}<2.63\times 10^{-10}$. Converting the quark-
scalar coupling to nucleon level coupling using Eq. (77) , we get
$y_{N}<2.49\times 10^{-9}$.

###  4.2 Experimental Constraints on $y_{\nu}$

Dirac $\nu-\phi$ coupling: The analysis in this case is similar as for the
$\phi-N$ coupling. We require that the scalar $\phi$ and right-handed
neutrinos should decouple from the plasma at a temperature above the QCD phase
transition temperature which will dilute the $\Delta N_{\rm
eff}=3+\frac{4}{7}\sim 3.57$ by the time BBN occurs, in agreement with the
currently allowed range from Planck Aghanim:2018eyx . Thus, requiring that the
interaction rate of processes like $\nu\bar{\nu}\rightarrow\phi$ should be
lower than the Hubble rate at $T=200$ MeV yields an upper bound of
$y_{\nu}\sim 2.6\times 10^{-10}$.

Majorana $\nu-\phi$ coupling: Presence of NSI can lead to re-thermalization of
the neutrinos, which otherwise decouple at $T\sim 1$ MeV in the standard
scenario. This can leave a signature in the cosmological observables. The
analysis in Ref. Forastieri:2019cuf  constrains the couplings in the secret
interaction of neutrinos with a light mediator. Assuming model independence,
we use the upper bound on coupling $y_{\nu}$ from Ref. Forastieri:2019cuf ,
which yields a stringent limit of $y_{\nu}<2\times 10^{-7}$.

##  5 Thermal mass of scalar $\phi$

Figure 2: Feynman diagram responsible for the thermal mass of the scalar
$\phi$.

If the interactions of the scalar $\phi$ with the medium are significant
enough, then it might get thermalized with the medium. Since the scalar field
in consideration is ultra-light, medium effects might lead to substantial
correction to the vacuum mass of the $\phi$. The medium induced mass at one-
loop is shown in the Feynman diagram in Fig. 2. The relevant contribution to
the mass of $\phi$ at finite density and temperature is given by:

$\mathcal{M}\ =\
4y_{f}^{2}\int\frac{d^{4}k}{(2\pi)^{4}}\left(k^{2}-\frac{p^{2}}{4}+m_{f}^{2}\right)\left[\frac{\Gamma(k+p/2)}{(k-p/2)^{2}-m_{f}^{2}}+\frac{\Gamma(k-p/2)}{(k+p/2)^{2}-m_{f}^{2}}\right]~{}.$
(82)

We refer the reader to Appendix C for the evaluation of the scalar mass
integral. In the limit $m_{\phi}\rightarrow 0$, the mass correction for scalar
is found to be:

$\Delta m_{\phi}^{2}\ =\
\frac{y_{f}^{2}}{\pi^{2}}\int_{m_{f}}^{\infty}dk_{0}\>n_{f}(k_{0})\sqrt{k_{0}^{2}-m_{f}^{2}}\
=\ \frac{y_{f}^{2}}{3}\left(\frac{\pi^{2}n_{f}}{12\zeta(3)}\right)^{2/3}\,,$
(83)

where we have used the result for the same integral appearing in Eq. (14) in
the high-temperature limit. Eq. (83) will be relevant to deriving neutrino
self-interaction limits from early universe cosmology. It has been found that
CMB anisotropy data constrain such interactions to be
$(y_{\nu}^{2}/m_{\phi}^{2})\leq(3~{}{\rm MeV})^{-2}$ (for $m_{\phi}>1$ keV)
Kreisch:2019yzn . If the scalar field indeed thermalizes with the medium,
which occurs for $y_{\nu}\geq 10^{-10}$ or so, then one should use the thermal
mass of $\phi$, Eq. (83) in this constraint, which can weaken the constraint
significantly. In cosmological simulations involving a light scalar, the
thermal mass effects of Eq. (83) should be included. Such interactions may be
testable in future CMB and large-scale structure observations through the
thermally induced mass in such settings.

In this limit for $\Delta m_{\phi}^{2}$, the scalar NSI expression Eq. (14)
takes a special form:

$\Delta m_{\nu,{\alpha\beta}}\ =\ \frac{y_{\alpha\beta}}{y_{f}}m_{f}~{}.$ (84)

Note that Eq. (84) is independent of the scalar mass $m_{\phi}$. This scenario
may be realized in supernovae, provided that $\phi$ has significant
interactions with matter. From discussions in Sec. 4.1.3, it is clear that for
high enough values of $y_{e}$ or $y_{N}$, the scalar gets trapped and thermal
correction to the mass should be taken into account. Thus, in case of
thermalization of the scalar, Eq. (84) should be used in lieu of Eqs. (15),
(16) and (17).

##  6 Quantum-mechanical bound on light scalar mass

Here we show that the uncertainty principle of quantum mechanics sets a lower
limit on the minimum $q^{2}$ that appears in neutrino forward scattering. This
limit applies to neutrino propagating through Earth, where it interacts either
with electrons in atoms, or with nucleons inside the nuclei. This limit is not
applicable to Sun or supernovae, where the temperature of the medium is too
high that there are no bound atoms/nuclei.

Consider $\nu_{\alpha}-e$ elastic scattering. Working in the rest frame of the
electron, the initial and final four-momenta of the electron can be written as

$p^{\mu}\ =\ (m_{e},\,0,\,0,\,0)\,,\qquad p^{\prime\mu}\ =\
(\sqrt{p_{e}^{2}+m_{e}^{2}},\,0,\,0,\,p_{e})\,,$ (85)

where $p$ is the recoil momentum of the electron. The $q^{2}$ related to
coherent forward scattering is then

$q^{2}\ =\ (p^{\prime}-p)^{2}\ =\ 2m_{e}(m_{e}-\sqrt{p_{e}^{2}+m_{e}^{2}})\
\simeq\ -p_{e}^{2},$ (86)

where in the second step $q^{2}\ll m_{e}^{2}$ is assumed.

Now, the recoil momentum of the electron is subject to the uncertainty
relation. Its position is not precisely known inside the atom, so we have

$\Delta p\,\Delta x\ \gtrsim\ \hbar~{}.$ (87)

When we set $q^{2}=0$ in the computation of forward scattering, we only know
this up to an uncertainty in $q^{2}$ given by

$q^{2}\ \simeq\ p_{e}^{2}\ \sim\ \frac{1}{(\Delta x)^{2}}~{}.$ (88)

Using $\Delta x=140\times 10^{-8}$ cm, which is the radius of 26Fe – the most
abundant element in Earth’s matter, one obtains for the uncertainty in $q^{2}$
to be

$q^{2}_{y_{e}}\ \approx\ (14~{}{\rm eV})^{2}~{}.$ (89)

Thus, when the mediator mass becomes much smaller than $14$ eV, one should use
this quantum mechanical cut-off in computing scalar NSI. Similarly for
coupling to nucleon, the cut-off would be given by the inverse of the nuclear
radius of 26Fe. Using nuclear diameter $\Delta x=9.6$ fm, we obtain

$q^{2}_{y_{N}}\ \approx\ (21~{}{\rm MeV})^{2}~{}.$ (90)

These rough quantum-mechanical bounds can be better motivated by using
atomic/nuclear form factors for coherent forward scattering. In Earth, the
expression for scalar NSI will get modified with the inclusion of a form
factor.

$\Delta m_{\nu,{\alpha\beta}}\ =\
\frac{y_{f}y_{\alpha\beta}N_{f}}{m_{\phi}^{2}-q^{2}}F(m_{\phi}^{2})\,,$ (91)

The original result in Eq. (15) was obtained by setting $q^{2}=0$ and
$F(m_{\phi}^{2})=1$, but if the mass of the scalar $m_{\phi}\rightarrow 0$
then the denominator is not well-defined. This is remedied by the atomic form
factor $F(m_{\phi}^{2})$ which is of the form hubbell1979relativistic :

$F(m^{2})\ =\ \frac{m^{2}}{m^{2}+q_{0}^{2}}\,,$ (92)

where $q_{0}=1/4\pi a_{0}$ and $a_{0}$ is the radius of the first orbit for
hydrogen-like atoms. Similar qualitative results should apply for the
outermost-orbit electrons in 26Fe. For high values of $m^{2}_{\phi}\gg
q_{0}^{2}$, $F(m_{\phi}^{2})\sim 1$ as expected. Thus, the vanishing $q^{2}$
limit is well-defined and yields the original result in Eq. (15). Difference
appears in the regime $m_{\phi}^{2}\ll q_{0}^{2}$, where $F(m_{\phi}^{2})\sim
m_{\phi}^{2}/q_{0}^{2}$. The form of Eq. (91) in the low $m_{\phi}$ limit and
with $q^{2}\to 0$ is thus given by:

$\Delta m_{\nu,{\alpha\beta}}\ =\
\frac{y_{f}y_{\alpha\beta}N_{f}}{q_{0}^{2}}~{},$ (93)

which is independent of $m_{\phi}$. This result agrees with the quantum-
mechanical bound discussed above based on the uncertainty principle.

When a scalar mediator couples to the electron, from fifth force constraints
either the mass of the mediator should be larger than a keV, or its coupling
to the electron should be extremely weak, of order $10^{-24}$. For such tiny
couplings, to generate scalar NSI in the observable range, one could naively
make the mediator mass of order $10^{-8}$ eV. In this case, the quantum-
mechanical intrinsic bound should be applied for computing forward scattering.
The result is that scalar NSI arising from coupling to electrons cannot be in
the observable range for neutrino propagation in Earth.

##  7 Results and discussion

Figure 3: Different experimental constraints on Yukawa coupling of scalar to
electron for the case of Dirac neutrinos. The shaded regions are excluded.
Some representative values of scalar NSI in Earth, Sun and supernova are also
shown.

Figure 4: Same as in Fig. 4, but for scalar coupling to nucleons.

Figure 5: Same as in Fig. 4, but for Majorana neutrinos.

Figure 6: Same as in Fig. 4, but for Majorana neutrinos.

We have discussed the calculation for scalar NSI and experimental constraints
in previous sections. The results for different cases with scalar coupling to
electron/nucleon and in case of either Dirac or Majorana neutrinos have been
presented in Figs. 4, 4, 6, 6. Here we have fixed the value of $y_{\nu}$ at
its maximum allowed value in each case, as discussed in Sec. 4.2, whereas the
other Yukawa coupling (either $y_{e}$ or $y_{N}$) is varied, along with the
scalar mass $m_{\phi}$.

###  7.1 Earth and Sun

In case of Earth and Sun, the background medium of electrons and nucleons are
non-relativistic. Therefore, the expression used for scalar NSI in these media
is given by Eq. (15) with $N_{\bar{f}}=0$:

$\Delta m_{\nu,{\alpha\beta}}\ =\
\frac{y_{f}y_{\alpha\beta}}{m_{\phi}^{2}}N_{f}~{}.$ (94)

From the discussion in Sec. 6, when the mediator mass becomes lower than the
quantum mechanical cut-off of $m_{0}\sim 14$ eV, $m_{0}^{2}$ should be used in
the denominator of Eq. (94) in lieu of $m_{\phi}^{2}$ for Earth. This leads to
the turning of the scalar NSI line in the plots for Earth. We have used
$N_{e}^{\text{Earth}}=5.4\>N_{A}\text{ cm}^{-3}$ Tanabashi:2018oca  and
$N_{N}^{\text{Earth}}=\frac{2.9}{m_{N}}\>\text{ g}\text{ cm}^{-3}$ PREM ,
where the nucleon mass $m_{N}=931.5$ MeV and the Avogadro number
$N_{A}=6.022\times 10^{23}$. As can be seen from the plots, there are no
prospects for observable scalar NSI to be detected on Earth in any of the four
cases (Dirac/Majorana and coupling to electrons/nucleons). It can be seen from
Table. 2 that highest allowed value of scalar NSI in case of Earth is around
$10^{-14}$ eV for the case of $\phi$ coupling to Majorana neutrinos and
electrons.

For the case of Sun, there will also be correction to the scalar NSI from
finite size of the medium in the case of light mediators masses
$m_{\phi}\simeq R_{\rm Sun}^{-1}$ as discussed in Sec. 3 and Ref. Smirnov_2019
. We calculate the form factor for Sun using Eq. (46) and the number density
of electrons/nucleons, which is obtained by fitting the known solar density
profile given in Refs. SolarM ; Bahcall:1998wm ; Bahcall_2005 . We have used
the following best fit to the number density profile for Sun:

$\displaystyle N(r)_{\text{e}}$ $\displaystyle\ =\
111.61N_{A}\;e^{-(4.81\;r+10.21\>r^{2})}\text{ cm}^{-3}\qquad$
$\displaystyle\text{(Electron)}\,,$ (95) $\displaystyle N(r)_{N}$
$\displaystyle\ =\ \frac{157.13}{m_{N}}\;e^{-(6.1\;r+5.2\>r^{2})}\text{
g}\text{ cm}^{-3}$ $\displaystyle\text{(Nucleon)}\,.$ (96)

As can be seen from the plots, the existing laboratory and astrophysical
constraints do allow for a sizable scalar NSI in the Sun, especially for
$m_{\phi}\lesssim 1\mu$eV where the NSI can be as large as $10^{5}$ eV for the
case of $\phi$ coupling to Dirac/Majorana neutrinos and electrons. However,
this will lead to a large correction term to the solar neutrino mass, which is
severely constrained by solar neutrino data. Using the $\chi^{2}$-analysis of
the Borexino data from Ref. Ge:2018uhz , we find a $3\sigma$ upper bound on
the scalar NSI in Sun: $\Delta m_{\rm Sun}\lesssim 7.4\times 10^{-3}$ eV, as
shown by the yellow shaded region in Figs. 4, 4, 6, and 6. This still leaves
some room for observable scalar NSI effects in future solar neutrino data,
especially for ultra-light scalar mediators. Note that very small coupling
values for which $y_{f}^{2}\lesssim Gm^{2}_{\nu}=(m_{\nu}/M_{\rm Pl})^{2}\sim
10^{-30}$ are disfavored by the weak gravity conjecture ArkaniHamed:2006dz
which suggests gravity as the weakest force in nature.

###  7.2 Supernovae

In the case of supernovae with a typical core temperature $T\sim 30$ MeV, the
electron background is relativistic while the nucleon background can be
essentially treated to be at rest. Thus, there are two different expressions
to be used [cf. Eqs. (15) and (16)]:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}\ =\ $
$\displaystyle\frac{y_{f}y_{\alpha\beta}}{m_{\phi}^{2}}N_{N}^{\rm SN}$ (for
nucleon) (97) $\displaystyle\Delta m_{\nu,{\alpha\beta}}\ =\ $
$\displaystyle\frac{y_{\alpha\beta}y_{f}}{m_{\phi}^{2}}\frac{m_{e}}{2}\left(\frac{3N_{e}^{\rm
SN}}{\pi}\right)^{\frac{2}{3}}$ (for electron) . (98)

Similar to the case in Sun, there will be correction to the scalar NSI in
supernova from the finite size of the medium. Therefore, we numerically
integrate Eq. (46) to obtain the form factor for a realistic supernova density
profile. We use the fiducial model parameters from Ref. Chang_2017  given
below:

$\displaystyle\rho(r)\ =\ \rho_{c}\times$ $\displaystyle
1+k_{\rho}(1-r/R_{c})\;$ ($r<R_{c}$) (99) $\displaystyle\rho(r)\ =\
\rho_{c}\times$ $\displaystyle(r/R_{c})^{-\eta}$ ($r\geq R_{c}$) (100)

where $\rho_{c}=3\times 10^{14}\text{ g}\text{ cm}^{-3}$ is the density at
core radius $R_{c}=10$ km , $k_{\rho}=0.2$ and $\eta=5$. Assuming the medium
to be electrically neutral and using a proton fraction $Y_{p}=0.3$, we can
obtain the number density for electrons from $\rho(r)$.

An interesting feature emerges for scalar NSI in supernova. Due to the high
temperature, a light scalar might develop a considerable thermal mass if it
has strong enough coupling to the background as discussed in Sec. 5. This
leads to Eq. (84) which is independent of $m_{\phi}$. Trapping leads to the
thermalization of the scalar in the medium. Thus, we have only plotted the
scalar NSI expression for supernova as long as it is not trapped inside.

Scalar NSI produced in supernova cannot be arbitrarily high. If it becomes too
large, then neutrino production would be affected in direct conflict with
observations from SN1987A. For typical supernova core temperature around
$T\simeq 30$ MeV, we constrain the scalar NSI to be less than $5$ MeV
Smirnov_2019 , so that 10-MeVish neutrinos could be detected on Earth from
SN1987A. In the plots, this bound is shown as a dashed line marked $\Delta
m_{\rm SN}>5$ MeV. In any case, we find that sizable scalar NSI can still be
observed in supernovae, while being consistent with all other constraints.

Case | Max. NSI (eV) | Scalar Mass Range (eV) | Range for $y_{f}$
---|---|---|---
Dirac $\nu$, $\phi-e$ |  |  |
Earth | $3\times 10^{-17}$ | 0.04 -14 | $\sim 7\times 10^{-16}$
Sun | $7\times 10^{-3}$ | $<10^{-11}$ | $3.3\times 10^{-34}\>-10^{-26}$
Supernova | $5\times 10^{6}$ | $10^{-11}\>-10^{-9}$ | $10^{-26}\>-1.8\times 10^{-23}$
Dirac $\nu$, $\phi-N$ |  |  |
Earth | $10^{-24}$ | $5.3\times 10^{3}\>-2.1\times 10^{7}$ | $\sim 2.4\times 10^{-10}$
Sun | $7\times 10^{-3}$ | $<3.3\times 10^{-13}$ | $2.4\times 10^{-34}\>-7.5\times 10^{-30}$
Supernova | $5\times 10^{6}$ | $3.3\times 10^{-13}\>-1.8\times 10^{-7}$ | $7.5\times 10^{-30}\>-4.9\times 10^{-22}$
Majorana $\nu$, $\phi-e$ |  |  |
Earth | $10^{-14}$ | 0.04 -14 | $\sim 6\times 10^{-16}$
Sun | $7\times 10^{-3}$ | $<10^{-11}$ | $4.4\times 10^{-37}\>-8.7\times 10^{-30}$
Supernova | $5\times 10^{6}$ | $10^{-11}\>-7\times 10^{-8}$ | $8.7\times 10^{-30}\>-9.3\times 10^{-23}$
Majorana $\nu$, $\phi-N$ |  |  |
Earth | $10^{-21}$ | $5.3\times 10^{3}\>-2.1\times 10^{7}$ | $\sim 2.1\times 10^{-10}$
Sun | $7\times 10^{-3}$ | $<3.5\times 10^{-13}$ | $3.1\times 10^{-37}\>-8.4\times 10^{-33}$
Supernova | $5\times 10^{6}$ | $3.5\times 10^{-13}\>-1.3\times 10^{-5}$ | $8.4\times 10^{-33}\>-2\times 10^{-21}$
Table 2: The maximum allowed value of scalar NSI in different cases and
domains with corresponding ranges for the scalar mass $\phi$ and the coupling
strength $y_{f}$.

##  8 UV-complete model for scalar NSI

In this section we sketch possible ultra-violet completions that would induce
interactions of neutrinos with a light scalar. This discussion is intended
only as a proof of principle. We focus on the case of Dirac neutrinos, with a
light scalar $\phi$ coupling to the neutrinos and the electron.

First we construct two effective operators that are invariant under the SM
gauge symmetry. One induces couplings of the scalar $\phi$ to neutrinos and
the other to the electron. These operators are

$(i)~{}~{}~{}~{}\overline{\psi}_{L}\widetilde{H}\nu_{R}\frac{\phi}{\Lambda_{\nu}},~{}~{}~{}~{}~{}~{}~{}(ii)~{}~{}~{}~{}\overline{\psi}_{L}He_{R}\frac{\phi^{2}}{\Lambda_{e}^{2}}~{}.$
(101)

Here $\phi$ is a real scalar field, which is a singlet under SM symmetry,
$H=\big{(}\begin{smallmatrix}H^{+}\\\ H^{0}\end{smallmatrix}\big{)}$ is the SM
Higgs doublet and $\psi_{L}=\big{(}\begin{smallmatrix}\nu\\\
e\end{smallmatrix}\big{)}_{L}$ is the left-handed lepton doublet. These
effective operators exhibit a $Z_{2}$ symmetry (apart from lepton number)
under which $\nu_{R}$ and $\phi$ are odd, with other fields being even. $\phi$
develops a vacuum expectation value, $\langle\phi\rangle=v_{\phi}\sim 10$ eV,
which breaks the $Z_{2}$ symmetry. The neutrino Yukawa coupling $y_{\nu}$ and
the electron Yukawa coupling $y_{e}$ with the $\phi$ field are respectively
given by

$y_{\nu}\ =\ \frac{v}{\Lambda_{\nu}}\,,\qquad y_{e}\ =\
\frac{2vv_{\phi}}{\Lambda_{e}^{2}}$ (102)

where $v=174$ GeV is the VEV of the SM Higgs doublet. Once $\phi$ acquires a
VEV, the operator $(i)$ generates a mass term for the neutrino given by

$m_{\nu}\ =\ \frac{v_{\phi}v}{\Lambda_{\nu}}~{}.$ (103)

While this may be the leading contribution, it is not required to be so, as
there could be other contributions as well. In any case, this would imply an
upper limit on $y_{\nu}$ given by

$y_{\nu}\ <\ \frac{m_{\nu}}{v_{\phi}}~{}.$ (104)

The cut-off scale $\Lambda_{e}$ is expected to be at least a hundred GeV,
while $\Lambda_{\nu}$ may be lower. Choosing $\Lambda_{e}\sim v$, we would
have $y_{e}\sim v_{\phi}/v$. For $y_{e}\sim 10^{-10}$, as our analysis
requires for observable scalar NSI, $v_{\phi}\sim 10$ eV is preferable. This
in turn implies from Eq. (104) that $y_{\nu}<5\times 10^{-3}$, using
$m_{\nu}\equiv\sqrt{\Delta m^{2}_{\rm atm}}\sim 0.05$ eV. $y_{\nu}$ of course
can be smaller than this value, which would be in the interesting range for
observable scalar NSI.

The operators in Eq. (101) can be generated by adding new vector-like fermions
to the SM. For example, operator $(i)$ can arise by the addition of SM singlet
fermions $N_{L,R}$ with a a lepton number preserving Dirac mass. The relevant
Lagrangian is given by

$\mathcal{L}\ \supset\
y_{N}\overline{\psi}_{L}\widetilde{H}N_{R}+M_{N}\overline{N}_{R}N_{L}+y^{\nu}_{\phi}\overline{N}_{L}\nu_{R}\phi+{\rm
H.c.}$ (105)

These interactions also preserve the $Z_{2}$ symmetry with $N_{L,R}$ being
even under it. The diagram generating operator $(i)$ is shown in Fig. 7, left
panel.

Operator $(ii)$ is induced by integrating out a pair of vector-like leptons,
$E,\,E^{\prime}$, both being singlets of $SU(2)_{L}$ and carrying hypercharge
$Y=-2$. Their interaction Lagrangian is given by

$\mathcal{L}\ \supset\
y_{E}\overline{\psi}_{L}{H}E_{R}+\mu_{E}\overline{E}_{R}E_{L}+y^{E}_{\phi}\overline{E}_{L}E^{{}^{\prime}}_{R}\phi+M_{E^{{}^{\prime}}}\overline{E}_{R}^{{}^{\prime}}E_{L}+y^{e}_{\phi}\overline{E}_{L}^{{}^{\prime}}e_{R}\phi+{\rm
H.c.}$ (106)

Here $E_{L,R}$ are even and $E_{L,R}^{{}^{\prime}}$ are odd under $Z_{2}$. The
effective operator involving electron and $\phi$ is generated by Fig. 7, right
panel.

Figure 7: Explicit models generating operators of Eq. (101).

Integrating out the heavy degrees of freedom we obtain the following effective
Lagrangian terms:

$(i)~{}~{}~{}\frac{y_{N}y^{\nu}_{\phi}}{M_{N}}\overline{\psi}_{L}\widetilde{H}\nu_{R}\phi,~{}~{}~{}~{}~{}~{}~{}(ii)~{}~{}~{}\frac{y_{E}y^{E}_{\phi}y^{e}_{\phi}}{M_{E}M_{E^{{}^{\prime}}}}\overline{\psi}_{L}{H}e_{R}\phi^{2}~{}.$
(107)

These expressions can be mapped to Eq. (101) to identify the cut-off scales
$\Lambda_{\nu}$ and $\Lambda_{e}$, and the constraints discussed in terms of
the cut-off scales will apply to them. We thus see broad consistency of the
model. In particular, the induced neutrino mass from these interactions is not
excessive and the vector-like leptons having mass of order few hundred GeV is
consistent with collider data. Note that breaking the $Z_{2}$ at a scale of
order 10 eV does not cause cosmological domain wall problem, since the energy
density carried by the walls is quite small. We have ignored here possible
mixing between the $\phi$ and $H$ fields since such mixing is small, of order
$v_{\phi}/v$ and is controlled by a new quartic coupling which may also be
small.

##  9 Conclusion

We have performed a systematic study of scalar NSI of neutrinos with matter
due to a light scalar mediator. First, a general field-theoretic derivation of
the scalar NSI formula is given, which is valid at arbitrary temperature and
density, and hence, applicable in widely different environments, such as
Earth, Sun, supernovae and early universe. We have also extended the analysis
of long-range force effects for all background media, including both
relativistic and non-relativistic limits. Using these results and applying
various experimental and astrophysical constraints, we find that observable
scalar NSI has been precluded in terrestrial experiments, primarily due to
atomic form factor effects, which can also be understood from simple quantum-
mechanical uncertainty principle. Nevertheless, sizable scalar NSI effects are
still possible in the Sun, supernovae and early universe environments, which
could be detected in future solar and supernova neutrino data, as well as in
the form of extra relativistic species ($\Delta N_{\rm eff}$) and neutrino
self-interactions in cosmological observations. We have also presented
examples of UV-complete models that could give rise to such scalar NSI
effects.

## Acknowledgments

We thank Mark Alford, Steven Harris, Ahmed Ismail and Pedro Machado for
discussions. This work was supported in part by US Department of Energy Grant
Numbers DE-SC 0016013 (K.S.B.) and DE-SC0017987 (B.D. and G. C.). This work
was also supported by the Neutrino Theory Network Program under Grant No. DE-
AC02-07CH11359. K.S.B. and B.D. thank the Fermilab Theory Group for warm
hospitality, where part of this work was done. B.D. and G.C. also thank the
Department of Physics at Oklahoma State University for warm hospitality during
the completion of this work.

## Appendices

##  Appendix A Limiting cases for scalar NSI expression

In this Appendix we evaluate the self-energy given in Eq. (14) corresponding
to the tadpole diagram of Fig. 1. We shall evaluate only the fermionic
contribution to Eq. (14), from which it is easy to read of the anti-fermionic
background contribution as well. We also provide an exact expression for the
medium-dependent neutrino mass, which can be evaluated numerically.

###  A.1 Case 1: $\mu>m_{f}\gg T$

Breaking the integration limits and expanding the occupation number as an
infinite series, we can write Eq. (14) as follows:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}$ $\displaystyle\ =\
\frac{m_{f}y_{\alpha\beta}y_{f}}{2\pi^{2}m_{\phi}^{2}}\left(\left[\mu\sqrt{\mu^{2}-m_{f}^{2}}+m_{f}^{2}\ln{\left(\frac{m_{f}}{\mu+\sqrt{\mu^{2}-m_{f}^{2}}}\right)}\right]\right.$
$\displaystyle+\left.\sum_{n=1}^{\infty}(-1)^{n}\left[\int_{m_{f}}^{\mu}dE\>e^{n(E-\mu)/T}\sqrt{E^{2}-m_{f}^{2}}+\int_{\mu}^{\infty}dE\>e^{-n(E-\mu)/T}\sqrt{E^{2}-m_{f}^{2}}\right]\right)~{}.$
(108)

As $T\rightarrow 0$, the first term in the series dominates the result. We
know that sum over all momentum states weighted by occupation number yields
the number density. Inverting the relation to obtain $\mu$, we get:

$\mu^{2}\ =\ (3\pi^{2}N_{f})^{\frac{2}{3}}+m_{f}^{2}\ \simeq\
(3\pi^{2}N_{f})^{\frac{2}{3}},$ (109)

where in the second relation we assumed $\mu^{2}\gg m_{f}^{2}$. Thus, for
$\mu\gg m_{f}$ we have

$\Delta m_{\nu,{\alpha\beta}}\ \simeq\
\frac{y_{\alpha\beta}y_{f}}{m_{\phi}^{2}}\frac{m_{f}}{2}\left(\frac{3N_{f}}{\pi}\right)^{\frac{2}{3}}~{},$
(110)

as given in Eq. (16).

###  A.2 Case 2: $T\ll\mu<m_{f}$

When $\mu<m_{f}$, the expression for $\Sigma$ of Eq. (14) can be written as a
weighted series of modified Bessel function of the second kind:

$\Delta m_{\nu,{\alpha\beta}}\ =\
\frac{m_{f}y_{\alpha\beta}y_{f}}{\pi^{2}m_{\phi}^{2}}\sum_{n=1}^{\infty}(-1)^{n+1}\>\frac{m_{f}T}{n}e^{n\mu/T}K_{1}\left(\frac{nm_{f}}{T}\right)\,.$
(111)

For $z\rightarrow\infty$, we can use the asymptotic form for
$K_{\nu}\left(z\right)$ :

$K_{\nu}\left(z\right)\ \simeq\
e^{-z}\sqrt{\frac{\pi}{2z}}\left(1+\frac{4\nu^{2}-1}{8z}+...\right)~{}.$ (112)

Due to the exponential suppression, the $n=1$ term in the sum will be dominant
in Eq. (111). This yields:

$\Delta m_{\nu,{\alpha\beta}}\ \simeq\
\frac{2y_{f}y_{\alpha\beta}}{m_{\phi}^{2}}\left(\frac{m_{f}T}{2\pi}\right)^{\frac{3}{2}}e^{-(m_{f}-\mu)/T}~{}.$
(113)

To relate the above function to the number density $N_{f}$, we use

$\displaystyle N_{f}$ $\displaystyle\ =\
2\int\frac{d^{3}k}{(2\pi)^{3}}\frac{1}{e^{(E-\mu)/T}+1}$ $\displaystyle\ =\
\frac{1}{\pi^{2}}\int_{m_{f}}^{\infty}dE\>\frac{E\sqrt{E^{2}-m_{f}^{2}}}{e^{(E-\mu)/T}+1}$
$\displaystyle\ =\
\frac{1}{\pi^{2}}\sum_{n=1}^{\infty}\int_{m_{f}}^{\infty}dE\>E\sqrt{E^{2}-m_{f}^{2}}\>e^{-n(E-\mu)/T}(-1)^{n+1}$
$\displaystyle\ =\
\frac{1}{\pi^{2}}\sum_{n=1}^{\infty}(-1)^{n+1}\>\frac{m_{f}^{2}T}{n}e^{n\mu/T}K_{2}\left(\frac{nm_{f}}{T}\right)~{}.$
(114)

Using Eq. (112) in the expression above and retaining only the dominant $n=1$
term, we have

$N_{f}\ \simeq\
2\left(\frac{m_{f}T}{2\pi}\right)^{\frac{3}{2}}e^{-(m_{f}-\mu)/T}~{}.$ (115)

Thus, the medium-induced neutrino mass in the limit $T\ll\mu<m_{f}$ evaluates
to:

$\Delta m_{\nu,{\alpha\beta}}\ \simeq\
\frac{y_{f}y_{\alpha\beta}}{m_{\phi}^{2}}N_{f}\,,$ (116)

as given in Eq. (15).

###  A.3 Case 3: $\mu<m_{f}\ll T$

For $z\rightarrow 0$, the asymptotic form for $K_{\nu}\left(z\right)$ is:

$K_{\nu}\left(z\right)\ \simeq\
\frac{\Gamma(\nu)}{2}\left(\frac{z}{2}\right)^{-\nu}~{}.$ (117)

Using the above in Eq. (111), we can write the mass correction as:

$\displaystyle\Delta m_{\nu,{\alpha\beta}}$ $\displaystyle\ \simeq\
\frac{m_{f}y_{\alpha\beta}y_{f}}{\pi^{2}m_{\phi}^{2}}\sum_{n=1}^{\infty}(-1)^{n+1}\>\frac{T^{2}}{n^{2}}e^{n\mu/T}$
(118) $\displaystyle\ =\
-\frac{m_{f}y_{\alpha\beta}y_{f}T^{2}}{\pi^{2}m_{\phi}^{2}}\text{Li}_{2}(-e^{\mu/T})\,,$
(119)

where $\text{Li}_{\nu}(z)$ is the polylogarithm. In the case $|z|\rightarrow
0$, $\text{Li}_{n}(-e^{z})\simeq-(1-2^{1-n})\zeta(n)$. Using this one obtains:

$\Delta m_{\nu,{\alpha\beta}}\ \simeq\
\frac{y_{f}y_{\alpha\beta}m_{f}T^{2}}{12m_{\phi}^{2}}~{}.$ (120)

Again using Eq. (117) in Eq. (A.2) and retaining only the $n=1$ term we get:

$\displaystyle N_{f}\ \simeq\
-\frac{2T^{3}}{\pi^{2}}\text{Li}_{3}(-e^{\mu/T})\ =\
\frac{3T^{3}}{2\pi^{2}}\zeta(3)~{}.$ (121)

Thus, the scalar NSI expression for $\mu<m_{f}\ll T$ evaluates to:

$\Delta m_{\nu,{\alpha\beta}}\ \simeq\
\frac{y_{\alpha\beta}y_{f}m_{f}}{3\>m_{\phi}^{2}}\left(\frac{\pi^{2}N_{f}}{12\>\zeta(3)}\right)^{\frac{2}{3}}~{},$
(122)

as given in Eq. (17).

##  Appendix B Calculation of self-energy in neutrino background

Here we evaluate the neutrino self-energy arising from a neutrino background
as given in Eq. (31). We can rewrite the delta function in Eq. (31) as
follows:

$\displaystyle\delta\left[{\left(k+\frac{p}{2}\right)^{2}-m_{\phi}^{2}}\right]$
$\displaystyle\ =\
\frac{1}{|\mathbb{k}||\mathbb{p}|}\delta(\cos{\theta}-\cos{\theta_{0}})\,,$
(123)

where

$\cos{\theta_{0}}\ =\
\frac{k_{0}^{2}-|\mathbb{k}|^{2}+\frac{p^{2}}{4}-m_{\phi}^{2}+k_{0}p_{0}}{|\mathbb{k}||\mathbb{p}|}~{}.$
(124)

Using kinematical arguments and $|\cos{\theta_{0}}|\leq 1$, we find the range
for $k_{0}$ and $|\mathbb{k}|^{2}$:

$k_{0}\ :\
\left\\{\frac{-p_{0}}{2}+m_{\nu},\infty\right\\},~{}~{}~{}~{}|\mathbb{k}|^{2}\
:\ \left\\{|\mathbb{k}|^{2}_{-}\;,\>|\mathbb{k}|^{2}_{+}\right\\}$ (125)

where

$|\mathbb{k}|^{2}_{\pm}\ =\
\frac{1}{4}\bigg{(}|(\mathbb{p}|\pm\sqrt{|(\mathbb{p}|)^{2}+4k_{0}p_{0}-4m_{\nu}^{2}+4k_{0}^{2}+p^{2}}\bigg{)}^{2}~{}.$
(126)

Changing the integration variables to spherical coordinates and integrating
over $\cos{\theta}$ we obtain:

$\Sigma_{\alpha\beta}^{\nu}\ =\
-\frac{y_{\alpha\gamma}y_{\gamma\beta}}{16\pi^{2}|\textbf{p}|}\int_{k_{0}^{\textbf{min}}}^{k_{0}^{\textbf{max}}}dk_{0}\int_{|\mathbb{k}|^{2}_{-}}^{|\mathbb{k}|^{2}_{+}}d|\textbf{k}|^{2}\>\frac{(\not{k}+\frac{\not{p}}{2}+m_{\nu})}{k_{0}^{2}-|\textbf{k}|^{2}+\frac{p^{2}}{4}-\frac{m_{\phi}^{2}+m_{\nu}^{2}}{2}}\>n_{\nu}\left(k_{0}+\frac{p_{0}}{2}\right)~{}.$
(127)

This contribution can be decomposed as given in Eq. (22). By defining

$\mathcal{I}\ =\
\int_{m_{\nu}}^{\infty}dk_{0}\;n_{\nu}(k_{0})\ln\left[\frac{k_{0}p_{0}-p^{2}+\frac{m_{\phi}^{2}-m_{\nu}^{2}}{2}+|\textbf{p}|\sqrt{k_{0}^{2}-m_{\nu}^{2}}}{k_{0}p_{0}-p^{2}+\frac{m_{\phi}^{2}-m_{\nu}^{2}}{2}-|\textbf{p}|\sqrt{k_{0}^{2}-m_{\nu}^{2}}}\right]\,,$
(128)

the quantities $J_{u},J_{m},J_{p}$ in Eq. (24) can be written succinctly as:

$\displaystyle J_{m}\ $ $\displaystyle=\ -2m_{\nu}\>\mathcal{I}\,,$ (129)
$\displaystyle J_{p}\ $ $\displaystyle=\
-(p^{2}+m_{\nu}^{2}-m_{\phi}^{2})\>\mathcal{I}-2|\textbf{p}|\int_{m_{\nu}}^{\infty}dk_{0}\;n_{\nu}(k_{0})\sqrt{k_{0}^{2}-m_{\nu}^{2}}\,,$
(130) $\displaystyle J_{u}\ $ $\displaystyle=\
-2\int_{m_{\nu}}^{\infty}dk_{0}\;k_{0}\;n_{\nu}(k_{0})\ln\left[\frac{k_{0}p_{0}-p^{2}+\frac{m_{\phi}^{2}-m_{\nu}^{2}}{2}+|\textbf{p}|\sqrt{k_{0}^{2}-m_{\nu}^{2}}}{k_{0}p_{0}-p^{2}+\frac{m_{\phi}^{2}-m_{\nu}^{2}}{2}-|\textbf{p}|\sqrt{k_{0}^{2}-m_{\nu}^{2}}}\right]~{}.$
(131)

These integrals $(J_{m},\,J_{p},\,J_{u})$ cannot be evaluated analytically in
general. However, they may be evaluated in the high temperature limit. For
this purpose we set $m_{\nu}$ to zero and assume the chemical potential $\mu$
is small. This condition should be realized when the results are applied to
early universe. The integrals in this limit are evaluated to be:

$\displaystyle J_{m}$ $\displaystyle\ \simeq\ $ $\displaystyle-2m_{\nu}T\,{\rm
ln}2\,\,{\rm ln}\left(\frac{2\sqrt{2}|{\bf p}|T}{m_{\phi}^{2}}\right)\,,$
(132) $\displaystyle J_{p}$ $\displaystyle\ \simeq\ $
$\displaystyle\frac{\pi^{2}T^{2}|{\bf p}|}{3}+|{\bf p}|^{2}T\,{\rm
ln}2\,\,{\rm ln}\left(\frac{2\sqrt{2}|{\bf p}|T}{m_{\phi}^{2}}\right)\,,$
(133) $\displaystyle J_{u}$ $\displaystyle\ \simeq\ $
$\displaystyle\frac{\pi^{2}T^{2}}{6}\left(12\zeta^{\prime}(-1)+{\rm
ln}\left(\frac{16\pi|{\bf p}|T}{m_{\phi}^{2}}\right)\right)~{}.$ (134)

These results have been applied to derive the energy shift for neutrinos and
antineutrinos in Sec. 2.2, see Eq. (33).

##  Appendix C Calculation of thermal mass for the scalar field

Here we carry out the evaluation of the self-energy diagram of $\phi$ to
calculate its thermal mass. As shown in Sec. 5, $\phi$ can develop a medium-
dependent mass, which is given by Eq. (82). This contribution can be written
as:

$\mathcal{M}\ =\ \mathcal{M}_{1}+\mathcal{M}_{2}\,,$ (135)

where

$\displaystyle\mathcal{M}_{1}\ $ $\displaystyle=\
4y_{f}^{2}\int\frac{d^{4}p}{(2\pi)^{4}}\left(k^{2}-\frac{p^{2}}{4}+m_{f}^{2}\right)\frac{\Gamma_{f}(k+p/2)}{(k-p/2)^{2}-m_{f}^{2}}\,,$
(136) $\displaystyle\mathcal{M}_{2}\ $ $\displaystyle=\
4y_{f}^{2}\int\frac{d^{4}p}{(2\pi)^{4}}\left(k^{2}-\frac{p^{2}}{4}+m_{f}^{2}\right)\frac{\Gamma_{f}(k-p/2)}{(k+p/2)^{2}-m_{f}^{2}}~{}.$
(137)

Since $\mathcal{M}_{1}\rightarrow\mathcal{M}_{2}$ with the replacement
$p\rightarrow-p$, we will focus only on simplifying the expression for
$\mathcal{M}_{1}$.

$\mathcal{M}_{1}\ =\
4y_{f}^{2}\int_{\frac{-p_{0}}{2}}^{\infty}dk_{0}\int\frac{d^{3}p}{(2\pi)^{3}}\left(k^{2}-\frac{p^{2}}{4}+m_{f}^{2}\right)\frac{\delta((k+p/2)^{2}-m_{f}^{2})}{(k-p/2)^{2}-m_{f}^{2}}n_{f}\left(k_{0}+\frac{p_{0}}{2}\right)\,.$
(138)

The delta function can be written as

$\displaystyle\delta\left[{\left(k+\frac{p}{2}\right)^{2}-m_{f}^{2}}\right]$
$\displaystyle\ =\
\frac{1}{|\mathbb{k}||\mathbb{p}|}\delta(\cos{\theta}-\cos{\theta_{0}})\,,$
(139)

where

$\cos{\theta_{0}}\ =\
\frac{k_{0}^{2}-|\mathbb{k}|^{2}+\frac{p^{2}}{4}-m_{f}^{2}+k_{0}p_{0}}{|\mathbb{k}||\mathbb{p}|}~{}.$
(140)

Using kinematical arguments and $|\cos{\theta_{0}}|\leq 1$, we find the range
for $k_{0}$ and $|\mathbb{k}|^{2}$:

$k_{0}\ :\
\left\\{\frac{-p_{0}}{2}+m_{f},\infty\right\\},~{}~{}~{}|\mathbb{k}|^{2}\ :\
\left\\{|\mathbb{k}|^{2}_{-}\;,\>|\mathbb{k}|^{2}_{+}\right\\}\,,$ (141)

where

$|\mathbb{k}|^{2}_{\pm}\ =\
\frac{1}{4}\bigg{(}|(\mathbb{p}|\pm\sqrt{|(\mathbb{p}|)^{2}+4k_{0}p_{0}-4m_{f}^{2}+4k_{0}^{2}+p^{2}}\bigg{)}^{2}~{}.$
(142)

Thus, changing the integration variables to spherical coordinates and
integrating over $\cos{\theta}$ we get:

$\mathcal{M}_{1}\ =\
-\frac{y_{f}^{2}}{4\pi^{2}|\textbf{p}|}\int_{k_{0}^{\textbf{min}}}^{k_{0}^{\textbf{max}}}dk_{0}\int_{|\mathbb{k}|^{2}_{-}}^{|\mathbb{k}|^{2}_{+}}d|\textbf{k}|^{2}\>\frac{k_{0}^{2}-|\textbf{k}|^{2}-\frac{p^{2}}{4}+m_{f}^{2}}{k_{0}^{2}-|\textbf{k}|^{2}+\frac{p^{2}}{4}+m_{f}^{2}}\>n_{f}\left(k_{0}+\frac{p_{0}}{2}\right)~{}.$
(143)

Integrating the above integral with respect to $|\textbf{k}|^{2}$ and adding
the contribution from both $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ yields:

$\displaystyle\mathcal{M}$ $\displaystyle\ =\
\frac{y_{f}^{2}}{\pi^{2}}\int_{m_{f}}^{\infty}dk_{0}\>n_{f}\left(k_{0}\right)\sqrt{k_{0}^{2}-m_{f}^{2}}$
$\displaystyle-\frac{y_{f}^{2}}{2\pi^{2}|\textbf{p}|}\left(m_{f}^{2}-\frac{m_{\phi}^{2}}{4}\right)\int_{m_{f}}^{\infty}dk_{0}\>n_{f}\left(k_{0}\right)\>\ln\left(\frac{\left(|\textbf{p}|\sqrt{k_{0}^{2}-m_{f}^{2}}-\frac{m_{\phi}^{2}}{2}\right)^{2}-k_{0}^{2}p_{0}^{2}}{\left(|\textbf{p}|\sqrt{k_{0}^{2}-m_{f}^{2}}+\frac{m_{\phi}^{2}}{2}\right)^{2}-k_{0}^{2}p_{0}^{2}}\right)~{}.$
(144)

In the limit $m_{\phi}\rightarrow 0$, the mass correction for scalar reduces
to Eq. (83).

## References

  * (1) Particle Data Group collaboration, _Review of Particle Physics_ , _Phys. Rev._ D98 (2018) 030001.
  * (2) L. Wolfenstein, _Neutrino Oscillations in Matter_ , _Phys. Rev._ D17 (1978) 2369.
  * (3) K. S. Babu, P. S. B. Dev, S. Jana and A. Thapa, _Non-Standard Interactions in Radiative Neutrino Mass Models_ , 1907.09498.
  * (4) Y. Farzan, _A model for large non-standard interactions of neutrinos leading to the LMA-Dark solution_ , _Phys. Lett._ B748 (2015) 311 [1505.06906].
  * (5) Y. Farzan and I. M. Shoemaker, _Lepton Flavor Violating Non-Standard Interactions via Light Mediators_ , _JHEP_ 07 (2016) 033 [1512.09147].
  * (6) K. S. Babu, A. Friedland, P. A. N. Machado and I. Mocioiu, _Flavor Gauge Models Below the Fermi Scale_ , _JHEP_ 12 (2017) 096 [1705.01822].
  * (7) P. S. B. Dev et al., _Neutrino Non-Standard Interactions: A Status Report_ , 1907.00991.
  * (8) P. Coloma, I. Esteban, M. C. Gonzalez-Garcia and M. Maltoni, _Improved global fit to Non-Standard neutrino Interactions using COHERENT energy and timing data_ , 1911.09109.
  * (9) S.-F. Ge and S. J. Parke, _Scalar Nonstandard Interactions in Neutrino Oscillation_ , _Phys. Rev. Lett._ 122 (2019) 211801 [1812.08376].
  * (10) A. K. Das, _Finite Temperature Field Theory_. World Scientific, New York, 1997.
  * (11) D. Notzold and G. Raffelt, _Neutrino Dispersion at Finite Temperature and Density_ , _Nucl. Phys._ B307 (1988) 924.
  * (12) P. B. Pal and T. N. Pham, _A Field Theoretic Derivation of Wolfenstein’s Matter Oscillation Formula_ , _Phys. Rev._ D40 (1989) 259.
  * (13) P. Langacker and J. Liu, _Standard Model contributions to the neutrino index of refraction in the early universe_ , _Phys. Rev._ D46 (1992) 4140 [hep-ph/9206209].
  * (14) M. B. Wise and Y. Zhang, _Lepton Flavorful Fifth Force and Depth-dependent Neutrino Matter Interactions_ , _JHEP_ 06 (2018) 053 [1803.00591].
  * (15) A. Y. Smirnov and X.-J. Xu, _Wolfenstein potentials for neutrinos induced by ultra-light mediators_ , _Journal of High Energy Physics_ 2019 (2019) .
  * (16) M. E. Peskin and D. V. Schroeder, _An Introduction to quantum field theory_. Addison-Wesley, Reading, USA, 1995.
  * (17) J. Liu, C. E. M. Wagner and X.-P. Wang, _A light complex scalar for the electron and muon anomalous magnetic moments_ , _JHEP_ 03 (2019) 008 [1810.11028].
  * (18) J. K. Hoskins, R. D. Newman, R. Spero and J. Schultz, _Experimental tests of the gravitational inverse-square law for mass separations from 2 to 105 cm_ , _Phys. Rev. D_ 32 (1985) 3084.
  * (19) J. C. Long, H. W. Chan, A. B. Churnside, E. A. Gulbis, M. C. M. Varney and J. C. Price, _Upper limits to submillimeter-range forces from extra space-time dimensions_ , hep-ph/0210004.
  * (20) D. J. Kapner, T. S. Cook, E. G. Adelberger, J. H. Gundlach, B. R. Heckel, C. D. Hoyle et al., _Tests of the gravitational inverse-square law below the dark-energy length scale_ , _Phys. Rev. Lett._ 98 (2007) 021101 [hep-ph/0611184].
  * (21) A. A. Geraci, S. J. Smullin, D. M. Weld, J. Chiaverini and A. Kapitulnik, _Improved constraints on non-Newtonian forces at 10 microns_ , _Phys. Rev._ D78 (2008) 022002 [0802.2350].
  * (22) S.-Q. Yang, B.-F. Zhan, Q.-L. Wang, C.-G. Shao, L.-C. Tu, W.-H. Tan et al., _Test of the gravitational inverse square law at millimeter ranges_ , _Phys. Rev. Lett._ 108 (2012) 081101.
  * (23) W.-H. Tan, S.-Q. Yang, C.-G. Shao, J. Li, A.-B. Du, B.-F. Zhan et al., _New test of the gravitational inverse-square law at the submillimeter range with dual modulation and compensation_ , _Phys. Rev. Lett._ 116 (2016) 131101.
  * (24) E. Adelberger, J. Gundlach, B. Heckel, S. Hoedl and S. Schlamminger, _Torsion balance experiments: A low-energy frontier of particle physics_ , _Progress in Particle and Nuclear Physics_ 62 (2009) 102 .
  * (25) E. Hardy and R. Lasenby, _Stellar cooling bounds on new light particles: plasma mixing effects_ , _JHEP_ 02 (2017) 033 [1611.05852].
  * (26) S. Knapen, T. Lin and K. M. Zurek, _Light Dark Matter: Models and Constraints_ , _Phys. Rev._ D96 (2017) 115021 [1709.07882].
  * (27) G. G. Raffelt, _Stars as laboratories for fundamental physics_. 1996.
  * (28) N. Ishizuka and M. Yoshimura, _Axion and Dilaton Emissivity From Nascent Neutron Stars_ , _Prog. Theor. Phys._ 84 (1990) 233.
  * (29) B. Batell, A. Freitas, A. Ismail and D. Mckeen, _Probing Light Dark Matter with a Hadrophilic Scalar Mediator_ , _Phys. Rev._ D100 (2019) 095020 [1812.05103].
  * (30) C. Alexandrou et al., _Nucleon scalar and tensor charges using lattice QCD simulations at the physical value of the pion mass_ , _Phys. Rev._ D95 (2017) 114514 [1703.08788].
  * (31) Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , 1807.06209.
  * (32) F. Forastieri, M. Lattanzi and P. Natoli, _Cosmological constraints on neutrino self-interactions with a light mediator_ , _Phys. Rev._ D100 (2019) 103526 [1904.07810].
  * (33) C. D. Kreisch, F.-Y. Cyr-Racine and O. Doré, _The Neutrino Puzzle: Anomalies, Interactions, and Cosmological Tensions_ , 1902.00534.
  * (34) J. H. Hubbell and Overboi, _Relativistic atomic form factors and photon coherent scattering cross sections_ , _Journal of Physical and Chemical Reference Data_ 8 (1979) 69.
  * (35) A. M. Dziewonski and D. L. Anderson, _Preliminary reference Earth model_ , _Physics of the Earth and Planetary Interiors_ 25 (1981) 297.
  * (36) “Standard solar model.” http://www.sns.ias.edu/~jnb/SNdata/solarmodels.html.
  * (37) J. N. Bahcall, S. Basu and M. H. Pinsonneault, _How uncertain are solar neutrino predictions?_ , _Phys. Lett._ B433 (1998) 1 [astro-ph/9805135].
  * (38) J. N. Bahcall, A. M. Serenelli and S. Basu, _New solar opacities, abundances, helioseismology, and neutrino fluxes_ , _The Astrophysical Journal_ 621 (2005) L85–L88.
  * (39) N. Arkani-Hamed, L. Motl, A. Nicolis and C. Vafa, _The String landscape, black holes and gravity as the weakest force_ , _JHEP_ 06 (2007) 060 [hep-th/0601001].
  * (40) J. H. Chang, R. Essig and S. D. McDermott, _Revisiting supernova 1987a constraints on dark photons_ , _Journal of High Energy Physics_ 2017 (2017) .