# Roadmap to Thermal Dark Matter Beyond the WIMP Unitarity Bound

Ronny Frumkin<EMAIL_ADDRESS>Eric Kuflik
<EMAIL_ADDRESS>Itay Lavie<EMAIL_ADDRESS>Tal
Silverwater<EMAIL_ADDRESS>Racah Institute of Physics,
Hebrew University of Jerusalem, Jerusalem 91904, Israel

(August 27, 2024)

###### Abstract

We study the general properties of the freezeout of a thermal relic. We give
analytic estimates of the relic abundance for an arbitrary freezeout process,
showing when instantaneous freezeout is appropriate and how it can be
corrected when freezeout is slow. This is used to generalize the relationship
between the dark mater mass and coupling that matches the observed abundance.
The result encompasses well-studied particular examples, such as WIMPs, SIMPs,
coannihilation, coscattering, inverse decays, and forbidden channels, and
generalizes beyond them. In turn, this gives an approximate perturbative
unitarity bound on the dark matter mass for an arbitrary thermal freezeout
process. We show that going beyond the maximal masses allowed for freezeout
via dark matter self-annihilations (WIMP-like, $m_{\rm
DM}\gg\mathcal{O}(100\leavevmode\nobreak\ \rm TeV)$) predicts that there are
nearly degenerate states with the dark matter and that the dark matter is
generically metastable. We show how freezeout of a thermal relic may allow for
dark matter masses up to the Planck scale.

##  I Introduction and General Idea

The identity of dark matter (DM) remains one of the biggest mysteries of the
universe. Although the evidence for its existence is overwhelming, we still do
not know much about it. All we know about DM’s existence and properties is
inferred from its gravitational effects. We still do not know its particle
nature: its mass, which forces it interacts with, the number of DM particles,
if it is elementary or composite, and perhaps most importantly, how we may
discover it.

It has long been known that the current abundance of DM can be related to its
microscopic properties, such as its mass and interaction cross-sections. This
is analogous, for instance, to how the abundance of the light nuclear elements
are determined Alpher:1948ve ; Cyburt:2015mya . There are many different types
of interactions and processes in the early universe that can determine a
particle’s abundance today, but only a small number of those can give rise to
the DM we observe.

The possibility that DM is a thermal relic—namely, that its interactions
brought it into thermal and chemical equilibrium with the big bang bath, and
that its abundance is determined by those interactions when it finally departs
equilibrium—is an exciting possibility Lee:1977ua ; Scherrer:1985zt . The
abundance today is then insensitive to initial conditions because thermal
equilibrium keeps no memory about the past. Additionally, since the relic
abundance would determine the size of the interactions, thermal relics are
very predictive. Some recent ideas for new thermal processes controlling the
DM abundance beyond ordinary DM-DM annihilations, include DEramo:2010keq ;
Hochberg:2014dra ; Kuflik:2015isi ; Cline:2017tka ; DAgnolo:2017dbv ;
Garny:2017rxs ; Smirnov:2020zwf ; Kramer:2020sbb ; Fitzpatrick:2020vba ;
Frumkin:2021zng .

The common lore on thermal relics within a standard cosmological history has
been that an elementary DM candidate has an upper bound on its mass
Griest:1989wd  (for heavy dark matter in other scenarios see the review in
Carney:2022gse  and references within). This result comes from considering
WIMP-like DM, whose abundance is determined by the freezeout of annihilations
of DM pairs. The rate for a particle $\chi$, representing a DM candidate, to
self-annihilate is

$\Gamma=n_{\chi}\left<\sigma v\right>_{\chi\chi}\,,$ (1)

where $n_{\chi}$ is the DM density and $\left<\sigma v\right>_{\chi\chi}$ is
the thermally averaged cross section. We assume that this process is fast
enough to keep the DM in chemical equilibrium with the bath prior to
freezeout. As the temperature drops below the mass of the dark matter
$m_{\chi}$, the density of DM becomes exponentially Boltzmann suppressed,

$n_{\chi}=n_{\chi}^{\rm eq}\propto e^{-m_{\chi}/T}\,.$ (2)

This implies that the rate of annihilation—the process which keeps the DM in
equilibrium and annihilates away its abundance—is also becoming exponentially
slower. Eventually the process and the number of DM particles freezes out,
leaving a relic abundance of DM particles. Additionally, since cross sections
are (generally) bounded above by perturbative unitarity, there is a limit on
the amount of time the DM can stay in equilibrium and reduce its abundance.
This results in a lower bound on the DM number density at freezeout. Since the
energy density (mass$\times$density) of DM is fixed to the observed value,
this gives an upper bound on the DM mass. This number is around

$m_{\rm WIMP}\lesssim\left({m_{\rm pl}T_{\rm eq}}\right)^{1/2}\sim
100\leavevmode\nobreak\ {\rm TeV}\,,$ (3)

where $T_{\rm eq}=0.8$ eV is the temperature at matter-radiation equality and
$m_{\rm pl}=2\times 10^{18}$ GeV is the reduced Planck mass Griest:1989wd .
Another way to state the WIMP unitarity bound is that if a stable relic has
mass above the bound, it will be overabundant if the freezeout process is
self-annihilations.

However, the unitarity bound can be weakened if DM can more efficiently
annihilate away and stay in equilibrium longer. The strong bound can be
traced, in part, to the exponentially decaying rates from Eqs. (1) and (2). If
there is less exponential suppression, or it is removed entirely, then the DM
density can be reduced more and the mass limit increased. For a given
exponential dependence on the annihilation rate,

$\Gamma\propto e^{-b\,m_{\chi}/T},$ (4)

for some value of $b$, one can recalculate the perturbative unitarity bound on
the DM mass. Provided that $b$ is not too close to zero, the unitarity bound
becomes

$m_{\chi}\lesssim\left({m_{\rm pl}T_{\rm eq}^{b}}\right)^{\frac{1}{1+b}}\,.$
(5)

We will derive the exact form of the above bound, and the case where $b\sim
0$, later in the text.

At the minimum, this would require that the DM annihilate with a particle
other than itself, and lighter in mass, in order for the rate to be less
exponentially suppressed. For instance, if the DM, $\chi$, annihilates with a
particle $\psi$ (via a process $\chi+\psi\to\rm\phi_{1}+\phi_{2}$), then the
rate of annihilation scales as

$\Gamma\propto e^{-m_{\psi}/T}=e^{-(m_{\psi}/m_{\chi})m_{\chi}/T}$ (6)

and $b=m_{\psi}/m_{\chi}<1$ implies that $m_{\psi}<m_{\chi}$. Additionally, to
go beyond the WIMP unitarity bound (Eq. (3)), the total mass in the final
state in the annihilation process should not be larger than the initial
particles’ masses (except for special circumstances; see inverse decays
later). Otherwise, the process may enter a forbidden regime Griest:1990kh ;
DAgnolo:2015ujb , where additional exponential suppression in the thermally
averaged cross section can come about, since only DM particles on the
exponential Boltzmann tail can annihilate.

Moreover, the DM should be heavier than the sum of all other particles in the
process, or else it will be too unstable to be DM (in the above process Eq.
(6) this corresponds to $m_{\chi}<m_{\psi}+m_{\phi_{1}}+m_{\phi_{2}}$).

Together, these constraints suggest that in order to go beyond the WIMP
unitarity bound, there will be particles degenerate with the DM. We quantify
this in the paper and show that, in fact,

$m_{\chi}\leq\sum_{j}m_{j}<3m_{\chi},$ (7)

where $\sum_{j}m_{j}$ is the sum of the mass of all other particles (initial
and final) in the annihilation process. Therefore, if the DM mass is beyond
the WIMP unitarity bound, there is another particle close to mass with it,
which we will call $\chi^{\prime}$. However, $\chi^{\prime}$ must have
abundance less than the DM (since it’s not the DM), but its mass is above the
WIMP unitarity bound. Its abundance will be too large, unless it has some
process that efficiently depletes its abundance, or it decays away. If it is
unstable, then it potentially causes the DM itself to be unstable. If it is
stable, then it must annihilated away via a non-WIMP like process, and Eq. (7)
applies to $\chi^{\prime}$. Then there is another particle close to mass with
it, $\chi^{\prime\prime}$, and the problem continues, leading to a chain of DM
particles and interactions. This can end when we reach a particle below the
WIMP unitarity bound or one that decays. Such cases are exemplified in
scenarios within a standard cosmological history Kim:2019udq ; Kramer:2020sbb
and within an early matter dominated scenario Berlin:2017ife .

This paper is organized as follows. In Section II we study thermal freezeout
for a general process, giving a relationship between the DM mass and coupling
that matches the observed abundance. We show when instantaneous freezeout is
appropriate and how it is corrected when freezeout is slow. Section III
derives the unitarity bound for an arbitrary process and shows what is needed
to in order to go beyond the WIMP unitarity bound. Explicit examples of DM
chains with thermal DM well above the WIMP unitarity bound are given in
Section IV. We conclude with a discussion in Section V.

##  II Thermal Freezeout

| Process | Masses | $b$ |  $\beta-1$ | Unitarity bound (GeV)
---|---|---|---|---|---
WIMP | $\chi\chi\to\phi\phi$ | $m_{\phi}\ll m_{\chi}$ | 1 | $\frac{1}{x_{d}}\log(1+x_{d})$ | $\mathcal{O}(10^{5})$
SIMP | $\chi\chi\chi\to\chi\chi$ |  | 2 | $\frac{1}{2x_{d}}\log(1+x_{d}/2)$ | $\mathcal{O}(1)$
Forbidden | $\chi\chi\to\psi\psi$ | $m_{\psi}=m_{\chi}(1+\Delta)$ | $2\Delta\leavevmode\nobreak\ \leavevmode\nobreak\ $ | $\frac{1}{x_{d}}\log(1+x_{d})$ | $\mathcal{O}(10^{5-14\frac{\Delta}{1+\Delta}})$
Coannihilations | $\begin{array}[]{l}\psi\psi\to\phi\phi\\\ \chi\chi\leftrightarrow\psi\psi\end{array}$ | $m_{\psi}=m_{\chi}(1+\Delta)$ | $1+\Delta$ | $\frac{1}{x_{d}}\log(1+x_{d})$ | $\mathcal{O}(10^{5-7\frac{\Delta}{1+\Delta/2}})$
Coscattering | $\chi\phi\to\psi\phi$ | $m_{\phi}\ll m_{\psi}=m_{\chi}(1+\Delta)$ | $\Delta$ | $\begin{cases}\frac{1}{\Delta x_{d}}&{\rm for\leavevmode\nobreak\ }c=1{\rm\leavevmode\nobreak\ or\leavevmode\nobreak\ }\Delta\gg 0\\\ \frac{1}{c-2}&{\rm for\leavevmode\nobreak\ }c>2{\rm\leavevmode\nobreak\ and\leavevmode\nobreak\ }\Delta=0\\\ \end{cases}$ | $\mathcal{O}(10^{9\frac{2-\Delta}{1+\Delta}})$
Zombie | $\chi\psi\to\psi\psi$ | $m_{\psi}=m_{\chi}(1-\Delta)$ | $1-\Delta$ | $\frac{1}{(1-\Delta)x_{d}}$ | $\mathcal{O}(10^{5+7\frac{\Delta}{1-\Delta/2}})$
Inverse decays | $\chi\phi\to\psi$ | $m_{\phi}\ll m_{\psi}=m_{\chi}(1+\Delta)$ | $\Delta$ | $\frac{1}{\Delta x_{d}}\left(1+\frac{1}{\Delta x_{d}}\right)$ | $\mathcal{O}(10^{9\frac{2-\Delta}{1+\Delta}})$
Table 1: Parameters and the maximum mass allowed by perturbative unitarity for
different freezeout processes. Here $\chi$ is the DM, $\psi$ is particle close
to mass with $\chi$, and $\phi$ is a particle much lighter than both $\chi$
and $\psi$. Both $\psi$ and $\phi$ are assumed to be in thermal and chemical
equilibrium throughout the freezeout process. In the $2\to 2$ freezeout cases,
we assumed s-wave annihilations for calculating $\beta$. For coannihilations,
the process $\psi\psi\to\phi\phi$ controls the abundance, where
$\chi\chi\leftrightarrow\psi\psi$ is assumed to be in equilibrium throughout
freezeout.

Consider a DM candidate in thermal equilibrium with the bath and in chemical
equilibrium via the $n\to m$ process of the form

$\underbrace{\chi\,\chi\,\cdots\chi}_{k}\,i_{1}\,i_{2}\cdots
i_{n-k}\to\underbrace{\chi\,\chi\,\cdots\chi}_{\ell}\,f_{1}\,f_{2}\cdots\,f_{m-\ell}.$
(8)

The general form of the Boltzmann equation for the $\chi$ number density is

$\dot{n}_{\chi}+3Hn_{\chi}=-\langle\sigma v\rangle(n_{i_{1}}^{\rm eq}\cdots
n_{i_{n}}^{\rm
eq})\left(n_{\chi}^{k}-n_{\chi}^{k-\ell}\left(n_{\chi}^{eq}\right)^{k}\right),$
(9)

neglecting Pauli-blocking and stimulated emission. Here $\langle\sigma
v\rangle$ is the thermally averaged cross-section of the $n\to m$ process in
Eq. (8). On dimensional grounds, this thermally averaged cross section will
have the form

$\langle\sigma v\rangle\equiv\frac{\alpha_{\rm
eff}^{m+n-2}}{m_{\chi}^{3n-4}}\times\frac{e^{-b^{\prime}x}}{x^{c^{\prime}}}\,.$
(10)

for some coupling $\alpha_{\rm eff}$, and parameters $b^{\prime}$ and
$c^{\prime}$. Here we define the dimensionless temperature,
$x={m_{\chi}}/{T}$. The power of $\alpha_{\rm eff}$ is chosen to be the form
expected if the amplitude is determined solely by cubic interactions, such as
gauge or Yukawa interactions, but this need not be the case. The thermally
averaged cross section can also have temperature dependence. Usually the
temperature dependence can be well approximated at freezeout by a power law
(for instance $c^{\prime}=1$ for p-wave annihilation) or exponential
dependence on the temperature, which occurs for instance for forbidden
channels.

Defining the yield $Y={n}/{s}$, the Boltzmann equation takes the familiar form

$\frac{\partial Y_{\chi}}{\partial x}=-\frac{\langle\sigma v\rangle
s^{k-1}}{xH}(n_{i_{1}}^{\rm eq}\cdots n_{i_{n}}^{\rm
eq})\left(Y_{\chi}^{k}-Y_{\chi}^{k-\ell}\left(Y_{\chi}^{eq}\right)^{k}\right)\,.$
(11)

Typically, this is the form of the Boltzmann equation that is solved
numerically for the asymptotic value of the DM abundance, $Y_{\chi}(\infty)$.
An estimate for the relic abundance of $\chi$ can be obtained by the
instantaneous freezeout approximation. In this approximation, it is assumed
that the DM density departs from equilibrium and instantly stops annihilating,
where its co-moving density freezes to the value when it departed equilibrium.
Departure from equilibrium, which defines the temperature $x_{d}$, occurs
roughly when

$\frac{\Gamma(x_{d})}{H(x_{d})x_{d}}=1.$ (12)

The annihilation rate is given as:

${\Gamma(x)}={(n_{i_{1}}^{\rm eq}\cdots n_{i_{n}}^{\rm eq})(n^{\rm
eq}_{\chi})^{k-1}}{}\langle\sigma v\rangle\equiv\alpha_{\rm
eff}^{m+n-2}m_{\chi}\frac{ae^{-bx}}{x^{c}}$ (13)

for some constants, $a$, $b$, and $c$, which are defined by this equation. The
values $b$ in different freezeout processes studied in the literature can be
seen in Table 1.

Requiring that the co-moving density at the time of freezeout matches the one
observed today by Planck Aghanim:2018eyx  yields the constraint

$Y_{\chi}(\infty)=0.55\times\frac{T_{\rm eq}}{m_{\chi}}.$ (14)

Assuming instantaneous freezeout, namely that the DM density abruptly changes
from an equilibrium density to a fixed a co-moving density at $x=x_{d}$, then

${Y_{\chi}(\infty)=Y^{\rm
eq}_{\chi}(x_{d})}=\frac{45g_{\chi}}{2^{5/2}\pi^{7/2}g_{\star
s}}x_{d}^{3/2}e^{-x_{d}}$ (15)

for a DM particles freezing out while non-relativistic. Using Eqs. (12) –
(15), one finds that the mass-coupling relationship that gives the observed
abundance is

$m_{\chi}\simeq\left(\frac{\alpha^{n+m-2}_{\rm eff}m_{\rm
pl}T_{eq}^{b}}{x_{d}^{c+\frac{3}{2}b-1}}\right)^{\frac{1}{1+b}}\,.$ (16)

where $x_{d}$ is the solution to

$x_{d}\simeq\frac{1}{b+1}\log\left(\frac{\alpha_{\rm eff}^{m+n-2}m_{\rm
pl}}{x_{d}^{c-\frac{5}{2}}T_{\rm eq}}\right).$ (17)

Overall numerical factors are left out for simplicity. These factors can be
recovered here and later with the replacement

$\displaystyle m_{\rm pl}$ $\displaystyle\to$ $\displaystyle
a\sqrt{\frac{90}{\pi^{2}g_{\star}(x_{d})}}m_{\rm pl}$ $\displaystyle T_{\rm
eq}$ $\displaystyle\to$ $\displaystyle{\frac{\pi^{2}(2\pi)^{3/2}g_{\star
s}(x_{d})}{30(1+\Omega_{\rm b}/\Omega_{\rm dm})g_{\chi}}}T_{\rm eq}\,.$ (18)

The factor is typically an $\mathcal{O}(1-10)$ correction to the mass-coupling
relationship.

The relic abundance is mostly determined by the mass, coupling, and the degree
of exponential suppression of the annihilation rate, $b$. From Eq. (16) one
can obtain as particular examples the mass-coupling relationships for well-
studied scenarios such as WIMPs Lee:1977ua ; Scherrer:1985zt , SIMPs
Hochberg:2014dra , forbidden Griest:1990kh ; DAgnolo:2015ujb , co-
annihilations Griest:1990kh , co-scattering DAgnolo:2017dbv ; Kim:2019udq ,
zombie Berlin:2017ife ; Kramer:2020sbb  and inverse decay processes
Garny:2017rxs ; Frumkin:2021zng , which are given in Table 1.

Figure 1:  Left: $m_{\chi}Y_{\chi}$ vs. $x$ for DM freezeout. The orange
curves shows the solution of the Boltzmann equation for a $m_{\chi}=1$ GeV
WIMP. The red point shows the value $x_{d}$, given by Eqs (17) and (24). After
decoupling, the abundance continues to reduce by a factor
$e^{-x_{d}(\beta-1)}$ until it completely freezes out to a final yield,
$m_{\chi}Y_{\chi}(\infty)=0.55T_{\rm eq}$. In gray we also show the abundance
curves for for a $m_{\chi}=1$ GeV DM candidate whose relic abundance is
determined by a coscattering or $3\to 2$ (SIMP) freezeout process. For
coscattering, decoupling happens very early, freezeout is slow and $\beta$ is
quite large. On the other hand, SIMP DM freezes out very fast and instanaeous
freezeout works as a better approximation of the final relic abundance. Right:
Mass vs. coupling relationship that matches the observed abundance. The solid
line shows the results from the numerical simulation, while the dashed line
shows the analytical solution in Eq. (23). The definitions of $\Delta$ for
zombie and forbidden annihilations are given in Table 1.

Instantaneous freezeout is a good approximation in many scenarios, because the
annihilation rate is dropping exponentially with time (or temperature).
However, when exponential suppression is weaker ($b<1$) instantaneous
freezeout may no longer be a good approximation for estimating the relic
abundance. This will often be the case when going beyond the WIMP unitarity
bound, which requires $b<1$. Indeed, it was shown in Ref. Kim:2019udq ;
Kramer:2020sbb  that freezeout can be very slow.

We define a parameter $\beta$ as a measure of the deviation from instantaneous
freezeout as

$Y_{\chi}(\infty)\equiv Y^{\rm eq}_{\chi}(x_{d})e^{-x_{d}(\beta-1)}$ (19)

For $\beta=1$, the instantaneous freezeout approximation is perfect. For $b\ll
1$, the correction can be large. Large deviations from instantaneous freezeout
are expected when the DM is annihilating with a lighter particle than itself.
Therefore, we consider here corrections to instantaneous freezeout when there
is only one DM particle in the initial state ($k=1$), although the method here
can be applied to any scenario. The Boltzmann equation (11) for $k=1$, using
the definitions above, can be written as

$\frac{\partial Y_{\chi}}{\partial x}=-\frac{\alpha_{\rm
eff}^{m+n-2}\frac{m_{\chi}}{x^{c}}e^{-bx}}{xH}\left(Y_{\chi}-Y_{\chi}^{eq}\right)\,.$
(20)

There exists an integral solution to this Boltzmann equation. However, an easy
way to get a good estimate for the correction to Eq. (16) is to evolve the
equation from the moment $\chi$ departs equilibrium, which we approximate as
occurring at $x=x_{d}$. For $x>x_{d}$ the inverse reaction can be dropped. The
solution is

$Y_{\chi}(\infty)=Y^{\rm
eq}_{\chi}(x_{d})\exp\left(-x_{d}\int_{1}^{\infty}x^{-c+1}e^{-bx_{d}(x-1)}\,dx\right)\,.$
(21)

Thus we see an additional possible exponential suppression compared to the
instantaneous freezeout case given in Eq. (15), where

$\beta=1+\int_{1}^{\infty}dx\,x^{-c+1}e^{-bx_{d}(x-1)}.$ (22)

In this form, the mass-coupling relationship in Eq. (16), can easily be found
as before, but using Eq. (19) instead of Eq. (15). One finds

$m_{\chi}\simeq\left(\frac{\alpha^{n+m-2}_{\rm eff}m_{\rm
pl}T_{eq}^{b/\beta}}{x_{d}^{c+\frac{3}{2}\frac{b}{\beta}-1}}\right)^{\frac{1}{1+b/\beta}}\,.$
(23)

where

$x_{d}\simeq\frac{1}{b+\beta}\log\left(\frac{\alpha_{\rm eff}^{m+n-2}m_{\rm
pl}}{x_{d}^{c-\frac{5}{2}}T_{\rm eq}}\right).$ (24)

Again, overall numerical factors can be restored with the replacement Eq.
(18). Notice that Eq. (22) and Eq. (24) still provide two equations needed to
solve for $x_{d}$ and $\beta$. However, the integral $\beta$ can usually be
done semi-analytically. For $b\gtrsim 1$, where instantaneous freezeout it
expected to be good, one finds that $\beta\simeq 1+1/(bx_{d})$. For $bx_{d}\ll
1$ (and further more for $c\to 0$), $\beta$ can be quite large. In the extreme
$b\to 0$ case, $\beta\simeq 1+1/(c-2)$. Note that for $b=0$ and $c\leq 2$,
freezeout never occurs, because the annihilation rate is decreasing slower
relative to the Hubble rate; this occurs for decay processes and can also
occur for co-scattering with a massless mediator. The value of $\beta$ for co-
scattering, inverse decays, and zombie processes are given in Table 1.111We
also include values of $\beta$ for WIMPs and SIMPs and forbidden, still
defined by Eq. (19). Here $\beta$ is calculated in the same way, but starting
with the Boltzmann equation Eq. (11), and dropping the inverse annihilation
term for $x>x_{d}$.

In left panel of Fig. 1 we show a schematic of freezeout, exemplifying the
meaning of the different parameters defined in this section. In the right
panel of Fig. 1 we show the mass vs coupling relationship that matches the
observed abundance. The solid line indicates the results from the numerical
simulation, while the dashed line indicates the analytical solution in Eq.
(23). The numerical and analytical results are in excellent agreement across
many different freezeout processes.

##  III Unitary bound

From Eq. (16), the perturbative unitarity bound on WIMP DM is immediately
apparent. For DM-DM annihilations into two light bath particles, we have $b=1$
and $m=n=2$, which gives

$m_{\rm WIMP}\simeq\alpha_{\rm eff}\left(m_{\rm pl}T_{\rm
eq}\right)^{1/2}\simeq\alpha_{\rm eff}\times 30\leavevmode\nobreak\ \rm TeV.$
(25)

For perturbative couplings, one sees that the mass of a WIMP DM candidate is
never greater than a few hundred TeV. However, going beyond this mass is still
possible for a thermal relic, provided that $b<1$—in other words, that the
annihilation rate is less exponentially suppressed as the DM density is
annihilating away.222There are another two interesting possibilities to get an
exponentially fast rate, but they require species to be out of chemical
equilibrium. One is to consider enhancement via stimulated emission if the DM
annihilates to a boson with large occupation number. Another possibility is if
DM annihilates with a particle that has already frozen out and therefore has a
large chemical potential. See Refs. Garny:2017rxs ; Dror:2016rxc ;
Kramer:2020sbb ; DAgnolo:2017dbv  for freezeout mechanisms that involve
particles chemically decoupled from the bath.  The unitarity bound for an
arbitrary interaction is given by

$m_{\chi}\lesssim m_{\rm pl}\left(\frac{T_{\rm eq}}{m_{\rm
pl}}\right)^{\frac{1}{1+\beta/b}}{x_{d}^{-\frac{c-1}{1+b/\beta}}}.$ (26)

Thus for sufficiently small $b$ and $c\leq 1$, one can expect that there might
be a perturbative thermal candidate for DM all the way up to the Planck scale.

As we will now show, going above the perturbative WIMP unitarity bound
generically implies that there are additional nearly degenerate particles to
the DM, and that the DM is metastable. To go beyond WIMP unitarity, first
consider the process

$\chi+i_{1}+\cdots i_{n}\to f_{1}+\cdots f_{m}\ $ (27)

setting the abundance.333We assume that this is the dominant process, which
itself can impose hierarchy on the masses of the particles involved. We need
only consider the case with one DM in the initial state, otherwise $b\geq 1$.
The rate for this process scales as

$\Gamma\propto\exp\left[-\begin{cases}{\frac{\sum_{i}m_{i}}{m_{\chi}}x}&\scriptstyle
m_{\chi}+\sum\limits_{i}m_{i}>\sum\limits_{f}m_{f}\\\
{\frac{\sum_{f}m_{f}-m_{\chi}}{m_{\chi}}x}&\scriptstyle
m_{\chi}+\sum\limits_{i}m_{i}\leq\sum\limits_{f}m_{f}\end{cases}\right]$ (28)

where the second case is for forbidden DM. Requiring $b<1$ imposes

$\sum_{i}m_{i}<m_{\chi}\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\
m_{\chi}+\sum_{i}m_{i}>\sum_{f}m_{f}$ (29)

or

$\sum_{i}m_{f}<2m_{\chi}\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\
m_{\chi}+\sum_{i}m_{i}\leq\sum_{f}m_{f}.$ (30)

Finally, we require that the DM is stable from the decay induced by moving all
initial particles to the final state

$\sum_{i}m_{i}+\sum_{f}m_{f}\geq m_{\chi}.$ (31)

Combining these conditions gives one of the main results of this paper,

$m_{\chi}\leq\sum_{i}m_{i}+\sum_{f}m_{f}<3m_{\chi}.$ (32)

Therefore, one of the particles, which we will call a partner, will always
have mass close to the DM particle.444It is possible that there is a very
large number of particles in the process, possibly avoiding a degeneracy.
However, this would give a large phase space suppression in the cross-section
which would ultimately reduce the DM mass. However, since the partner is close
to the mass of the DM, it must have a process with $b<1$ in order to reduce
its abundance for the same reasons. Therefore this degeneracy repeats itself,
leading to a potential chain of DM reactions. The chain will end either when
the mass of a partner goes below the WIMP unitarity bound, or when it reaches
a particle that decays in equilibrium with the bath.

Eq. (32) gives a constraint on the masses of the particles involved in the
freezeout process, when we require that the DM mass is above the WIMP
unitarity bound ($b=1$). A stronger constraint can be obtained if we require
that the DM mass be even heavier. For instance, we can impose that the DM is
heavier than the mass given by the unitarity limit in Eq. (26) for arbitrary
$b$. We repeat the same steps as above, but instead of requiring that the rate
in Eq. (28) is faster than the WIMP, we require it to be faster than freezeout
with an exponential suppression ${\Gamma\propto e^{-bx}}$. Eqs. (29) and (30)
become

$\sum_{i}m_{i}<bm_{\chi}\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\
m_{\chi}+\sum_{i}m_{i}>\sum_{f}m_{f}$ (33)

or

$\sum_{i}m_{f}<(1+b)m_{\chi}\leavevmode\nobreak\ {\rm and}\leavevmode\nobreak\
m_{\chi}+\sum_{i}m_{i}\leq\sum_{f}m_{f}.$ (34)

Therefore, the constraint on the masses for DM above the unitarity bound value
found in Eq. (26) is

$m_{\chi}\leq\sum_{i}m_{i}+\sum_{f}m_{f}<(1+2b)m_{\chi}.$ (35)

##  IV Freezeout chains

Having established a roadmap to going beyond the WIMP unitarity limit for
thermal DM, we now turn to new mechanisms that predict super heavy DM. We will
focus on two chain mechanisms allowing for orders of magnitude larger mass
than the WIMP unitarity bound. First we discuss a zombie chain, in which
absolutely stable DM is possible well beyond the WIMP unitarity bound with
only a small chain. Next we study an inverse decay chain, where DM is
metastable, but whose mass may be much larger.

Figure 2: Thermal evolution of the comoving energy density
$m_{\chi_{i}}Y_{\chi_{i}}$ vs. temperature $x=m_{\chi_{i}}/T$ for the zombie
chain (left) and inverse decay chain (right). The solid lines shows the
numerical abundance for each particle, while the dashed lines represent
equilibrium abundance curves. In each case, we indicate the DM mass, mass-
splitting, effective coupling and number of particles in the chain on the
plot.

###  IV.1 The Zombie Chain

Next we study a chain based on zombie type-interactions. It was shown in
Kramer:2020sbb  that a zombie process can allow for a heavy DM candidate
without a chain, but with an unstable DM particle. Here we introduce a chain
of zombie interactions, and show that it supports exponentially larger mass
for each additional particle in the chain, up to the bound set by Eq. (23),
while allowing for completely stable DM.

The zombie chain consists of $N$ DM particles, $\chi_{i}$ ($i=1...N$) with
zombie-type nearest neighbor interactions

$\chi_{i}+\chi_{i+1}\leftrightarrow\chi_{i+1}+\chi_{i+1}\,.$ (36)

For simplicity, we will assume that each chain interaction has the same
strength and that the mass hierarchy is constant going from the heaviest to
lightest. Namely, we take $\left<\sigma
v\right>_{i,{i+1}\to{i+1},{i+1}}\equiv\alpha^{2}/m_{\chi_{i}}^{2}$ and
$m_{i+1}/m_{i}$ independent of $i$. The last particle in the chain,
$\chi_{N}$, is assumed to also have direct annihilations into the SM bath
particles

$\chi_{N}+\chi_{N}\leftrightarrow{\rm SM}+{\rm SM}\,.$ (37)

Since the abundance of $\chi_{N}$ is determined by standard self-
annihilations, and is assumed to be stable, its mass must be less that the
WIMP unitarity bound, i.e., $m_{\chi_{N}}\lesssim 100$ TeV.

For the zombie interactions, the annihilation rate is

$\Gamma_{i,{i+1}\to{i+1},{i+1}}=n_{{i+1}}^{\rm eq}\left<\sigma
v\right>_{i,{i+1}\to{i+1},{i+1}}\propto
e^{-\frac{m_{{i+1}}}{m_{i}}\frac{m_{{i}}}{T}}$ (38)

and here $b_{\rm zombie}=m_{i+1}/m_{i}$. Following Eq. (32) and Eq. (35), we
see that in order to go beyond the WIMP unitarity bound, each link in the
chain must satisfy

$m_{i}\leq 3m_{i+1}=3b\,m_{i}$ (39)

The highest mass is obtained for the smallest allowed value $b=1/3$; below
this value, the DM is unstable. Furthermore, for the interactions above, each
$\chi_{i}$ is absolutely stable, although each will have an exponentially
smaller relic abundance than the previous.

Precise solutions to the relic abundance can be obtained by solving the
coupled Boltzmann equations. However, the freezeout occurs sequentially: first
$\chi_{1}$ freezes out from $\chi_{2}$, then $\chi_{2}$ freezes out from
$\chi_{3}$, and so on. At each stage each particle freezes out with its
neighbors still in chemical equilibrium. For each $\chi_{i}$, we can simply
use Eqs. (16) and (17) to determine the abundance.

For $b=1/3$, the maximum allowed mass for $\chi_{1}$ is given by Eq. (26),
which gives

$m_{\chi_{1}}\lesssim 10^{11}\leavevmode\nobreak\ \rm GeV$ (40)

(for $\alpha_{\rm eff}=10$). In terms of the last particle in the chain, the
${\chi_{1}}$ mass is $m_{\chi_{1}}=3^{N-1}m_{\chi_{N}}$. At least $N=13$
particles are required in the chain to achieve $m_{\chi_{1}}\simeq
10^{11}\leavevmode\nobreak\ \rm GeV$, if $m_{\chi_{N}}\simeq
100\leavevmode\nobreak\ \rm TeV$.

In the left panel of Fig. 2, we show the numerical solution for $N=18$,
$b=0.337$, $\alpha_{\rm eff}=1$, and ${m_{\chi_{1}}=3.6\times 10^{9}}$ GeV,
which exceeds the WIMP unitarity bound Eq. (3). The thermal evolution of the
comoving density is shown for each $\chi_{i}$. The abundances of $\chi_{1}$
through $\chi_{N-1}$ freezeout via zombie processes and have decreasing relic
abundance $m_{\chi_{i}}Y_{\chi_{i}}/m_{\chi_{i+1}}Y_{\chi_{i+1}}\simeq 70$.
The abundance of $\chi_{N}$ freezes out via self annihilations Eq. (37), and
therefore has a relatively large relic abundance, but much smaller than the
${\chi_{1}}$ abundance.

###  IV.2 The Inverse Decay Chain

Consider freezeout of $\chi$ via inverse decay (INDY) process Garny:2017rxs ;
Frumkin:2021zng

$\chi+\gamma^{\prime}\to\psi,$ (41)

where $\gamma^{\prime}$ is a light particle in equilibrium with the bath. It
is easy to see from detailed balance that the rate for depleting $\chi$ is

$\Gamma_{\chi+\gamma^{\prime}\to\psi}=\frac{n_{\psi}^{\rm eq}}{n_{\chi}^{\rm
eq}}\Gamma_{\psi\to\chi+\gamma^{\prime}}\propto
e^{-\frac{m_{\psi}-m_{\chi}}{m_{\chi}}x}$ (42)

where $\Gamma_{\psi\to\chi+\gamma^{\prime}}$ is the partial width of $\psi$.
For inverse decays, in terms of the definition of the rate in Eq. (13),

$b=\frac{m_{\psi}-m_{\chi}}{m_{\chi}}\equiv\Delta,\qquad c=0,$ (43)

and

$\alpha_{\rm
eff}=\frac{\Gamma_{\psi\to\chi+\gamma^{\prime}}}{m_{\chi}}(1+\Delta)^{5/2}\,.$
(44)

The suppression of the rate in temperature is only due to the mass-splitting
between $\psi$ and $\chi$. Therefore, the rate only decreases when the mass
splitting is of order the temperature, i.e. $x\Delta>1$. For sufficiently
small $\Delta$, decoupling and freezeout can occur arbitrarily late. For these
reasons it is expected that arbitrarily high DM masses are possible.

However, for the reasons discussed in the previous section, if $m_{\chi}$ is
well above the WIMP unitarity bound, then there needs to be a way for $\psi$
to deplete in order to not be over abundant. If it is unstable while decaying
in equilibrium, then certainly $\chi$ will be too short lived to be DM. Thus,
we need to consider $\psi$ depleting by inverse decays as well. To this end,
we consider a chain of inverse decays in order to achieve arbitrarily high DM
mass. The chain must be long enough so that the phase space in the $\chi$
decay is small enough to stabilize the DM on cosmological scales.

Consider the chain of inverse decay processes

$\chi_{i}+\gamma^{\prime}\to\chi_{i+1},$ (45)

for $i=1\ldots N$, and the decay for $\chi_{N}$ directly into the Standard
Model bath

$\chi_{N}\to{\rm sm}+{\rm sm}.$ (46)

The decay rate $\chi_{N}$ is taken to be fast, so that it is an equilibrium
process, i.e., $\Gamma_{\chi_{N}\to{\rm sm}+{\rm sm}}>H(m_{\chi_{N}})$. For
simplicity we consider the case where the decay width for each process in the
chain is constant
${\Gamma\equiv\Gamma_{\chi_{i+1}\to\chi_{i}+\gamma^{\prime}}}$ and a constant
mass splitting ${m_{\chi_{i+1}}=m_{\chi_{i}}(1+\Delta)}$.

The relic abundance can be found by solving the $N$ coupled Boltzmann
equations:

$\displaystyle xH\frac{\partial Y_{1}}{\partial x}$ $\displaystyle=$
$\displaystyle{\Gamma}{}\left(Y_{2}-Y_{1}\frac{n_{{2}}^{\rm eq}}{n_{1}^{\rm
eq}}\right)$ $\displaystyle\vdots$ $\displaystyle xH\frac{\partial
Y_{i}}{\partial x}$ $\displaystyle=$
$\displaystyle{\Gamma}{}\left(Y_{{i+1}}-Y_{i}\frac{n_{{i+1}}^{\rm
eq}}{n_{i}^{\rm eq}}\right)+{\Gamma}{}\left(Y_{{i-1}}\frac{n_{{i}}^{\rm
eq}}{n_{{i-1}}^{\rm eq}}-Y_{{i}}\right)$ $\displaystyle\vdots$ $\displaystyle
xH\frac{\partial Y_{N}}{\partial x}$ $\displaystyle=$
$\displaystyle{\Gamma}{}\left(Y_{{N-1}}\frac{n_{{N}}^{\rm eq}}{n_{{N-1}}^{\rm
eq}}-Y_{{N}}\right)-{\Gamma_{N\to{\rm sm}}}\left(Y_{N}-Y_{N}^{\rm
eq}\right)\,,$

where $x=m_{\chi_{1}}/T$.

We show a numerical solution in the right panel of Fig. 2 for $N=30$,
$\Delta=0.01$, $\alpha_{\rm eff}=1$, and $m_{\chi_{1}}=7.35\times 10^{11}$
GeV, which exceeds the WIMP unitarity bound, Eq. (3). One can see that the
analytic mass-coupling relationship derived in Eq. (23) does not well-describe
the results here. Similar to what was observed for the coscattering chain in
Ref. Kim:2019udq , the chain can weaken the strength of the DM depletion. For
scattering, it was shown that the effective rate for annihilation is
suppressed by $N^{2}$. In other words, one should replace $\alpha_{\rm
eff}\to\alpha_{\rm eff}/N^{2}$ in order to obtain an estimate for correct
coupling needed to match the observed abundance. Exact analytic solutions to
the coupled set of equations are difficult to obtain, and we leave this to
future work inversechain .

Finally, we comment on the stability of the DM for the example shown. The DM
candidate, $\chi_{1}$, can decay to $N$ SM particles via $N-1$ off-shell
$\chi_{i}$ particles. We estimate the phase-space and choose $N=30$ as a safe
value that guarantees that the phase space in the $\chi$ decay is small enough
to stabilize the DM on cosmological scales.

##  V Outlook

In this paper we have derived general formula for the thermal relic abundance
of dark matter for an arbitrary freezeout process, and have outlined the
roadmap to obtain a thermal relic with mass above the standard WIMP
perturbative unitarity bound of a few hundred TeV.

Thermal freezeout provides a well-motivated framework for studying the relic
abundance of DM. The relic abundance can be determined in the instantaneous
freezeout approximation for most cases with Eq. (16), or with Eq. (23) which
corrects for when freezeout is slow. Depending on the process controlling
freezeout, there is a different bound on the dark matter mass from
perturbative unitarity, which can be substantially higher than that of the
WIMP standard lore. To go beyond the standard WIMP unitarity bound, we find
that there are particles degenerate with the DM, possibly constituting a long
chain of DM interactions.

Much work remains to be done to better understand heavy dark matter in context
of a chain and to populate the model parameter space that achieves such heavy
thermal relics. We leave this to future work.

Acknowledgments. We thank Hyun Min Lee and Josh Ruderman for useful
discussions. We are especially grateful to Yonit Hochberg and Nadav
Outmezguine for useful discussions and comments on the manuscript. The works
of RF, EK, IL, and TS are supported by the US-Israeli Binational Science
Foundation (grants 2016153 and 2020220), the Israel Science Foundation (grant
No. 1111/17), and by the I-CORE Program of the Planning and Budgeting
Committee (grant No. 1937/12).

## References

  * (1) R.A. Alpher, H. Bethe and G. Gamow, _The origin of chemical elements_ , _Phys. Rev._ 73 (1948) 803.
  * (2) R.H. Cyburt, B.D. Fields, K.A. Olive and T.-H. Yeh, _Big Bang Nucleosynthesis: 2015_ , _Rev. Mod. Phys._ 88 (2016) 015004 [1505.01076].
  * (3) B.W. Lee and S. Weinberg, _Cosmological Lower Bound on Heavy Neutrino Masses_ , _Phys. Rev. Lett._ 39 (1977) 165.
  * (4) R.J. Scherrer and M.S. Turner, _On the Relic, Cosmic Abundance of Stable Weakly Interacting Massive Particles_ , _Phys. Rev. D_ 33 (1986) 1585.
  * (5) F. D’Eramo and J. Thaler, _Semi-Annihilation of Dark Matter_ , _JHEP_ 06 (2010) 109 [1003.5912].
  * (6) Y. Hochberg, E. Kuflik, T. Volansky and J.G. Wacker, _Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles_ , _Phys. Rev. Lett._ 113 (2014) 171301 [1402.5143].
  * (7) E. Kuflik, M. Perelstein, N.R.-L. Lorier and Y.-D. Tsai, _Elastically Decoupling Dark Matter_ , _Phys. Rev. Lett._ 116 (2016) 221302 [1512.04545].
  * (8) J.M. Cline, H. Liu, T. Slatyer and W. Xue, _Enabling Forbidden Dark Matter_ , _Phys. Rev. D_ 96 (2017) 083521 [1702.07716].
  * (9) R.T. D’Agnolo, D. Pappadopulo and J.T. Ruderman, _Fourth Exception in the Calculation of Relic Abundances_ , _Phys. Rev. Lett._ 119 (2017) 061102 [1705.08450].
  * (10) M. Garny, J. Heisig, B. Lülf and S. Vogl, _Coannihilation without chemical equilibrium_ , _Phys. Rev. D_ 96 (2017) 103521 [1705.09292].
  * (11) J. Smirnov and J.F. Beacom, _New Freezeout Mechanism for Strongly Interacting Dark Matter_ , _Phys. Rev. Lett._ 125 (2020) 131301 [2002.04038].
  * (12) E.D. Kramer, E. Kuflik, N. Levi, N.J. Outmezguine and J.T. Ruderman, _Heavy Thermal Dark Matter from a New Collision Mechanism_ , _Phys. Rev. Lett._ 126 (2021) 081802 [2003.04900].
  * (13) P.J. Fitzpatrick, H. Liu, T.R. Slatyer and Y.-D. Tsai, _New Pathways to the Relic Abundance of Vector-Portal Dark Matter_ , 2011.01240.
  * (14) R. Frumkin, Y. Hochberg, E. Kuflik and H. Murayama, _Thermal Dark Matter from Freezeout of Inverse Decays_ , 2111.14857.
  * (15) K. Griest and M. Kamionkowski, _Unitarity Limits on the Mass and Radius of Dark Matter Particles_ , _Phys. Rev. Lett._ 64 (1990) 615.
  * (16) D. Carney et al., _Snowmass2021 Cosmic Frontier White Paper: Ultraheavy particle dark matter_ , 2203.06508.
  * (17) K. Griest and D. Seckel, _Three exceptions in the calculation of relic abundances_ , _Phys. Rev._ D43 (1991) 3191.
  * (18) R.T. D’Agnolo and J.T. Ruderman, _Light Dark Matter from Forbidden Channels_ , _Phys. Rev. Lett._ 115 (2015) 061301 [1505.07107].
  * (19) H. Kim and E. Kuflik, _Superheavy Thermal Dark Matter_ , _Phys. Rev. Lett._ 123 (2019) 191801 [1906.00981].
  * (20) A. Berlin, _Wimps with Guts: Dark Matter Coannihilation with a Lighter Species_ , _Phys. Rev. Lett._ 119 (2017) 121801 [1704.08256].
  * (21) Planck collaboration, _Planck 2018 Results. Vi. Cosmological Parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * (22) J.A. Dror, E. Kuflik and W.H. Ng, _Codecaying Dark Matter_ , _Phys. Rev. Lett._ 117 (2016) 211801 [1607.03110].
  * (23) R. Frumkin, Y. Hochberg, E. Kuflik, I. Lavie and T. Silverwater work in progress.


11institutetext: School of Physical Science and Technology, Xinjiang
University, Urumqi 830046, PR China
11email<EMAIL_ADDRESS>22institutetext: Center for Theoretical
Physics, Xinjiang University, Urumqi 830046, PR China
33institutetext: College of Mechanical and Electronic Engineering, Tarim
University, Alar 843300, PR China
44institutetext: CAS Key Laboratory of Optical Astronomy, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
55institutetext: School of Astronomy and Space Science, University of Chinese
Academy of Sciences, Beijing 100101, China

# Li-rich and super Li-rich giants produced by element diffusion††thanks: The
full version of Table LABEL:tab1 only available in electronic form at the CDS
via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5)

Jun Gao 1122 Chunhua Zhu 1122 Jinlong Yu 33 Helei Liu 1122 Xizhen Lu 1122
Jianrong Shi 4455 Guoliang Lü 112*2*

###### Abstract

Context. About 0.2-2% of giant stars are Li-rich, whose lithium abundance
(A(Li)) is higher than 1.5 dex. Among them, near 6% are super Li-rich with
A(Li) exceeding 3.2 dex. Meanwhile, the formation mechanism of these Li-rich
and super Li-rich giants is still under debate.

Aims. Considering the compact He core of red giants, attention is paid to the
effect of element diffusion on A(Li). In particular, when the He core flash
occurs, the element diffusion makes the thermohaline mixing zone extend inward
and connect to the inner convection region of stars. Then, a large amount of
7Be produced by the He flash can be transferred to stellar surface, finally
turning into 7Li. Thus, the goal of this work is to propose the mechanism of
A(Li) enrichment and achieve the consistency between the theoretical and
observation data.

Methods. Using the Modules for Experiments in Stellar Astrophysics (MESA), we
simulate the evolution of low-mass stars, with considering the effects of
element diffusion on the Li abundances. The timescale ratio of Li-rich giants
to normal giants is estimated by population synthesis method. Then we get the
theoretical value of A(Li) and make a comparison with observations.

Results. Considering the influence of element diffusion in the model results
in the increase of lithium abundance up to about 1.8 dex, which can reveal Li-
rich giants. Simultaneously, introducing high constant diffusive mixing
coefficients ($D_{\rm mix}$) with the values from $10^{11}$ to $10^{15}$cm2
s-1 in the model allows A(Li) to increase from 2.4 to 4.5 dex, which can
explain the most of Li-rich and super Li-rich giant stars. The population
synthesis method reveals that the amount of Li-rich giants among giants is
about 0.2-2%, which is consistent with observation estimated levels.

Conclusions. In our model, the element diffusion mainly triggered by the
gravity field changes the mean molecular weight at the junction zone between
the stellar envelope and the He core, which makes the thermohaline mixing
region expanding to the inner convection region of stars. A transport channel,
efficiently transporting 7Be in the hydrogen burning region of the star to the
convective envelope where 7Be decays into 7Li, is formed. Combing a high
constant diffusive mixing coefficients, the transport channel can explain the
origin of Li-rich and super Li-rich giant, even the most super Li-rich giants.

###### Key Words.:

stars: evolution; standards; stars: low-mass; diffusion; stars: abundances

##  1 Introduction

Lithium (Li) is one of the important elements to study the origin of the
universe. In the evolution of low-mass stars, Li begins to deplete in the main
sequence (MS) stage. This process experienced the first dredge-up and some
deep mixing, and most of the Li will be consumed. The phenomenon has been
predicted by a standard staller evolution model(Deepak & Reddy 2019), and
confirmed by numerous observations (Brown et al. 1989; Lind et al. 2009; Liu
et al. 2014; Kirby et al. 2016).

Since the first giant star with a high Li abundance was discovered by
Wallerstein & Sneden (1982), the development of the related research field
challenged the traditional stellar evolution mechanism. These stars were
called Li-rich giants, with a classic definition of A(Li) $\geq$ 1.5 dex(Iben
1967a, b). 111Here, A(Li) is the Li abundance expressed as A(Li) = log [n (Li)
/ n (H)] + 12, where n (Li) and n (H) is the number density of Li and
hydrogen, respectively. It is noteworthy that Li-rich giants are extremely
rare, accounting about 1%-2% or even less among all normal giants(Brown et al.
1989; Kumar et al. 2011; Liu et al. 2014; Casey et al. 2016; Kirby et al.
2016; Monaco et al. 2011; Li et al. 2018; Gao et al. 2019). In particular, the
ratios estimated from some large survey programs are $\sim$0.9% from Gaia-ESO
survey (Casey et al. 2016; Smiljanic et al. 2018), $\sim$0.8% from RAVE
(Ruchti et al. 2011) and $\sim$0.2%-0.3% from SDSS and GALAH data (Martell &
Shetrone 2013; Deepak & Reddy 2019).

In the past 40 years, various Li-rich giant stars have been successively
identified. Based on the Large Sky Area Multi-Object Fiber Spectroscopic
Telescope (LAMOST) low resolution spectra acquired in China: in 2018, a star
with the highest Li abundance was found,with A(Li) $\sim$ 4.5 dex(Yan et al.
2018). The Li-rich giants with A(Li) higher than 3.2 dex are called super Li-
rich giants, whose amount among all Li-rich giants is about 6%(Singh et al.
2021). In the period from October 2011 to June 2019, a total of 10,535 Li-rich
giants with A(Li) $\geq$ 1.5 dex were screened from the LAMOST low-resolution
spectra, which allowed one to expand the existing observation sample database
by about 5 times and greatly enriches research samples.

There are two main hypotheses about the origin of Li-rich giant stars: one is
that Li comes from the outside of stars, such as planetary engulfment or
pollution of binary companion stars (Stephan et al. 2018; Lodders 2019;
Holanda et al. 2020). The other is through
3He($\alpha$,$\gamma$)7Be(e,$\nu$)7Li, also known as Cameron-Fowler (CF)
mechanism (Cameron 1955; Cameron & Fowler 1971). In the latter process, 7Be
produced at a high-temperature region of H burning has to be quickly carried
away to a low-temperature region, such as the convective envelope, where it
will decay into 7Li, and this fresh 7Li will survive. As is predicted by
Schwab (2020) that there is enough 7Be in the H-burning shell prior to the
first helium (He) subflash.

Combining the Kepler and spectroscopy information, Yan et al. (2021) confirmed
that most Li-rich stars are red clump (RC), while a few of them are red giant
branch (RGB). In order to explain the Li enhancement in RC stars, Mori et al.
(2021) introduces a neutrino magnetic moment (NMM), which shows that 7Be
production becomes more active owing to the fact that the delay of He flash
makes themohaline mixing more effective when the NMM is excited. Thus Li was
produced at the tip of the red giant branch (TRGB). However, they did not
attempt to explain the origin of super Li-rich giants. At the same time, Kumar
has also used the GALAH DR2(Buder et al. 2018) and Gaia DR2(Gaia Collaboration
2018) data to confirm that the Li abundance of RC stars is 40 times higher
than those at TRGB(Kumar et al. 2020). Perhaps, this abnormality of Li
abundance might be due to a complex interaction between TRGB and RC phase. The
most significant event between these two stages is the occurrence of a He
flash, especially the first, strongest He subflash.

Recently, Schwab (2020) has reported that the He flash induced mixing links
the H-burning shell and the convective zone, and plenty of 7Be circulate
toward the cooler convective zone to turn into 7Li through the CF mechanism.
The model proposed by Schwab (2020) implies there is the possibility of very
high Li abundances (see Fig.4), although the author notes that these high Li
abundances are quickly depleted in their model. Therefore, searching for a
physical mechanism of Li enrichment (including super Li-rich giant stars) to
obtain the consistency between observation and theory is of great significance
for fundamental and applied research. It is well known that the surface
chemical abundance of a star can be affected by many factors such as
convection, thermohaline mixing or element diffusion(e. g., Kippenhahn et al.
1980; Dupuis et al. 1992; Zhu et al. 2021). The effects of convection and
thermohaline mixing on Li abundances of RGB stars have been investigated by
Yan et al. (2018) and Martell et al. (2021). However, the element diffusion is
seldom considered. Element diffusion is a dynamic process that changes the
distribution of chemical elements in stars. It is mainly the result of the
joint action of pressure, temperature, material concentration, and other
factors. Element diffusion plays a very important role in the stellar
evolution, especially in the chemical element distribution on the stellar
surface(Semenova et al. 2020).

In this paper, we consider the effects of the element diffusion impact on the
Li abundance for stars in the RC phase. In particular, Sect. 2 describes the
details of the stellar model and element diffusion. Sect. 3 presents the Li
abundances predicted by our model and the comparative analysis of the
observation and theoretical results. The summary is given in Sect. 4.

##  2 Models

We use Modules for Experiments in Stellar Astrophysics (MESA, [rev. 12778];
Paxton et al. (2011, 2013, 2015, 2018, 2019)) to construct one-dimensional
low-mass stellar models. MESA adopts the equation of state of Rogers &
Nayfonov (2002) and Timmes & Swesty (2000) and the opacity of Iglesias &
Rogers (1996, 1993) and Ferguson et al. (2005).

Our model use the standard MESA $pp\\_and\\_cno\\_extras$ nuclear network,
which includes 25 species and the reactions covering the pp-chains and CNO-
cycles. We adopt nuclear reaction rates compiled by JINA REACLIB (Cyburt et
al. 2010). Treatment of electron screening is based on Alastuey & Jancovici
(1978) and Itoh et al. (1979). The mass loss formula in Reimers (1975) is
adopted. We use the electron-capture rate on 7Be from Simonucci et al. (2013),
as made available in machine-readable form by Vescovi et al. (2019). Models
are initialized on the pre-main sequence with the Asplund et al. (2009) solar
abundance pattern and Z = 0.014. This initializes Li to the meteoritic
abundance A(Li) = 3.26 dex.

The size of convective zone depends on mixing length parameter. We adopt a
mixing length of 1.8 times the pressure scale height. These models include
thermohaline mixing($\alpha_{\rm th}$) using the Kippenhahn et al. (1980)
prescription with an efficiency of $\alpha_{\rm th}$ = 100. This gives the
deep mixing necessary to destroy 7Li on the first ascent giant branch.

Elemental diffusion in stars is mainly driven by a combination of pressure
gradients (or gravity), temperature gradients, compositional gradients, and
radiation pressures. The main driving factor of element diffusion is gravity
sedimentation(Paxton et al. 2015, 2018). In the model, we input a mixing
diffusive coefficient $D_{\rm mix}$$\sim$($\Delta$R)2/($\Delta$t) to show the
efficiency of element diffusion in different regions. Standard stellar
evolution theory points out that the H-burning core in the MS stage
continuously generates helium elements, which are deposited into the stellar
interior, so as to form a helium core. The He core reaches a certain mass, it
will begin to burn. We speculate that the influence of element diffusion will
affect the element abundance on the star surface, which will affect the
formation process of helium core, the occurrence time of the first He flash,
and then stimulate the thermohline mixing to make the mixing process more
sufficient.

By solving the Burgers equation (Burgers 1969), Thoul et al. (1994) proposed a
general method to arrange the whole set of equations into a single matrix
equation, that is to input Burgers equation into the matrix structure without
readjusting any number. There is no approximation of the relative
concentrations of various species, nor is there any limitation on the number
of elements to be considered. Therefore, this method is suitable for a wide
variety of astrophysical problems. Using the method of Thoul et al. (1994),
MESA can calculate the diffusion of chemical elements in stellar interior
(Paxton et al. 2015, 2018).

The inputs provided by the MESA model are the number densities $n_{\rm s}$,
temperature $T$, the gradients of these quantities $d$ ln $n_{\rm s}$/$dr$ and
$d$ ln $T$/$dr$, species mass in atomic units $A_{\rm s}$, species mean charge
as an average ionization state $\overline{Z}_{\rm s}$, and the resistance
coefficients $K_{\rm st}$, $z_{\rm st}$, $z_{\rm st}$′ and $z_{\rm st}$′′, as
defined by Equation (86) in Paxton et al. (2015). In our model, diffusion
coefficients $D_{\rm mix}$ derived from Paquette et al. (1986) and updated by
Stanton & Murillo (2016). By calculating the characteristic duration of the
first He flash and the characteristic length scale of the mixing region, this
suggest a effective mixing diffusion coefficient requires $D_{\rm mix}$ > 1010
cm 2 s-1(Schwab 2020). Together with the mean ionization states, these are key
parts of the input physics that determine the diffusion of all ions. The
additional acceleration terms $g_{\rm rad,s}$ of radiation suspension is set
to zero by default.

##  3 Result

Figure 1: The evolution of Li on the stellar surface for 1 $M_{\odot}$ model.
The x axis is the Li abundance. The y axis represents the logarithm of
luminosity. The lines indicates the evolution process of Li abundance.

According to Yan et al. (2021), the Li-rich RGB and RC stars have different
mass distributions with the peaks at about 1.7 $M_{\odot}$ and 1.2
$M_{\odot}$, respectively. In general, when the stellar mass at a zero age
main sequence (ZAMS) is lower than 0.9 $M_{\odot}$, the star hardly evolves
into a giant phase. Simultaneously, when the stellar mass is larger than 1.8
$M_{\odot}$, the temperature in the envelope of the star in the giant phase is
sufficiently high. As a result, beryllium (Be) elements produced by the
H-burning shell are quickly destroyed and cannot be brought to the surface of
the star. Therefore, we calculated the evolution of Li abundance for the
models with masses of 0.9, 1.0, 1.2, 1.4, 1.6, and 1.8 $M_{\odot}$. For
simplicity, a 1.0 $M_{\odot}$ model was taken as a sample.

###  3.1 Effect of element diffusion

The solid line in Fig. 1 shows the standard evolutionary track. It shows that
low mass star samples usually start from the MS turnoff with A(Li) $\sim$ 3.2
dex and suffer depletion through the first dredge-up. Then the star reaches
the RGB bump with a luminosity of 101.6 $L_{\odot}$, where it begins to
consume Li rapidly again(Iben 1967a; Charbonnel & Zahn 2007; Lattanzio et al.
2015). The black solid line represents the stellar evolution without element
diffusion at Stage II. It shows the luminosity drops sharply to the level of
RC, whilst maintaining the RGB tip A(Li)(Kumar et al. 2020). The black dotted
line denotes the model with element diffusion. A(Li) successfully increases
from -1 to 1.8 dex in Stage II. As a result, the element diffusion improves
the efficiency of thermohaline mixing, and its activity range gets connected
with the convective zone of stellar interior. Also, the 7Be in the H-burning
shell get transferred to the convective envelope, where they decay into 7Li
through the CF mechanism. This figure reveals that the increase of Li
abundance up to 1.8 dex can be realized by considering the element diffusion.

Almost all explanations proposed for Li-rich giants involve the He flash(Kumar
et al. 2020; Mori et al. 2021; Schwab 2020). Based on the standard model of
stellar structure and evolution, the He-core burning occurs when the mass of
the He core increases to $M_{\rm He}$ $\approx$ 0.45 $M_{\odot}$(Thomas 1967;
Bildsten et al. 2012). Therefore, the stellar evolution in this work was
divided into two stages: Stage I is from the MS stage to $M_{\rm He}$ = 0.45
$M_{\odot}$. The next range until the RC stage is called Stage II.

Figure 1 displays the evolution trajectory of Li abundance on the surface of a
1$M_{\odot}$ star without and with the element diffusion. As shown by the
dark-blue solid line, the evolution of Li abundance in both models was similar
during Stage I. In the MS phase, the Li abundance is kept constant. When a
low-mass star evolves into a red giant, it undergoes the first dredge-up.
After the first dredge-up, Li starts to rapidly decrease because the Li
elements are mixed up the stellar interior material by the envelope deepen
process and then are diluted. This depletion is due to dilution as the
envelope deepens and mixes up interior material heavily depleted in Li. Before
the RGB tip, the Li abundance continued to drop due to the thermohaline
mixing. From the ZAMS to the RGB tip, the Li abundance on the stellar surface
decreased by about 4 orders of magnitudes, which was consistent with the
results of Kumar et al. (2020).

Figure 2:  Profiles of element abundances (top left corner) and diffusion
coefficient (top right corner) ($D_{\rm mix}$) on the first He flash. The
solid lines are for model with element diffusion, and the dashed lines are for
models without element diffusion. The left and right panels at the bottom show
the relative changes between the mean molecular weights ($\mu_{\rm diff}$ and
$\mu$) with and without element diffusion, and between the mean molecular
weight gradients ($\nabla\mu_{\rm diff}$ and $\nabla\mu$) with and without
element diffusion, respectively.

The He flash emerges on the RGB tip. The Li abundance in the model without
element diffusion was falling during Stage II. However, A(Li) in the case of
element diffusion increased from -1 to 1.8 dex. It indicates that element
diffusion can enhance the Li abundance on the stellar surface. The main reason
for this phenomenon is depicted in Fig. 2. Especialy, a compact He core is
formed on the RGB tip, whose strong gravity can produce the efficient element
diffusion, resulting in the expansion of the thermohaline mixing zone. As
shown in the top-right panel of Fig. 2, for the model with element diffusion,
the thermohaline mixing zone connected to the convection zone in the stellar
interior in which the He flash occurred, that is, it extended to the deeper
interior of the hydrogen burning region. Therefore, the large amounts of 7Be,
produced by the H-burning shell, could be transferred to the stellar surface,
finally turning into 7Li. On the other hand, for the model without element
diffusion, 7Be could not or rarely be brought to the stellar envelope due to
the disconnection of the stellar interior convection zone and the thermohaline
mixing zone.

Thermohaline convection is a turbulent mixing process that can occur in
stellar radiative regions whenever the mean molecular weight increases with
radius. In some cases, it can have a significant observable impact on stellar
structure and evolution(Ulrich 1972; Kippenhahn et al. 1980; Brown et al.
2013; Garaud 2018). The left and right panels at the bottom of the Fig. 2 show
the relative changes of mean molecular weights ($\mu_{\rm diff}$ and $\mu$)
and mean molecular weight gradients ($\nabla\mu_{\rm diff}$ and $\nabla\mu$)
respectively at the junction zone between the stellar envelope and the He
core. It can be seen that due to the influence of element diffusion, the
$\mu_{\rm diff}$ and $\nabla\mu_{\rm diff}$ have decreased significantly. The
local decrease of mean molecular weight can drive a more efficient
thermohaline mixing. It expands to the inner convection region of stars, which
is shown by pink lines in the top-right panel. These indicate that element
diffusion can greatly affect the mean molecular weight and the mean molecular
weight gradient, which leads to the expansion of the mixing region of
thermohaline convection.

The essence of this phenomenon is that the element diffusion mainly triggered
by the gravity field suppress the the mean molecular weight and changes the
element concentration gradient at the junction zone between the stellar
envelope and the He core. It is one of the most important factors which affect
the thermohaline mixing. It is well known that the mean molecular weight
greatly affects the efficiency and range of thermohaline mixing, and make the
thermohaline mixing region expanding to the inner convection region of stars.
That is, the thermohaline mixing can occur in the internal area of hydrogen
burning, and bring products and by-products of nuclear reactions to the
surface. Thus, a transport channel, efficiently transporting 7Be in the
hydrogen burning region of the star to the convective envelope where 7Be
decays into 7Li, is formed.

###  3.2 Effect of element diffusion with constant diffusive coefficients

Although A(Li) on the surface of the star predicted by the model with element
diffusion can increase up to about 1.8 dex, it cannot explain the formation of
the super Li-rich giants.

Figure 3: Similar to Fig. 1, but for the evolutional tracks of Li abundances of four models. |
---|---
|
Figure 4: Similar to Fiure 1, but for the models with different masses and
constant diffusive mixing coefficient of $D_{\rm mix}$ to $10^{11}$ cm2 s-1,
$10^{12}$ cm2 s-1, $10^{13}$ cm2 s-1 and $10^{15}$ cm2 s-1. Different values
are displayed in the upper left corner of each subgraph. The solid and dashed
lines show the evolutional tracks at Stages I and II, respectively. The red
and blue circles are the Li-rich giant stars and RC stars listed in Table 1,
respectively. The star represent the most Li-rich giant star, TYC 429-2097-1
observed by Yan et al. (2018). Figure 5: The evolution of A(Li) after the
first He flash with time. The solid lines represents the increase produced
with element diffusion and different constant diffusive coefficient( $D_{\rm
mix}$). The dotted line indicates the A(Li) $\sim$ 1.5 dex and A(Li) $\sim$
3.2 dex.

According to Fig. 2, the diffusive coefficient is a very important factor for
the formation of the Li-rich giant. Very recently, in order to produce the Li-
rich giant, Schwab (2020) has considered that turbulent convective motions can
excite internal gravity waves, and a chemical mixing occurs when the
luminosity of the He flash ($L_{\rm He}$) is higher than $10^{4}$ L⊙. In
particular, the effective diffusive mixing coefficient was estimated to be
about $10^{11}$ cm2 s-1. Given the constant diffusive mixing coefficients
$D_{\rm mix}$ of $10^{10}$, $10^{12}$ and $10^{14}$ cm2 s-1 from different
models, the maximum value of A(Li) calculated by Schwab (2020) was about 3.6
dex. Meanwhile, the model proposed by Schwab (2020) fails to explain the
origin of super Li-rich giants.

In this section, combining the element diffusion and constant diffusive mixing
coefficients, we calculate the evolution of A(Li) on the stellar surface.
Following Schwab (2020), we assume a constant diffusive mixing coefficient for
the mixing when $L_{\rm He}$ ¿ $10^{4}$ L⊙. According to Fig. 2, $D_{\rm mix}$
during the He flash could reach $10^{15}$ cm2 s-1. Therefore, in this work,
$D_{\rm mix}=10^{11}$, $10^{12}$, $10^{13}$ and $10^{15}$ cm2 s-1 were adopted
in different models.

Figure 3 displays the A(Li) evolution on the surface of a star with the mass
of 1 $M_{\odot}$. Using the model with a constant diffusive mixing coefficient
of $10^{12}$ cm2 s-1 but without element diffusion could enhance the Li
abundance up to about 1.0 dex, which agreed with the result of Schwab (2020).
In the model with element diffusion, the Li abundance could be increased to
1.8 dex. At this time, the inner convection region was connected with the
thermohaline mixing zone to form a channel for transporting 7Be elements,
which greatly increased the 7Li on the stellar surface. Surprisingly, the
model combining the element diffusion and the constant diffusive mixing
coefficient exhibited an increase of A(Li) up to 3.4 dex. The reason for this
is that the diffusive mixing coefficient improves the mixing efficiency of the
channel excited by the element diffusion. Therefore, A(Li) in the model with
the element diffusion and the constant diffusive mixing coefficient can keep a
constant value above 3.0 dex.

###  3.3 Li-rich giants and super Li-rich giants

In recent years, many large survey programs have revealed the existence of
numerous Li-rich giants. Combining the astrometric data from the Gaia
satellite (Gaia Collaboration et al. 2016) with spectroscopic abundance
surveys (such as GALAH survey, LAMOST survey), twenty Li-rich abundances could
be identified. Based on GALAH DR2 and DR3 surveys, Deepak & Reddy (2019),
Deepak et al. (2020), Kumar et al. (2020) and Martell et al. (2021) measured
the Li abundances of 1872 giant stars. According to LAMOST survey, Singh et
al. (2019) and Yan et al. (2021) explored the Li abundances of 456 giant stars
which are in the Kepler field. In order to compare with the theoretical
results with observation samples in this work, we selected 351 published Li-
rich giants with precise values of luminosity, temperature and Li abundance as
our samples (see Table LABEL:tab1).

Figure 4 depicts the observed data on Li-rich giants and the theoretical
results for the models with different $D_{\rm mix}$s and masses. In this
study, an increase in $D_{\rm mix}$ from $10^{11}$ to $10^{15}$ cm2 s-1 led to
a rise in A(Li) from 2.4 to 4.5 dex. For the models with element diffusion and
$D_{\rm mix}>10^{12}$ cm2 s-1, the evolutionary tracks passed through most of
observed samples of super Li-rich giants (A(Li) $\geq$ 3.2 dex), Li-rich
giants (A(Li $\geq$ 1.5 dex) to the normal giants. Especially, the value of
A(Li), calculated in the model with element diffusion and $D_{\rm
mix}=10^{15}$ cm2 s-1, reached 4.5 dex, which could explain the Li abundance
of the most super Li-rich giants.

###  3.4 Population synthesis for Li-rich giants

As mentioned in the Introduction, the Li-rich giants among giants are scarce
(about 0.2-2%). Based on the models in this work, we estimate the theoretical
ratio by the population synthesis method which is used in the previous
investigations by our group(Lü et al. 2006, 2009, 2013, 2020; Yu et al. 2019,
2021; Zhu et al. 2021).

In the population synthesis method for single-star systems, the initial mass
function (IMF) is the most important input parameter. The IMF used in the
present research was derived from the stellar distribution toward both
Galactic poles as well as that within 5.2 pc of the Sun by Kroupa et al.
(1993). Based on this IMF, 106 stars were produced by Monte Carlo calculation.
In order to estimate their percentage, the lifetimes of Li-rich giants and
giants were afterward estimated.

Figure 5 displays the evolution of A(Li) in all models. After the first He
flash, the element diffusion firstly initiates and enhances the Li abundance.
According to Schwab (2020), the constant diffusive mixing coefficient can work
when $L_{\rm He}>$ 104 $L_{\odot}$. The Li abundance greatly increased after
about 0.2 Myr of the first He flash, and remained constant within several Myr
due to a constant $D_{\rm mix}$. Of course, A(Li) decreases when the stellar
luminosity is lower than 104 $L_{\odot}$. When the mass is greater than about
1.9 M⊙, the element diffusion and constant diffusive mixing coefficient are
not excited because the temperature within the stellar envelope is too high.
Simultaneously, the lifetimes of the Li-rich giants decrease with the increase
of $D_{\rm mix}$. The main reason for this is that the high diffusive
coefficient accelerates the circulation process of elements, so that Be or Li
elements can be quickly carried to a high-temperature zone and destroyed.

In this study, MESA was applied to calculate the evolution of stars with
initial masses of 0.9 M⊙, 1.0 M⊙, 1.2 M⊙, 1.4 M⊙, 1.6 M⊙ and 1.8 M⊙. Through a
linear interpolation method, the lifetimes of Li-rich giants and giants for
these 106 stars were estimated. Then, the percentage of Li-rich giants among
giants was assessed. In particular, the ratio values of models with $D_{\rm
mix}=10^{11}$, $10^{12}$, $10^{13}$ and $10^{15}$ cm2 s-1 were 0.5, 1.2, 1.1
and 0.2%, respectively, which was consistent with the observational estimates.

Very recently, Zhang et al. (2021) has shown that the deterioration of Li in
the RC stage is not so obvious and the low ratio (such as 0.2%) may be an
anomaly. If it is true, the high diffusive coefficient (such as $10^{15}$ cm2
s-1) may be undesirable. It means that the model proposed in this study fails
to produce all super Li-rich giants, especially the most Li-rich giant star
TYC 429-2097-1. Thus, the binary merging model proposed by Zhang et al. (2020)
may be competitive.

##  4 Conclusions

Considering the element diffusion, we used MESA to calculate the evolution of
Li abundance. The element diffusion mainly triggered by the gravity field
suppress the the mean molecular weight and changes the element concentration
gradient at the junction zone between the stellar envelope and the He core.
The local decrease of mean molecular weight greatly affects the efficiency and
range of thermohaline mixing, and make the thermohaline mixing region
expanding to the inner convection region of stars. A transport channel,
efficiently transporting 7Be in the hydrogen burning region of the star to the
convective envelope where 7Be decays into 7Li, is formed. Therefore, a large
amount of 7Be, produced by the He flash, could be transferred to the stellar
surface, finally turning into 7Li. However, the value of A(Li) could be
increased up to only 1.8 dex, which was insufficient to produce super Li-rich
giant stars.

In turn, combining the element diffusion and constant diffusive mixing
coefficients enabled one to increase the theoretical A(Li) values from 2.4 to
4.5 dex by increasing $D_{\rm mix}$ from $10^{11}$ to $10^{15}$ cm2 s-1. This
means that the element diffusion in the proposed model can result in the
extension of the thermohaline mixing zone and its connection with the stellar
interior convection zone. Then, 7Be produced by He burning can be mixed in the
stellar envelope. The high diffusive mixing coefficient can improve the
efficiency of 7Be transfer to the stellar surface. Therefore, our model can
produce the Li-rich giants, and even the most of super Li-rich giants. The
results provided by population synthesis method were also consistent with the
observations, which confirmed the feasibility of this mechanism. However, the
accuracy of the results in the model under consideration may be affected by
the uncertain input parameter, $D_{\rm mix}$. Since calculating an accurate
$D_{\rm mix}$ is beyond the scope of this work, attention is rather paid to
the diffusive mixing coefficients.

###### Acknowledgements.

We are grateful to anonymous referee for careful reading of the paper and
constructive criticism. This work received the generous support of the
National Natural Science Foundation of China, project Nos. U2031204, 11863005,
12163005, and 12090044, the science research grants from the China Manned
Space Project with NO. CMS-CSST-2021-A10, and the Natural Science Foundation
of Xinjiang No.2021D01C075.

Table 1: From about 11000 observational samples, the 351 Li-rich giant stars whose luminosities are measured are selected in this work. The observational data come from the references listed in the last column. (The full version of this Table is available at the CDS.) Object ID | Teff (K) | log g (dex) | [Fe/H] (dex) | log(L/$L_{\odot}$) | A(Li) | Reference
---|---|---|---|---|---|---
HD 8676 | 4860 | 2.95 | 0.02 | 1.68 | 3.55 | Kumar et al.(2011)
HD 10437 | 4830 | 2.85 | 0.1 | 1.77 | 3.48 | Kumar et al.(2011)
HD 12203 | 4870 | 2.65 | -0.27 | 1.69 | 2.08 | Kumar et al.(2011)
HD 37719 | 4650 | 2.4 | 0.09 | 1.76 | 2.71 | Kumar et al.(2011)
HD 40168 | 4800 | 2.5 | 0.1 | 2.1 | 1.7 | Kumar et al.(2011)
HD 51367 | 4650 | 2.55 | 0.2 | 1.59 | 2.6 | Kumar et al.(2011)
HD 77361 | 4580 | 2.35 | -0.02 | 1.66 | 3.8 | Kumar et al.(2011)
HD 88476 | 5100 | 3.1 | -0.01 | 1.87 | 2.21 | Kumar et al.(2011)
HD 107484 | 4640 | 2.5 | 0.18 | 1.78 | 2.14 | Kumar et al.(2011)
HD 118319 | 4700 | 2.2 | -0.25 | 1.68 | 2.02 | Kumar et al.(2011)
HD 133086 | 4940 | 2.98 | 0.02 | 1.7 | 2.14 | Kumar et al.(2011)
HD 145457 | 4850 | 2.75 | -0.08 | 1.61 | 2.49 | Kumar et al.(2011)
HD 150902 | 4690 | 2.55 | 0.09 | 1.83 | 2.65 | Kumar et al.(2011)
HD 167304 | 4860 | 2.95 | 0.18 | 1.93 | 2.85 | Kumar et al.(2011)
HD 170527 | 4810 | 2.85 | -0.1 | 1.69 | 3.12 | Kumar et al.(2011)
TYC 429-2097-1 | 4696 | 2.25 | -0.36 | 1.95 | 4.51 | Yan et al. (2018)
Gaia DR2 6423511482552457344 | 4828.68 | 2.84 | 0.18 | 1.56 | 3.54 | Deepak and Reddy(2019)
Gaia DR2 6216747182780840576 | 4773.08 | 2.69 | 0.12 | 1.54 | 3.41 | Deepak and Reddy(2019)
Gaia DR2 3080569351805501824 | 4995.53 | 2.6 | 0.03 | 1.71 | 3.41 | Deepak and Reddy(2019)
Gaia DR2 5920543908525756800 | 4815.52 | 2.68 | 0.14 | 1.54 | 3.39 | Deepak and Reddy(2019)
Gaia DR2 5676420200792553600 | 4854.1 | 2.31 | -0.11 | 1.83 | 3.38 | Deepak and Reddy(2019)
Gaia DR2 6721793108675117440 | 4911.04 | 2.45 | -0.04 | 1.64 | 3.33 | Deepak and Reddy(2019)
Gaia DR2 4488063566731544960 | 4778.94 | 2.37 | -0.02 | 1.52 | 3.27 | Deepak and Reddy(2019)
Gaia DR2 2939800046333110272 | 4985.2 | 2.56 | -0.13 | 1.55 | 3.26 | Deepak and Reddy(2019)
Gaia DR2 4168437628181576192 | 4749.88 | 2.71 | 0.19 | 1.21 | 3.26 | Deepak and Reddy(2019)
Gaia DR2 5229729170925959552 | 5038.38 | 2.79 | -0.15 | 1.65 | 3.24 | Deepak and Reddy(2019)
Gaia DR2 5293680581122445184 | 4832.25 | 2.5 | -0.01 | 1.65 | 3.23 | Deepak and Reddy(2019)
Gaia DR2 5242382659974594688 | 4786.39 | 2.8 | 0.28 | 1.6 | 3.23 | Deepak and Reddy(2019)
Gaia DR2 5628302754467688576 | 4868.94 | 2.37 | -0.2 | 1.57 | 3.21 | Deepak and Reddy(2019)
Gaia DR2 3202012502737830784 | 4906.32 | 2.42 | -0.14 | 1.77 | 3.21 | Deepak and Reddy(2019)
Gaia DR2 5460011229840058880 | 4813.14 | 2.63 | 0.16 | 1.7 | 3.21 | Deepak and Reddy(2019)
Gaia DR2 5452473905831060480 | 4711.6 | 2.17 | -0.3 | 1.84 | 3.2 | Deepak and Reddy(2019)
Gaia DR2 6162898261508964992 | 4835.99 | 2.58 | 0.01 | 1.7 | 3.2 | Deepak and Reddy(2019)
Gaia DR2 6779302244026689920 | 4541.3 | 2.17 | -0.48 | 1.99 | 3.2 | Deepak and Reddy(2019)
Gaia DR2 3496188144418768640 | 4776.14 | 2.59 | 0.04 | 1.63 | 3.2 | Deepak and Reddy(2019)
… | … | … | … | … | … | …

## References

  * Alastuey & Jancovici (1978) Alastuey, A. & Jancovici, B. 1978, ApJ, 226, 1034
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Bildsten et al. (2012) Bildsten, L., Paxton, B., Moore, K., & Macias, P. J. 2012, ApJ, 744, L6
  * Brown et al. (1989) Brown, J. A., Sneden, C., Lambert, D. L., & Dutchover, Edward, J. 1989, ApJS, 71, 293
  * Brown et al. (2013) Brown, J. M., Garaud, P., & Stellmach, S. 2013, The Astrophysical Journal, 768, 34
  * Buder et al. (2018) Buder, S., Asplund, M., Duong, L., et al. 2018, MNRAS, 478, 4513
  * Burgers (1969) Burgers, J. M. 1969, Flow Equations for Composite Gases
  * Cameron (1955) Cameron, A. G. W. 1955, ApJ, 121, 144
  * Cameron & Fowler (1971) Cameron, A. G. W. & Fowler, W. A. 1971, ApJ, 164, 111
  * Casey et al. (2016) Casey, A. R., Ruchti, G., Masseron, T., et al. 2016, MNRAS, 461, 3336
  * Charbonnel & Zahn (2007) Charbonnel, C. & Zahn, J. P. 2007, A&A, 467, L15
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240
  * Deepak et al. (2020) Deepak, Lambert, D. L., & Reddy, B. E. 2020, MNRAS, 494, 1348
  * Deepak & Reddy (2019) Deepak & Reddy, B. E. 2019, Monthly Notices of the Royal Astronomical Society, 484, 2000
  * Dupuis et al. (1992) Dupuis, J., Fontaine, G., Pelletier, C., & Wesemael, F. 1992, ApJS, 82, 505
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  * Gaia Collaboration (2018) Gaia Collaboration. 2018, VizieR Online Data Catalog, I/345
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., Brown, A. G. A., & Vallenari, A. 2016, A&A, 595, A1
  * Gao et al. (2019) Gao, Q., Shi, J.-R., Yan, H.-L., et al. 2019, The Astrophysical Journal Supplement Series, 245, 33
  * Garaud (2018) Garaud, P. 2018, Annual Review of Fluid Mechanics, 50, 275
  * Holanda et al. (2020) Holanda, N., Drake, N. A., & Pereira, C. B. 2020, AJ, 159, 9
  * Iben (1967a) Iben, Icko, J. 1967a, ApJ, 147, 650
  * Iben (1967b) Iben, Icko, J. 1967b, ApJ, 147, 624
  * Iglesias & Rogers (1993) Iglesias, C. A. & Rogers, F. J. 1993, ApJ, 412, 752
  * Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
  * Itoh et al. (1979) Itoh, N., Totsuji, H., Ichimaru, S., & Dewitt, H. E. 1979, ApJ, 234, 1079
  * Kippenhahn et al. (1980) Kippenhahn, R., Ruschenplatt, G., & Thomas, H. C. 1980, A&A, 91, 175
  * Kirby et al. (2016) Kirby, E. N., Guhathakurta, P., Zhang, A. J., et al. 2016, VizieR Online Data Catalog, J/ApJ/819/135
  * Kroupa et al. (1993) Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545
  * Kumar et al. (2020) Kumar, Y. B., Reddy, B. E., Campbell, S. W., et al. 2020, Nature Astronomy, 4, 1059
  * Kumar et al. (2011) Kumar, Y. B., Reddy, B. E., & Lambert, D. L. 2011, ApJ, 730, L12
  * Lattanzio et al. (2015) Lattanzio, J. C., Siess, L., Church, R. P., et al. 2015, MNRAS, 446, 2673
  * Li et al. (2018) Li, H., Aoki, W., Matsuno, T., et al. 2018, ApJ, 852, L31
  * Lind et al. (2009) Lind, K., Primas, F., Charbonnel, C., Grundahl, F., & Asplund, M. 2009, A&A, 503, 545
  * Liu et al. (2014) Liu, Y. J., Tan, K. F., Wang, L., et al. 2014, ApJ, 785, 94
  * Lodders (2019) Lodders, K. 2019, in Nuclei in the Cosmos XV, Vol. 219, 165–170
  * Lü et al. (2006) Lü, G., Yungelson, L., & Han, Z. 2006, MNRAS, 372, 1389
  * Lü et al. (2013) Lü, G., Zhu, C., & Podsiadlowski, P. 2013, ApJ, 768, 193
  * Lü et al. (2020) Lü, G., Zhu, C., Wang, Z., et al. 2020, ApJ, 890, 69
  * Lü et al. (2009) Lü, G., Zhu, C., Wang, Z., & Wang, N. 2009, MNRAS, 396, 1086
  * Martell & Shetrone (2013) Martell, S. L. & Shetrone, M. D. 2013, MNRAS, 430, 611
  * Martell et al. (2021) Martell, S. L., Simpson, J. D., Balasubramaniam, A. G., et al. 2021, MNRAS, 505, 5340
  * Monaco et al. (2011) Monaco, L., Villanova, S., Moni Bidin, C., et al. 2011, A&A, 529, A90
  * Mori et al. (2021) Mori, K., Kusakabe, M., Balantekin, A. B., Kajino, T., & Famiano, M. A. 2021, MNRAS, 503, 2746
  * Paquette et al. (1986) Paquette, C., Pelletier, C., Fontaine, G., & Michaud, G. 1986, ApJS, 61, 177
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Reimers (1975) Reimers, D. 1975, Memoires of the Societe Royale des Sciences de Liege, 8, 369
  * Rogers & Nayfonov (2002) Rogers, F. J. & Nayfonov, A. 2002, ApJ, 576, 1064
  * Ruchti et al. (2011) Ruchti, G. R., Fulbright, J. P., Wyse, R. F. G., et al. 2011, ApJ, 743, 107
  * Schwab (2020) Schwab, J. 2020, ApJ, 901, L18
  * Semenova et al. (2020) Semenova, E., Bergemann, M., Deal, M., et al. 2020, A&A, 643, A164
  * Simonucci et al. (2013) Simonucci, S., Taioli, S., Palmerini, S., & Busso, M. 2013, ApJ, 764, 118
  * Singh et al. (2019) Singh, R., Reddy, B. E., Bharat Kumar, Y., & Antia, H. M. 2019, ApJ, 878, L21
  * Singh et al. (2021) Singh, R., Reddy, B. E., Campbell, S. W., Kumar, Y. B., & Vrard, M. 2021, ApJ, 913, L4
  * Smiljanic et al. (2018) Smiljanic, R., Franciosini, E., Bragaglia, A., et al. 2018, A&A, 617, A4
  * Stanton & Murillo (2016) Stanton, L. G. & Murillo, M. S. 2016, Phys. Rev. E, 93, 043203
  * Stephan et al. (2018) Stephan, A. P., Naoz, S., & Gaudi, B. S. 2018, AJ, 156, 128
  * Thomas (1967) Thomas, H. C. 1967, in Late-Type Stars, ed. M. Hack, 395
  * Thoul et al. (1994) Thoul, A. A., Bahcall, J. N., & Loeb, A. 1994, ApJ, 421, 828
  * Timmes & Swesty (2000) Timmes, F. X. & Swesty, F. D. 2000, ApJS, 126, 501
  * Ulrich (1972) Ulrich, R. K. 1972, The Astrophysical Journal, 172, 165
  * Vescovi et al. (2019) Vescovi, D., Piersanti, L., Cristallo, S., et al. 2019, A&A, 623, A126
  * Wallerstein & Sneden (1982) Wallerstein, G. & Sneden, C. 1982, ApJ, 255, 577
  * Yan et al. (2018) Yan, H.-L., Shi, J.-R., Zhou, Y.-T., et al. 2018, Nature Astronomy, 2, 790
  * Yan et al. (2021) Yan, H.-L., Zhou, Y.-T., Zhang, X., et al. 2021, Nature Astronomy, 5, 86
  * Yu et al. (2019) Yu, J., Li, Z., Zhu, C., et al. 2019, ApJ, 885, 20
  * Yu et al. (2021) Yu, J., Zhang, X., & Lü, G. 2021, MNRAS, 504, 2670
  * Zhang et al. (2021) Zhang, J., Shi, J.-R., Yan, H.-L., et al. 2021, ApJ, 919, L3
  * Zhang et al. (2020) Zhang, X., Jeffery, C. S., Li, Y., & Bi, S. 2020, ApJ, 889, 33
  * Zhu et al. (2021) Zhu, C., Liu, H., Wang, Z., & Lü, G. 2021, A&A, 654, A57