text:
\section{Introduction}

Until popularized by Morris and Thorne \cite{1} the idea of spacetime
wormhole introduced in the 1930s by Einstein and Rosen \cite{2} was
considered non-physical and largely as a fantasy. Although all kinds of
spherically / cylindrically symmetric metrics known to date were tried the
wormhole concept was shadowed by the required negative total energy. While
it was easy to resort to quantum field theoretical negative energy for
remedy the absence of large scale quantum systems persisted as another
serious handicap. For such reasons relying on classical physics and
searching for support within this context seems indispensable. Even to
minimize the negative (exotic) energy the idea of thin-shell wormhole (TSW)
was developed (See \cite{3} and those cited therein). By construction in all
these studies the otherwise non-traversable wormhole throat that connects
two different universes has circular topology.

In this study we add oblateness as a new degree of freedom represented by
the parameter $\delta $ ($-\infty <\delta <\infty $) and show that for
certain range of $\delta $ total energy becomes positive to avoid exotic
sources. This happens in the Zipoy-Voorhees-Weyl (ZVW) spacetime \cite{4}
with a quadrupole moment $Q=\frac{1}{3}M^{3}\delta \left( 1-\delta
^{2}\right) $, where $M=$mass of the bumpy object (or the black hole).
Naturally for $\delta =1$ one recovers the spherical Schwarzschild geometry.
It should be supplemented that integrability and chaotic behaviours of the
ZVW spacetime still is not well understood \cite{5}. Asymptotically flat,
rotating ZVW metric was discovered by Tomimatsu and Sato (TS) \cite{6},
which similar to its static predecessor remains from physics stand point yet
unclear. Once the problem of gravitational wave detection is overcome we
expect that non-Kerr (i.e., $\delta \neq 1$) multipoles of the entire TS
family can be detected. We note also that non-asymptotically flat extension
of ZVW metric is also available whose physics is yet to be understood \cit
{7}. Herein we concentrate on the static ZVW metric in general relativity
and construct thin-shell wormhole (TSW) in this spacetime.

We should add that in the context of the modified theories of gravity,
previously, there were some attempts to introduce thin-shell wormhole
supported by positive / normal matter \cite{8}. From this token it was
realized that the normal matter is possible only for the exotic branch
solution of the Einstein-Gauss-Bonnet field equation \cite{8}.

The paper is organized as follows. Construction of TSW in ZVW spacetime is
carried out in Section II. Integration of the total energy is achieved in
Section III. Stability analysis follows in Section IV and Conclusion in
Section V completes the paper.

\section{ZVW thin-shell wormhole (TSW)}

The two parameter ZVW spacetime in the prolate spheroidal coordinates is
described by the line elemen
\begin{multline}
ds^{2}=-A\left( x\right) dt^{2}+ \\
B\left( x,y\right) dx^{2}+C\left( x,y\right) dy^{2}+F\left( x,y\right)
d\varphi ^{2}
\end{multline
in which
\begin{eqnarray}
A &=&\left( \frac{x-1}{x+1}\right) ^{\delta },\text{ \ }B=\frac{k^{2}}{A
\left( \frac{x^{2}-1}{x^{2}-y^{2}}\right) ^{\left( \delta ^{2}-1\right) } \\
C &=&B\left( \frac{x^{2}-1}{1-y^{2}}\right) ,\text{ \ }F=\frac{k^{2}}{A
\left( x^{2}-1\right) \left( 1-y^{2}\right) ,  \notag
\end{eqnarray
where $k=\frac{M}{\delta }$ and the range of coordinates are $1<x$, $0\leq
y^{2}\leq 1$, $\varphi \in \left[ 0,2\pi \right] $ and $-\infty <t<\infty $.
We note that $-\infty <\delta <\infty $ such that $\delta =0$ corresponds to
a flat spacetime and with $\delta =1$ one finds the Schwarzschild black hole
solution with the horizon located at $x=1$. For the case $\delta \neq 1$ the
hypersurface $x=1$ is a true curvature singularity (naked singularity) \cit
{9}. As we shall see, the asymptotic behaviour of the ZV spacetime for
x\rightarrow \infty $ and $\delta >1$ is of our interest. It can be seen
from (1) that, in the limit $x\rightarrow \infty ,$ it become
\begin{equation}
ds^{2}=-dt^{2}+k^{2}dx^{2}+k^{2}x^{2}\left( \frac{dy^{2}}{1-y^{2}}+\left(
1-y^{2}\right) d\varphi ^{2}\right)
\end{equation
which after a redefinition of $kx=r$ and $y=\cos \theta ,$ the line element
become
\begin{equation}
ds^{2}=-dt^{2}+dr^{2}+r^{2}\left( d\theta ^{2}+\sin ^{2}\theta d\varphi
^{2}\right)
\end{equation
which is flat.

Construction of thin-shell wormhole (TSW) follows the standard procedure of
cutting and pasting \cite{8}. We consider two copies of ZVW spacetime and we
remove from each
\begin{equation}
\mathcal{M}^{\pm }=\left\{ x^{\pm }<a,\text{ }1<a\right\}
\end{equation
in which $a=$ constant outside the singularities / horizons. We should
comment here that the minimality conditions of Hochberg and Visser \cite{10}
which is also known as \textit{the generalized flare-out conditions} in
static wormholes do not apply in the present case of TSW \cite{11}. At the
throat two spacetimes are identified to make a complete manifold. We
introduce next the induced coordinates $\xi ^{i}=\left( \tau ,y,\phi \right)
$ on the wormhole's throat with its proper time $\tau $. The two coordinates
are related b
\begin{equation}
g_{ij}=\frac{\partial x^{\alpha }}{\partial \xi ^{i}}\frac{\partial x^{\beta
}}{\partial \xi ^{j}}g_{\alpha \beta }
\end{equation
so that induced metric on the throat $\Sigma $ reads
\begin{equation}
g_{ij}=diag\left[ -1,C\left( a\left( \tau \right) ,y\right) ,F\left( a\left(
\tau \right) ,y\right) \right] .
\end{equation
The Israel junction conditions \cite{12} on $\Sigma $ take the form ($c=8\pi
G=1$)
\begin{equation}
\left\langle K_{i}^{j}\right\rangle -\left\langle K\right\rangle \delta
_{i}^{j}=-S_{i}^{j},
\end{equation
in which $\left\langle .\right\rangle $ stands for a jump across the
hypersurface. $K_{i}^{j}$ is the extrinsic curvature defined by
\begin{equation}
K_{ij}^{\left( \pm \right) }=-n_{\gamma }^{\left( \pm \right) }\left( \frac
\partial ^{2}x^{\gamma }}{\partial \xi ^{i}\partial \xi ^{j}}+\Gamma
_{\alpha \beta }^{\gamma }\frac{\partial x^{\alpha }}{\partial \xi ^{i}
\frac{\partial x^{\beta }}{\partial \xi ^{j}}\right) _{\Sigma }
\end{equation
with the normal unit vector
\begin{equation}
n_{\gamma }^{\left( \pm \right) }=\left( \pm \left\vert g^{\alpha \beta
\frac{\partial \mathcal{H}}{\partial x^{\alpha }}\frac{\partial \mathcal{H}}
\partial x^{\beta }}\right\vert ^{-1/2}\frac{\partial \mathcal{H}}{\partial
x^{\gamma }}\right) .
\end{equation
Note that $\left\langle K\right\rangle =Trace\left\langle
K_{i}^{j}\right\rangle $ and $S_{i}^{j}=$diag$\left( -\sigma ,P_{y},P_{\phi
}\right) $ is the energy momentum tensor on the thin-shell. The parametric
equation of the hypersurface $\Sigma $ is given b
\begin{equation}
\mathcal{H}\left( x,a\left( \tau \right) \right) =x-a\left( \tau \right) =0.
\end{equation
The normal unit vectors to $\mathcal{M}_{\pm }$ is found to b
\begin{equation}
n_{\gamma }^{\left( \pm \right) }=\pm \left( -\sqrt{AB}\dot{a},B\sqrt{\Delta
},0,0\right) _{\Sigma }
\end{equation
with $\Delta =\frac{1}{B}+\dot{a}^{2}$ and $\dot{a}=\frac{da}{d\tau }$. The
resulting extrinsic curvature components are
\begin{eqnarray}
K_{\tau }^{\tau \left( \pm \right) } &=&\pm \frac{\ddot{a}+\left( \frac{B_{a
}{B}+\frac{A_{a}}{A}\right) \frac{\dot{a}^{2}}{2}+\frac{A_{a}}{2AB}}{\sqrt
\Delta }} \\
K_{y}^{y\left( \pm \right) } &=&\pm \frac{C_{a}}{2C}\sqrt{\Delta },  \notag
\\
K_{\varphi }^{\varphi \left( \pm \right) } &=&\pm \frac{F_{a}}{2F}\sqrt
\Delta },  \notag
\end{eqnarray
in which a subscript $a$ stands for $\frac{\partial }{\partial a}$. The
surface energy momentum tensor has components defined as
\begin{equation}
\sigma =-\left( \frac{C_{a}}{C}+\frac{F_{a}}{F}\right) \sqrt{\Delta }
\end{equation
\begin{equation}
P_{y}=\frac{2\ddot{a}+\left( \frac{B_{a}}{B}+\frac{A_{a}}{A}\right) \dot{a
^{2}+\frac{A_{a}}{AB}}{\sqrt{\Delta }}+\frac{\sqrt{\Delta }F_{a}}{F}
\end{equation
\begin{equation}
P_{\varphi }=\frac{2\ddot{a}+\left( \frac{B_{a}}{B}+\frac{A_{a}}{A}\right)
\dot{a}^{2}+\frac{A_{a}}{AB}}{\sqrt{\Delta }}+\frac{\sqrt{\Delta }C_{a}}{C}.
\end{equation}

\begin{figure}[tbp]
\includegraphics[width=80mm,scale=0.7]{Figure01}
\captionsetup{justification=raggedright, singlelinecheck=false}
\caption{A $3D$ plot of the positive part of $\protect\sigma _{0}$ in terms
of $a_{0}$ and $y$ with $\protect\delta =2.0$. We see that although $\protec
\sigma _{0}$ gets positive values for some interval of $y$ but it is not
positive everywhere on $y$. When the value of $\protect\delta $ decreases
the interval of $y$ on which $\protect\delta $ is positive gets smaller and
ultimately for $\protect\delta \leq 1$ the interval disappears so that
\protect\sigma _{0}$ gets only negative values. (Note that $\protect\sigma
_{0}$ is an even function with respect to $y$ and only a section has been
plotted.)}
\end{figure}
\begin{figure}[tbp]
\includegraphics[width=80mm,scale=0.7]{Figure02}
\captionsetup{justification=raggedright, singlelinecheck=false}
\caption{A $3D$ plot of the positive amount of energy $\frac{\Omega }{2}$
versus $a$ and $\protect\delta $ with constant mass parameter i.e. $M=1.$
This figure shows clearly that for $\protect\delta >2$ there exist $a_{c}$
in which with $1<a<a_{c}$ the total energy is positive and therefore the
resultant thin shell wormhole is supported by ordinary / normal matter. It
should be added that for a given $\protect\delta >2$, there exists
1<a<a_{c} $ which leads to a physically accepted TSW.}
\end{figure}

\begin{figure}[tbp]
\includegraphics[width=80mm,scale=0.7]{Figure03}
\captionsetup{justification=raggedright, singlelinecheck=false}
\caption{$\frac{\Omega }{2}$ versus $\protect\delta $ for different values
of $a=1.0001,1.0010,1.0100,1.1000,1.5000$ and $M=1.$ As we commented in Fig.
2, here it is more clear that when the value of $a\rightarrow 1$, the
threshold $\protect\delta $ which admits a positive total energy approaches
to $\protect\delta =2.$ This means that $\protect\delta =2$ is a critical
value to have a thin shell wormhole supported by ordinary matter. For larger
$a$ the threshold $\protect\delta $ gets larger values.}
\end{figure}

\section{Positive Matter Sources}

The energy-momentum components at equilibrium condition, i.e. $a=a_{0}=
constant with $\dot{a}=\ddot{a}=0$ yield
\begin{eqnarray}
\sigma _{0} &=&\left. -\left( \frac{C_{a}}{C}+\frac{F_{a}}{F}\right) \frac{
}{\sqrt{B}}\right\vert _{a=a_{0}}  \notag \\
P_{y0} &=&\left. \left( \frac{A_{a}}{A}+\frac{F_{a}}{F}\right) \frac{1}
\sqrt{B}}\right\vert _{a=a_{0}}, \\
P_{\phi 0} &=&\left. \left( \frac{A_{a}}{A}+\frac{C_{a}}{C}\right) \frac{1}
\sqrt{B}}\right\vert _{a=a_{0}}.  \notag
\end{eqnarray
The explicit form of the energy momentum components are then found to b
\begin{eqnarray}
\sigma _{0} &=&\frac{2\left[ 2\delta \left( a_{0}^{2}-y^{2}\right)
-2a_{0}^{3}+a_{0}\left( 1+y^{2}\right) +a_{0}\delta ^{2}\left(
y^{2}-1\right) \right] }{\left( a_{0}^{2}-1\right) \left(
a_{0}^{2}-y^{2}\right) },  \notag \\
P_{y0} &=&\frac{2a_{0}}{a_{0}^{2}-1}, \\
P_{\phi 0} &=&\frac{2a_{0}\left[ a_{0}^{2}-1+\left( 1-y^{2}\right) \delta
^{2}\right] }{\left( a_{0}^{2}-1\right) \left( a_{0}^{2}-y^{2}\right) }.
\notag
\end{eqnarray
Note that since $1<a_{0}$ there are no singularities in the foregoing
expressions. From the conditions $a_{0}>1$ and $y<1$ one observes that
P_{y0}$ and $P_{\phi 0}$ are both positive while $\sigma _{0}$ may be
positive, negative or zero. Fig. 1 displays $\sigma _{0}$ in terms of $a_{0}$
and $y$. We observe that for $\delta >1$ there exists regions that $\sigma
_{0}$ becomes positive. This is seen clearly in Fig. 1. Also in the interval
on which $\sigma _{0}\geq 0$ the weak and strong energy conditions are
satisfied.

We note that, in Fig. 1 the energy density $\sigma _{0}$ is shown in terms
of $y$ and $a_{0}$ but only $y$ is variable and $a_{0}$ is fixed for a
specific TSW. This means that at the throat $x=a_{0}$ and only $y$ and $\phi
$ are variable. As one sees from (18), there is angular symmetry and as a
result, not only $\sigma _{0}$ but also $P_{y0}$ and $P_{\phi 0}$ are only
functions of $y.$ In addition, $P_{y0}$ is a constant function of $y$.
Therefore once we set the radius of the throat i.e., $a_{0},$ our energy
momentum tensor's components are left with the only variable $y.$ Hence,
depending on $y$, the energy density of the TSW i.e. $\sigma _{0}$ is
locally positive or negative. This is what we see in Fig. 1. Of course, the
situation is completely different for $P_{y0}$ and $P_{\phi 0}$ which are
positive for entire domain of $y.$

In addition to the energy conditions we are mainly interested in the total
energy supporting the TSW given by
\begin{equation}
\Omega =2\int\nolimits_{0}^{2\pi
}\int\nolimits_{0}^{1}\int\nolimits_{1}^{\infty }\sigma _{0}\mathcal{\delta
\left( x-a_{0}\right) \sqrt{-g}dxdyd\phi ,
\end{equation
which simplifies to
\begin{equation}
\Omega =4\pi \int\nolimits_{0}^{1}\sigma _{0}\sqrt{-g_{0}}dy.
\end{equation
In Eq. (19) $\mathrm{\delta }\left( x-a_{0}\right) $ is the Dirac delta
function. In Fig. 2 we plot $\Omega $ versus $a$ and $\delta $ with fixed
value of mass $M=1.$ Fig. 3 reveals more details. These plots overall show
that TSW\ supported by normal matter is possible provided $\delta >2$. That
explains also why given ordinary matter alone in Schwarzschild spacetime
with $\delta =1,$ there was no such traversable wormhole. In this regard let
us add that even an arbitrarily small energy condition violation is
considered respectable \cite{13}.

\section{Stability Analysis}

In this section we apply the small velocity perturbation on the shell under
the condition that the EoS of the TSW after the perturbation is same as its
EoS at its static equilibrium. This is possible if the perturbation process
occurs slow enough in which all the intermediate states can also be
considered as static equilibrium points. Therefore the EoS after the
perturbation read
\begin{equation}
\frac{P_{y}}{\sigma }=-\frac{\frac{A_{a}}{A}+\frac{F_{a}}{F}}{\frac{C_{a}}{C
+\frac{F_{a}}{F}}
\end{equation
an
\begin{equation}
\frac{P_{\varphi }}{\sigma }=-\frac{\frac{A_{a}}{A}+\frac{C_{a}}{C}}{\frac
C_{a}}{C}+\frac{F_{a}}{F}}.
\end{equation
These in explicit form amount t
\begin{equation}
2\ddot{a}+\left( \frac{B_{a}}{B}\right) \dot{a}^{2}=0
\end{equation
which upon integration yield
\begin{equation}
\dot{a}=\dot{a}_{0}\sqrt{\frac{B_{0}}{B}}=\dot{a}_{0}\sqrt{\frac{\left(
\frac{a_{0}+1}{a_{0}-1}\right) ^{\delta }\left( \frac{a_{0}^{2}-1}
a_{0}^{2}-y^{2}}\right) ^{\left( \delta ^{2}-1\right) }}{\left( \frac{a+1}
a-1}\right) ^{\delta }\left( \frac{a^{2}-1}{a^{2}-y^{2}}\right) ^{\left(
\delta ^{2}-1\right) }}}.\text{ \ }
\end{equation
In Figs. 4 and 5 we plot $\dot{a}$ in terms of $a$ and $y$ for the initial
values $a_{0}=1.2$ and $\dot{a}_{0}=\pm 0.1,$ respectively, and $\delta =2.5
. As one can see in both cases the velocity does not get zero, which means
that the throat does not go back to its initial position.

\begin{figure}[tbp]
\includegraphics[width=80mm,scale=0.7]{Figure04}
\captionsetup{justification=raggedright, singlelinecheck=false}
\caption{$\dot{a}$ versus $a$ and $y$ with $\protect\delta =2.5,$ $\dot{a
_{0}=0.1$ and $a_{0}=1.2.$ The velocity never vanishes which is an
indication of instability of the throat under the small velocity
perturbation.}
\end{figure}
\begin{figure}[tbp]
\includegraphics[width=80mm,scale=0.7]{Figure05}
\captionsetup{justification=raggedright, singlelinecheck=false}
\caption{A $3D$ plot of $\dot{a}$ versus $a$ and $y$ with $\protect\delta
=2.5,$ $\dot{a}_{0}=-0.1$ and $a_{0}=1.2.$ The velocity fails to vanish
which amounts to instability of the throat under the small velocity
perturbation.}
\end{figure}

\section{Conclusion}

We considered the possibility of having physical wormhole solution in
Einstein's general relativity which is supported by normal matter and the
energy conditions are satisfied. Our emphasis is on Einstein's general
relativity instead the modified theories such as the Lovelock theory. In
such theories there were attempts to find such physical wormholes \cite{8}
with non-physical solution and it was shown that a TSW supported by positive
energy was possible. In this context we have shown that TSWs with oblate
sources can be employed to admit an overall physical (i.e. non-exotic)
matter even in Einstein's general relativity. As we have depicted in
figures, for $\delta >2$ not only the total energy of the wormhole is
positive but also the WECs and SECs are satisfied for limited $y$-interval,
which increases for larger $\delta .$ Certain range of deviation from
spherical symmetry can be chosen from the total energy integral to render
this possible. Locally it can easily be checked from Eq. (16) for $y=0$, for
instance, that we have a negative energy density, however, this is
compensated within the total energy expression. It is also expected that
once the metric and the surface energy-momentum become time dependent energy
conservation on the thin shell will not be valid any more. Another important
aspect concerning TSWs which has not been discussed here is their stability
against perturbations. It is shown that small velocity perturbations in the
x-$direction leads to an unstable wormhole throat. Finally, from curiosity
we wish to ask: \textit{does the deformation parameter }$\delta $\textit{\
saves wormholes other than TSWs in Einstein's theory?}. It remains to be
seen.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

Since $^7$Li is a fragile nucleus which is easily destroyed by the proton capture reaction, its surface abundance reflects detailed stellar structure. In low-mass giants, stellar models predict surface lithium depletion \citep{1967ApJ...147..624I}. However, spectroscopic surveys have shown that $\sim1\%$ of giant stars have the lithium abundance as high as $A(\mathrm{Li})=\log(\mathrm{Li/H})+12>1.5$ \citep[e.g.][]{2016MNRAS.461.3336C,2018NatAs...2..790Y,2018A&A...617A...4S,2019MNRAS.484.2000D}. This is a long-standing problem in our understanding of low-mass stars \citep{1982ApJ...255..577W,1989ApJS...71..293B}.

Stars in the red giant (RG) branch and the red clump (RC) have the similar luminosity and the effective temperature, so the boundary between them is ambiguous in the Hertzsprung-Russell diagram. Some authors have suggested that a  part of the lithium-rich giants are RC stars \citep{2014ApJ...784L..16S,2014A&A...564L...6M}. Recent works  \citep{2019MNRAS.482.3822S,2019ApJ...878L..21S,2020NatAs.tmp..139K}  distinguished RC stars from RGs in data of spectroscopic surveys with the help of asteroseismological data \citep{2011Natur.471..608B,2016A&A...588A..87V}. They concluded that all of RC stars have the  lithium abundances of $A(\mathrm{Li})>-0.9$, which are higher than the predicted values by stellar models. This implies that a ubiquitous process produces $^7$Li during or before central helium burning.

The mechanism of this lithium enhancement is under debate. Some authors suggest engulfment of substellar objects which keep high lithium abundances \citep[e.g.][]{1999MNRAS.308.1133S,2012A&A...538A..36L,2016ApJ...833L..24A}. Others discuss in situ production by the Cameron-Fowler (CF) mechanism \citep{1971ApJ...164..111C}. In the hydrogen burning shell, $^7$Be is produced via the $^3$He$(\alpha,\;\gamma)^7$Be reaction. The produced $^7$Be is conveyed to the stellar surface and decays to $^7$Li by the electron capture. In the standard model, the CF mechanism is insufficient to reproduce the abundance in lithium-rich giants. However, \citet{2016MNRAS.461.3336C} point out that extra mixing induced by the tidal interaction with a binary companion can drive the lithium production. Also,  a recent study \citep{2020arXiv200901248S} focused on mixing induced by the helium flash and claimed that it can naturally solve the lithium problem in RC stars proposed by \citet{2020NatAs.tmp..139K}.

In order to explain the ubiquitous enhancement of lithium in RC stars, we introduce the additional energy losses induced by physics beyond the Standard Model. In this {paper}, as an example we focus on the neutrino magnetic moment (NMM), which is denoted as $\mu_\nu$. The existence of neutrino masses was established by the detection of the neutrino oscillations, starting with atmospheric neutrinos  \citep{1998PhRvL..81.1562F} and later by solar and reactor neutrinos. In the Standard Model of particle physics  massive neutrinos have magnetic moments which are too small to be detected by present and near-term future experiments  \citep{1980PhRvL..45..963F,1982NuPhB.206..359S,2015RvMP...87..531G,2018ARNPS..68..313B}. The current best experimental limit $\mu_\nu<2.9\times10^{-11}\mu_\mathrm{B}$, where $\mu_\mathrm{B}$ is the Bohr magneton, comes from the GEMMA experiment \citep{2013PPNL...10..139B}, which measured the scattering cross sections of target electrons and reactor anti-electron neutrinos.

The NMM induces the additional energy loss and affects various stages of
stellar evolution for wider mass-range of stars.  The effect on the
evolution of intermediate-mass stars has recently been studied in detail
\citep{2020arXiv200808393M} and it was found that in the presence of a sufficiently large NMM the duration of blue giants is
shorter and the blue loops are eliminated.  Stellar plasma of low-mass
stars also is effected by NMM and the helium flash delays \citep{1994ApJ...425..222H}. As a result, a heavier inert helium core is formed and the luminosity of the tip of the RG branch (TRGB) increases \citep[e.g.][]{1996slfp.book.....R}. This enables one to use low-mass stars in globular clusters to give tighter constraints of $\mu_\nu<1.5\times10^{-12}\mu_\mathrm{B}$ \citep{2020arXiv200703694C} and $\mu_\nu<2.2\times10^{-12}\mu_\mathrm{B}$ \citep{2015APh....70....1A}. Also, the delayed helium flash may result in activation of the CF mechanism induced by thermohaline mixing \citep{1999ApJ...510..217S,2015MNRAS.446.2673L}.
In this paper we show that NMM needed to enhance $A(\mathrm{Li})$ in RC stars and  reduce the discrepancy between the observations and the theory is slightly higher than current astrophysical limits, so other  channels may be in the play.

Section 2 describes the stellar models and the treatment of the NMM. Section 3 shows the results of our  calculations and compares them with the observational data. In Section 4, we summarize our results and discuss the future perspective.

\section{Methods}

We use Modules for Experiments in Stellar Astrophysics \citep[\texttt{MESA};][]{MESA1,MESA2,MESA3,MESA4,MESA5} version 10398 to construct one-dimensional low-mass stellar models. \texttt{MESA} adopts the equation of state of \citet{Rogers2002} and \citet{Timmes2000} and the opacity of \cite{Iglesias1996,Iglesias1993} and \cite{Ferguson2005}. The electron conductivity is from \citep{2007ApJ...661.1094C}.  We adopt nuclear reaction rates compiled by  NACRE \citep{1999NuPhA.656....3A} and \citet{1988ADNDT..40..283C}. If a reaction rate appears in both, the one tabulated in NACRE is adopted. {In particular}, the rates for the triple-$\alpha$ reaction and $^{14}$N$(p,\;\gamma)^{15}$O are from NACRE.  The adopted nuclear reaction network is \texttt{pp\_and\_cno\_extras.net}, which includes $^{1,2}$H, $^{3,4}$He, $^7$Li, $^7$Be, $^8$B, $^{12}$C, $^{14}$N, $^{14,16}$O, $^{19}$F, $^{18,19,20}$Ne, and $^{22,24}$Mg.  Treatment of  electron screening is based on for the strong regime and \citet{1973ApJ...181..457G} for the weak regime.  The mass loss formula in \citet{1975MSRSL...8..369R} with the scaling factor $\eta=0.3$ is adopted.

\texttt{MESA} makes use of the mixing length theory \citep{1968pss..book.....C} to calculate the convective luminosity from the temperature gradient. We adopt the Ledoux criterion \citep{1947ApJ...105..305L} to calculate the convective instability.  Our model considers thermohaline mixing as well because it affects the lithium abundance \citep{1999ApJ...510..217S,2015MNRAS.446.2673L}. Thermohaline mixing is treated as a diffusive process with the diffusion coefficient \citep{MESA2}
\begin{equation}
    D_\mathrm{thm}=\alpha_\mathrm{thm}\frac{3K}{2\rho C_\mathrm{P}}\frac{B}{\nabla_\mathrm{T}-\nabla_\mathrm{ad}},
\end{equation}
where $\alpha_\mathrm{thm}$ is a free parameter, $\rho$ is the density, $C_\mathrm{P}$ is the specific heat, $B$ is the Ledoux term \citep{1989nos..book.....U}, $\nabla_\mathrm{T}$ is the actual temperature gradient, and $\nabla_\mathrm{ad}$ is the adiabatic temperature gradient. $K$ is the thermal conductivity written as
\begin{equation}
        K=\frac{4acT^3}{3\kappa\rho},\label{K}
\end{equation}
where $a$ is the radiation density constant, $c$ is the speed of light, $T$ is the temperature, and $\kappa$ is the opacity.

The parameters in our models follow those in \citet{2020NatAs.tmp..139K}. The initial mass is fixed to $1M_{\odot}$ and the initial metallicity is fixed to be solar: $Z=0.0148$ \citep{2019arXiv191200844L}. However, the initial lithium abundance in the pre-main sequence is set to $A(\mathrm{Li})=2.8$. Although there is no observational evidence, this value is used in previous works \citep{2020NatAs.tmp..139K,2020arXiv200901248S}. We follow their prescription to compare our result with theirs.  The mixing length is $\alpha=1.6$ and the thermohaline coefficient is $\alpha_\mathrm{thm}=50$, 100 and 1000, because $A(\mathrm{Li})$ after the RG branch bump is sensitive to thermohaline mixing.

We consider the plasmon decay and neutrino pair production as the additional energy loss induced by the NMM. The energy loss rate due to the plasmon decay is given by \citep{1994ApJ...425..222H,2009ApJ...696..608H}
\begin{equation}
\epsilon^\mu_\mathrm{plas}=0.318\left(\frac{\omega_\mathrm{pl}}{10\;\mathrm{keV}}\right)^{-2}\left(\frac{\mu_\nu}{10^{-12}\mu_\mathrm{B}}\right)^2\epsilon_\mathrm{plas},\label{plas}
\end{equation}
where $\epsilon_\mathrm{plas}$ is the standard plasmon decay  rate \citep{1996ApJS..102..411I} and $\omega_\mathrm{pl}$ is the plasma frequency \citep{1996slfp.book.....R}
\begin{equation}
    \omega_\mathrm{pl}=28.7\;\mathrm{eV}\frac{(Y_\mathrm{e}\rho)^\frac{1}{2}}{(1+(1.019\times10^{-6}Y_\mathrm{e}\rho)^\frac{2}{3})^\frac{1}{4}},
\end{equation}
where $Y_\mathrm{e}$ is the electron mole fraction and $\rho$ is the density in units of g cm$^{-3}$. The energy loss rate due to the pair production is given by \citep{2009ApJ...696..608H}
\begin{equation}
    \epsilon^\mu_\mathrm{pair}=1.6\times10^{11}\;\mathrm{erg\;g^{-1}\;s^{-1}}\left(\frac{\mu_\nu}{10^{-10}\mu_\mathrm{B}}\right)^2\frac{e^{-\frac{118.5}{T_8}}}{\rho_4},\label{pair}
\end{equation}
where $T_8$ is the temperature in units of 0.1 GK and $\rho_4=\rho/(10^4\;\mathrm{g\; cm^{-3}})$.
\begin{figure}
\centering
\includegraphics[width=8.5cm]{energy_neu.pdf}
\caption{The energy loss rates induced by the neutrino emission with $\mu_\nu=5\times10^{-12}\mu_\mathrm{B}$. The density of $\rho=10^6$ g cm$^{-3}$ is assumed. The black line shows the total standard rate and the other lines show the rates induced} by the NMM. The total energy loss rate is given by $\epsilon=\epsilon_\mathrm{std}+\epsilon^\mu_\mathrm{plas}+\epsilon^\mu_\mathrm{pair}$. \label{energy}
\end{figure}

Fig. \ref{energy} shows the energy loss rates with $\mu_\nu=5\times10^{-12}\mu_\mathrm{B}$ at $\rho=10^6$ g cm$^{-3}$, which is the typical central density at the helium flash. At the temperature of $\sim10^8$ K, the additional energy loss rate is comparable with the standard rate. Also, it is seen that the pair production is negligible in the temperature range of interest.

\section{Result}
\subsection{Evolution of the Fiducial Models}
\label{result1}
In this Section, we describe the fiducial model with $\mu_\nu=0$. Fig. \ref{fig:lumi} shows the evolution of the stellar models in the $L-A(\mathrm{Li})$ plane, where $L$ is the luminosity. The solid lines show the evolution of the fiducial model. The upper panel adopts $\alpha_\mathrm{thm}=50$, the middle panel adopts $\alpha_\mathrm{thm}=100$, and the lower panel adopts $\alpha_\mathrm{thm}=1000$.

The evolution starts from a low luminosity (the lower-left side of Fig. \ref{fig:lumi}). The lithium abundance $A(\mathrm{Li})$ stays constant during the main sequence. When the star reaches the main-sequence turnoff at $\log(L/L_\odot)=0.4$, $A(\mathrm{Li})$ starts to decrease from 2.6 to 0.8. This is because surface lithium is conveyed to the stellar interior due to the first dredge-up \citep{1967ApJ...147..624I}, and is destroyed by the proton capture. The lithium depletion becomes slower as the star evolves, but $A(\mathrm{Li})$ starts to decrease again when the star reaches the RG branch bump at $\log(L/L_\odot)=1.5$. At this point, the star develops thermohaline mixing between the convective envelope and the hydrogen burning shell \citep{2007A&A...467L..15C,2015MNRAS.446.2673L}. This happens because the mean molecular weight is inverted by $^3$He($^3$He, $2p)^4$He. Because of thermohaline mixing, lithium in the envelope is conveyed to the inner hot region and destroyed. One can see that $A$(Li) after the RG branch bump is smaller when a larger $\alpha_\mathrm{thm}$ is adopted. {In particular}, when $\alpha_\mathrm{thm}=1000$, lithium is depleted so much that the model cannot reproduce the observed samples. \citet{2015MNRAS.446.2673L} calculated $A$(Li) in RGs with $C_t=10^2$, $10^3$, and $10^4$, where $C_t=3\alpha_\mathrm{thm}/2$. {Although the value of A(Li) varies by changing a stellar evolution code adopted}, they also reported that a larger $\alpha_\mathrm{thm}$ leads to a larger variation of $A$(Li) during the evolution along the RG branch.

The decrease of $A(\mathrm{Li})$ stops when $\log(L/L_\odot)=3.2$ and it starts increasing. This is because thermohaline mixing becomes more efficient as the star expands \citep{2015MNRAS.446.2673L}. The efficient mixing helps the CF mechanism work and hence increases $A(\mathrm{Li})$. After the TRGB, the core becomes non-degenerate because of the helium flash and $L$ decreases suddenly.  As a result, core helium burning begins and a RC star is formed.
\begin{figure*}
\centering
\includegraphics[width=12cm]{L_Li.png}
\includegraphics[width=12cm]{L_Li_thm100.png}
\includegraphics[width=12cm]{L_Li_thm1000.png}
\caption{The lines show the evolution of our models with $\mu_{12}=0-5$ in the $L-A$(Li) plane. The upper panel adopts $\alpha_\mathrm{thm}=50$, the middle panel adopts $\alpha_\mathrm{thm}=100$, and the lower panel adopts $\alpha_\mathrm{thm}=1000$. The grey dots are  GALAH DR2 samples \citep{2018MNRAS.478.4513B} with reliable lithium abundances and the red dots are RC samples selected by \citet{2020NatAs.tmp..139K}. }
\label{fig:lumi}
\end{figure*}

In order to understand the evolution of $A$(Li), we introduce three timescales. We denote the timescale for the electron capture of $^7$Be as $t_\mathrm{prod}$, the timescale for $^7$Li$(p,\;\alpha)^4$He as $t_\mathrm{dest}$, and the timescale on which the material at a given radius is conveyed to the bottom of the convective envelope by thermohaline mixing as $t_\mathrm{mix}$.
\begin{figure}
\centering
\includegraphics[width=8.5cm]{D_26-crop.pdf}
\includegraphics[width=8.5cm]{timescale_26-crop.pdf}
\includegraphics[width=8.5cm]{composition_26-crop.pdf}
\caption{The profile of (upper) $D_\mathrm{thm}$, (middle) the mixing and nuclear reaction timescales, and (lower) the abundances of $^7$Be and $^7$Li. The NMM is not adopted, and $\alpha_\mathrm{thm}=100$ is assumed. The solid lines show the model when $M_\mathrm{He}=0.45M_\odot$ and the broken lines show the model when $M_\mathrm{He}=0.26M_\odot$. \label{fiducial}}
\end{figure}

Fig. \ref{fiducial} shows the profile of the $D_\mathrm{thm}$, the mixing and nuclear reaction timescales, and the $^7$Be and $^7$Li abundances in the region where $^7$Be is produced. The solid lines indicate the time when $M_\mathrm{He}=0.45M_\odot$ and the broken lines indicate the time when $M_\mathrm{He}=0.26M_\odot$. {$^7$Be is produced by $^3$He$(\alpha,\;\gamma)^7$Be at $R\sim10^{-1.3}R_\odot$.} It is seen that thermohaline mixing becomes more effective as the core grows. This is because the density in this region decreases as a function of time. The lower density leads to a larger $K$ as we can see from Eq. (\ref{K}) and thus a larger $D_\mathrm{thm}$ \citep{2015MNRAS.446.2673L}.

This efficient mixing below the convective envelope explains why $A$(Li) decreases after reaching the RG branch bump and increases near the TRGB.  In order for the CF mechanism to work, mixing in this region should be efficient enough to convey $^7$Be to the bottom of the convective envelope. When $M_\mathrm{He}=0.26M_\odot$, $^7$Be decays into $^7$Li and is destroyed by the proton capture before being conveyed to the convective envelope because $t_\mathrm{mix}$ is much longer than {$t_\mathrm{dest}$} in $R<10^{-0.2}R_\odot$. Rather, since $^7$Li in the envelope diffuses into the thermohaline region, $A$(Li) decreases in time. On the other hand, when $M=0.45M_\odot$, $t_\mathrm{mix}$ becomes shorter than before. As a result, a part of $^7$Be is conveyed to the envelope before it decays and $A$(Li) increases.

\subsection{Dependence on the NMM}
We perform stellar evolution calculations with $\mu_{12}=1-5$, where $\mu_{12}=\mu_\nu/(10^{-12}\mu_\mathrm{B})$. The adopted parameters and results are summarized in Table 1. $A(\mathrm{Li})$ and $L$  at the TRGB increase when a larger value of $\mu_\nu$ is adopted. The argument that TRGB stars become luminous when $\mu_\nu$ is adopted has been used to constrain $\mu_\nu$ \citep[e.g.][]{1996slfp.book.....R,2015APh....70....1A}.
\begin{table}
\centering
 \begin{tabular}{cc|ccc}
 $\alpha_\mathrm{thm}$&$\mu_{12}$&$M_\mathrm{He,TRGB}/M_\odot$&$\log (L/L_\odot)_\mathrm{TRGB}$&$A(\mathrm{Li})_\mathrm{RC}$\\\hline\hline
 1000&0&0.467&3.39&$-1.74$\\
1000&1&0.471&3.41&$-1.60$\\
1000&2&0.480&3.46&$-1.26$\\
1000&3&0.490&3.51&$-0.91$\\
1000&4&0.500&3.56&$-0.56$\\
1000&5&0.509&3.61&$-0.28$\\\hline
100&0&0.467&3.39&$-1.08$\\
100&1&0.471&3.41&$-0.98$\\
100&2&0.480&3.46&$-0.75$\\
100&3&0.490&3.51&$-0.46$\\
100&4&0.500&3.56&$-0.18$\\
100&5&0.509&3.61&0.07\\\hline
50&0&0.467&3.39&$-0.73$\\
50&1&0.471&3.41&$-0.70$\\
50&2&0.480&3.46&$-0.57$\\
50&3&0.490&3.51&$-0.36$\\
50&4&0.500&3.56&$-0.13$\\
50&5&0.509&3.61&0.10\\
\end{tabular}
\caption{The parameters of the models. The zero age main sequence mass and the initial metallicity are fixed to $M_\mathrm{ZAMS} =1 M_\odot$ and $Z=0.0148$, respectively. $\alpha_\mathrm{thm}$ is the thermohaline coefficient, $\mu_{12}$ is the NMM, $M_\mathrm{He,TRGB}$ is the mass of the helium core at the TRGB, $L_\mathrm{TRGB}$ is the luminosity at the TRGB, and $A(\mathrm{Li})_\mathrm{RC}$ is the lithium abundance of RC stars. Since $A(\mathrm{Li})$ decreases during the evolution of RC stars, the initial values just after the helium flash are shown.}
\end{table}

The temperature in the core gradually increases as the star evolves because the gravitational energy is released. The energy released near the TRGB is estimated to be $\sim100\;\mathrm{erg\;g^{-1}\;s^{-1}}$ \citep[e.g. Section 2.5.2 in][]{1996slfp.book.....R}. Enhanced energy loss rates in the center
can result in a more rapid evolution to advanced
stellar burning and growth of the helium core.
If we assume $\mu_{12}=5$, the additional energy loss rate induced by the NMM at the stellar center is $\epsilon^\mu_\mathrm{plas}\approx6\;\mathrm{erg\;g^{-1}\;s^{-1}}$ {just before the onset of the helium-core flash}.  Since this additional energy loss is not negligible compared with the released gravitational energy, the helium flash would be delayed.

The additional energy loss affects the structure of the helium core just before the helium flash. Fig. \ref{flash} shows the temperature profile in the helium core when the helium luminosity $L_\mathrm{He}$ reaches $100L_\odot$. The solid line shows the fiducial model with $\mu_{12}=0$ and the other lines show the models with $\mu_{12}=2$ and 5. {Because neutrinos carry away energies from the center, the temperature inversion is developed in the core. As a result,  $T$ is slightly higher in an off-center shell as seen in Fig. \ref{flash}.}
Since the core is degenerate and helium burning is highly sensitive to the temperature, the helium flash is ignited in the hottest shell. When the NMM is adopted, the location of the helium flash shifts outward because the energy loss  is higher \citep{1978ApJS...36..405S,1992A&A...264..536R}.

The previous works \cite[e.g.][]{1996slfp.book.....R} showed that core mass at the TRGB is increased by the additional energy loss. Fig. \ref{coremass1} shows the increase of helium core mass
\begin{equation}
    \delta M_\mathrm{He}=M_\mathrm{He}(\mu_{12})-M_\mathrm{He}(\mu_{12}=0),\label{deltam}
\end{equation}
where $M_{\rm He}(\mu_{12})$ is the core mass when the helium flash occurs. The solid line shows our model and the red crosses are $1M_\odot$ models with $Z=0.01$ in a previous study \citep{2015APh....70....1A}.  Although physical inputs to the models are different, our result is similar to the previous calculation.

\begin{figure}
\centering
\includegraphics[width=8.5cm]{flash_temp.pdf}
\caption{The temperature profile in the  core when the helium luminosity reaches $L_\mathrm{He}=100L_\odot$. The solid line shows the model without the NMM and the other lines adopt $\mu_{12}=2$ and 5.  \label{flash}}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=8.5cm]{core_mass-crop.pdf}
\caption{The core mass as a function of the NMM. The definition of $\delta M_\mathrm{He}$ is shown in Eq. (\ref{deltam}). The red crosses are $1M_\odot$ models with $Z=0.01$ calculated by \citet{2015APh....70....1A}. \label{coremass1}}
\end{figure}

Fig. \ref{coremass} shows the evolution of $A(\mathrm{Li})$ as a function of the helium core mass $M_\mathrm{He}$ when $\alpha_\mathrm{thm}=100$. The peaks around $M_\mathrm{He}\sim0.5M_\odot$ correspond to the helium flash. One can confirm that $A(\mathrm{Li})$ at the TRGB is higher when $\mu_\nu$ is larger.

The physical mechanisms of the lithium enhancement are twofold. When $\mu_\nu>0$ is adopted, the helium core on the ignition of the helium flash becomes heavier because of the additional energy loss. Because the flash is delayed, the CF mechanism can continue to produce more $^7$Li and $A(\mathrm{Li})$ in RC stars becomes higher. The more massive core leads to a smaller density above the hydrogen burning shell and hence a larger thermal conductivity \citep{2015MNRAS.446.2673L}. Since $D_\mathrm{thm}$ increases as a fuction of conductivity, thermohaline mixing at the TRGB becomes more efficient with a larger NMM. Fig. \ref{fig:_TRGB} shows the profile of $D_\mathrm{thm}$ and the mass fractions of $^7$Be and $^7$Li at the TRGB. It is seen that $D_\mathrm{thm}$ is larger if the NMM is larger.  Therefore the CF mechanism can convey $^7$Be to the convective envelope more effectively and  $A(\mathrm{Li})$ in RC stars becomes higher.
\begin{figure}
\centering
\includegraphics[width=8.5cm]{HeCore_Li.pdf}
\caption{The evolution of the surface lithium abundance as a function of the helium core mass. The thermohaline coefficient is fixed to $\alpha_\mathrm{thm}=100$. \label{coremass}}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=8.5cm]{D_TRGB.pdf}
\includegraphics[width=8.5cm]{composition_TRGB.pdf}
\caption{Structure of our model at the TRGB with $\mu_{12}=0$, 2, and 5. The upper panel shows the thermohaline diffusion coefficient and the lower panel shows the mass fractions of $^7$Li and $^7$Be as a function of the radius. The thermohaline coefficient is fixed to $\alpha_\mathrm{thm}=100$. }
\label{fig:_TRGB}
\end{figure}

It is seen from Fig. \ref{coremass} that $A$(Li) starts to deviate from the standard model even before the helium flash. This is explained by changes of stellar structure induced by the NMM. Fig. \ref{structure} shows the thermohaline diffusion coefficient $D_\mathrm{thm}$ and the mass fractions of $^7$Li and $^7$Be for the models with $\mu_{12}=0$, 2, and 5 in the region where  thermohaline mixing is effective. In this figure, the helium core mass is fixed to $M_\mathrm{He}=0.45M_\odot$. When a larger NMM is adopted, the radius of the helium core becomes smaller and the density in the envelope decreases.
Since $^7$Be produced via $^3$He($\alpha$,$\gamma$) is conveyed to the outer region by thermohaline mixing, the more efficient mixing leads to a larger $A$(Li).
\begin{figure}
\centering
\includegraphics[width=8.5cm]{D_thm.pdf}
\includegraphics[width=8.5cm]{massfrac.pdf}
\caption{Structure of our model when $\mu_{12}=0$, 2, and 5 and  $M_\mathrm{He}=0.45M_\odot$. The upper panel shows the thermohaline diffusion coefficient and the lower panel shows the mass fractions of $^7$Li and $^7$Be as a function of the radius. The thermohaline coefficient is fixed to $\alpha_\mathrm{thm}=100$. \label{structure}}
\end{figure}

Recently, the GALAH (Galactic Archaeology with HERMES) survey second data release (DR2) provided spectroscopic data of 342,682  stars in the Milky Way \citep{2018MNRAS.478.4513B}. \citet{2020NatAs.tmp..139K} selected stars with $\log L/L_\odot\in[1.55,\;1.85]$ and the effective temperature $T_\mathrm{eff}\in[4650,\;4900]$ K from the GALAH DR2 samples and identified them as RC stars. \citet{2020NatAs.tmp..139K} used GALAH samples that overlap with an astroseismic catalog \citep{2018ApJ...858L...7T} to distinguish RC and RG stars. They concluded that the contamination of RGs in their RC samples accounts for only $\sim10\%$.

\citet{2020NatAs.tmp..139K} found that the lithium abundance in RC stars is distributed around $A(\mathrm{Li})\sim0.71\pm0.39$. This ubiquitous enhancement of lithium has not been predicted by stellar models. When the NMM is not adopted in our model, the lithium abundance in RC stars is only $A(\mathrm{Li})=-1.08\;(-0.73)$ when $\alpha_\mathrm{thm}=100$ (50). We find that, if $\mu_{12}=5$ is adopted, $A(\mathrm{Li})$ reaches 0.07 (0.10), which is closer to the observed $A$(Li). When $\alpha_\mathrm{thm}$ is larger, $^7$Li is destroyed to a greater extent after the RG branch bump and thus $A$(Li) in RC stars becomes smaller. Although $A$(Li) is not sufficiently large when $\mu_{12}=2-4$ in both cases, the discrepancy in $A$(Li) becomes smaller if the NMM is adopted. The additional energy loss induced by the NMM is thus a candidate of a ubiquitous mechanism of the high $A$(Li) in RC stars.

Traditionally, giants with $A(\mathrm{Li})>1.5$ have been called lithium-rich giants \citep{1989ApJS...71..293B}. It is difficult to explain such extremely high lithium abundances with the NMM only. \citet{2020NatAs.tmp..139K} point out that lithium-rich giants with $A(\mathrm{Li})>1.5$ account only for $\sim3.0\%$ of RC stars. The rare population implies another mechanism which works only in a certain kind of stars.
\subsection{Dependence on Stellar Mass}
In Section \ref{result1}, we investigated $1M_\odot$ models in detail because the stellar samples \citep{2020NatAs.tmp..139K} are distributed around $(1.0\pm0.4)M_\odot$. In this section, we study $1.5M_\odot$ and $0.8M_\odot$ models to see dependence on stellar masses. Fig. \ref{mass} shows the evolutionary tracks in the plane of the luminosity $L$ and $A$(Li) in the models with $\mu_{12}=0$, 2, and 4. Efficiency of thermohaline mixing is fixed to $\alpha_\mathrm{thm}=50$. In the $1.5M_\odot$ models, the lithium enhancement near the TRGB is not seen. Also, the effect of the NMM is smaller than in the $1M_\odot$ case. {This is because the lithium depletion due to the first dredge-up is less significant in heavier stars. As a result, $A$(Li) in a more massive RG is higher. Although thermohaline mixing works even in the $1.5M_\odot$ model, the relative change in $A$(Li) during the evolution along the RG branch is smaller and the contribution of the $^7$Be transport to $A$(Li) is negligible. Such larger values of $A$(Li) in heavier RGs are reported also in previous works \citep[e.g.][]{2020arXiv200901248S}.}

On the other hand, the $0.8M_\odot$ models show large enhancements of $A$(Li) even if $\mu_{12}=0$. The lithium enhancement becomes larger when the NMM is considered as we saw in the $1M_\odot$ models. The model with $\mu_{12}=4$ does not experience the helium flash. This is  because the larger luminosity near the TRGB results in a higher mass loss rate and thus the model loses the hydrogen envelope before the helium flash. {Such a star would become a massive helium white dwarf. Although the lack of the flash would contradict observations of globular clusters, it depends on mass loss rates, which are still uncertain. We performed additional calculations of $0.8M_\odot$ models with $(\eta,\;Z)=(0.3,\;0.007)$, (0.2, 0.0148), and (0.2, 0.007) to check the dependence on the mass loss. We found that all of the additional models undergo the flash because of smaller mass loss rates. Therefore non-existence of single massive helium white dwarfs does not necessarily constrain the NMM.}
\begin{figure}
\centering
\includegraphics[width=8.5cm]{L_Li_mass.pdf}
\caption{The evolution of our $1.5M_\odot$ and $0.8M_\odot$ models with $\mu_{12}=0$, 2, and 4 and $Z=0.148$ in the $L-A$(Li) plane. Efficiency of thermohaline mixing is fixed to $\alpha_\mathrm{thm}=50$.  \label{mass}}
\end{figure}
\subsection{Dependence on Metallicity}
In Section \ref{result1}, we investigated Solar metallicity models in detail because the stellar samples are distributed around $\mathrm{[Fe/H]}=(-0.1\pm0.2)$. In this section, we study $Z=0.00136$ and 0.007 models to see dependence on metallicities. The metallicity $Z=0.00136$ corresponds to the metal-poor environment in the globular cluster M5 \citep{2013A&A...558A..12V} and $Z=0.007$ is between the Solar and M5 metallicities. Fig. \ref{metal} shows the evolutionary tracks in the plane of $L$ and $A$(Li) in $Z=0.00136$ and 0.007 models. Efficiency of thermohaline mixing is fixed to $\alpha_\mathrm{thm}=50$. It is seen that, if $\mu_{12}$ and $\alpha_\mathrm{thm}$ are fixed, $A$(Li) in RC stars is larger when $Z$ is lower. In the metal-poor models, the lithium enhancement near the TRGB is not seen if $\mu_{12}\leq2$. In the $Z=0.00136$ models, $A$(Li) is not enhanced even if $\mu_{12}=4$. In the $Z=0.007$ model, $A$(Li) slightly increases near the TRGB when $\mu_{12}=4$, although the effect of the NMM is smaller than in the $Z=Z_\odot$ case.
\begin{figure}
\centering
\includegraphics[width=8.5cm]{L_Li_M5.pdf}
\caption{The evolution of our $1M_\odot$ models with $\mu_{12}=0$, 2, and 4 and $Z=0.00136$ and 0.007 in the $L-A$(Li) plane. Efficiency of thermohaline mixing is fixed to $\alpha_\mathrm{thm}=50$.  \label{metal}}
\end{figure}
\section{Discussion}
In this paper, we discussed the effects of additional energy loss channels on $A$(Li) in RC stars. In the special case of NMM we found that the production of $^7$Li near the TRGB is activated when $\mu_{12}=1-5$ is adopted.  Although this value of the neutrino magnetic moment is smaller than current limits obtained from the reactor experiments, it exceeds the recent astrophysical upper limits. It is therefore difficult to solve the lithium problem only with the NMM. However, the additional energy loss can be induced by other physics like extra dimensions \citep{2000PhLB..481..323C} and axion-like particles \citep{1987PhRvD..36.2211R,2014PhRvL.113s1302A}. Since they are expected to result in the similar enhancement of $A$(Li), they can be a candidate of the mechanism of the lithium enhancement.

The destruction and production of $^7$Li are dependent on deep mixing including thermohaline mixing \citep{2007A&A...467L..15C,2015MNRAS.446.2673L} and magnetic buoyancy \citep{2007ApJ...671..802B}. It is desirable to investigate these mechanisms in detail.

The enhancement of energy loss rate from neutrino emission affects Li abundances on stellar surfaces through a change in stellar structure, including He core mass  and the mixing time scale. This characterizes the current theoretical prediction distinguished from other possibilities. For example, in addition to the $^7$Be production via the $^3$He($\alpha$,$\gamma$) reaction operating deep inside the stars, stellar surfaces can be polluted from outside by accretion of companion stellar ejecta or nucleosynthesis via flare-accelerated nuclei on stellar surfaces. If the observed high abundances of Li originate from nucleosynthesis in companion asymptotic giant branch stars \citep{2010MNRAS.402L..72V}, observed stars can have enhanced abundances of carbon and $s$-nuclei. On the other hand, if nuclear reactions of flare-accelerated nuclei \citep{2007A&A...469..265T} are providing $^{6,7}$Li, the isotopic fraction of $^6$Li is expected to be high. Furthermore, the observed Li-rich stars must be associated with very strong flare activities and simultaneous production of Be and B. In this way, respective possibilities are associated with different astronomical observables to be measured in future.

\begin{figure}
\centering
\includegraphics[width=8.5cm]{time_Li.pdf}
\caption{The evolution of $A$(Li) as a function of time. The black and red lines show $A$(Li) in our model with $\mu_{12}=0$ and 5, respectively. The blue line shows the $1M_\odot$ model with helium-flash induced mixing \citep{2020arXiv200901248S}. Efficiency of thermohaline mixing is fixed to $\alpha_\mathrm{thm}=50$ in our model.  \label{timescale_comparison}}
\end{figure}
A recent study \citep{2020arXiv200901248S} showed that the lithium problem in RC stars can be solved if mixing is induced by the helium flash. He assumed that mixing with the diffusion coefficient $D=10^{12}\;\mathrm{cm^2\;s^{-1}}$ is induced when the helium luminosity $L_\mathrm{He}$ exceeds $10^4L_\odot$. A major difference between  his and our models is the timescale of the surface lithium enhancement. Fig. \ref{timescale_comparison} shows $A$(Li) as a function of time. The black and red lines show our model while the blue line shows the $1M_\odot$ model in \citet{2020arXiv200901248S}. It is seen that $A$(Li) is enhanced in $\sim100$ yr during the helium flash in
his model, while it occurs $\sim1$ Myr before the flash in our model. Although the time evolution of $A$(Li) near the TRGB is dependent on model parameters including $\alpha_\mathrm{thm}$, it is robust that the lithium enhancement in our model is much slower than the model in  \citet{2020arXiv200901248S}. Hence it might be possible to distinguish the two models by searching for lithium-rich RGs with $A(\mathrm{Li})\sim 0$ near the TRGB. The model in \citet{2020arXiv200901248S} predicts that it would be impossible to find such stars, while our model predicts that such stars would be found if $A$(Li) values of a sufficient number of stars are observed.
\section*{Acknowledgements}
We thank Yerra Bharat Kumar and Gang Zhao for providing observational data. K.M. is supported
by JSPS
KAKENHI Grant Number JP19J12892.
T.K. is supported in part by Grants-in-Aid for Scientific Research of JSPS (17K05459, 20K03958).  A.B.B. is supported in part by the U.S. National Science Foundation Grants No. PHY-1806368 and
PHYS-2020275.  M.A.F. is supported by National Science Foundation Grant No. PHY-1712832 and by NASA Grant No. 80NSSC20K0498.  M.A.F. and A.B.B. acknowledge support from the NAOJ Visiting Professor program.

\section*{Data Availability}

The code and the inlist used in this paper are available at https://doi.org/10.5281/zenodo.4281702.


text:
\section{Introduction}

The initial state that leads to a type-II (or type-Ib or type-Ic)
supernova is essentially an iron white dwarf (the endpoint of nuclear
fusion reactions) in the center of a massive star; when it reaches the
Chandrasekhar mass, this core will collapse, along with the rest of
the star. The direct microphysical messengers of the gravitational
explosion are the neutrinos emitted from the hot and dense
proto-neutron star.  The detectable (few tens of MeV) neutrinos are
emitted on a timescale of several seconds, corresponding to their
diffusion outward.  These neutrinos from a specific type-II supernova
have only been detected once, from SN 1987A~\citep{Hirata,Bionta}.

The initial state that leads to a type-Ia supernova seems to be a
carbon-oxygen white dwarf (the endpoint of nuclear fusion reactions in
a low-mass star) that accretes material from a bound binary companion;
when it reaches the Chandrasekhar mass, it will experience runaway
nuclear burning.  The direct microphysical messengers of the
thermonuclear explosion are the gamma rays from the nuclear decays of
freshly produced $^{56}$Ni and $^{56}$Co.  The detectable (few MeV)
gamma rays are emitted on a timescale of up to several months,
corresponding to the nuclear lifetimes.  These gamma rays have never
been convincingly detected from a specific type-Ia supernova, though
there have been three cases in which restrictive limits were placed
(e.g.~\cite{Milne}).

The difficulty is that nearby supernovae are rare.  Although gamma
rays are of course easier to detect, type-Ia supernovae are more rare,
and the number of emitted particles per supernova is less.  In both
neutrino astrophysics and gamma-ray astrophysics, the detection of a
nearby supernova is a high priority.  (Proposed neutrino detectors
should be able to reach even to several Mpc, if neutrinos can be
detected one or two at a time using a coincidence technique
\citep{Ando}.)  It should be noted immediately that the optical
emission from both types of supernovae, while readily detected even at
great distances, does not directly or faithfully reveal the details of
the explosions or the attendant extreme conditions that allow new
tests of particle physics.

I probably don't need to further convince anyone that detecting these
direct messengers of supernova explosions would be good.  But how?
What I always say is ``{\bf Everyone complains about the supernova
rate, but nobody does anything about it.}"  While we can't wish up a
nearby supernova, we can attempt to detect the diffuse glows of
neutrinos and gamma rays made by all past type-II and type-Ia
supernovae, respectively (practically speaking, redshifts less than
about $z = 1$ are the most important for detection).  I will argue
that while this is very challenging, the signal is always there, and
the prospects are quite encouraging.  The parallel nature of the above
discussion now develops an interesting split.  The neutrino background
has never been detected, but it is believed to be dominantly produced
by type-II supernovae.  On the other hand, the gamma-ray background
has been detected, but it is now believed that the contribution from
type-Ia supernovae is subdominant.

\begin{figure}
\centerline{\includegraphics[width=2.85in]{malek}
\includegraphics[width=2.7in]{speclin}}
\caption{\label{fig:malek} {\bf Left figure:} The measured spectrum
(points with error bars) of $e^-$/$e^+$-like events in
Super-Kamiokande.  The solid line is the sum of the two atmospheric
neutrino background components, shown as lines below the solid line.
The dashed line {\it above} the solid line shows the largest allowed
fit for the detector background plus the type-II supernova neutrino
signal.  {\bf Right figure:} The projected spectrum if dissolved
gadolinium is added to Super-Kamiokande.  The shaded band is the
allowed range of supernova neutrino signals; our more recent
work~\citep{Concordance} favors the largest values.  The atmospheric
neutrino backgrounds have here been reduced by a factor $\sim 5$,
which would allow a much lower threshold, more like 10 MeV, where the
background from reactor neutrinos becomes overwhelming.  In the figure
on the right, the counts are per 1-MeV bin, whereas on the left, they
are per 4-MeV bin.  [Left figure from \cite{Malek}.  Right figure
adapted from \cite{GADZOOKS!}; the arXiv version defines
``GADZOOKS!".]}
\end{figure}

\cite{ClaytonSilk} proposed long ago that the prospects for detecting
the gamma-ray background from type-Ia supernovae were quite promising.
It wasn't until about 15 years later that people even started
considering the neutrino background from type-II supernovae; for early
work, see \cite{Seidov}, \cite{KGS}, and \cite{Woosley}.  Until very
recently, the prospects for detecting the neutrino background from
type-II supernova were thought to be at best exceedingly unlikely.
However, in 2003, the Super-Kamiokande collaboration published a flux
limit~\citep{Malek} that was about 200 times more restrictive than the
earlier limit from Kamiokande~\citep{Zhang}, and in fact in the range
of realistic models.  One surprising point that I will argue is that
this limit on the neutrino background strongly constrains the star
formation rate;\ with reasonable assumptions, this then constrains the
type-Ia supernova rate.  Using this limit, or going directly from the
measured type-Ia supernova rate data, one finds that the corresponding
gamma-ray background is well below the measured data.  Thus it now
appears {\it easier} to detect the neutrino background from type-II
supernovae than the gamma-ray background from type-Ia supernovae,
which is certainly not what anyone expected.

I will take the Super-Kamiokande neutrino limit as a starting point,
and then briefly review these constraints,
following~\cite{Concordance}.  The prospects for very significantly
improving the Super-Kamiokande sensitivity are discussed
in~\cite{GADZOOKS!}.  In brief, it was proposed that adding about
0.2\% dissolved gadolinium trichloride to Super-Kamiokande would allow
the detection of neutrons by their radiative captures on gadolinium.
This would then allow a coincidence detection of the positron and the
neutron in $\bar{\nu}_e + p \rightarrow e^+ + n$, which would greatly
reduce backgrounds which can mimic the positron by itself.  This
proposal is undergoing extensive research and development testing by
Vagins, and the prospects continue to be quite promising.

\begin{figure}
\centerline{\includegraphics[width=4in]{models}}
\caption{\label{fig:models} Recent measurements of the cosmic star
formation rate history.  The entire shaded region is consistent with
the results of the 2dF and SDSS cosmic optical spectrum.  The upper
cross-hatched region is ruled out by the limit on the type-II
supernova neutrino flux, while the lower shaded region is allowed.
Three recent (dust corrected by those authors) results are also shown
in the middle of the region; in one case, the result without dust
correction is shown for illustration.  The heavy solid line is the
result from GALEX~\citep{GALEX}.  The concordance region is defined by
the tension between the neutrino bound and the astronomical data, and
is thus concentrated at the upper edge of the lower band.  [Figure
taken from \cite{Concordance}.]}
\end{figure}

\section{The Cosmic Stellar Birth and Death Rates}

Recently, the measurements of the star formation rate history have
improved dramatically (e.g., see the summary in~\cite{Hopkins}).  An
important aspect of this is the larger and better-understood
corrections needed to compensate for obscuration by dust.  As a very
recent example of a new measurement, the GALEX results~\citep{GALEX}
are especially noteworthy.  Using the new Super-Kamiokande neutrino
flux limit~\citep{Malek}, an upper limit can be placed on the
normalization of the star formation rate history.  This requires
making some reasonable assumptions about the neutrino emission per
supernova.

\begin{figure}
\centerline{\includegraphics[width=3.5in]{SNrates}}
\caption{\label{fig:SNrates} The corresponding ranges for the type-II
and type-Ia supernova rates, following from the allowed lower band
from the previous figure.  [Figure taken from \cite{Concordance}.]}
\end{figure}

In \cite{Concordance}, we showed how the neutrino flux limit
significantly constrains the allowed range of realistic star formation
rate histories.  Interestingly, the dust-corrected results from GALEX,
for example, come in just above our deduced limit.  Very likely, this
means that the neutrino emission per supernova is somewhat less than
assumed.  On the other hand, it also means that the true flux is
probably close to the Super-Kamiokande limit, indicating that with
improved sensitivity~\citep{GADZOOKS!}, a discovery could soon be
made.  These results for the neutrino background are similar to other
recent calculations, e.g.,
\cite{Fukugita,Strigari,AndoSato,Lunardini,Daigne}, up to some
variations in the chosen inputs.  Since the astronomical factors are
already well known (and their precision is rapidly improving), a
measurement of the supernova neutrino background would alternatively
allow a new direct measurement of the supernova neutrino emission
parameters~\citep{Yuksel}, which would help resolve the lingering
mysteries of the SN 1987A data (e.g., as displayed by~\cite{Mirizzi}).

We defined a ``Concordance Model" by the tension between the upper
limit from the neutrino data, and the lower limit from the
astronomical data (i.e., these results with a somewhat smaller
correction for dust).  To validate this model, we considered the
type-II and type-Ia supernova rate histories.  Matching the type-II
supernova rate history depends just on the assumed stellar initial
mass function, whereas matching the type-Ia supernova rate history
also requires making some reasonable assumptions about the formation
efficiency and typical time delay between formation and
explosion~\citep{Watanabe}.  Using these results, which were also
confirmed by other indicators, we predicted the gamma-ray background
from type-Ia supernovae, showing that it must be far below the
measured data, which is surprising (see also similar results from
\cite{Ahn}).

\begin{figure}
\centerline{\includegraphics[width=3.5in]{MeVSpectrum}}
\caption{\label{fig:MeVSpectrum} Using the previous results, the
prediction for the type-Ia supernova gamma-ray background; the shaded
band has the same meaning as the previous figures.  [Figure taken from
\cite{Concordance}.]}
\end{figure}

\section{Conclusions}

Why would it be so important to detect the diffuse glows of neutrinos
and gamma rays from past supernovae?  That depends on your
perspective, but I should begin by reminding you that, except for the
$\simeq 20$ neutrinos from SN 1987A, and the at-most weak hints for
the gamma rays from $^{56}$Ni/$^{56}$Co from three nearby type-Ia
supernovae, these emissions have never been directly detected.  That's
strange, and frustrating, since we know that they {\it must} be
produced, and that they arise from the same physics that sources the
explosions.  For type-II supernovae, the neutrino emission can be
directly estimated from the final mass and radius of the neutron star,
which determines the gravitational energy release of $\sim 10^{53}$
erg.  For type-Ia supernovae, the gamma ray emission can be directly
estimated from the amount of synthesized $^{56}$Ni, which decays to
$^{56}$Co and then stable $^{56}$Fe.  These decays also power the
light curve, and the deduced thermonuclear energy release is $\sim
10^{51}$ erg.

For the neutrino background, the existing limit is very constraining,
allowing only something close to the most recent predictions for the
contribution from type-II supernovae.  If any other source contributed
significantly, it would thus be a big surprise; the present upper
limit on the flux of course applies also to any new sources.  Once the
neutrino flux is discovered, the spectrum shape can be tested, and the
shape corresponding to type-II supernovae is quite distinctive.  The
primary test of new physics is thus connected with the received
neutrino flux per type-II supernova.  This may be affected by neutrino
mixing among the known active flavors, possible mixing with postulated
sterile neutrinos, neutrino decay en route, modified neutrino emission
because of significant emission of new particles, etc.  The neutrino
emission per supernova can also be tested by effects on
nucleosynthetic yields, e.g., \cite{Yoshida}.

For the gamma-ray background, the situation is reversed, in that the
observed background is a more sensitive test of new physics than the
gamma ray emission per type-Ia supernova.  There are measurements of
the gamma-ray background, albeit with large error bars, over a broad
range of energies.  Just judging by the spectrum shape, there is no
clear indication of the characteristic shape corresponding to the
contribution from type-Ia supernovae near 1 MeV, and there are clearly
other sources at lower and higher energies.  An even stronger limit on
the contribution from type-Ia supernovae can be set using the star
formation and type-Ia supernova rates.  The properties of photons are
known, and it is believed that the physical conditions of the
explosion do not generally permit significant perturbations due to new
physics.  However, models of new physics, including the decay or
annihilation of dark matter (or other particles arising in
extra-dimensional models), can contribute significantly to the
observed gamma-ray background.  Once the measurements improve, and the
astrophysical components are better understood, they can be
subtracted, leading to more stringent limits on the contributions due
to new physics.  If the type-Ia supernova contribution can somehow be
isolated~\citep{CGB}, it will provide a new test of the gamma-ray
emission per supernova and the evolution of the supernova rate, both
of which could help improve the understanding required for using
type-Ia supernovae as standard candles to measure dark energy.

To summarize, a better understanding of the diffuse neutrino and
gamma-ray backgrounds, and specifically the portions associated with
supernovae, would have important implications for several fields:
\begin{itemize}
\item {\bf Particle physics:}
Neutrino properties; novel energy-loss channels in type-II supernovae; etc.
\item {\bf Nuclear physics:} Production of the light and heavy
elements; neutron star equation of state; etc.
\item {\bf Astrophysics:} Cycle of stellar birth, life, and death;
constraints on new sources; etc.
\item {\bf Cosmology:} Supernova distance indicators and dark energy;
dark matter decay or annihilation; etc.
\end{itemize}
Thus the faint fossil records of supernovae, revealed by their diffuse
neutrino and gamma ray backgrounds, offer a challenging but very
promising route to better understanding the physics of their
explosions.

This work was supported by The Ohio State University.  I thank the
organizers for the invitation to this very nice meeting, and also
especially Eli Dwek, Peter H\"oflich, and Stan Woosley for their
interesting and helpful questions and comments on my talk.  I thank
Shin'ichiro Ando, Louie Strigari, Mark Vagins, Terry Walker, Hasan
Y\"uksel, and Pengjie Zhang for very enjoyable collaborations on this
and closely related work.