# The Potential to Probe Solar Neutrino Physics with LiCl Water Solution

Wenhui Shao  Department of Engineering Physics, Tsinghua University, Beijing
100084, China  Key Laboratory of Particle & Radiation Imaging (Tsinghua
University), Ministry of Education, Beijing 100084, China  Center for High
Energy Physics, Tsinghua University, Beijing 100084, China  Weiran Xu111Now at
Laboratory for Nuclear Science, Massachusetts Institute of Technology, MA
02139, USA Department of Engineering Physics, Tsinghua University, Beijing
100084, China  Key Laboratory of Particle & Radiation Imaging (Tsinghua
University), Ministry of Education, Beijing 100084, China  Center for High
Energy Physics, Tsinghua University, Beijing 100084, China  Ye Liang
Department of Engineering Physics, Tsinghua University, Beijing 100084, China
Key Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry
of Education, Beijing 100084, China  Center for High Energy Physics, Tsinghua
University, Beijing 100084, China  Wentai Luo  School of Physical Sciences,
University of Chinese Academy of Sciences, 100049, China  Tong Xu  Department
of Engineering Physics, Tsinghua University, Beijing 100084, China  Key
Laboratory of Particle & Radiation Imaging (Tsinghua University), Ministry of
Education, Beijing 100084, China  Center for High Energy Physics, Tsinghua
University, Beijing 100084, China  Ming Qi  School of Physics, Nanjing
University, Nanjing 210093, China  Jialiang Zhang  School of Physics, Nanjing
University, Nanjing 210093, China  Benda Xu  Department of Engineering
Physics, Tsinghua University, Beijing 100084, China  Key Laboratory of
Particle & Radiation Imaging (Tsinghua University), Ministry of Education,
Beijing 100084, China  Center for High Energy Physics, Tsinghua University,
Beijing 100084, China  Zhe Wang222Corresponding author: wangzhe-
<EMAIL_ADDRESS>Department of Engineering Physics, Tsinghua
University, Beijing 100084, China  Key Laboratory of Particle & Radiation
Imaging (Tsinghua University), Ministry of Education, Beijing 100084, China
Center for High Energy Physics, Tsinghua University, Beijing 100084, China
Shaomin Chen  Department of Engineering Physics, Tsinghua University, Beijing
100084, China  Key Laboratory of Particle & Radiation Imaging (Tsinghua
University), Ministry of Education, Beijing 100084, China  Center for High
Energy Physics, Tsinghua University, Beijing 100084, China

###### Abstract

Lithium chloride water solution is a good option for solar neutrino detection.
The $\nu_{e}$ charged-current (CC) interaction cross-section on
$\rm{{}^{7}Li}$ is evaluated with new B(GT) experimental measurements. The
total CC interaction cross-section weighted by the solar 8B electron neutrino
spectrum is $3.759\times 10^{-42}~{}\rm{cm}^{2}$, which is about 60 times that
of the neutrino-electron elastic scattering process. The final state effective
kinetic energy after the CC interaction on $\rm{{}^{7}Li}$ directly reflects
the neutrino energy, which stands in sharp contrast to the plateau structure
of recoil electrons of the elastic scattering. With the high solubility of
LiCl of 74.5 g/100 g water at 10∘C and the high natural abundance of 92.41%,
the molarity of $\rm{{}^{7}Li}$ in water can reach 11 mol/L for safe operation
at room temperature. The CC event rate of $\nu_{e}$ on $\rm{{}^{7}Li}$ in the
LiCl water solution is comparable to that of neutrino-electron elastic
scattering. In addition, the $\nu_{e}$ CC interaction with the contained
$\rm{{}^{37}Cl}$ also contributes a few percent of the total CC event rate.
The contained $\rm{{}^{35}Cl}$ and $\rm{{}^{6}Li}$ also make a delay-
coincidence detection for electron antineutrinos possible. The
recrystallization method is found to be applicable for LiCl sample
purification. The measured attenuation length of $11\pm 1$ m at 430 nm shows
that the LiCl solution is practicable for a 10-m diameter detector for solar
neutrino detection. Clear advantages are found in studying the upturn effect
of solar neutrino oscillation, light sterile neutrinos, and Earth matter
effect. The sensitivities in discovering solar neutrino upturn and light
sterile neutrinos are shown.

Lithium-7, Lithium chloride, LiCl water solution, Solar neutrino upturn,
Sterile neutrino, Earth effect

##  I Introduction

The propagation of solar neutrinos has several special features. At high
energy, the $\nu_{e}$ survival probability is low and dominated by the matter
effect, i.e., the MSW effect Wolfenstein (1978); Mikheyev and Smirnov (1985),
while at low energy, the probability is high, and the flavor change occurs as
in vacuum. Between the high- and low-energy regions, there is a smooth
“upturn” of the survival probability. When arriving at the Earth, they are
decoherent mass eigenstates. The survival probability of $\nu_{e}$ arriving at
a terrestrial experiment is further modulated according to their path in the
Earth, and in the first order, it shows a day-night asymmetry. The upturn and
Earth matter effects are poorly constrained by the SNO Aharmim  _et al._
(2013); Tesic (2008) and Super-Kamiokande experiments Abe  _et al._ (2016);
SKT . In addition to confirming these theoretical predictions, future precise
solar neutrino experiments Abe  _et al._ (2018); Wurm  _et al._ (2012); Askins
_et al._ (2020); Beacom  _et al._ (2017) are expected to probe new physics.
The weakly mixed light sterile neutrino model can influence the upturn curve
of the $\nu_{e}$ survival probability and make it “dip and wiggle” in the
expected upturn region de Holanda and Smirnov (2004, 2011). Nonstandard
interaction (NSI) Friedland  _et al._ (2004) and light dark matter Lopes
(2020) are also interesting solar neutrino physics topics to investigate.

The charged-current (CC) interaction of $\nu_{e}$ on nuclei is most favorable
for such types of physics studies because the recoil electron energy is
strongly correlated with the incident neutrino energy. In radiochemical
neutrino experiments with the CC reactions, such as the Homestake experiment
Cleveland  _et al._ (1998) with $\rm{{}^{37}Cl}$, GALLEX/GNO Hampel  _et al._
(1999); Altmann  _et al._ (2005), and SAGE Abdurashitov  _et al._ (2009)
experiments with $\rm{{}^{71}Ga}$, the energy of the recoil electrons,
however, is not measured. The SNO experiment Aharmim  _et al._ (2013) measured
the solar neutrino oscillation with heavy water, and it is the only experiment
by now to give a real-time energy measurement of the CC interactions on
nuclei, i.e., deuterium. However, the CC signals are obscured by the products
of neutrino neutral-current (NC) interactions on deuterium because of their
similar energy, position, and angular distributions. Other experimental
attempts can also be found for $\rm{{}^{11}B}$ Kovacs  _et al._ (1990) and
$\rm{{}^{115}In}$ Pfeiffer  _et al._ (1978), as well as a few recent ideas for
solar neutrino detection methods with $\rm{{}^{116}Cd}$ Zuber (2003),
$\rm{{}^{71}Ga}$ Wang  _et al._ (2021), $\rm{{}^{131}Xe}$, $\rm{{}^{136}Xe}$
Haselschwardt  _et al._ (2020), and others Fujita  _et al._ (2021).

The feasibility of using lithium-7 as a target for neutrino detection was
recognized in references Bahcall (1964); Reines and Woods (1965); Kuzmin and
Zatsepin (1965); Bahcall (1969) in the 1960s. Experimentally, researchers have
proposed to detect the recoil electron signals in an electronic detector with
a water solution of lithium chloride Peak (1980) or to extract the final state
7Be signal in radiochemical methods Rowley (1978); Kopylov  _et al._ (2009).
Recently, the large cross-section is updated and highlighted again Haxton
(1996), and by a solar angle cut, the interaction products of $\nu_{e}$ with
$\rm{{}^{7}Li}$ can be distinguished from the neutrino elastic scattering on
electrons Alonso  _et al._ (2014); Aharmim  _et al._ (2013).

In this work, we discuss that the detection approach with LiCl water solution
is practical and efficient to explore energy-dependent solar neutrino physics.
In Sec. II, we present the detection channels of neutrinos in LiCl water
solution and, in particular, the CC process on $\rm{{}^{7}Li}$. We present a
cross-section estimation of $\nu_{e}$ CC interaction $\rm{{}^{7}Li}$ with new
experimental transition matrix element inputs. In Sec. III, we provide our
recent investigation of the properties of LiCl water solution and a compact
detector proposal with estimated detector performance. In section IV, we
clarify the advantage of LiCl water solution in measuring the solar neutrino
upturn effect, the search for light sterile neutrinos, and the study of the
Earth matter effect. The sensitivity of probing the upturn and sterile
neutrinos with the detector proposal is presented in section V. The paper
concludes in section VI.

##  II MeV neutrino detection in LiCl water solution

In this section, we introduce the detection channels of MeV neutrinos in LiCl
water solution. The detection of $\nu_{e}$ is discussed first and then
followed by the cross-section calculation of the CC interactions of $\nu_{e}$
on $\rm{{}^{7}Li}$. The result is compared with the CC interactions of
$\nu_{e}$ on $\rm{{}^{37}Cl}$ and the elastic scattering of $\nu_{x}$ on
electron, where $\nu_{x}$ represents $\nu_{e}$, $\nu_{\mu}$, and $\nu_{\tau}$.
The $\bar{\nu}_{e}$ detection is discussed in the end.

###  II.1 Detection of neutrinos

The dominant interactions of neutrinos of 1-20 MeV in LiCl water solution are
(1) the CC process of $\nu_{e}$ on $\rm{{}^{7}Li}$ (LiCC), (2) the CC process
of $\nu_{e}$ on $\rm{{}^{37}Cl}$ (ClCC), (3) the neutral-current process of
$\nu_{x}$ on $\rm{{}^{7}Li}$ (LiNC), and (4) the elastic scatter of $\nu_{x}$
on $e^{-}$ (Elas).

The LiCC process, as shown in Fig. 1, is

$\displaystyle\nu_{e}+{\rm{{}^{7}Li}}\to{\rm{{}^{7}Be}}+e^{-}(+\gamma).$ (1)

The interaction can go through the ground state of $\rm{{}^{7}Be}$ with a
threshold of 0.862 MeV, in which both the Fermi and superallowed Gamow-Teller
(GT) transitions are possible Fujita (2019); Fujita  _et al._ (2021). The
interaction can also go through the first excited state of $\rm{{}^{7}Be}$
with a threshold of 1.291 MeV, producing an extra 0.429 MeV deexcitation
$\gamma$, and the transition is also a superallowed GT transition Fujita
(2019); Fujita  _et al._ (2021).

Figure 1:  Energy levels of $\rm{{}^{7}Li}$ and $\rm{{}^{7}Be}$ involved in
the neutrino charged current and neutral current interactions on
$\rm{{}^{7}Li}$ Fujita (2019); Fujita  _et al._ (2021).

The neutrino energy threshold of the ClCC process is 0.814 MeV, which is very
close to the LiCC threshold. More details of the ClCC process can be found
elsewhere Bahcall (1989); Shimbara  _et al._ (2012).

The LiNC process happens as

$\displaystyle\nu_{x}+{\rm{{}^{7}Li}}\to\nu_{x}+{\rm{{}^{7}Li}}+\gamma.$ (2)

The energy of the emitted $\gamma$ is 0.478 MeV Fujita (2019); Fujita  _et
al._ (2021). Because it is very low, we skip the discussion involving the NC
process in this paper.

The neutrino-electron elastic scattering is

$\displaystyle\nu_{x}+e^{-}\to\nu_{x}+e^{-},$ (3)

which is one important solar neutrino detection channel Agostini  _et al._
(2020); Hosaka  _et al._ (2006).

The final state effective kinetic energy, $T$, of the LiCC process includes
the final state $e^{-}$ and deexcitation $\gamma$ kinetic energy. It is
calculated with the neutrino energy, $E_{\nu}$, as

$\displaystyle T=E_{\nu}-{\rm{0.862~{}MeV}}.$ (4)

Similarly, the final state effective kinetic energy of the ClCC process can be
calculated as

$\displaystyle T=E_{\nu}-{\rm{0.814~{}MeV}}.$ (5)

The $T$ of the CC processes represents the neutrino energy well, and this is
critical for energy-dependent physics studies, such as the upturn, sterile
neutrinos, and NSI effects.

For the Elas process in Eq. 3, the kinetic energy, $T$, of the recoil electron
shows a plateau structure Bahcall (1987), which smooths out many energy-
dependent features.

The $T$ of the LiCC, ClCC, and Elas processes is called kinetic energy in the
rest of the paper, and it is thought not to cause any confusion.

The angular distribution of the recoiling electrons of the LiCC and ClCC
processes, as pointed out in Bahcall (1989, 1964); Haxton (1996), is close to
a uniform distribution with respect to the incident neutrino direction. With
the contamination of the deexcitation gamma(s) of final states, the
reconstructed direction distribution can be even more uniform in a real
experiment. In contrast, the recoil electrons from the Elas process favor the
forward direction. Figure 2 shows the distribution of the reconstructed solar
angle, $\theta_{\rm{Sun}}$, which is the angle between the reconstructed final
particle direction and the solar direction calculated with the position of the
Sun. A uniform distribution is assumed for the CC processes, and a real
distribution for the Elas process is extracted from the Super-Kamiokande
Hosaka  _et al._ (2006) and SNO+ results SNO . Their difference can be used to
separate these two types of signals Aharmim  _et al._ (2013); Alonso  _et al._
(2014).

Figure 2:  Distribution of the reconstructed solar angles,
$\theta_{\rm{Sun}}$, with $T>$5 MeV. The distribution of $\nu_{x}$-elastic
scattering is extracted from the Super-Kamiokande Hosaka  _et al._ (2006) and
SNO+ SNO  results. The distribution of $\nu_{e}$-$\rm{{}^{7}Li}$ and
$\nu_{e}$-$\rm{{}^{37}Cl}$ CC processes is assumed to be uniform. The ratio of
the number of CC events to that of Elas events is set according to the last
column of Tab. 3.

###  II.2 $\rm{{}^{7}Li}$ charged-current cross-section estimation

The LiCC cross-section, $\sigma$, for a specific $E_{\nu}$ and an energy level
of $\rm{{}^{7}Be}$ is calculated according to Bahcall (1989, 1978, 1966);
Barinov  _et al._ (2018)

$\displaystyle\sigma$
$\displaystyle=\sigma_{0}\frac{\omega_{e}p_{e}}{2\pi\alpha Z}F(\omega_{e},Z),$
(6)

where $m_{e}$ is the electron mass, $\omega_{e}$ and $p_{e}$ are the electron
energy and momentum in unit of $m_{e}$, respectively, $\alpha$ is the fine
structure constant, $Z$ is the atomic number of $\rm{{}^{7}Be}$,
$F(\omega_{e},Z)$ is the Fermi function Bahcall (1966), and $\sigma_{0}$ is

$\displaystyle\sigma_{0}$ $\displaystyle=\frac{2\alpha
Zm_{e}^{2}}{\hbar^{4}}(G_{V}^{2}\langle
1\rangle^{2}+G_{A}^{2}\langle\sigma\rangle^{2}),$ (7)

where $G_{V}$ and $G_{A}$ are the vector and axial coupling constants,
respectively, and $\langle 1\rangle^{2}$ and $\langle\sigma\rangle^{2}$ are
the corresponding squares of the Fermi and Gamow-Teller transition matrix
elements, represented by B(F) and B(GT) in experimental measurements. The
B(GT) values for the ground and first excited states are measured to be 1.19
and 1.06, respectively Fujita (2019); Fujita  _et al._ (2021). The energy
$\omega_{e}$ is determined by

$\displaystyle\omega_{e}=\frac{E_{\nu}+[M(A,~{}Z-1)-M(A,~{}Z)]+m_{e}-\bar{E}_{ex}}{m_{e}},$
(8)

which depends on $E_{\nu}$, $m_{e}$, the final atomic mass, $M(A,~{}Z)$, the
initial atomic mass, $M(A,~{}Z-1)$, and the average excitation energy of the
final atom, $\bar{E}_{ex}$ Bahcall (1978). Then, the result is corrected for
the screening effect of atomic electrons Bahcall (1966). The cross-section is
calculated for the $\rm{{}^{7}Be}$ ground state and the first excited state.
The result as a function of $E_{\nu}$ is shown in Fig. 3, and the total cross-
section weighted by an undistorted 8B spectrum Bahcall  _et al._ (1996) is
$3.759\times 10^{-42}$ cm2. These quantities for the LiCC cross-section
calculation are tabulated in Tab. 1.

Figure 3:  Shown is the cross-section of $\nu_{e}$ charged-current interaction on $\rm{{}^{7}Li}$ as a function of $\nu_{e}$ energy, where the total cross-section (red solid), the contribution from the ground state of $\rm{{}^{7}Be}$ (red dashed), and the contribution from the first excited state (red dot-dashed) are shown separately. The total cross-section of $\nu_{e}$ charged-current interaction on $\rm{{}^{37}Cl}$ is also shown (blue solid). The sharp increase at about 5.8 MeV is caused by a Fermi transition to an excited level of $\rm{{}^{37}Ar}$ at 5 MeV. The elastic scattering cross-sections of $\nu_{e}$ on $e^{-}$ (black dashed) and $\nu_{\nu,\tau}$ (black dot-dashed) on $e^{-}$ are also overlaid for comparison. Table 1: Listed are the energy levels, $E_{l}$, of $\rm{{}^{7}Be}$ of the ground level to ground level (gs.-gs.) and ground level to first excited level (gs.-ex.) transitions for the $\nu_{e}$-$\rm{{}^{7}Li}$ CC process and the corresponding B(GT) and B(F) valuesFujita (2019); Fujita  _et al._ (2021) and the cross-sections weighted by an undistorted 8B neutrino spectrum Bahcall  _et al._ (1996). The total CC cross-section of $\nu_{e}$-$\rm{{}^{7}Li}$ and $\nu_{e}$-$\rm{{}^{37}Cl}$ and the $\nu_{e}$-$e^{-}$ elastic scattering cross-section are also shown for comparison. Channel |  $E_{l}$ (MeV) | B(GT) | B(F) | $\sigma(\rm{{}^{8}B})(10^{-42}~{}\rm{cm}^{2})$
---|---|---|---|---
$\rm{{}^{7}Li}$ gs.-gs. | 0 | 1.19 | 1 | 2.470
$\rm{{}^{7}Li}$ gs.-ex. | 0.429 | 1.06 |  | 1.289
$\rm{{}^{7}Li}$ total |  |  |  | 3.759
$\rm{{}^{37}Cl}$ total |  |  |  | 1.069
$e^{-}$ |  |  |  | 0.061

Following the same calculation procedure and with the BT strength input from
Shimbara  _et al._ (2012), the ClCC process cross-section is repeated. The
differential cross-section is shown in Fig. 3, and the total cross-section
weighted by an undistorted 8B spectrum Bahcall  _et al._ (1996) is
$1.069\times 10^{-42}$ cm2. The information is also tabulated in Tab. 1.

Our result for $\rm{{}^{37}Cl}$ is consistent with the result of $1.08\times
10^{-42}$ cm2 in Shimbara  _et al._ (2012) within 1%, where the same procedure
and input parameters are taken. However, our new result for $\rm{{}^{7}Li}$ is
about 7.5% higher than those calculated in Haxton (1996) and the difference
comes from B(GT) inputs. The $\nu_{e}-e$ and $\nu_{\mu,\tau}-e$ elastic cross-
sections Bahcall (1987) are also shown in Fig. 3 and Tab. 1 for comparison.
For the 8B neutrinos, the LiCC cross-section is about 3.5 times that of the
ClCC and about 60 times that of the Elas process.

A thorough uncertainty analysis is not carried out here, but considering the
discrepancy in the neutrino-gallium cross-section validation experiments
Hampel  _et al._ (1998); Abdurashitov  _et al._ (1999), we think, in the
future, a 2% uncertainty in the total LiCC cross-section is realistic.

###  II.3 Detection of $\bar{\nu}_{e}$

LiCl water solution is also convenient for $\bar{\nu}_{e}$ detection. The
water solution contains many hydrogens, i.e., free protons, which are the
target of the inverse-beta-decay process of $\bar{\nu}_{e}$,

$\displaystyle\bar{\nu}_{e}+p\to n+e^{+}.$ (9)

The neutron can be captured on $\rm{{}^{6}Li}$ or $\rm{{}^{35}Cl}$ to form a
delayed signal. The capture cross-sections are 940 barns and 44 barns for
$\rm{{}^{6}Li}$ and $\rm{{}^{35}Cl}$, respectively. The delayed signal is a
triton and an $\alpha$ for the capture on $\rm{{}^{6}Li}$ Ashenfelter  _et
al._ (2019) and several gammas with a total energy of 8.6 MeV for the capture
on $\rm{{}^{35}Cl}$ Aharmim  _et al._ (2013). With the natural abundance input
of the Li and Cl isotopes, the neutron capture probability on $\rm{{}^{35}Cl}$
is about 30%, and the capture probability on $\rm{{}^{6}Li}$ is about 70%. The
delayed coincidence is excellent in extracting the $\bar{\nu}_{e}$ signals and
suppressing backgrounds.

##  III Detector with LiCl water solution

In this section, we report the $\rm{{}^{7}Li}$, $\rm{{}^{37}Cl}$, and electron
molarities in LiCl water solution and the related LiCl purification and the
attenuation length of LiCl water solution. Then, we present a compact detector
setup and the corresponding properties, such as the energy resolution, angular
resolution, detection threshold, and fiducial volume.

###  III.1 $\rm{{}^{7}Li}$, $\rm{{}^{37}Cl}$, and electron molarities in LiCl
water solution

A LiCl water solution with high $\rm{{}^{7}Li}$ molarity can be easily
achieved at room temperature. The density of a saturated LiCl water solution
at room temperature is measured to be about 1.2 g/cm3 YeL . The solubility of
LiCl in water is rather high, i.e., 74.5 g/100 g water at 10∘C, and Tab. 2
shows the solubility at several temperatures che . The natural abundance of
$\rm{{}^{7}Li}$ is 92.41% ens , and the natural abundance of $\rm{{}^{37}Cl}$
is only 24.24%. Considering a detector running at room temperature of 20∘C, to
avoid precipitation, we can set the actual LiCl concentration to be the value
at 10∘C for safe operation. In such a LiCl water solution, the molarities of
$\rm{{}^{7}Li}$, $\rm{{}^{37}Cl}$, and electron are 11, 2.9, and 610 mol/L,
respectively, as shown in Tab. 3.

Table 2: LiCl solubility in water at several temperatures che . Temperature (∘C) | Solubility (g/100 g water)
---|---
0 | 68.3
10 | 74.5
20 | 83.2
40 | 89.4
60 | 98.8
80 | 112.3
Table 3: Molarity of $\rm{{}^{7}Li}$, $\rm{{}^{37}Cl}$, and $e^{-}$ in LiCl solution, in which the LiCl concentration is 74.5 g/100 g water, i.e., the saturation solubility at $10^{\circ}$C. The event rates for the charged-current interactions of $\nu_{e}$ on $\rm{{}^{7}Li}$, $\rm{{}^{37}Cl}$, and the elastic scatterings of $\nu_{x}$ on $e^{-}$ are also shown, where they are calculated with the undistorted $\rm{{}^{8}B}$ $\nu_{e}$ spectrum Bahcall  _et al._ (1996), oscillated spectrum and oscillated spectrum plus a $T>5$ MeV cut. | Molarity | Event rate | Event rate | Event rate
---|---|---|---|---
|  | No osci. | Osci. | Osci. & $>$5 MeV
| (mol/L) | (/100 ton-year) | (/100 ton-year) | (/100 ton-year)
$\rm{{}^{7}Li}$ | 11 | 305 | 101 | 87.3
$\rm{{}^{37}Cl}$ | 2.9 | 22.7 | 7.28 | 7.17
All CC |  | 328 | 108 | 94.4
$e^{-}$ | 610 | 271 | 124 | 34.5

###  III.2 LiCl purification

For the application in a neutrino detector, the purification of LiCl is a key
question. For a market sample, usually with a purity of $\geq$99%, filtration
is the first essential procedure to remove the dominant impurity. Observing
that the solubility of LiCl in water varies with temperature, we find that
LiCl can be recrystallized and purified by adjusting its saturated solution
temperature. After a round of filtration with a 0.2 $\mu$m membrane and a
round of recrystallization, the potassium, uranium, and thorium concentrations
can be suppressed. More details are given in YeL .

###  III.3 Attenuation length of saturated LiCl water solution

Our initial test results of the attenuation length and spectrum of a saturated
LiCl water solution are satisfactory for a compact neutrino detector. The
attenuation length is measured with an 80 cm long tube and an LED. The
emission spectrum of the LED peaks at 430 nm and spans from 375 to 550 nm. The
attenuation length of the sample is measured to be 11$\pm$1 m. More details of
the measurement, the emission spectrum, and the absorption spectrum are given
in YeL . We expect that this result can be further improved when higher
purification is achieved.

###  III.4 Neutrino detector proposal and property

With the above information, we see that LiCl water solution can be used for a
low background and compact detector. The detector structure is similar to the
SNO (SNO+) experiment Aharmim  _et al._ (2013). It has a spherical array of
photomultipliers (PMTs), a water buffer, an acrylic vessel, and LiCl water
solution in the center as the neutrino detection medium. The diameter of the
acrylic vessel is 10 meters, which matches the measured attenuation length. To
run the detector safely at 20∘C, we assume that the LiCl concentration is 74.5
g/100 g water (see the discussion in Sec. III.1). We assume a photocathode
coverage of 50% and a PMT photon detection efficiency (quantum efficiency
times collection efficiency) of 30%. Note that LiCl water solution is
corrosive to metal. Glass, acrylic, Teflon, or Teflon-lined containers or
tools are necessary. In this work, we focus on the neutrino signals from 8B
neutrinos. With the experience from SNO+ SNO  and Super-Kamiokande Abe  _et
al._ (2016) experiments, the energy resolution, angular resolution, fiducial
volume, and detection threshold are discussed below.

The Cherenkov light signals from charged particles are detected with the PMTs.
The energy resolution for particles follows the Poisson uncertainty of the
total number of detected photoelectrons (PEs). With the higher photocathode
coverage and PMT photon detection efficiency, the total photon detection
efficiency is about twice that of the Super-Kamiokande experiment Abe  _et
al._ (2016), so that a light yield of approximately 20 PE/MeV is expected.

Direction reconstruction is performed with the Cherenkov light. A resolution
of about 35 degrees (68% C.L.) has been achieved for electrons with energy
greater than 5 MeV at the SNO and Super-Kamiokande experiments Abe  _et al._
(2016); Aharmim  _et al._ (2013). However, for the solar neutrino Elas event
study, because the distribution of the reconstructed solar angle,
$\theta_{\rm{Sun}}$, is not a simple gauss distribution, the actual
distribution from the two experiments is extracted and used in this study,
which is shown in Fig. 2.

To reduce the radiative background from PMTs and detector structures, only a
central fiducial volume is available for physics studies. Because we focus on
signals with energy greater than 4 or 5 MeV, the central spherical volume with
an 8 m diameter is considered as the fiducial volume for the proposed 10-m
diameter acrylic vessel. The fiducial region is 268 m3 or 320 tons. The total
amount of $\rm{{}^{7}Li}$ in the fiducial volume is 2.95$\times 10^{6}$ moles
or 20.6 tons.

A 5 MeV detection threshold is realistic for a background free solar neutrino
study, as shown in the SNO+ result SNO . Given the better energy resolution in
this assumed detector to suppress low energy radioactive or instrumental
background, a 4 MeV threshold may also be possible.

We are also interested in a solution with some scintillation component, i.e.,
a water-based liquid scintillator, in which the scintillation light yield is
comparable with Cherenkov light Alonso  _et al._ (2014); Biller  _et al._
(2020); Guo  _et al._ (2019). This further enhances the light yield and energy
resolution and maintains a similar performance for the direction
reconstruction Land  _et al._ (2021); Luo .

##  IV Advantage in solar neutrino physics study

In this section, we demonstrate that the LiCl-based detection strategy has
advantages for energy-dependent neutrino physics studies. The candidate event
rates and spectra of the LiCC, ClCC, and Elas processes are explained first,
and then we look at the solar neutrino upturn issue, the light sterile
neutrino search, and the Earth effect.

###  IV.1 Candidate event rates

We focus on the 8B neutrino studies with a 4- or 5-MeV detection threshold, as
explained in section III.4. In this work, the total 8B flux is assumed to be
$4.59\times 10^{-14}/({\rm{cm^{2}s}})$ for the AGS09 low metallicity
prediction Serenelli  _et al._ (2011), so that it gives a conservative
statistical estimation for sensitivity studies. We adopt the 8B neutrino
spectrum prediction, $\Phi_{\nu}(E_{\nu})$, in Bahcall  _et al._ (1996), which
is shown in Fig. 4.

Figure 4:  Solar 8B electron neutrino energy spectrum with different
oscillation configurations (upper panel) and their corresponding survival
probability curves (lower panel). For the sterile neutrino mixing, $\alpha$ =
0.021 and $\Delta m^{2}_{01}$ = $1.56\times 10^{-5}$ eV2.

The differential energy spectra of all neutrino flavors at a terrestrial
detector are evaluated numerically. The calculation starts with the generation
of 8B neutrinos according to their spatial probability density function Num .
Then, their propagation and oscillation probability are estimated. Their
survival or appearance probability of $\nu_{e}$, $\nu_{\mu}$ and $\nu_{\tau}$
at a terrestrial detector is represented with $P_{ee}(E_{\nu})$,
$P_{e\mu}(E_{\nu})$ and $P_{e\tau}(E_{\nu})$, respectively. The calculation
methods for the three-active-neutrino propagation in the Sun, the case with
one sterile neutrino (3+1), and the neutrino propagation in the Earth are not
identical. Next, the three-active-neutrino propagation in the Sun is
introduced first. More detailed explanations for the (3+1) case and the
propagation in the Earth are given in the relevant sections.

For the three-active-neutrino propagation in the Sun, the following
calculation is carried out according to the MSW theory. The 8B neutrino
generation zone ($r<0.135R_{\odot}$, where $R_{\odot}$ is the radius of the
Sun) is divided into 60 shells. Neutrinos are generated according to the
predicted probability density in each shell Num . The initial fusion produced
$\nu_{e}$ flux is decomposed into mass eigenstates $\nu_{1}$, $\nu_{2}$, and
$\nu_{3}$ according to the local number density of electrons Num . With the
adiabatic assumption Wolfenstein (1978); Mikheyev and Smirnov (1985), these
mass eigenstates are recombined into flavor eigenstates at the solar surface
under vacuum condition. The neutrino oscillation parameters used are
$\theta_{12}=0.587$ Abe  _et al._ (2016), $\theta_{13}=0.148$ Adey  _et al._
(2018), $\theta_{23}=0.849$ Zyla  _et al._ (2020), $\Delta
m^{2}_{21}=7.49\times 10^{-5}~{}{\rm{eV}^{2}}$ Abe  _et al._ (2016), $\Delta
m^{2}_{31}=2.53\times 10^{-3}~{}{\rm{eV}^{2}}$ Zyla  _et al._ (2020). The
oscillated 8B neutrino energy spectrum and $P_{ee}(E_{\nu})$ are also shown in
Fig. 4.

The neutrino event rates and spectra of LiCC and ClCC are estimated according
to the cross-sections and the molarities of $\rm{{}^{7}Li}$ and
$\rm{{}^{37}Cl}$ in the LiCl solution. The event rate, $N_{\rm{LiCC}}$, as a
function of the kinetic energy (defined in Eq. 4) is calculated as

$\begin{split}N_{\rm{LiCC}}(T)=tN_{\rm{Li}}\sum_{i}^{\rm{levels}}\Phi_{\nu}(E_{\nu})\sigma_{\rm{Li-i}}(E_{\nu},T)P_{ee}(E_{\nu}),\end{split}$
(10)

where $\sigma_{\rm{Li-i}}(E_{\nu},T)$ is calculated in Eq. 4, 6, 7, and 8,
including all final state levels of $\rm{{}^{7}Be}$, $t$ is the data-taking
time, and $N_{\rm{Li}}$ is the number of target $\rm{{}^{7}Li}$ per unit LiCl
water solution as described in Sec. III.4. Similarly, we obtain the event rate
and spectrum, $N_{\rm{ClCC}}(T)$, for the ClCC process. The energy spectrum of
all CC events on nuclei is

$\begin{split}N_{\rm{CC}}(T)=N_{\rm{LiCC}}(T)+N_{\rm{ClCC}}(T).\end{split}$
(11)

Both the oscillated and undistorted $T$ spectra of 8B neutrino CC events are
shown in Fig. 5. The integrated LiCC, ClCC, and all CC rates with undistorted
8B neutrino spectrum, oscillated spectrum, and oscillated spectrum plus a
$T>5$ MeV cut are calculated and tabulated in Tab. 3. The ClCC event rate is
7% of the LiCC process.

Figure 5:  The final state effective kinetic energy spectra of the LiCC and
ClCC signal events of 8B neutrinos with different oscillation configurations
(upper panel) and the ratios of the oscillated spectra to the undisturbed
spectrum (lower panel). For the sterile neutrino mixing, $\alpha$ = 0.021 and
$\Delta m^{2}_{01}$ = $1.56\times 10^{-5}$ eV2. The exposure is set to 100-ton
LiCl solution $\times$ 1 data-taking year, and the LiCl concentration in the
water solution is assumed to be 74.5 g/100 g water.

The surviving $\nu_{e}$, $\nu_{\mu}$, and $\nu_{\tau}$ all scatter on
electrons. The kinetic energy spectrum, $N_{\rm{Elas}}(T)$, of recoil
electrons contains all three contributions, and it is

$\begin{split}N_{\rm{El}}(T)=&tN_{e}\int
dE_{\nu}\Phi_{\nu}(E_{\nu})\\{\sigma_{e}(E_{\nu},T)P_{ee}(E_{\nu})\\\
&+\sigma_{\mu,\tau}(E_{\nu},T)\left[1-P_{ee}(E_{\nu})\right]\\},\end{split}$
(12)

where $\sigma_{e}(E_{\nu},T)$ and $\sigma_{\mu,\tau}(E_{\nu},T)$ are the
differential scattering cross sections as a function of electron kinetic
energy for $\nu_{e}$ and $\nu_{\mu,\tau}$ Bahcall (1987), respectively, and
$N_{e}$ is the total number of target electrons per unit LiCl water solution
as described in Sec. III.4. Both the oscillated and undistorted $T$ spectra of
the 8B neutrinos of the Elas process are shown in Fig. 6 The integrated Elas
event rates with the undistorted 8B neutrino spectrum, oscillated spectrum,
and oscillated spectrum plus a $T>5$ MeV cut are calculated and tabulated in
Tab. 3. The ratio of the CC event rate to that of Elas events is 108:124 for
the oscillated spectrum, and it is further enhanced to 94:35 with the $T>5$
MeV cut. In Fig. 2, correspondingly, the ratio of the number of CC events to
the number of Elas events is set according to 94:35.

Figure 6:  The recoil electron kinetic energy spectra of the Elas signal
events of 8B neutrinos with different oscillation configurations (upper panel)
and the ratios of the oscillated spectra to the undisturbed spectrum (lower
panel). For the sterile neutrino mixing, $\alpha$ = 0.021 and $\Delta
m^{2}_{01}$ = $1.56\times 10^{-5}$ eV2. The exposure is set to 100-ton LiCl
solution $\times$ 1 data-taking year, and the LiCl concentration in the water
solution is assumed to be 74.5 g/100 g water.

In summary, as seen in Tab. 3, the total event rates of CC, mainly LiCC, and
Elas are comparable in the proposed LiCl water solution. With the $T>5$ MeV
cut, the CC events are dominant.

###  IV.2 Upturn effect

In the following discussion, we compare the signal strengths of the upturn
effect of the CC and Elas processes.

In Fig. 4 of the $P_{ee}(E_{\nu})$ of the MSW solution, a clear upturn can be
seen from the high energy to low energy. To obtain a specific quantitative
size of upturn, the relative difference between the survival probability at
4.862 MeV and 10,862 MeV is calculated as a figure of merit, and the result is

$\displaystyle\frac{P_{ee}(4.862~{}{\rm{MeV}})-P_{ee}(10.862~{}{\rm{MeV}})}{P_{ee}(10.862~{}{\rm{MeV}})}=23\%.$
(13)

Subtracting the $\rm{{}^{7}Li}$ interaction threshold (Eq. 4), they correspond
to 4 and 10 MeV for the LiCC final state effective kinetic energy $T$.

As a reference, a flat survival probability $P_{ee}$ is plotted in Fig. 4, in
which $P_{ee}$ is set to a constant of 0.33, and $P_{e\mu}+P_{e\tau}$ = 0.67.

With the CC interactions on $\rm{{}^{7}Li}$ or $\rm{{}^{37}Cl}$, the final
state effective kinetic energy spectrum is calculated with Eq. 11 and shown in
Fig. 5. We define $R_{\rm{CC}}(T)$ as

$\displaystyle
R_{\rm{CC}}(T)=\frac{N_{\rm{CC}}(T)~{}|~{}\rm{Osci}}{N_{\rm{CC}}(T)~{}|~{}\rm{No~{}Osci}},$
(14)

which is the ratio of the kinetic energy spectra under the oscillation
condition (Osci) to no oscillation (No Osci). The spectrum of $R_{\rm{CC}}(T)$
is plotted in the lower panel of Fig. 5. The relative difference between
$R_{\rm{CC}}(4~{}{\rm{MeV}})$ and $R_{\rm{CC}}(10~{}{\rm{MeV}})$ for the MSW
oscillation study is

$\displaystyle\rm{MSW:}~{}\frac{R_{\rm{CC}}(4~{}{\rm{MeV}})-R_{\rm{CC}}(10~{}{\rm{MeV}})}{R_{\rm{CC}}(10~{}{\rm{MeV}})}=23\%.$
(15)

This CC channel result has a consistent signal strength with the original
result in Eq. 13. For the flat $P_{ee}$ configuration, there is no interesting
feature in the $R_{\rm{CC}}(T)$ spectrum, as plotted in Fig. 5.

For the Elas process, the kinetic energy spectrum of the recoil electron is
calculated according to Eq. 12 and shown in Fig. 6. A similar ratio of
$R_{\rm{El}}(T)$ is defined as

$\displaystyle
R_{\rm{El}}(T)=\frac{N_{\rm{El}}(T)~{}|~{}\rm{Osci}}{N_{\rm{El}}(T)~{}|~{}\rm{No~{}Osci}},$
(16)

which is the ratio of the oscillated kinetic energy spectrum to the
undistorted spectrum. The spectrum of $R_{\rm{El}}(T)$ is plotted in the lower
panel of Fig. 6. The relative difference between $R_{\rm{El}}(4~{}{\rm{MeV}})$
and $R_{\rm{El}}(10~{}{\rm{MeV}})$ for the MSW oscillation study is

$\displaystyle{\rm{MSW}:}~{}\frac{R_{\rm{El}}(4~{}{\rm{MeV}})-R_{\rm{El}}(10~{}{\rm{MeV}})}{R_{\rm{El}}(10~{}{\rm{MeV}})}=9.6\%.$
(17)

However, we notice that in the Elas process, even for the flat survival
probability, there is a minor upturn in the ratio plot, which is shown in the
lower panel of Fig. 6. The relative difference

$\displaystyle{\rm{flat}:}~{}\frac{R_{\rm{El}}(4~{}{\rm{MeV}})-R_{\rm{El}}(10~{}{\rm{MeV}})}{R_{\rm{El}}(10~{}{\rm{MeV}})}=2.7\%.$
(18)

This is caused by the difference in the differential cross-section in
$\nu_{e}$-e and $\nu_{\mu,\tau}$-e. An extra $\nu_{\mu,\tau}$ contribution
appears to the low energy part of the recoil electron spectrum. Therefore, the
net signal strength with the Elas process is only 6.9%, which is much smaller
than the 23% of Eq. 13 and Eq. 15.

In summary, in such a LiCl detector, as described in section III, the strength
of the upturn signal of the CC processes on LiCl is much larger than that of
the Elas process.

###  IV.3 Sterile neutrino

The neutrino transition probability $P_{ee}(E_{\nu})$, $P_{e\mu}(E_{\nu})$,
and $P_{e\tau}(E_{\nu})$ calculation with the (3+1) situation de Holanda and
Smirnov (2004, 2011) is done in the following way.

We adopt the parameter convention of sterile neutrinos in de Holanda and
Smirnov (2004, 2011). One parameter, mixing angle $\alpha$, describes the
mixing between the sterile neutrino, $\nu_{s}$, and active neutrinos, and the
other parameter, mass squared difference $\Delta m^{2}_{01}$, is the mass
squared difference between $\nu_{0}$ and $\nu_{1}$, in which $\nu_{0}$ is
introduced along with $\nu_{s}$. With a small mixing angle $\alpha$, $\nu_{s}$
mixes weakly with active neutrinos, and $\nu_{s}$ almost coincides with
$\nu_{0}$.

The nonadiabatic situation must be considered for the propagation of (3+1)
neutrinos in the Sun for some parameter settings of $\alpha$ and $\Delta
m^{2}_{01}$. The 8B neutrino generation zone $r<0.135R_{\odot}$ Num  is
divided into $60\times 60\times 60$ ($x\times y\times z$) small cells, and 8B
neutrinos are generated according to the probability density function Num . We
use two different numerical methods to calculate the flavor transition
probability $P_{ee}(E_{\nu})$, $P_{e\mu}(E_{\nu})$, and $P_{e\tau}(E_{\nu})$.
One is the multislab method. The neutrino outgoing path is divided into many
slabs with a step size of $R_{\odot}/1000$. Besides the local number density
of electrons, the local number density of neutrons Num  is also considered for
the important neutral current process de Holanda and Smirnov (2004, 2011).
Each slab is assumed to have a uniform material with constant number
densities. Flavor and mass eigenstate transfer occurs at each slab interface.
The other method is a 4th-order Runge-Kutta method, which is used to solve the
propagation differential equations with $10^{6}$ steps. The probabilities
$P_{ee}(E_{\nu})$, $P_{e\mu}(E_{\nu})$, and $P_{e\tau}(E_{\nu})$ calculated by
the two methods are in good agreement in our test, as in Long (2008). The
multislab method is easier and faster to obtain a stable solution for this
calculation. No Earth matter effect is considered since the difference
introduced by the Earth effect is not significant for this study Tesic (2008).

Taking one set of sterile neutrino mixing parameters ($\alpha$ = 0.021 and
$\Delta m^{2}_{01}$ = $-1.56\times 10^{-5}$ eV2) as an example, the oscillated
neutrino spectrum and the $\nu_{e}$ survival probability $P_{ee}(E_{\nu})$ are
shown in Fig. 4. The kinetic energy spectrum of the CC processes and its ratio
to the undistorted spectrum are both shown in Fig. 5. The recoil electron
kinetic energy spectrum of the Elas process and its ratio to the undistorted
spectrum are both shown in Fig. 6.

With these comparisons, we see that the rich structure information in Fig. 4
is reserved in the CC process, as shown in Fig. 5; however, it is almost
smeared out in Fig. 6 with the Elas process. To distinguish the sterile
neutrinos with the CC channels, fewer signal statistics are needed.

###  IV.4 Earth matter effect

In this section, we follow the three-active-neutrino oscillation calculation
in Sec. 3 and then continue with the oscillation calculation of the Earth
matter effect. A multilayer Earth model Giunti  _et al._ (1998) is adopted,
and a multislab numerical calculation (Sec. IV.3) is implemented.

In the upper panel of Fig. 7, we show the $\nu_{e}$ survival probability
$P_{ee}$ as a function of the neutrino energy, $E_{\nu}$, and the cosine of
the solar angle, $\cos(\theta_{\rm{Sun}})$. For the neutrino going through the
Earth, $\cos(\theta_{\rm{Sun}})$ is negative. The pattern can be understood as
in Akhmedov  _et al._ (2004).

Figure 7:  Electron neutrino survival probability $P_{ee}$ as a function of
the neutrino energy and the cosine of the solar angle (upper). For the
neutrino going through the Earth, $\cos(\theta_{\rm{Sun}})$ is negative. Ratio
of the oscillated LiCl charged-current event rate to the undistorted LiCl
charged-current event rate, i.e., $R_{\rm{CC}}$, in Eq. 14 as a function of
the kinetic energy of the final state particles and the cosine of the solar
angle (middle). Ratio of the oscillated Elas event rate to the undistorted
Elas event rate, i.e., $R_{\rm{El}}$ in Eq. 16 as a function of the recoil
electron kinetic energy and the cosine of the solar angle (lower).

With the CC processes on LiCl, the kinetic energy spectrum is obtained for
each solar angle as done with Eq. 11. In the middle panel of Fig. 7, we show
the $R_{\rm{CC}}$ (defined in Eq. 14) as a function of the kinetic energy,
$T$, and $\cos(\theta_{\rm{Sun}})$.

The lower panel of Fig. 7 is the plot for $R_{\rm{El}}$ (defined in Eq. 16)
for the Elas process.

First, the original rich pattern in the neutrino plot (upper panel of Fig. 7)
is well repeated in the LiCl CC plot (middle panel) but is almost smeared out
in the Elas plot (lower panel). Second, the structures are mostly seen in the
4-12 MeV region of the neutrino plot (upper and middle panels), where the CC
kinetic energy spectrum has the most statistics, as seen in Fig. 5. There are
some residual structures in the Elas plot higher than 10 MeV (lower panel),
but the statistics of the Elas signals is low, as seen in Fig. 6. In
conclusion, the CC process is most helpful in distinguishing the Earth matter
effect.

##  V Sensitivity study for upturn and sterile neutrino

In this section, more realistic detector effects and signal selection criteria
are considered. The sensitivity of the upturn study and sterile neutrino
search is presented. The Earth matter effect study requires a different
detector setup, such as a higher energy resolution and higher signal
statistics, and the required detector configuration and sensitivity are given
in a separate paper.

###  V.1 Detector effects and signal selection criteria

With the detector proposal and expected property described in Sec. III.4, more
realistic predictions are made, and based on the predictions, many random
samples are generated for the following sensitivity studies.

A gauss energy smearing with an energy resolution of 20 PE/MeV,
$R(T,E_{\rm{rec}})$, is applied to the kinetic energy spectrum of
$N_{\rm{CC}}(T)$ (Eq. 11) and $N_{\rm{El}}(T)$ (Eq. 12).

$\begin{split}N_{\rm{CC}}(E_{\rm{rec}})=N_{\rm{CC}}(T)\otimes
R(T,E_{\rm{rec}}),\\\ N_{\rm{El}}(E_{\rm{rec}})=N_{\rm{El}}(T)\otimes
R(T,E_{\rm{rec}}),\end{split}$ (19)

where $E_{\rm{rec}}$ is the reconstructed energy, and the corresponding energy
spectra $N_{\rm{CC}}(E_{\rm{rec}})$ and $N_{\rm{El}}(E_{\rm{rec}})$ are
obtained for the CC and Elas processes, respectively. For the limited data-
taking time and target mass and a binned fitting later, each $E_{\rm{rec}}$
spectrum is divided into several 1 MeV bins, and the bin content in each bin
is denoted with $N_{\rm{CC},i}$ or $N_{\rm{El},i}$, where $i$ is the bin
number.

The CC and Elas events are the background to each other. As discussed in Sec.
III.4, we expect that there is no other background with an $E_{\rm{rec}}>$5
MeV cut or a more aggressive $E_{\rm{rec}}>$4 MeV cut.

The Elas and CC signals must be separated with a reconstructed solar angle cut
(see Fig. 2). Applying a solar angle cut at 60 degrees, i.e.,
$\cos(\theta_{\rm{Sun}})=0.5$, a CC-rich sample with 10% of the Elas events
and 75% of the CC events and an El-rich sample with 90% of the Elas events and
25% of the CC events are obtained. They are represented by

$\begin{split}N_{{\rm{CC\text{-}rich}},i}=N_{{\rm{CC}},i}75\%+N_{{\rm{El}},i}10\%,\\\
N_{{\rm{El\text{-}rich}},i}=N_{{\rm{CC}},i}25\%+N_{{\rm{El}},i}90\%.\end{split}$
(20)

The predicted $N_{\rm{CC\text{-}rich},i}$ and $N_{\rm{El\text{-}rich},i}$ are
shown in Fig. 8, and the LiCC, ClCC, and Elas components are also shown.

Figure 8:  Predicted CC-rich (upper) and El-rich (lower) samples as a function
of $E_{\rm{rec}}$ are shown. The contained LiCC, ClCC, and Elas components are
also overlaid. The exposure is set to 320-ton LiCl solution $\times$ 10 data-
taking years, and the LiCl concentration in the water solution is assumed to
be 74.5 g/100 g water. The detector resolution and signal selection
information can be seen in Sec. V.1.

###  V.2 Upturn effect

Poisson random sampling is performed according to the predictions in Eq. 20
for an exposure of 320$\times$5 ton-year. A total of 1000 random samples are
generated, and for each sample, there are two sub datasets,
$D_{{\rm{CC\text{-}rich}},i}$ and $D_{{\rm{El\text{-}rich}},i}$, for the CC-
rich and El-rich, respectively.

Each random sample is fitted with the following $\chi^{2}$

$\displaystyle\begin{split}\chi^{2}=&\sum_{i}(D_{{\rm{CC\text{-}rich}},i}-P_{{\rm{CC\text{-}rich}},i})^{2}/D_{{\rm{CC\text{-}rich}},i}\\\
&+\sum_{i}(D_{{\rm{El\text{-}rich}},i}-P_{{\rm{El\text{-}rich}},i})^{2}/E_{{\rm{El\text{-}rich}},i}\\\
&+{\rm{Pull}},\end{split}$ (21)

in which the predictions are $P_{{\rm{CC\text{-}rich}},i}$ and
$P_{{\rm{El\text{-}rich}},i}$ for the CC-rich and El-rich samples,
respectively. With respect to Eq. 20, they further consider the detector
systematic uncertainties as follows for the fit.

$\displaystyle\begin{split}P_{{\rm{CC\text{-}rich}},i}=&[N_{{\rm{CC}},i}75\%(1+\eta_{\sigma})+N_{{\rm{El}},i}10\%(1+\eta_{\epsilon})]\\\
&\times(1+\eta_{\rm{Norm}}),\\\
P_{{\rm{El\text{-}rich}},i}=&[N_{{\rm{CC}},i}25\%(1+\eta_{\sigma})+N_{{\rm{El}},i}(1-10\%(1+\eta_{\epsilon}))]\\\
&\times(1+\eta_{\rm{Norm}}),\\\ \end{split}$ (22)

where the LiCC cross-section uncertainty, $\eta_{\sigma}$, the selection
efficiency uncertainty of the Elas process, $\eta_{\epsilon}$, and the
normalization uncertainty, $\eta_{Norm}$ are taken into account. The 8B
neutrino spectrum theoretical error in the energy region of interest Roger
_et al._ (2012) is very small compared to the expected statistical error, and
the theoretical spectrum uncertainty is ignored in this study. The uncertainty
of the LiCC cross-section is assumed to be 2%, as suggested in Sec. II.2. The
Elas process selection efficiency uncertainty requires a dedicated detector
calibration and is assigned as 10%. The normalization uncertainty includes the
impact of the high-low metallicity ambiguity Serenelli  _et al._ (2011) and
fiducial mass and is assigned as 10%. The fit basically only considers the
shape of the observed spectrum. Correspondingly, the pull term is

$\displaystyle\begin{split}{\rm{Pull}}=(\eta_{\sigma}/2\%)^{2}+(\eta_{\epsilon}/10\%)^{2}+(\eta_{\rm{Norm}}/10\%)^{2}.\end{split}$
(23)

For the upturn study, we adopt two kinds of simplified $P_{ee}(E_{\nu})$ for
Eq. 22. The first one is a quadratic function of $E_{\nu}$ as used in the SNO
Aharmim  _et al._ (2013) and Super-Kamiokande Abe  _et al._ (2016)
experiments. Using the 10 MeV survival probability as the reference point, the
survival probability is

$\displaystyle\begin{split}P_{ee}(E_{\nu})&=c_{0}+c_{1}(E_{\nu}-10)+c_{2}(E_{\nu}-10)^{2}.\end{split}$
(24)

The term $c_{0}$ is the best fit of the average survival probability at 10
MeV, which is not the topic of this paper. To focus on the upturn issue, the
best fit and rms spread of $c_{1}(E_{\nu}-10)+c_{2}(E_{\nu}-10)^{2}$ for the
LiCl water solution detector with the exposure of 320$\times$5 ton-year are
shown in Fig. 9, where both the 5 MeV cut and 4 MeV cut results are analyzed.
The result from the SNO experiment is also overlaid for comparison. A
significant improvement with the compact LiCl detector can be expected.

Figure 9:  The best fit and rms spread of
$c_{1}(E_{\nu}-10)+c_{2}(E_{\nu}-10)^{2}$ of $P_{ee}(E_{\nu})$ in Eq. 24. The
upper panel is the fit result with a 5 MeV cut and the lower panel is with a 4
MeV cut. They present the evolution trend of $\nu_{e}$ survival probability as
a function of $E_{\nu}$ in the upturn sensitive region. The result from the
SNO experiment Aharmim  _et al._ (2013) is overlaid for comparison. The LiCl
water solution detector has an exposure of 320$\times$5 ton-year, where the
LiCl concentration is assumed to be 74.5 g/100 g water.

The sensitivity of observing the upturn is scanned with exposures of 320
tons$\times$1, 3, 5, 10, 15, and 20 years. One thousand data samples are
generated with each exposure setting. Because of the joint contribution of
$c_{1}$ and $c_{2}$, the upward or downward slope is difficult to quantify.
The $P_{ee}(E_{\nu})$ in Eq. 24 is replaced by a simple linear function as
below

$\displaystyle\begin{split}P_{ee}(E_{\nu})&=c_{0}+c_{1}(E_{\nu}-10).\end{split}$
(25)

If a statistically significant negative fit result of $c_{1}$ is obtained, an
upturn is observed. On the contrary, a zero or positive fit result of $c_{1}$
corresponds to no upturn effect, i.e., the null assumption. The 1000 $c_{1}$
fit results of each exposure setting are plotted. They are gauss distributed,
and the mean $\mu_{c1}$ and resolution $\sigma_{c1}$ are obtained. For our
proposed detector setup and exposures, $\mu_{c1}$ is always negative, i.e.,
the upturn is favored. Therefore, we define $\mu_{c1}/\sigma_{c1}$ as the
sensitivity for rejecting the null assumption. The sensitivity versus the
exposures is shown in Fig. 10, and the results with the 5 MeV and 4 MeV cuts
are both shown.

Figure 10:  The sensitivity of rejecting no upturn effect versus the number of
data-taking years for a 320-ton LiCl water solution. More detector resolution
and signal selection information can be seen in Sec. V.1.

In summary, we expect an improvement in studying the upturn effect with the
LiCl water solution detector.

###  V.3 Sterile neutrino

We use the Feldman-Cousin method Feldman and Cousins (1998) to determine the
exclusion sensitivity Agostini and Neumair (2020) of a detector with the
exposure of 320$\times$5 ton-year. The $\alpha$ and $\Delta m^{2}_{01}$
parameter space, $\log_{10}(\sin^{2}2\alpha)\in$ $[-5$, $-2]$ and $\Delta
m^{2}_{01}\in$ $[0$, $25\times 10^{-6}]~{}eV^{2}$, is split into 40$\times$40
grids. For each grid, 1000 statistically random samples are generated
according to the exposure configuration. The $P_{ee}(E_{\nu})$ calculation
procedure is described in Sec. IV.3, and the detector effect is added as in
Eq. 19 and 20 of Sec. V.1. With the $\chi^{2}$ definition in Eq. 21, 22 and
23, the exclusion sensitivity is estimated. Figure. 11 shows the sensitivity
contours with a 5 MeV cut. Besides the statistical power, the proposed
detector and settings are not sensitive to some sterile parameter regions due
to the 5 MeV cut and the systematic uncertainties.

Figure 11:  Exclusion sensitivity of sterile neutrinos using a LiCl detector
with an exposure of 320$\times$5 ton-year. More detector resolution and signal
selection information can be seen in Sec. V.1.

##  VI Conclusion

In this work, we study MeV neutrino detection in LiCl water solution. We
reevaluate the $\nu_{e}$ charged-current interaction cross-section on
$\rm{{}^{7}Li}$ with new B(GT) experimental inputs. The total CC interaction
cross-section weighted by the solar 8B electron neutrino spectrum is
$3.759\times 10^{-42}~{}\rm{cm}^{2}$, which is about 60 times the neutrino-
electron elastic scattering process. In addition, $\rm{{}^{37}Cl}$ also
contributes about seven percent to the CC event rate. The contained
$\rm{{}^{35}Cl}$ and $\rm{{}^{6}Li}$ also make the delay-coincidence detection
for electron antineutrinos possible. The detector with LiCl water solution is
basically a MeV-scale $\nu_{e}$ and $\bar{\nu}_{e}$ spectrometer.

We investigated the physical properties of LiCl and its water solution. LiCl
can be purified by recrystallization. The water solution has an attenuation
length of $11\pm 1$ m. A very high molarity of $\rm{{}^{7}Li}$ of 11 mol/L can
be achieved for operations at room temperature. The event rate of $\nu_{e}$ on
$\rm{{}^{7}Li}$ in a LiCl water solution is comparable to the elastic
scattering. A compact detector proposal is made, and its energy and angular
resolutions are estimated. The CC signals can be well separated from the
elastic signals by a solar angle cut.

The kinetic energy after the interaction on $\rm{{}^{7}Li}$ directly reflects
the neutrino energy. It demonstrates clear advantages in studying the solar
neutrino upturn effect, the light sterile neutrinos, and the Earth matter
effect. The sensitivities in studying the solar neutrino upturn and light
sterile neutrinos are reported for a detector with a 320-ton fiducial mass.

##  VII Acknowledgement

This work is supported in part by the National Natural Science Foundation of
China (Nos. 12141503 and 11620101004), the Ministry of Science and Technology
of China (No. 2018YFA0404102), the Key Laboratory of Particle & Radiation
Imaging (Tsinghua University), and the CAS Center for Excellence in Particle
Physics (CCEPP).

## References

  * Wolfenstein (1978) L. Wolfenstein, Phys. Rev. D 17, 2369 (1978).
  * Mikheyev and Smirnov (1985) S. P. Mikheyev and A. Y. Smirnov, Sov. J. Nucl. Phys. 42, 913 (1985).
  * Aharmim  _et al._ (2013) B. Aharmim _et al._ (SNO), Phys. Rev. C 88, 025501 (2013).
  * Tesic (2008) G. Tesic,  _Extraction of active and sterile neutrino mixing parameters with the Sudbury Neutrino Observatory_ , Doctor thesis, Carleton University (2008).
  * Abe  _et al._ (2016) K. Abe _et al._ (Super-Kamiokande), Phys. Rev. D 94, 052010 (2016).
  * (6) Recent results and future prospects from Super-Kamiokande, Presentation at the Neutrino 2020 conference.
  * Abe  _et al._ (2018) K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 [physics.ins-det] .
  * Wurm  _et al._ (2012) M. Wurm _et al._ (LENA), Astropart. Phys. 35, 685 (2012).
  * Askins  _et al._ (2020) M. Askins _et al._ (Theia), Eur. Phys. J. C 80, 416 (2020).
  * Beacom  _et al._ (2017) J. F. Beacom _et al._ (Jinping), Chin. Phys. C 41, 023002 (2017).
  * de Holanda and Smirnov (2004) P. C. de Holanda and A. Yu. Smirnov, Phys. Rev. D69, 113002 (2004).
  * de Holanda and Smirnov (2011) P. C. de Holanda and A. Yu. Smirnov, Phys. Rev. D83, 113011 (2011).
  * Friedland  _et al._ (2004) A. Friedland, C. Lunardini, and C. Pena-Garay, Phys. Lett. B 594, 347 (2004).
  * Lopes (2020) I. Lopes, Astrophys. J. 905, 22 (2020).
  * Cleveland  _et al._ (1998) B. T. Cleveland, T. Daily, R. Davis, Jr., J. R. Distel, K. Lande, C. K. Lee, P. S. Wildenhain, and J. Ullman, Astrophys. J. 496, 505 (1998).
  * Hampel  _et al._ (1999) W. Hampel _et al._ (GALLEX), Phys. Lett. B 447, 127 (1999).
  * Altmann  _et al._ (2005) M. Altmann _et al._ (GNO), Phys. Lett. B 616, 174 (2005).
  * Abdurashitov  _et al._ (2009) J. N. Abdurashitov _et al._ (SAGE), Phys. Rev. C 80, 015807 (2009).
  * Kovacs  _et al._ (1990) T. Kovacs _et al._ , Solar Phys. 128, 61 (1990).
  * Pfeiffer  _et al._ (1978) L. Pfeiffer, A. P. Mills, R. S. Raghavan, and E. A. Chandross, Phys. Rev. Lett. 41, 63 (1978).
  * Zuber (2003) K. Zuber, Phys. Lett. B571, 148 (2003).
  * Wang  _et al._ (2021) Z. Wang, B. Xu, and S. Chen, Astropart. Phys. 126, 102519 (2021).
  * Haselschwardt  _et al._ (2020) S. Haselschwardt, B. Lenardo, P. Pirinen, and J. Suhonen, Phys. Rev. D 102, 072009 (2020).
  * Fujita  _et al._ (2021) Y. Fujita, K. Zuber, and H. Fujita, Phys. Rev. D 104, 013004 (2021).
  * Bahcall (1964) J. N. Bahcall, Phys. Lett. 13, 332 (1964).
  * Reines and Woods (1965) F. Reines and R. M. Woods, Phys. Rev. Lett. 14, 20 (1965).
  * Kuzmin and Zatsepin (1965) V. A. Kuzmin and G. T. Zatsepin, Proceedings of the 9th International Cosmic Ray Conference 1, 1023 (1965).
  * Bahcall (1969) J. N. Bahcall, Phys. Rev. Lett. 23, 251 (1969).
  * Peak (1980) L. S. Peak, Austral. J. Phys. 33, 821 (1980).
  * Rowley (1978) J. K. Rowley, “The ${}^{7}\mathrm{Li}$-${}^{7}\mathrm{Be}$ experiment,” (Proceedings informal conference on status and future of solar neutrino research, NY, US, 1978).
  * Kopylov  _et al._ (2009) A. V. Kopylov, I. V. Orekhov, V. V. Petukhov, and A. E. Solomatin, “A lithium-beryllium method for the detection of solar neutrinos,” (2009), arXiv:0910.3889 [nucl-ex] .
  * Haxton (1996) W. C. Haxton, Phys. Rev. Lett. 76, 1562 (1996), [Erratum: Phys. Rev. Lett.77,1662(1996)].
  * Alonso  _et al._ (2014) J. R. Alonso _et al._ , “Advanced scintillator detector concept (asdc): A concept paper on the physics potential of water-based liquid scintillator,” (2014), arXiv:1409.5864 [physics.ins-det] .
  * Fujita (2019) Y. Fujita, “Nuclei as neutrino detectors,” (Proceedings, 5th International Solar Neutrino Conference, Dresden, Germany, 2019).
  * Bahcall (1989) J. N. Bahcall,  _Neutrino Astrophysics_ (Cambridge University Press, 1989).
  * Shimbara  _et al._ (2012) Y. Shimbara _et al._ , Phys. Rev. C 86, 024312 (2012), [Erratum: Phys.Rev.C 88, 059901 (2013)].
  * Agostini  _et al._ (2020) M. Agostini _et al._ (Borexino), Phys. Rev. D 101, 062001 (2020).
  * Hosaka  _et al._ (2006) J. Hosaka _et al._ (Super-Kamiokande), Phys. Rev. D 73, 112001 (2006).
  * Bahcall (1987) J. N. Bahcall, Rev. Mod. Phys. 59, 505 (1987).
  * (40) KamLAND-Zen and SNO+, Presentation at the Neutrino 2020 conference.
  * Bahcall (1978) J. N. Bahcall, Rev. Mod. Phys. 50, 881 (1978).
  * Bahcall (1966) J. N. Bahcall, Nuclear Physics 75, 10 (1966).
  * Barinov  _et al._ (2018) V. Barinov, B. Cleveland, V. Gavrin, D. Gorbunov, and T. Ibragimova, Phys. Rev. D97, 073001 (2018).
  * Bahcall  _et al._ (1996) J. N. Bahcall _et al._ , Phys. Rev. C 54, 411 (1996).
  * Hampel  _et al._ (1998) W. Hampel _et al._ , Physics Letters B 420, 114 (1998).
  * Abdurashitov  _et al._ (1999) J. N. Abdurashitov _et al._ (The SAGE Collaboration), Phys. Rev. C 59, 2246 (1999).
  * Ashenfelter  _et al._ (2019) J. Ashenfelter _et al._ (PROSPECT), Nucl. Instrum. Meth. A 922, 287 (2019).
  * (48) Ye Liang and others, to be published.
  * (49) From Chemister database as of Aug. 20, 2020. Version available at http://www.chemister.ru/Database/search-en.php.
  * (50) From ENSDF database as of Nov. 18, 2019. Version available at http://www.nndc.bnl.gov/ensarchivals/.
  * Biller  _et al._ (2020) S. D. Biller, E. J. Leming, and J. L. Paton, Nucl. Instrum. Meth. A 972, 164106 (2020).
  * Guo  _et al._ (2019) Z. Guo _et al._ , Astropart. Phys. 109, 33 (2019).
  * Land  _et al._ (2021) B. J. Land _et al._ , Phys. Rev. D 103, 052004 (2021).
  * (54) Wentai Luo and others, to be published.
  * Serenelli  _et al._ (2011) A. M. Serenelli, W. C. Haxton, and C. Pena-Garay, Astrophys. J. 743, 24 (2011).
  * (56) From John Bahcall’s webpage. Available at http://www.sns.ias.edu/ jnb/SNdata/sndata.html.
  * Adey  _et al._ (2018) D. Adey _et al._ (Daya Bay), Phys. Rev. Lett. 121, 241805 (2018).
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Long (2008) H. Long,  _Phenomenological studies of Active-Sterile Neutrino Oscillations_ , Doctor thesis, University of Science and Technology of China (2008).
  * Giunti  _et al._ (1998) C. Giunti, C. W. Kim, and M. Monteno, Nucl. Phys. B 521, 3 (1998).
  * Akhmedov  _et al._ (2004) E. K. Akhmedov, M. A. Tortola, and J. W. F. Valle, JHEP 05, 057 (2004).
  * Roger  _et al._ (2012) T. Roger _et al._ , Phys. Rev. Lett. 108, 162502 (2012).
  * Feldman and Cousins (1998) G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998).
  * Agostini and Neumair (2020) M. Agostini and B. Neumair, Eur. Phys. J. C 80, 750 (2020).


# Convective Penetration in Early-Type Stars

Adam S. Jermyn Center for Computational Astrophysics, Flatiron Institute, New
York, NY 10010, USA  Evan H. Anders CIERA, Northwestern University, Evanston,
IL 60201, USA  Daniel Lecoanet CIERA, Northwestern University, Evanston, IL
60201, USA  Department of Engineering Sciences and Applied Mathematics,
Northwestern University, Evanston IL 60208, USA  Matteo Cantiello Center for
Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA
Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA

###### Abstract

Observations indicate that the convective cores of stars must ingest a
substantial amount of material from the overlying radiative zone, but the
extent of this mixing and the mechanism which causes it remain uncertain.
Recently, Anders et al. (2021) developed a theory of convective penetration
and calibrated it with 3D numerical hydrodynamics simulations. Here we employ
that theory to predict the extent of convective boundary mixing in early-type
main-sequence stars. We find that convective penetration produces enough
mixing to explain core masses inferred from asteroseismology and eclipsing
binary studies, and matches observed trends in mass and age. While there are
remaining uncertainties in the theory, this agreement suggests that most
convective boundary mixing in early-type main-sequence stars arises from
convective penetration. Finally, we provide a fitting formula for the extent
of core convective penetration for main-sequence stars in the mass range from
$1.1-60M_{\odot}$.

Stellar physics (1621); Stellar evolutionary models (2046); Stellar convection
zones (301)

††journal: ApJ††software:  MESA (Paxton et al., 2011, 2013, 2015, 2018, 2019,
http://mesa.sourceforge.net), MESASDK 20190830 (Townsend, 2019a, b),
matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011)

##  1 Introduction

Observations tell us that stars have more mixing near the boundaries of
convection zones than theoretical models predict. For instance
asteroseismology of B stars reveals that their convective cores are larger
than predicted by 1D evolutionary models by tens of percent (Pedersen et al.,
2021), and a similar conclusion comes from modelling of early type eclipsing
binary systems (Claret & Torres, 2019). Likewise observations of lithium
depletion in F stars hint at extra mixing beneath surface convection zones
(Balachandran, 1995), as do measurements of the location of the base of the
solar convection zone from helioseismology (Basu, 2016, Sct. 7.2.1).

The extent and profile of this extra mixing matters critically for our
understanding of stellar evolution. For instance, the lifetimes and end points
of massive star evolution depend directly on how much fresh fuel is mixed into
their cores (Brott et al., 2011; Salaris & Cassisi, 2017; Farmer et al.,
2019). Likewise the efficiency of processes like dredge up (Serrano, 1986) and
the position of the Red Giant Branch luminosity bump (Khan et al., 2018)
depend sensitively on convective boundary mixing (CBM).

These challenges and opportunities have motivated a great deal of theoretical
work identifying and characterizing mechanisms for CBM. These include
convective overshooting and penetration (Zahn, 1991; Brummell et al., 2002;
Korre et al., 2019; Shaviv & Salpeter, 1973; Herwig, 2000) as well as mixing
by internal gravity waves launched by convection (Garcia Lopez & Spruit, 1991;
Rogers & McElwaine, 2017).

Recently, Anders et al. (2021) proposed a theoretical model of convective
penetration and calibrated this to 3D hydrodynamical simulations of Boussinesq
convection. Here we study the implications of the resulting calibrated model
in early-type stars. We begin in Section 2 by reviewing the theory of
convective overshoot and penetration and presenting the model. We then explain
how we implement this model in Section 3. We show predictions of this model
for the cores of early-type stars in Section 4 and compare these with
observations. These predictions were produced by _post-processing_ stellar
evolution models run with no convective penetration, to better show the
instantaneous extent of this mixing rather than its cumulative effect over
evolutionary time-scales. We find that convective penetration produces enough
mixing to explain core masses inferred from observations, and matches the
observed trends in mass and age. We then discuss uncertainties and future work
in Section 5, and conclude in Section 6.

Figure 1: The structure of a convective boundary as found by Anders et al.
(2021) is shown schematically. The top of the domain is a radiative zone,
where $\nabla=\nabla_{\rm rad}$. Below this is the penetration zone, which has
$\nabla\approx\nabla_{\rm ad}$. Motions in the penetration zone which hit the
boundary with the radiative zone experience a large buoyant restoring force,
and so this mechanical overshooting only reaches a small distance into the
radiative zone. Beneath the penetration zone is the convection zone, which is
slightly superadiabatic. Because the penetration zone is nearly adiabatic and
hence nearly neutrally buoyant, motions which start in the convection zone
experience a very weak restoring force in traversing the penetration zone.
This allows the penetration zone to remain thermally well-mixed.

##  2 Theory

In 1D stellar evolution calculations, the boundary of a convection zone (CZ)
is typically set at the point where the Brunt-Väisälä frequency $N$ vanishes.
At this point local fluid motions feel no net gravitational acceleration. In a
system with uniform composition this is the Schwarzschild boundary, where

$\displaystyle\nabla_{\rm rad}=\nabla_{\rm ad}.$ (1)

Here

$\displaystyle\nabla_{\rm ad}\equiv\left.\frac{\partial\ln T}{\partial\ln
P}\right|_{s}$ (2)

is the adiabatic temperature gradient and $s$ is the entropy,

$\displaystyle\nabla_{\rm rad}=\frac{3\kappa\rho LP}{64\pi GM\sigma T^{4}}$
(3)

is the radiative temperature gradient in the diffusive approximation, $\kappa$
is the opacity, $\rho$ is the density, $L$ is the luminosity, $P$ is the
pressure, $T$ is the temperature, $G$ is the gravitational constant, $M$ is
the mass below the point of interest, and $\sigma$ is the Stefan-Boltzmann
constant. In the presence of composition gradients, the boundary is instead
where the Ledoux criterion is marginally satisfied, such that

$\displaystyle\nabla_{\rm rad}=\nabla_{\rm ad}+\nabla_{\mu}$ (4)

and for an ideal gas

$\displaystyle\nabla_{\mu}\equiv\frac{d\ln\mu}{d\ln P},$ (5)

where $\mu$ is the mean molecular weight.

Convective motions carry momentum, however, and so can extend past the surface
of neutral buoyancy, producing overshoot. After moving beyond the convective
boundary, eddies experience a buoyant restoring acceleration $\sim hN^{2}$,
where $N$ is evaluated in the stable layer. Here

$\displaystyle h\equiv-\frac{dr}{d\ln P}$ (6)

is the pressure scale height.

An eddy with velocity $v_{c}$ can then move past the convective boundary for a
time $\Delta t\sim v_{c}/hN^{2}$, during which time it moves a distance

$\displaystyle\Delta r\sim v_{c}\Delta t\sim\frac{v_{c}^{2}}{hN^{2}},$ (7)

where $N$ is evaluated in the stable layer, typically a distance $\sim h$ away
from the convective boundary, $v_{c}$ is the typical velocity of convective
motions, and $r$ is the radial coordinate.

The overshooting given by equation (7) is called mechanical overshoot and is
often very small, typically $\Delta r/h\lesssim 10^{-6}$ for convective cores
on the main sequence. More sophisticated models based on the same momentum
picture provide larger values (Korre et al., 2019), but these are generally
still much too small to match the extent of convective boundary mixing
inferred from observations (Pedersen et al., 2021).

In addition to carrying momentum, though, convective motions also carry heat.
This allows convective overshooting to alter the thermal structure of the star
near the point of neutral buoyancy and ultimately form a _layer_ of neutral
buoyancy, known as a penetration zone (PZ). The PZ is an adiabat, with
temperature gradient $\nabla\approx\nabla_{\rm ad}$, and so material can flow
freely through it unhindered by buoyancy. Overshooting still occurs, but
begins at the outer edge of the PZ, rather than the edge of the Schwarzschild-
unstable region (Figure 1). This process is known as convective penetration,
and can result in more extended mixing layers than mechanical overshooting.

Building on Zahn (1991) and Roxburgh (1989), Anders et al. (2021) developed a
theory of convective penetration to predict the extent of the PZ and
calibrated that theory to 3D cartesian hydrodynamical simulations of
Boussinesq convection111These simulations were conducted over a thermal time-
scale ($\sim 10^{4}$ turnover times) because that is the time over which a
statistically stationary penetration zone develops. This theory is based on an
exact rewriting of the time-and-spatially averaged energy equation, accounting
for all energy fluxes, sources, and sinks. Several of the terms arising can be
extracted from 1D stellar models, while others must be calibrated from
numerical simulations. All told, the extent of the penetration zone is then
given by the integral constraint

$-\frac{\int_{\rm{PZ}}L_{\rm{conv}}\,dr}{\int_{\rm{CZ}}L_{\rm{conv}}\,dr}+f\xi\frac{V_{\rm{PZ}}}{V_{\rm{CZ}}}=(1-f),$
(8)

where

$f=0.86\qquad\rm{and}\qquad\xi=0.6$ (9)

are calibrated from numerical simulations, $L_{\rm conv}$ is the convective
luminosity, $V_{\rm PZ}$ is the volume of the penetration zone, and $V_{\rm
CZ}$ is the volume of the convection zone. The convective luminosity in the PZ
is given in thermal equilibrium by

$\displaystyle L_{\rm conv,PZ}=-L\left(\frac{\nabla_{\rm ad}-\nabla_{\rm
rad}}{\nabla_{\rm rad}}\right).$ (10)

We emphasize that the calibration in equations (9) and our theory more
generally are only valid in the Boussinesq limit, and more general models must
be used in strongly stratified systems (Roxburgh, 1989).

While Anders et al. (2021) did not include the effects of composition
gradients, we do not expect these to change the equilibrium extent of the
penetration zone. Overshooting motions still occur when there is a composition
gradient at the convective boundary. These motions serve to erase the
composition gradient by mixing material in the overshooting region. Over time
this allows overshooting to erase the composition gradient and build a
penetration zone in its wake.

Note that in terms of compositional mixing convective penetration most
naturally matches the step overshooting prescription used in stellar evolution
software instruments, because motions in the penetration zone are nearly as
fast as those in the convection zone, and the boundary between the penetration
zone and the overlying stable layer is sharp (Anders et al., 2021). In
addition, however, convective penetration mixes entropy to produce an
adiabatic layer. This entropy mixing is not typically implemented as part of
step overshooting and is not standard in any stellar evolution software
instrument of which we are aware, though it has been implemented in individual
calculations (Pedersen et al., 2021).

##  3 Methods

We calculated stellar evolutionary tracks for stars ranging from
$1.1-60M_{\odot}$ using revision 15140 of the Modules for Experiments in
Stellar Astrophysics software instrument (MESA Paxton et al., 2011, 2013,
2015, 2018, 2019). Details on the MESA microphysics inputs are provided in
Appendix A. Our models use the Milky Way metallicity $Z=0.014$ and determine
convective stability using the Ledoux criterion. We further employ convective
premixing (Paxton et al., 2019, Section 5.2) to more accurately model the
evolution of the convective boundary in the presence of composition gradients.

Using equations (8) and (10) we calculated the extent of the penetration zone
just outside the convective core. We emphasize that this calculation is done
in post-processing: all stellar evolution calculations were performed with no
convective penetration or additional boundary mixing. This likely makes a
difference to the evolution of a model: as convective penetration grows the
convective core it changes the temperature profile, which in turn can change
the extent of convective penetration, particularly if the edge of the core
gets pushed into a region with a very different opacity profile.

Our motivation in doing this calculation in post-processing was just to
understand the rough scale and trends involved in convective penetration. To
enable precision comparisons with observations, future work should perform
this calculation self-consistently with stellar evolution.

Along similar lines, because our models have a composition gradient in the
vicinity of the convective boundary, and because that gradient would be erased
in the presence of a penetration zone, we calculate $\nabla_{\rm rad}$ and
$\nabla_{\rm ad}$ for purposes of equations (8) and (10) using the _core
composition_ everywhere. We likewise use these constant-composition
temperature gradients to determine the location of the nominal convective
boundary (i.e. the base of the penetration zone), which may differ from that
reported by MESA due to e.g. convective premixing and compositional gradients.

Our approach of using the core composition for the penetration zone is not
quite right, because penetrative motions mix material of a slightly different
composition into the core, but the amount of this mixing is small enough that
to first order we may use the core composition without any penetration to
determine the extent of the penetration zone. Improving on this approximation
again requires performing convective penetration self-consistently with
stellar evolution, which we leave for the future.

##  4 Results

###  4.1 Radial Extent

Figure 2: The radial extent of the penetration zone is shown on a Hertzsprung-
Russel diagram in terms of $\log T_{\rm eff}$ and $\log L$ for stellar models
ranging from $1.1-60M_{\odot}$ with Milky Way metallicity $Z=0.014$. (Left)
Normalized to the pressure scale height at the nominal convective boundary
reported by MESA. (Right) Normalized to the radius at the same. Values below
the lower edge of the scale occur at masses $M\lesssim 1.5M_{\odot}$ as the
core convection zone is just forming.

Figure 2 shows the radial extent of the penetration zone on a Hertzsprung-
Russel diagram, normalized to the pressure scale height $h$ (left) and the
radius (right), both at the nominal boundary of the convective core. The
quantity reported on the left is equivalent to the commonly-used step
overshooting parameter $\alpha_{\rm ov}$, which gives the step overshooting
distance as a fraction of a pressure scale-height.

We see clear trends in mass for $\alpha_{\rm ov}$ (left). Higher mass stars
show much larger penetration zones relative to the pressure scale height. We
see mixing of $\alpha_{\rm ov}\sim 0.1$ for A-stars, $\alpha_{\rm ov}\sim 0.3$
for B-stars, and up to $\alpha_{\rm ov}\sim 0.45$ for the most massive
O-stars.

By contrast, when we measure the penetration zone in terms of $\Delta R_{\rm
core}/R_{\rm core}$ (Figure 2, right) the trend with mass almost vanishes.
This tells us that the trend in $\alpha_{\rm ov}$ is being driven almost
entirely by the changing volume of the convection zone. That is, the ratio
$R_{\rm core}/h$ is varying (Figure 3), and that in turn causes the extent of
the penetration zone to vary via the term $f\xi V_{\rm PZ}/V_{\rm CZ}$
appearing in equation (8).

The fact that the volume of the convection zone sets the trends in Figure 2 is
supported by Figure 4, which shows mid-main-sequence (central hydrogen
abundance $X_{c}=0.36$) profiles of $\nabla_{\rm rad}/\nabla_{\rm ad}$ as a
function of mass coordinate normalized to the nominal convective boundary.
These profiles are calculated assuming a fixed composition matching that of
the core, and so are the ones entering into equation (10). They are _not_ the
profiles actually used in the stellar evolution.

Except for the $1.5M_{\odot}$ model, the profiles are extremely similar,
meaning that the ratios of convective to total luminosity in equation (10) are
similar, and hence so are the ratio of luminosity integrals in equation (8).
This leaves only the volume of the convection zone to set the scale of the
penetration zone.

The $1.5M_{\odot}$ model, and those at lower masses too, show a somewhat
different profile, which is driven by a composition discontinuity at the
convective boundary. Even though we use a fixed composition in computing these
profiles, we use the actual density profiles from our stellar models, and
those can exhibit a discontinuity at the nominal convective boundary if there
is a composition discontinuity at the same. Despite this jump, we see just a
small difference in $\Delta R/R_{\rm core}$ in Figure 2, and so believe that
even at low masses the trend is dominated by the volume of the convection zone
rather than any composition effects.

Figure 3: The ratio $R/h$ at the outer boundary of the nominal convective core
is shown on a Hertzsprung-Russel diagram in terms of $\log T_{\rm eff}$ and
$\log L$ for stellar models ranging from $1.1-60M_{\odot}$ with Milky Way
metallicity $Z=0.014$. Note that even though this ratio becomes large, the
Boussinesq approximation remains valid over much of the domain, as evidenced
by the density ratio in Figure 6. Figure 4: Profiles of $\nabla_{\rm
rad}/\nabla_{\rm ad}$ are shown for selected mid-main-sequence stellar models
as a function of mass coordinate normalized to the nominal convective
boundary. These profiles are calculated assuming a fixed composition matching
that of the core, and so are the ones entering into equation (10). They are
_not_ the profiles actually used in the stellar evolution.

We likewise see clear trends in age for $\alpha_{\rm ov}$. At the high mass
end ($M\gtrsim 30M_{\odot}$) the penetration zone shrinks by 20-30% with age
from the Zero-Age Main Sequence (ZAMS) to the Terminal-Age Main Sequence
(TAMS), while at the low mass end it actually expands slightly with age
($M\lesssim 3M_{\odot}$). This may be seen more clearly in Figure 5, which
shows our predicted $\alpha_{\rm ov}$ as a function of mass at the ZAMS, the
TAMS, and the halfway point between these (mid-MS). As with mass, this trend
almost entirely disappears for $\Delta R_{\rm core}/R_{\rm core}$, again
because the volume of the convection zone dominates the trends of $\alpha_{\rm
ov}$.

The residual trends we see in $\Delta R_{\rm core}/R_{\rm core}$ are
principally at the low-mass end ($M\lesssim 2M_{\odot}$). There, the core mass
increases with age as the nuclear luminosity of the CNO cycle increases,
increasing $\nabla_{\rm rad}$ and changing its profile (Aerts et al., 2010,
see figure 3.6). Later, once the CNO cycle is dominant, $\nabla_{\rm rad}$
falls and the core mass declines222This decline is responsible for the jump in
$\nabla_{\rm rad}/\nabla_{\rm ad}$ in our $1.5M_{\odot}$ model in Figure 4..
These effects are also all functions of mass, and so give rise to both a mass
and age dependence in $\Delta R_{\rm core}/R_{\rm core}$.

Figure 5: The $\alpha_{\rm ov}$ step overshooting parameter needed to match
the predicted extent of the penetration zone is shown as a function of mass
for three different points in time, the ZAMS, the TAMS, and the midpoint
between them as measured by central ${}^{1}{\rm H}$ abundance. The dashed line
shows the fit to the mid-aged models given by equation (11).

Figure 5 also shows good agreement with the trends inferred by Claret & Torres
(2019) from eclipsing binary stars, at least over the range from
$1-3M_{\odot}$. Both our predictions and their inferences yield $\alpha_{\rm
ov}$ increasing rapidly with increasing mass, up to around $2M_{\odot}$, and
then a much weaker trend beyond that. We further predict that this trend of
increasing $\alpha_{\rm ov}$ with increasing mass continues beyond the edge of
their data ($\approx 3M_{\odot}$) up to quite large masses ($60M_{\odot}$). We
emphasize that the trends we see in $\alpha_{\rm ov}$ are dominated by the
volume of the nominal convection zone and not by the temperature gradient
profiles, as the latter are quite uniform in mass and age.

We additionally have reasonable quantitative agreement on the magnitude of
$\alpha_{\rm ov}$: they infer a plateau around $2M_{\odot}$ of $f_{\rm
ov}\approx 0.015$ and say that the equivalent $\alpha_{\rm ov}$ is $11.4$
times larger, giving $\alpha_{\rm ov}\approx 0.17$, which is a little higher
than our prediction of $\alpha_{\rm ov}\approx 0.12$ for the mid-MS but well
within the uncertainties of their inferences and of e.g. converting between
step-overshooting ($\alpha_{\rm ov}$) and exponential overshooting ($f_{\rm
ov}$).

By contrast Deheuvels et al. (2010) inferred $\alpha_{\rm ov}\approx 0.17$ for
the solar-like star HD 203608, which is quite a bit larger than our predicted
$\alpha_{\rm ov}\sim 0.05$ for low-mass ZAMS stars, though that object is also
at sub-solar metallicity, which could change the extent of mixing.

While we recommend using equation (8) to predict the extent of the penetration
zone in stellar models, there are instances where a fitting formula is more
convenient. We find a good fit to our mid-main-sequence results with

$\displaystyle\alpha_{\rm
ov}=\tanh(m-1.1)\sqrt{m-1.1}\frac{a+bm^{2}+cm^{3}+dm^{5}}{e+m^{5}}$ (11)

where $m\equiv M_{\star}/M_{\odot}$ and

$\displaystyle a$ $\displaystyle=2.47109,$ (12) $\displaystyle b$
$\displaystyle=-1.19087,$ (13) $\displaystyle c$ $\displaystyle=0.724183,$
(14) $\displaystyle d$ $\displaystyle=0.0470249,$ (15) and $\displaystyle e$
$\displaystyle=0.560757$ (16)

are fitting coefficients. Equation (11) produces a relative error of at most
21% and an absolute error of at most $0.02$. These errors are smaller than the
difference between $\alpha_{\rm ov}$ on the ZAMS and TAMS.

###  4.2 Mass Extent

We can also view these predictions in mass-coordinate. The left panel of
Figure 7 shows the mass of the penetration zone relative to the mass of the
nominal convective core. We see that the additional fractional mass which
penetration mixes into the core declines with increasing mass, though does so
gradually, adding at most $30\%$ to the mass of the core and so extending the
lifetime of the star by a similar amount.

Figure 6: The number of pressure scale heights in the nominal convective core
is shown on a Hertzsprung-Russel diagram in terms of $\log T_{\rm eff}$ and
$\log L$ for stellar models ranging from $1.1-60M_{\odot}$ with Milky Way
metallicity $Z=0.014$.

This trend could be due in part to our formulation (equation 8) and
calibrating simulations neglecting density stratification. As the mass of the
star increases the core grows and the number of scale heights it contains
increases (Figure 6). Thus a fixed relative increase in the _radial_ extent of
the core, as seen in Figure 2, translates to a smaller increase in the mass
extent because of the increasing density contrast between the penetration zone
and the center of the star. We intend to examine the effects of density
stratification in future work.

Figure 7: (Left) The mass of the penetration zone is shown on a Hertzsprung-
Russel diagram in terms of $\log T_{\rm eff}$ and $\log L$ for stellar models
ranging from $1.1-60M_{\odot}$ with Milky Way metallicity $Z=0.014$,
normalized to the mass coordinate at the nominal convective boundary reported
by MESA. Values below the lower edge of the scale occur at masses $M\lesssim
1.5M_{\odot}$ as the core convection zone is just forming. (Right) The
fractional mass of the star inside the outer edge of the penetration zone is
shown on the same diagram for the same models. Figure 8: The fractional mass
of the star inside the outer edge of the penetration zone is shown as a
function of mass for three different points in time, the ZAMS, the TAMS, and
the midpoint between them as measured by central ${}^{1}{\rm H}$ abundance.
Dashed lines indicate the nominal core mass with no penetration zone and solid
lines show the core mass including the penetration zone. Light horizontal
lines indicate the approximate measurements on the ZAMS and TAMS from Pedersen
et al. (2021).

Note that the fraction of a star’s total mass in the core including the
penetration zone shows somewhat different trends. The right panel of Figure 7
shows the fraction of the total mass which resides in the well-mixed core. The
relative core mass increases monotonically with total mass and decreases
monotonically with age. This may also be seen in Figure 8, which shows our
predicted relative core mass as a function of total mass at the ZAMS, the
TAMS, and the halfway point between these. The nominal core mass increases
rapidly with stellar mass and decreases rapidly with age, and these effects
suffice to overcome trends in the relative mass of the penetration zone.

With Figure 8 we can also compare our predictions to the trends inferred from
asteroseismology of B stars. Pedersen et al. (2021) find approximately $M_{\rm
core}/M_{\star}\approx 0.25$ on the ZAMS and $\approx 0.1$ on the TAMS for
stars ranging from $4-8M_{\odot}$ (see their Figure 6). This is in
quantitative agreement with our models which include convective penetration
(solid lines) and significantly greater than if we exclude the penetration
zone (dashed lines). Pedersen et al. (2021) also see a weak trend on the TAMS
towards increasing fractional core mass with increasing stellar mass, of order
an increase of $\Delta(M_{\rm core}/M_{\star})\approx 0.1$ from
$4-8M_{\odot}$, which agrees roughly with our predictions.

##  5 Discussion

There are several uncertainties which require further study to provide a
complete theory of convective penetration. For instance Anders et al. (2021)
only examined cartesian domains in the Boussinesq limit. For stellar cores a
spherical geometry is more appropriate, and may change the obtained values of
$f$ and $\xi$ in equation (8), though we do not expect this to change the
orders of magnitude or the trends we see.

More importantly, at higher masses we see increasing density stratification in
stellar cores, and convective envelopes are even more strongly density-
stratified, so extending this work to higher mass main-sequence stars or Red
Giants dredge up requires moving beyond the Boussinesq approximation. Roxburgh
(1989) suggests a way to adapt this theory for stratified systems. This
approach was applied to stellar models by neglecting dissipation, producing an
upper bound on the amount of penetration (Roxburgh, 1992; Dowler & Vandenberg,
1993; VandenBerg & Stetson, 2004; VandenBerg et al., 2006), but so far as we
are aware this has not yet been done including dissipation and calibrating
with simulations. We know, however, that the existing formalism of Anders et
al. (2021) is insufficient because equation (8), applied to the Red Giant
Branch, predicts at points that the whole star is well-mixed, which is clearly
ruled out by observations.

Additionally the effects of rotation and magnetism on convective penetration
remain poorly understood, though there has been recent progress on the former
(Augustson & Mathis, 2018). Incorporating these effects may prove important to
explain the diversity of mixing profiles inferred from observations (Pedersen
et al., 2021).

##  6 Conclusions

We have applied a recent theory of convective penetration by Anders et al.
(2021) to predict the extent of convective core boundary mixing in early-type
stars owing to this mechanism. We find good agreement both in magnitude and in
trends with inferred core masses from asteroseismology and eclipsing binary
studies. This suggests that convective penetration may be responsible for most
of the convective boundary mixing which occurs in these stars, and motivates
the inclusion of this process in stellar evolution software instruments.
Importantly, while existing instruments support including the effects of
convective penetration on compositional mixing via a step overshoot
prescription, they do not generally feature a corresponding adiabatic layer,
which matters for making precise predictions of asteroseismic signatures
(Pedersen et al., 2021).

The Flatiron Institute is supported by the Simons Foundation. This research
was supported in part by the National Science Foundation under Grant No.
PHY-1748958. We are grateful to May Gade Pedersen and Ben Brown for productive
conversations on these topics.

##  Appendix A MESA

The MESA EOS is a blend of the OPAL (Rogers & Nayfonov, 2002), SCVH (Saumon et
al., 1995), FreeEOS (Irwin, 2004), HELM (Timmes & Swesty, 2000), and PC
(Potekhin & Chabrier, 2010) EOSes.

Radiative opacities are primarily from OPAL (Iglesias & Rogers, 1993, 1996),
with low-temperature data from Ferguson et al. (2005) and the high-
temperature, Compton-scattering dominated regime by Buchler & Yueh (1976).
Electron conduction opacities are from Cassisi et al. (2007).

Nuclear reaction rates are from JINA REACLIB (Cyburt et al., 2010) plus
additional tabulated weak reaction rates Fuller et al. (1985); Oda et al.
(1994); Langanke & Martínez-Pinedo (2000). Screening is included via the
prescription of Chugunov et al. (2007). Thermal neutrino loss rates are from
Itoh et al. (1996).

Models were constructed on the pre-main sequence with $Y=0.024+2Z$ and
$X=1-Y-Z$ and evolved from there.

Convective penetration is calculated by re-casting equation (8) in the form

$\displaystyle\int_{\rm{PZ}}-L_{\rm{conv}}+f\xi 4\pi r^{2}F_{\rm
avg}\,dr=(1-f)V_{\rm CZ}F_{\rm avg},$ (A1)

where

$\displaystyle F_{\rm avg}\equiv\frac{1}{V_{\rm CZ}}\int_{\rm CZ}L_{\rm
conv}dr.$ (A2)

and where $L_{\rm conv}$ is evaluated in the penetration zone using equation
(10). We then integrate the left-hand side of equation (A1) away from the
nominal convective boundary returned by MESA until we obtain the right-hand
side.

##  Appendix B Data Availability

The plotting scripts used in this work may be found in this GitHub repository
in the commit with short-sha 2fc2b48. The inlists and run scripts used in
producing the HR diagrams in this work are available in this other GitHub
repository on the main branch in the commit with short-sha 964e07d. The data
those scripts produced are available in Jermyn et al. (2022).

## References

  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology
  * Anders et al. (2021) Anders, E. H., Jermyn, A. S., Lecoanet, D., & Brown, B. P. 2021, arXiv e-prints, arXiv:2110.11356.  https://arxiv.org/abs/2110.11356
  * Augustson & Mathis (2018) Augustson, K. C., & Mathis, S. 2018, in SF2A-2018: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. P. Di Matteo, F. Billebaud, F. Herpin, N. Lagarde, J. B. Marquette, A. Robin, & O. Venot, Di.  https://arxiv.org/abs/1811.04413
  * Balachandran (1995) Balachandran, S. 1995, ApJ, 446, 203, doi: 10.1086/175779
  * Basu (2016) Basu, S. 2016, Living Reviews in Solar Physics, 13, 2, doi: 10.1007/s41116-016-0003-4
  * Brott et al. (2011) Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115, doi: 10.1051/0004-6361/201016113
  * Brummell et al. (2002) Brummell, N. H., Clune, T. L., & Toomre, J. 2002, ApJ, 570, 825, doi: 10.1086/339626
  * Buchler & Yueh (1976) Buchler, J. R., & Yueh, W. R. 1976, ApJ, 210, 440, doi: 10.1086/154847
  * Cassisi et al. (2007) Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094, doi: 10.1086/516819
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028, doi: 10.1103/PhysRevD.76.025028
  * Claret & Torres (2019) Claret, A., & Torres, G. 2019, ApJ, 876, 134, doi: 10.3847/1538-4357/ab1589
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240, doi: 10.1088/0067-0049/189/1/240
  * Deheuvels et al. (2010) Deheuvels, S., Bruntt, H., Michel, E., et al. 2010, A&A, 515, A87, doi: 10.1051/0004-6361/200913490
  * Dowler & Vandenberg (1993) Dowler, P. D., & Vandenberg, D. A. 1993, JRASC, 87, 185
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53, doi: 10.3847/1538-4357/ab518b
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585, doi: 10.1086/428642
  * Fuller et al. (1985) Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1, doi: 10.1086/163208
  * Garcia Lopez & Spruit (1991) Garcia Lopez, R. J., & Spruit, H. C. 1991, ApJ, 377, 268, doi: 10.1086/170356
  * Herwig (2000) Herwig, F. 2000, A&A, 360, 952.  https://arxiv.org/abs/astro-ph/0007139
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Iglesias & Rogers (1993) Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752, doi: 10.1086/172958
  * Iglesias & Rogers (1996) —. 1996, ApJ, 464, 943, doi: 10.1086/177381
  * Irwin (2004) Irwin, A. W. 2004, The FreeEOS Code for Calculating the Equation of State for Stellar Interiors.  http://freeeos.sourceforge.net/
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411, doi: 10.1086/192264
  * Jermyn et al. (2022) Jermyn, A., Anders, E., & Cantiello, M. 2022, Supporting information for v2 of “A Transparent Window into Early-Type Stellar Variability”, v2, Zenodo, doi: 10.5281/zenodo.5878965
  * Khan et al. (2018) Khan, S., Hall, O. J., Miglio, A., et al. 2018, ApJ, 859, 156, doi: 10.3847/1538-4357/aabf90
  * Korre et al. (2019) Korre, L., Garaud, P., & Brummell, N. H. 2019, Monthly Notices of the Royal Astronomical Society, 484, 1220, doi: 10.1093/mnras/stz047
  * Langanke & Martínez-Pinedo (2000) Langanke, K., & Martínez-Pinedo, G. 2000, Nuclear Physics A, 673, 481, doi: 10.1016/S0375-9474(00)00131-7
  * Oda et al. (1994) Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data and Nuclear Data Tables, 56, 231, doi: 10.1006/adnd.1994.1007
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pedersen et al. (2021) Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nature Astronomy, 5, 715, doi: 10.1038/s41550-021-01351-x
  * Potekhin & Chabrier (2010) Potekhin, A. Y., & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82, doi: 10.1002/ctpp.201010017
  * Rogers & Nayfonov (2002) Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064, doi: 10.1086/341894
  * Rogers & McElwaine (2017) Rogers, T. M., & McElwaine, J. N. 2017, ApJ, 848, L1, doi: 10.3847/2041-8213/aa8d13
  * Roxburgh (1989) Roxburgh, I. W. 1989, A&A, 211, 361
  * Roxburgh (1992) —. 1992, A&A, 266, 291
  * Salaris & Cassisi (2017) Salaris, M., & Cassisi, S. 2017, Royal Society Open Science, 4, 170192, doi: 10.1098/rsos.170192
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713, doi: 10.1086/192204
  * Serrano (1986) Serrano, A. 1986, PASP, 98, 1066, doi: 10.1086/131874
  * Shaviv & Salpeter (1973) Shaviv, G., & Salpeter, E. E. 1973, ApJ, 184, 191, doi: 10.1086/152318
  * Timmes & Swesty (2000) Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501, doi: 10.1086/313304
  * Townsend (2019a) Townsend, R. H. D. 2019a, MESA SDK for Mac OS, 20190503, Zenodo, doi: 10.5281/zenodo.2669543
  * Townsend (2019b) —. 2019b, MESA SDK for Linux, 20190503, Zenodo, doi: 10.5281/zenodo.2669541
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * VandenBerg et al. (2006) VandenBerg, D. A., Bergbusch, P. A., & Dowler, P. D. 2006, The Astrophysical Journal Supplement Series, 162, 375, doi: 10.1086/498451
  * VandenBerg & Stetson (2004) VandenBerg, D. A., & Stetson, P. B. 2004, PASP, 116, 997, doi: 10.1086/426340
  * Zahn (1991) Zahn, J. P. 1991, A&A, 252, 179