††thanks: Shared first authorship: These authors contributed equally to this
work.††thanks: Shared first authorship: These authors contributed equally to
this work.

# Effects of pair freeze-out on photon distributions in BBN epoch

Jeongyoon Choi  Department of Physics, Ulsan National Institute of Science and
Technology (UNIST), Ulsan 44919, Republic of Korea  Dukjae Jang  School of
Physics, and International Research Center for Big-Bang Cosmology and Element
Genesis, Beihang University, Beijing 100083, P. R. China  Youngshin Kwon
Corresponding author<EMAIL_ADDRESS>Department of Science
Education, Ewha Womans University, Seoul 03760, Republic of Korea  Gwangeon
Seong  Department of Physics, Ulsan National Institute of Science and
Technology (UNIST), Ulsan 44919, Republic of Korea  Myeong Hwan Mun
Department of Physics and OMEG Institute, Soongsil University, Seoul 06978,
Republic of Korea  Young-Min Kim  Korea Astronomy and Space Science Institute,
Daejeon 34055, Republic of Korea  Kyujin Kwak  Department of Physics, Ulsan
National Institute of Science and Technology (UNIST), Ulsan 44919, Republic of
Korea  Myung-Ki Cheoun  Department of Physics and OMEG Institute, Soongsil
University, Seoul 06978, Republic of Korea

###### Abstract

We investigate the evolution of non-extensivity in the photon distribution
during the Big Bang Nucleosynthesis (BBN) epoch using Tsallis statistics.
Assuming a minimal deviation from the Planck distribution, we construct the
perturbed Boltzmann equation for photons, including the collision terms for
pair creation and annihilation processes. We analyze the possibility that
these collisions could cause a slight increase in the number of high-frequency
photons within the BBN era, and consequently, the primordial plasma might be
temporarily placed in a state of chemical non-equilibrium. We also discuss the
restoration of the photon distribution to an equilibrium state as the Universe
enters the matter-dominated era. These findings, which suggest possible
changes in the photon distribution during the epoch between the BBN and the
recombination, offer insights that support the previously proposed ansatz
solution to the primordial lithium problem in Ref. [6].

Big bang nucleosynthesis, primordial lithium problem, non-extensive statistics

###### pacs:

##  I Introduction

Standard Big Bang Nucleosynthesis (SBBN) refers to the prevailing theoretical
framework that describes the production of light elements in the first few
minutes after the Big Bang. It provides a detailed account of the formation of
light elements such as deuterium ($\mathrm{D}$), helium isotopes
(${}^{3}\mathrm{He}$ and ${}^{4}\mathrm{He}$), and lithium
(${}^{7}\mathrm{Li}$). Remarkably, the primordial abundances of deuterium and
helium isotopes are well explained by SBBN. The success of SBBN in predicting
the observed abundances of these light elements, along with the cosmic
microwave background (CMB) observations, lends support to the hot Big Bang
theory and sheds light on the conditions of the early Universe. However, the
so-called primordial lithium problem, which indicates a discrepancy between
the theoretical prediction and the observational data of primordial
${}^{7}\mathrm{Li}$ abundance, suggests that the SBBN model may need to be
revised or updated.

To address this conundrum, various avenues have been explored in prior
research. These include the investigation of the potential role of long-lived
exotic particles as a solution [1] and the examination of possible
observational biases that could affect the determination of the primordial
lithium abundance [2]. Moreover, the implications of adopting Tsallis
statistics have been examined, leading to discussions about changes in the
nuclei distribution function [3, 4, 5]. By introducing the entropic index $q$,
the sole additional parameter for Tsallis statistics, it becomes feasible to
adjust the high-energy tail of the nuclei distribution functions and the
related reaction rates. However, it has been found that the application of
Tsallis statistics to nuclei distribution does not fully resolve the
primordial lithium problem [5].

Recently, some of the current authors [6] proposed that taking nonideal
properties of the plasma into account could be a clue to elucidate this
problem of BBN. SBBN assumes that BBN occurs within a thermal plasma, and
describes how the constituents of primordial plasma are distributed in terms
of their velocity or energy by using Boltzmann-Gibbs (BG) statistics. Because
this assumption was supported by rapid interactions between the plasma
constituents, which have characteristic time-scales shorter than the cosmic
expansion, the primordial plasma in an equilibrium state has won widespread
acceptance. However, we took note of the electron chemical potential
increasing dramatically at $T\lesssim 10^{9}\,\mathrm{K}$, which might be able
to lead the primordial plasma to a weakly nonideal state. Then, we adopted the
Tsallis statistics for describing photon distribution to reflect the nonideal
effect, which appears more appropriate than the standard BG statistics. (For
more details on the model construction, refer to Ref. [6].) We found that the
abundances of the light elements up to ${}^{7}\mathrm{Li}$ calculated with the
Tsallis distribution of photons can explain the observational data within
errors when a proper value of $q$ is chosen. Although this solution is
effective for the primordial lithium problem, it raises a crucial question:
whether or how using Tsallis statistics is physically valid for the photon
distribution in the BBN era. The present work aims at finding the answer to
that question by constructing and solving the Boltzmann equation for photons
with relevant collision terms.

When the temperature of the Universe reached about $10^{9}\,\mathrm{K}$, the
energy of most photons would no longer be sufficient for electron-positron
pair production ($2\gamma\to e^{+}e^{-}$). It disrupted the detailed balance
between the pair production and annihilation. While the pair annihilation
($e^{+}e^{-}\to 2\gamma$) continued until consuming all positrons and leaving
slightly excess electrons, the pair production could not replenish them.
Accordingly, the chemical equilibrium between the pairs and photons is not
established, at least for a while, and we hypothesize that therein lies the
last piece of the puzzle for primordial abundances. In this study, we examine
if the broken detailed balance between the pair production and annihilation
processes can induce a deviation from the Planck distribution for photons by
adding those processes into the collision term of the Boltzmann equation.

To achieve a comprehensive understanding of the modifications in photon
distribution during the BBN epoch, it is essential to solve the Boltzmann
equation by taking into account all pertinent constituents, namely photons,
electrons, positrons, and nuclei. This task would be intrinsically complex
because of the intricate nature of the system and the associated coupled set
of partial differential equations. To mitigate this complexity, we make some
approximation. We assume that only photons deviate from equilibrium to a small
extent, i.e., all particles except photons are in equilibrium and a change in
$q$, which describes the Tsallis statistics of photons, is small. With this
assumption, we expand the generalized Planck distribution utilizing Tsallis
statistics to first order, subsequently deriving the first-order Boltzmann
equation. Our analysis based on the perturbed Boltzmann equation suggests that
the pair annihilation process during the BBN epoch can potentially augment
$q$. Furthermore, we provide a discussion for the distorted photon
distribution restoring to a Planck distribution observed during the matter-
dominated era.

The remainder of this paper is organized as follows. In Sect. II, we present
the formalism utilized to describe the evolution of the parameter $q$ through
the perturbed Boltzmann equation. In Sect. III, we focus on the evaluation of
the collision term, particularly emphasizing key collision processes such as
pair creation and annihilation. In Sect. IV, we present our results regarding
the impact of Tsallis statistics on the photon distribution during the BBN
epoch. In Sect. V, we discuss the potential of restoring the distorted photon
spectrum to a Planck distribution to achieve consistency with the observed
CMB. Finally, Sect. VI provides a summary and conclusion of this paper. Note
that the natural unit of $\hbar=c=1$ is used throughout this paper.

##  II Boltzmann equation

The evolution of a particle distribution is generally governed by the
Boltzmann equation, represented as $\mathbb{L}[f]=\mathbb{C}[f]$. The
Liouville operator, $\mathbb{L}$, is associated with the streaming part of the
system in the absence of collisions. It captures the system’s behavior as it
approaches equilibrium or deviates from it by directing the time evolution of
the distribution function. On the other hand, the collision operator,
$\mathbb{C}$, accounts for interactions between particles in the system. The
streaming part for the evolution of photon distribution, $f_{\gamma}$, is
detailed in the relativistic Boltzmann equation as follows:

$E_{\gamma}\frac{\partial f_{\gamma}}{\partial
t}-H|\mathbf{p}_{\gamma}|^{2}\frac{\partial f_{\gamma}}{\partial
E_{\gamma}}=\mathbb{C}[f_{\gamma}(\mathbf{p_{\gamma}})]~{},$ (1)

where $H$ denotes the cosmic expansion rate.

###  II.1 Perturbed photon distribution

The model involved in this study accounts for the possibility of the deviation
of photon distribution, which can then be written and generalized according to
the Tsallis statistics [8, 7]:

$f_{\gamma}(E_{\gamma},T_{q};q)=\left[\left(1-(1-q)\frac{E_{\gamma}}{T_{q}}\right)^{\frac{1}{q-1}}-1\right]^{-1}~{},$
(2)

as a function of energy ($E^{2}_{\gamma}=\mathbf{p}^{2}_{\gamma}$ for photons)
and temperature. Here, we introduce a $q$-dependent effective temperature,
$T_{q}$, to distinguish it from the equilibrium temperature, $T$. The entropic
index $q$ quantifies the deviation from equilibrium. In the limit of $q\to 1$,
Eq. (2) converges to the Planck distribution with $T$. It is worth noting that
the solution obtained in Ref. [6] for the primordial lithium problem was
$q\approx 1.027$. Considering such a slight deviation of $q$ from unity, it
stands to reason that the photons are assumed to be nearly in a global
equilibrium. Consequently, Eq. (2) can be expanded in terms of perturbations
around the equilibrium distribution:

$f_{\gamma}(E_{\gamma},T_{q};Q)\simeq
f_{\gamma}^{(0)}(E_{\gamma},T)+f_{\gamma}^{(1)}(E_{\gamma},T_{q};Q)~{}.$ (3)

Here, the entropic index, $Q\equiv 1-q$, is redefined for ease of use as a
perturbative parameter. The first term on the right-hand side of Eq. (3)
stands for nothing but the Planckian form:

$f_{\gamma}^{(0)}(E_{\gamma},T)=f_{\gamma}(E_{\gamma},T_{q};Q\to
0)=\frac{1}{e^{E_{\gamma}/T}-1}~{},$ (4)

and the second term is specified by the perturbation truncated at first order:

$\displaystyle f_{\gamma}^{(1)}(E_{\gamma},T_{q};Q)=Q\left[\frac{\partial
f_{\gamma}(E_{\gamma},T_{q};Q)}{\partial Q}\right]_{Q\to 0}$
$\displaystyle=-\frac{Q}{2T^{2}}\,\frac{E_{\gamma}\,e^{E_{\gamma}/T}}{(e^{E_{\gamma}/T}-1)^{2}}\left(E_{\gamma}-2\left[\frac{dT_{q}}{dQ}\right]_{Q\to
0}\right)$
$\displaystyle=-\frac{Q}{2T^{2}}\,\frac{E_{\gamma}\,e^{E_{\gamma}/T}}{(e^{E_{\gamma}/T}-1)^{2}}\left(E_{\gamma}-\frac{450\zeta(5)T}{\pi^{4}}\right)~{},$
(5)

where $\zeta(x)$ is the Riemann zeta function.

###  II.2 Photon number density

By inserting Eq. (3) into Eq. (1) and performing an integration over the
momentum space, we can transform the Boltzmann equation into an equation
governing the evolution of the number density:

$\frac{d\left(n_{\gamma}^{(0)}+n_{\gamma}^{(1)}\right)}{dt}+3H\left(n_{\gamma}^{(0)}+n_{\gamma}^{(1)}\right)=\frac{g_{\gamma}}{(2\pi)^{3}}\int\mathbb{C}[f^{(0)}_{\gamma}+f^{(1)}_{\gamma}]\,\frac{d\mathbf{p}_{\gamma}}{E_{\gamma}}~{},$
(6)

where $g_{\gamma}$ represents the statistical degrees of freedom for photons
and is given as $g_{\gamma}=2$. The zeroth- and first-order number densities
are written, respectively, as follows:

$\displaystyle n_{\gamma}^{(0)}$ $\displaystyle=\frac{2}{(2\pi)^{3}}\int
d\mathbf{p}_{\gamma}\,f_{\gamma}^{(0)}(E_{\gamma},T;Q)=\frac{2\zeta(3)T^{3}}{\pi^{2}}~{},$
(7) $\displaystyle n_{\gamma}^{(1)}$ $\displaystyle=\frac{2}{(2\pi)^{3}}\int
d\mathbf{p}_{\gamma}\,f_{\gamma}^{(1)}(E_{\gamma},T_{q};Q)=-\mathcal{A}T^{3}Q~{},$
(8)

where

$\mathcal{A}=\frac{2\pi^{2}}{15}-\frac{1350\zeta(5)\zeta(3)}{\pi^{6}}~{}.$ (9)

###  II.3 Linear collision term

The general form of the collision integral $\mathbb{C}[f_{\gamma}]$ for
photon-induced processes is given in Ref. [9]:

$\frac{g_{\gamma}}{(2\pi)^{3}}\int\mathbb{C}[f_{\gamma}]\,\frac{d\mathbf{p}_{\gamma}}{E_{\gamma}}=-\int\left(\prod_{i=\gamma,1,2,3}d\Pi_{i}\right)\left|\mathcal{M}\right|^{2}(2\pi)^{4}\,\delta^{4}(p_{\gamma}+p_{1}-p_{2}-p_{3})\,F(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})$
(10)

with

$d\Pi_{i}=\frac{g_{i}}{(2\pi)^{3}}\frac{d\mathbf{p}_{i}}{2E_{i}}~{}.$ (11)

Here $g_{i}$ is the statistical degrees of freedom for species $i$.
$F(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})$ is
defined as:

$F(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})\equiv
f_{2}(\mathbf{p}_{2})f_{3}(\mathbf{p}_{3})\Big{(}1+f_{\gamma}(\mathbf{p}_{\gamma})\Big{)}\Big{(}1\pm
f_{1}(\mathbf{p}_{1})\Big{)}-f_{\gamma}(\mathbf{p}_{\gamma})f_{1}(\mathbf{p}_{1})\Big{(}1\pm
f_{2}(\mathbf{p}_{2})\Big{)}\Big{(}1\pm f_{3}(\mathbf{p}_{3})\Big{)}~{},$ (12)

where the positive and negative signs correspond to bosons and fermions,
respectively.

By substituting Eq. (3) into Eq. (12), one can find that the collision term is
linear with respect to each order, that is,
$\mathbb{C}[f_{\gamma}^{(0)}+f_{\gamma}^{(1)}]=\mathbb{C}[f_{\gamma}^{(0)}]+\mathbb{C}[f_{\gamma}^{(1)}]$.
The linearity of the collision term significantly simplifies the problem,
enabling us to consider the effects of each order of the distribution function
independently. Therefore, under the assumption that the zeroth order remains
in equilibrium, such that
$\mathbb{L}[f_{\gamma}^{(0)}]=\mathbb{C}[f_{\gamma}^{(0)}]$, Eq. (6) can be
rewritten as follows:

$\frac{dn_{\gamma}^{(1)}}{dt}+3Hn_{\gamma}^{(1)}=\frac{g_{\gamma}}{(2\pi)^{3}}\int\mathbb{C}[f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})]\,\frac{d\mathbf{p}_{\gamma}}{E_{\gamma}}~{}.$
(13)

###  II.4 Liouville term

The evolution of $n_{\gamma}^{(1)}$ in Eq. (13) can be converted into a
differential equation for $Q$. This subsection aims to reformulate the
Liouville term in Eq. (13) into an expression involving $Q$. To achieve this,
it is essential to establish the time-temperature relation, which can be
derived from the Friedmann equation. Owing to energy conservation governing
the photon energy density, $H$ is expressed in its standard form:

$H^{2}=\left(\frac{\dot{a}}{a}\right)^{2}=\frac{8\pi
G}{3}\rho_{\text{total}}~{},$ (14)

where $a$ represents the scale factor, and $\dot{a}$ is its time derivative.
By using the relation between the energy density and temperature in the
radiation-dominated era, it can be rewritten by

$H=\frac{T^{2}}{m_{\text{pl}}}\sqrt{\frac{8\pi^{3}g_{\ast}}{90}}\equiv\mathcal{B}T^{2}~{},$
(15)

where $g_{\ast}$ represents the total degrees of freedom for the effectively
relativistic particles and $m_{\text{pl}}$ stands for the Planck mass. Since
the scale factor follows $a\propto t^{1/2}$ in the radiation-dominated era,
the time-temperature relation is simply as follows:

$t=\frac{1}{2H}=\frac{1}{2\mathcal{B}T^{2}}~{},$ (16)

which gives rise to the time-derivative of $T$:

$\frac{dT}{dt}=-\mathcal{B}T^{3}~{}.$ (17)

By using the time-temperature relations of Eqs. (16) and (17), the Liouville
term in Eq. (13) takes the form of

$\displaystyle\frac{dn^{(1)}_{\gamma}}{dt}+3Hn_{\gamma}^{(1)}=\frac{\partial
n_{\gamma}^{(1)}}{\partial T}\frac{dT}{dt}+3Hn_{\gamma}^{(1)}$
$\displaystyle=\left(-3\mathcal{A}T^{2}Q-\mathcal{A}T^{3}\frac{dQ}{dT}\right)\left(-\mathcal{B}T^{3}\right)+3\left(\mathcal{B}T^{2}\right)\left(-\mathcal{A}T^{3}Q\right)$
$\displaystyle=\mathcal{A}\mathcal{B}T^{6}\frac{dQ}{dT}~{}.$ (18)

##  III Collision term in focus

In this section, we turn our attention to the evaluation of the collision
term, focusing primarily on key collision processes such as Compton scattering
and pair creation and annihilation. However, it should be noted that our
methodology effectively omits the Compton process within the non-relativistic
limit ($T<0.5\,\mathrm{MeV}$). This exclusion is well-founded, as the momentum
change in electrons due to interaction with photons is negligible in this low-
energy regime (see Appendix for further details). Additionally, higher-order
processes, such as double Compton scattering and bremsstrahlung, can be safely
neglected, given the significant reduction in electron number density at lower
temperatures. As a result, our analysis is specifically tailored to examine
pair creation and annihilation as the dominant collisional mechanism.

To apply Eq. (10) for pair creation and annihilation, that is,
$\gamma(E_{\gamma},\mathbf{p}_{\gamma})+\gamma(E_{1},\mathbf{p}_{1})\rightleftharpoons
e^{-}(E_{2},\mathbf{p}_{2})+e^{+}(E_{3},\mathbf{p}_{3})$, the transition
amplitude is calculated as follows:

$|\mathcal{M}|^{2}_{\text{pair}}=2e^{4}\left[\frac{p_{\gamma}\cdot
p_{3}}{p_{\gamma}\cdot p_{2}}+\frac{p_{\gamma}\cdot p_{2}}{p_{\gamma}\cdot
p_{3}}+2m_{e}^{2}\left(\frac{1}{p_{\gamma}\cdot
p_{2}}+\frac{1}{p_{\gamma}\cdot
p_{3}}\right)-m_{e}^{4}\left(\frac{1}{p_{\gamma}\cdot
p_{2}}+\frac{1}{p_{\gamma}\cdot p_{3}}\right)^{2}\right]~{},$ (19)

and

$F(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})=f_{e^{-}}(\mathbf{p}_{2})f_{e^{+}}(\mathbf{p}_{3})\Big{(}1+f_{\gamma}(\mathbf{p}_{\gamma})\Big{)}\Big{(}1+f_{\gamma}(\mathbf{p}_{1})\Big{)}-f_{\gamma}(\mathbf{p}_{\gamma})f_{\gamma}(\mathbf{p}_{1})\Big{(}1-f_{e^{-}}(\mathbf{p}_{2})\Big{)}\Big{(}1-f_{e^{+}}(\mathbf{p}_{3})\Big{)}~{}.$
(20)

By setting electrons and positrons to remain in equilibrium – a reasonable
assumption given the small momentum transfer, as supported by Ref. [10] – we
can rewrite Eq. (20) based on Eq. (3) as follows:

$F(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})=F^{(0)}(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})+QF^{(1)}(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})+\mathcal{O}(Q^{2})~{},$
(21)

where each term is expressed as:

$\displaystyle
F^{(0)}(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})$
$\displaystyle=f_{e^{-}}(\mathbf{p}_{2})f_{e^{+}}(\mathbf{p}_{3})\Big{(}1+f_{\gamma}^{(0)}(\mathbf{p}_{\gamma})\Big{)}\Big{(}1+f_{\gamma}^{(0)}(\mathbf{p}_{1})\Big{)}$
$\displaystyle\quad-
f_{\gamma}^{(0)}(\mathbf{p}_{\gamma})f_{\gamma}^{(0)}(\mathbf{p}_{1})\Big{(}1-f_{e^{-}}(\mathbf{p}_{2})\Big{)}\Big{(}1-f_{e^{+}}(\mathbf{p}_{3})\Big{)}$
(22)

and

$\displaystyle
F^{(1)}(\mathbf{p}_{\gamma},\mathbf{p}_{1};\mathbf{p}_{2},\mathbf{p}_{3})$
$\displaystyle=f_{e^{-}}(\mathbf{p}_{2})f_{e^{+}}(\mathbf{p}_{3})\left[\left(1+f_{\gamma}^{(0)}(\mathbf{p}_{\gamma})\right)f_{\gamma}^{(1)}(\mathbf{p}_{1})+f_{\gamma}^{(1)}(\mathbf{p}_{\gamma})\left(1+f_{\gamma}^{(0)}(\mathbf{p}_{1})\right)\right]$
$\displaystyle\quad-\left[f^{(0)}_{\gamma}(\mathbf{p})f^{(1)}_{\gamma}(\mathbf{p}_{1})+f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})f^{(0)}_{\gamma}(\mathbf{p}_{1})\right]\left(1-f_{e^{-}}(\mathbf{p}_{2})\right)\left(1-f_{e^{+}}(\mathbf{p}_{3})\right)~{}.$
(23)

Given that $Q$ does not depend on momentum integration, we can formulate the
collision integral for the first-order perturbation as follows:

$\frac{g_{\gamma}}{(2\pi)^{3}}\int\mathbb{C}[f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})]\frac{d\mathbf{p}_{\gamma}}{p_{\gamma}}=QM(T)~{}.$
(24)

Here, $M(T)$ is defined by the following expression:

$M(T)\equiv-\int\left(\prod_{i=\gamma,1,2,3}d\Pi_{i}\right)|\mathcal{M}|^{2}_{\text{pair}}F^{(1)}(\mathbf{p}_{\gamma},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})\,(2\pi)^{4}\delta^{3}(\mathbf{p}_{\gamma}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})\,\delta(E_{\gamma}+E_{1}-E_{2}-E_{3})~{}.$
(25)

The evaluation of this collision integral employs Gauss-Legendre and Gauss-
Laguerre quadrature methods, following the approach described in Ref. [11].
The resulting values of $M(T)$ are plotted in Fig. 1.

By applying Eqs. (18) and (24), we rewrite the perturbed Boltzmann equation
given in Eq. (13) as shown in below:

$\frac{dQ}{dT}=\frac{M(T)}{\mathcal{A}\mathcal{B}T^{6}}Q~{}.$ (26)

Figure 1:  Logarithmic plot of the collision term $M(T)$ as a function of
temperature, where $M(T)$ is evaluated in $\mathrm{MeV}^{4}$.

##  IV Results

The solution of Eq. (26) is given as

$Q(T)=Q_{0}\exp\left[\int_{T_{0}}^{T}\frac{M(T^{\prime})}{\mathcal{A}\mathcal{B}{T^{\prime}}^{6}}dT^{\prime}\right]~{},$
(27)

where a constant $Q_{0}$ can be positive or negative and will be determined by
an initial condition $Q(T_{0})$. Setting $Q_{0}=0$ at a specific $T_{0}$ leads
to a trivial solution, $Q=0$, as driving the system towards perfect
equilibrium by rapid collisions in the BBN epoch. This trivial solution is
equivalent to the assumption of SBBN. The cosmological principle, which posits
that the Universe is homogeneous and isotropic on large scales, underlies this
assumption.

Figure 2: Solution of $q(T)=1-Q(T)$ as a function of temperature with the
initial condition of a minuscule deviation $Q_{0}=-10^{-37}$ at
$T_{0}=0.0235\,\mathrm{MeV}$.

However, a perfect thermal equilibrium in the early Universe is unlikely. Even
negligibly small, non-zero fluctuations could emerge because of various
instabilities and other unknown factors. Minor seed fluctuations in
temperature could significantly impact the anisotropic evolution of the
Universe (akin to phenomena observed in the CMB). In this context, we
investigate the implications of a non-zero initial condition $Q_{0}$ in
solving Eq. (26). Fig. 2 presents, as an example, the numerical solution
obtained with an initial condition of $Q_{0}=-10^{-37}$ at
$T_{0}=0.0235\,\mathrm{MeV}$. This selected initial condition leads to
behavior that the original entropic index, $q$, exceeds unity and eventually
stabilizes at the solution of $q\simeq 1.027$ reported in Ref. [6], resulting
in a slight enhancement of the high-energy tail of the photon distribution. It
is not the unique condition under which $q$ is saturated at such a solution.
The initial condition could, in principle, be further optimized to a more
minuscule $Q_{0}$ at a higher $T_{0}$. As indicated in Fig. 2, the initial
small deviation in $q(T)$ evolves rapidly as the temperature cools. However,
this increasing trend in $q(T)$ does not persist because the collision term
$M(T)$ also swiftly decreases in magnitude, as shown in Fig. 1. Consequently,
the Boltzmann equation transitions to a collisionless regime at low
temperatures, where $M(T)\rightarrow 0$. In this regime, Eq. (26) simplifies
to $dQ/dT\rightarrow 0$, leading to a constant $Q(T)$. Hence, in Fig. 2,
$Q(T)$ becomes constant at $T<0.02\,\mathrm{MeV}$.

The freeze-out of electron-positron pairs is attributed to the disruption of
the detailed balance between pair creation and annihilation. This process
already begins during the early stages of BBN when the thermal energy per
particle becomes insufficient to create such pairs. The extremely small
deviation in $Q_{0}$ might represent a vestige of these early-stage BBN
conditions. Given that the magnitude of deviation $|Q|$ is on the order of the
temperature fluctuation $\delta T/T$, setting $|Q_{0}|=10^{-37}$ is
effectively negligible, especially when compared to the CMB fluctuation
$\delta T/T\sim 10^{-5}$. Nonetheless, the pair freeze-out process could cause
this small fluctuation to evolve and become pronounced enough to contribute to
BBN abundances, particularly affecting the primordial lithium abundance. This
would be one possible interpretation of our calculation in this work.

##  V Discussion

###  V.1 Reequilibration

In the previous section, we demonstrated that a seed fluctuation in the photon
distribution could grow up through the freeze-out of electron-positron pairs
during the BBN epoch. However, CMB observations suggest that the photon
distribution closely approximates the Planck distribution [12]. Therefore, it
is crucial for the altered $Q$ to revert to the Planck distribution prior to
the recombination era. While a detailed exploration of this recovery mechanism
would necessitate the consideration of specific micro-processes affecting
photons in the evolving Universe, we provide a rudimentary estimation based on
the Boltzmann equation for $Q$.

With the cooling of the Universe, the radiation energy density diminishes
rapidly, leading to a transition from a radiation-dominated era to a matter-
dominated era. In the matter-dominated era, where the scale factor follows
$a\propto t^{2/3}$, the time-temperature relation can be expressed as
$dt=-2H^{\prime}dT/(3H^{2})$, where $H^{\prime}=dH/dT$, and the annihilation
rate becomes zero, namely, $M(T)=0$. In this scenario, Eq. (13) becomes:

$\frac{dn^{(1)}_{\gamma}}{dt}+3Hn_{\gamma}^{(1)}=\frac{\partial
n^{(1)}_{\gamma}}{\partial T}\frac{dT}{dt}+3Hn_{\gamma}^{(1)}=0~{},$ (28)

yielding

$\frac{dQ}{dT}=Q\left(-\frac{3}{T}+\frac{2H^{\prime}}{H}\right)~{}.$ (29)

The solution to this differential equation is given by:

$Q(T)=\frac{C_{0}H^{2}}{T^{3}}~{},$ (30)

where the constant $C_{0}$ remains to be determined. In the matter-dominated
era, $H^{2}$ is proportional to the energy density of matters, and for non-
relativistic species, this energy density exhibits an exponential dependence
on temperature. Consequently, the $Q$ value in Eq. (30) decays rapidly as the
temperature lowers. This leads to a convergence of the deviated $Q$ value
toward zero in the matter-dominated era. A more detailed study of this
recovery mechanism, taking into account the relevant micro-processes for
photons, is left for future work.

###  V.2 Remarks on prior study

As mentioned earlier, Ref. [6] addressed the primordial lithium problem by
using Tsallis statistics to parameterize a distorted photon distribution
during the BBN epoch. It led to a specific value of $q\simeq 1.027$ at
$T\simeq 4\times 10^{8}\,\mathrm{K}$, indicating a deviation from the Planck
distribution. Such a deviation was a key to resolving the primordial lithium
problem. Our current work not only verifies these findings but also
significantly extends them. While the resulting $q(T)$ in this study is
markedly dependent on the chosen initial conditions, certain choices yield
results that align well with our previous study. Importantly, our studies
challenge the SBBN model by suggesting that the photon distribution can
deviate from the Planck distribution. It could be a crucial point that refines
the SBBN and opens up new avenues for understanding the early Universe and the
primordial abundances.

While the choice of initial conditions continues to be a subject for further
exploration, the consistency in the qualitative behavior across both studies
demonstrates the idea that deviations from the Planck distribution of photons
could be intrinsic features of the early Universe.

##  VI Conclusion

In this study, we investigated the temporal evolution of the entropic index
$Q$ in Tsallis statistics, employing the first-order Boltzmann equation with
pertinent collision terms. Our findings indicate that during the BBN epoch,
collisional processes for the freeze-out of $e^{+}e^{-}$ pairs can yield a
decrement in $Q$, thereby causing a distortion in the high-energy tail of the
photon distribution. To achieve a physically meaningful solution for $Q$ that
resolves the primordial lithium problem, we optimized the initial condition of
the Boltzmann equation, letting the initial fluctuation be minuscule but non-
vanishing. As a consequence, the decline in $Q$ ceases around $T\leq
0.02\,\mathrm{MeV}$, attributable to the decay of the collision term,
culminating in the stabilization of $Q$ at the terminal stage of the BBN
epoch.

Furthermore, we discussed the long-term behavior of $Q$ as the Universe
transitions into the matter-dominated era. Preliminary analysis suggests that
the influence of $Q$ on the cosmic expansion rate diminishes during this
period, hinting at a potential convergence of the deviated $Q$ towards zero.
This aspect requires further exploration, particularly to elucidate the micro-
processes contributing to this recovery mechanism in the evolution of the
Universe. Our results imply that the variations in $Q$ predicted by our model
could be detectable exclusively through BBN observations. A more comprehensive
analysis, encompassing the evolution of $Q$ throughout the early cosmic
history – including both the BBN epoch and subsequent eras – is essential for
a deeper understanding of the photon distribution and the fundamental physics
governing the early Universe.

It is noteworthy that observational evidence for tiny fluctuations during the
BBN era is not directly available, unlike the CMB, which offers a clear
snapshot of temperature fluctuations at the time of recombination. The
temperature fluctuations during BBN may have been precursors to the density
and temperature fluctuations observable in the CMB. Exploring the link between
these two epochs can yield valuable insights into the evolution of the early
Universe, particularly regarding inflationary theory and the development of
cosmic structures.

##  Appendix A Collision integral for Compton scattering

For the Compton scattering of
$\gamma(E_{\gamma},\mathbf{p}_{\gamma})+e^{-}(E_{1},\mathbf{p}_{1})\rightleftharpoons\gamma(E_{2},\mathbf{p}_{2})+e^{-}(E_{3},\mathbf{p}_{3})$,
Eq. (10) can be written as

$\frac{g_{\gamma}}{(2\pi)^{3}}\int\mathbb{C}_{\text{comp}}[f_{\gamma}(\mathbf{p}_{\gamma})]\frac{d\mathbf{p}_{\gamma}}{E_{\gamma}}=-\int\left(\prod_{i=\gamma,1,2,3}d\Pi_{i}\right)(2\pi)^{4}\delta^{(3)}(\mathbf{p}_{\gamma}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})\,\delta(E_{\gamma}+E_{1}-E_{2}-E_{3})\\\
\times|\mathcal{M}|^{2}_{\text{comp}}\left[f_{\gamma}(\mathbf{p}_{2})f_{e^{-}}(\mathbf{p}_{3})\left(1+f_{\gamma}(\mathbf{p}_{\gamma})\right)\left(1-f_{e^{-}}(\mathbf{p}_{1})\right)-f_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{1})\left(1+f_{\gamma}(\mathbf{p}_{2})\right)\left(1-f_{e^{-}}(\mathbf{p}_{3})\right)\right]~{},$
(31)

In the non-relativistic limit, the squared matrix elements
$|\mathcal{M}|^{2}_{\text{comp}}$ can be approximated as a constant $M$. Pauli
blocking and Bose-Einstein condensation effects can be neglected, simplifying
the terms $1-f_{e^{-}}$ and $1+f_{\gamma}$ to approximately one,

$f_{\gamma}(\mathbf{p}_{2})f_{e^{-}}(\mathbf{p}_{3})\Big{(}1+f_{\gamma}(\mathbf{p}_{\gamma})\Big{)}\Big{(}1-f_{e^{-}}(\mathbf{p}_{1})\Big{)}-f_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{1})\Big{(}1+f_{\gamma}(\mathbf{p}_{2})\Big{)}\Big{(}1-f_{e^{-}}(\mathbf{p}_{3})\Big{)}\\\
\simeq
f^{(1)}_{\gamma}(\mathbf{p}_{2})f_{e^{-}}(\mathbf{p}_{3})-f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{1})~{}.$
(32)

Then, the integral becomes

$\displaystyle\frac{g_{\gamma}}{(2\pi)^{3}}\int\mathbb{C}_{\text{comp}}[f(\mathbf{p}_{\gamma})]\frac{d\mathbf{p}_{\gamma}}{E_{\gamma}}$
$\displaystyle\simeq\frac{M}{(2\pi)^{8}m_{e}^{2}}\int\frac{d\mathbf{p}_{\gamma}}{p_{\gamma}}\int
d\mathbf{p}_{1}\int\frac{d\mathbf{p}_{2}}{p_{2}}\int d\mathbf{p}_{3}$
$\displaystyle\qquad\times\delta^{(3)}(\mathbf{p}_{\gamma}+\mathbf{p}_{1}-\mathbf{p}_{2}-\mathbf{p}_{3})\,\delta(E_{\gamma}+E_{1}-E_{2}-E_{3})$
$\displaystyle\qquad\times\left[f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{1})-f^{(1)}_{\gamma}(\mathbf{p}_{2})f_{e^{-}}(\mathbf{p}_{3})\right]$
(33)
$\displaystyle\simeq\frac{M}{(2\pi)^{8}m_{e}^{2}}\int\frac{d\mathbf{p}_{\gamma}}{p_{\gamma}}\int\frac{d\mathbf{p}_{2}}{p_{2}}\int
d\mathbf{p}_{3}$
$\displaystyle\qquad\times\delta\left(p_{\gamma}+\frac{\left|\mathbf{p}_{2}+\mathbf{p}_{3}-\mathbf{p}_{\gamma}\right|^{2}-\mathbf{p}^{2}_{3}}{2m_{e}}-p_{2}\right)$
$\displaystyle\qquad\times\left[f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{2}+\mathbf{p}_{3}-\mathbf{p}_{\gamma})-f^{(1)}_{\gamma}(\mathbf{p}_{2})f_{e^{-}}(\mathbf{p}_{3})\right]~{}.$
(34)

The momenta could be approximated as $\mathbf{p}_{\gamma}\sim\mathbf{p}_{2}$
and $\mathbf{p}_{1}\sim\mathbf{p}_{3}$ for small momentum transfer via Compton
processes, which is relevant when the mass is heavier than the photon energy
at low temperature. Therefore,

$\displaystyle\frac{g_{\gamma}}{(2\pi)^{3}}\int\mathbb{C}_{\text{comp}}[f(\mathbf{p}_{\gamma})]\frac{d\mathbf{p}_{\gamma}}{E_{\gamma}}$
$\displaystyle=\frac{M}{(2\pi)^{8}m_{e}^{2}}\int\frac{d\mathbf{p}_{\gamma}}{p_{\gamma}}\int\frac{d\mathbf{p}_{2}}{p_{2}}\int
d\mathbf{p}_{3}\,\delta(p_{\gamma}-p_{2})\,\left[f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{3})-f^{(1)}_{\gamma}(\mathbf{p}_{2})f_{e^{-}}(\mathbf{p}_{3})\right]$
$\displaystyle=\frac{M}{(2\pi)^{8}m_{e}^{2}}\int\frac{d\mathbf{p}_{\gamma}}{p_{\gamma}}\int
4\pi p_{\gamma}dp_{\gamma}\int 4\pi
p^{2}_{1}dp_{1}\,\left[f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{3})-f^{(1)}_{\gamma}(\mathbf{p}_{\gamma})f_{e^{-}}(\mathbf{p}_{3})\right]$
$\displaystyle=0~{}.$ (35)

Hence, the Compton scattering does not contribute to the collision integral
for $f^{(1)}_{\gamma}$.

###### Acknowledgements.

D.J. was partly supported by National Natural Science Foundation of China (No.
12335009), and the National Key R$\&$D Program of China (2022YFA1602401).
M.H.M. was supported by the National Research Foundation of Korea (NRF) grant
funded by the Korea government (Grant No. NRF- 2021R1F1A1060066). M.K.C. was
supported by NRF under grant numbers, NRF-2021R1A6A1A03043957 and
NRF-2021R1F1A1046575. M.H.M. and M.K.C. were also supported by
NRF-2020R1A2C3006177. K.K. was supported by NRF-2022R1F1A1073890 and 2023 UBSI
Innovation Project of UNIST. Y.M.K. was supported by NRF-2022R1C1C2012226.

## References

  * [1] M. Kusakabe, K. S. Kim, M. K. Cheoun, T. Kajino, Y. Kino and G. J. Mathews, Astrophys. J. Suppl. 214, 5 (2014).
  * [2] J. C. Howk, N. Lehner, B. D. Fields and G. J. Mathews, Nature 489, 121 (2012).
  * [3] C. A. Bertulani, J. Fuqua and M. S. Hussein, Astrophys. J. 767, 67 (2013).
  * [4] S. Q. Hou, J. J. He, A. Parikh, D. Kahl, C. A. Bertulani, T. Kajino, G. J. Mathews and G. Zhao, Astrophys. J. 834, no. 2, 165 (2017)
  * [5] M. Kusakabe, T. Kajino, G. J. Mathews and Y. Luo, Phys. Rev. D 99, no.4, 043505 (2019).
  * [6] D. Jang, Y. Kwon, K. Kwak and M. K. Cheoun, Astron. Astrophys. 650, A121 (2021).
  * [7] U. Tirnakli, F. Büyükkiliç and D. Demirhan, Physica A 240, 657 (1997).
  * [8] C. Tsallis, J. Statist. Phys. 52, 479-487 (1988).
  * [9] E. W. Kolb and M. S. Turner, Front. Phys. 69, 1-547 (1990).
  * [10] S. D. McDermott and M. S. Turner, arXiv:1811.04932 [hep-ph].
  * [11] L. C. Thomas, T. Dezen, E. B. Grohs and C. T. Kishimoto, Phys. Rev. D 101, no.6, 063507 (2020).
  * [12] Y. Akrami et al. [Planck], Astron. Astrophys. 641, A7 (2020).


11institutetext: 1Departament de Física, Universitat de les Illes Balears
(UIB), E-07122, Spain
Institute of Applied Computing & Community Code (IAC3), UIB, Spain
11email<EMAIL_ADDRESS>
2Departament de Ciències Matemàtiques i Informàtica, Universitat de les Illes
Balears (UIB), E-07122, Spain
Institute of Applied Computing & Community Code (IAC3), UIB, Spain

# One dimensional prominence threads

I. Equilibrium models

J. Terradas 11 M. Luna 11 R. Soler 11 R. Oliver 11 M. Carbonell 22 J. L.
Ballester 11

###### Abstract

Context. Threads are the building blocks of solar prominences and very often
show longitudinal oscillatory motions that are strongly attenuated with time.
The damping mechanism responsible for the reported oscillations is not fully
understood yet.

Aims. To understand the oscillations and damping of prominence threads it is
mandatory to investigate first the nature of the equilibrium solutions that
arise under static conditions and under the presence of radiative losses,
thermal conduction and background heating. This provides the basis to
calculate the eigenmodes of the thread models.

Methods. The nonlinear ordinary differential equations for hydrostatic and
thermal equilibrium under the presence of gravity are solved using standard
numerical techniques and simple analytical expressions are derived under
certain approximations. The solutions to the equations represent a prominence
thread, i.e., a dense and cold plasma region of a certain length that connects
with the corona through a prominence corona transition region (PCTR). The
solutions can also match with a chromospheric-like layer if a spatially
dependent heating function localised around the footpoints is considered.

Results. We have obtained static solutions representing prominence threads and
have investigated in detail the dependence of these solutions on the different
parameters of the model. Among other results, we have shown that multiple
condensations along a magnetic field line are possible, and that the effect of
partial ionisation in the model can significantly modify the thermal balance
in the thread and therefore their length. This last parameter is also shown to
be comparable to that reported in the observations when the radiative losses
are reduced for typical thread temperatures.

###### Key Words.:

Magnetohydrodynamics (MHD) — waves — Sun: magnetic fields

##  1 Introduction

A recent survey on longitudinal oscillations in solar filaments carried out by
Luna et al. (2018) has provided interesting results about the temporal
attenuation of the oscillatory motions. A measure of the attenuation is the
ratio between the damping time, assuming an exponential decay, and the period
of the oscillation. The mean value of this parameter over 106 small amplitude
events (with velocities below 10 km $\rm s^{-1}$) is 1.75, while for large
amplitude oscillations (above 10 km $\rm s^{-1}$) the 96 events give a mean
value of 1.25. This means that typically the oscillations do not last more
than 2 periods. The question that arises is the mechanism that produces such
strong damping, and several approaches can be used to investigate this
problem.

The first approach is to consider that the system is near equilibrium
satisfying the energetic balance between radiation losses, thermal conduction
and heating. The model may include the effect of the magnetic field and also
the gravitational force. Once the equilibrium is calculated, and this is the
main motivation of the present work (Paper I), the problem of linear and
nonlinear waves on these equilibrium configurations can be then addressed
(Paper II). Only a few works have focused on the determination of a static
equilibrium under thermal balance: Degenhardt & Deinzer (1993) modelled a
quiescent prominence assuming balance between heating and radiative losses but
ignoring heat conduction. These authors found reasonable values for prominence
temperatures and densities but significantly shorter threads were obtained (a
similar result was obtained by Milne et al., 1979). The connection with the
chromosphere was not included in the model. Later, Dahlburg et al. (1998)
demonstrated the role of a dipped geometry to support a prominence
condensation against gravity and included a localised footpoint heating to
match their solution with a chromospheric layer. Regarding the oscillations
Schmitt & Degenhardt (1995) and Rempel et al. (1999) performed a stability
analysis of a flux tube model based on the work of Degenhardt & Deinzer (1993)
under line-tying conditions. However, in these studies the perturbations were
assumed to be adiabatic and therefore no attenuation was reported, although
some hints of instability were found. The inclusion of nonadiabatic effects
under different boundary conditions will be addressed in Paper II. In
addition, the effect of partial ionisation, effective for plasma temperatures
below 10,000 K, has not been taken into account in the above mentioned
studies.

The second approach to investigate thread or prominence oscillations is based
on the analysis of a dynamical system that undergoes thermal nonequilibrium.
Oscillations are naturally produced in the system and their investigation
provides details about the damping processes. Along this line, Luna & Karpen
(2012) investigated the oscillations of multiple threads formed in long,
dipped flux tubes through the thermal nonequilibrium process previously
simulated (Luna et al., 2012). These authors found that the oscillation
properties predicted by their simulations are in agreement with the observed
behaviour, and that the main restoring force is the projected gravity along
the tube where the threads oscillate. Zhang et al. (2012, 2013) modelled
impulsive heating at one leg of the loop and an impulsive momentum deposition,
which cause the prominence to oscillate. The oscillation damps with time under
the presence of non-adiabatic processes and these authors concluded that
radiative cooling is the dominant factor leading to damping.

In the present paper we follow the first approach. We construct improved
static prominence/thread models by considering different forms of the
radiative losses under the presence of gravity. The basic conditions to have a
temperature minimum at the thread centre are derived, and the possibility of
having several condensations along the field line is discussed. The connection
of the thread solution with the prominence corona transition region (PCTR) and
eventually with a chromospheric layer near the footpoints is also
investigated. Partial ionisation is included in the model and its effect on
the obtained solutions is studied. The work presented here is the basis of the
analysis of the eigenmodes described in the forthcoming Paper II.

##  2 Basic model

We assume that the magnetic field dominates the plasma and reduce the problem
to a one-dimensional thread in which the magnetic field determines the
geometry but is not modified by the presence of the thread. The plasma
quantities change along the thread and therefore we concentrate on the forces
and equilibrium along the field only. This assumption is valid for prominences
embedded in sufficiently strong magnetic fields (e.g. Fiedler & Hood, 1992;
Hillier & van Ballegooijen, 2013; Jenkins et al., 2019). For a model with
equilibrium normal to the magnetic field see for example, Ballester & Priest
(1989). For simplicity here we concentrate on a 1D problem and the tube cross-
section is assumed to be constant. We adopt a symmetric geometry about the
midpoint with a dipped region in the central part as in Antiochos & Klimchuk
(1991), Dahlburg et al. (1998) and Zhou et al. (2014). The definition of the
field geometry is given by Eq. (6) of Dahlburg et al. (1998). We denote by
$g(s)$ the gravity acceleration parallel to the field line and $s$ the
distance along this line (starting at the midpoint of the structure, i.e., the
centre of the thread). Three different examples of field geometry investigated
in the present work are represented in Fig. 1.

Figure 1: Sketch of the assumed magnetic field geometry. The different
parameters allow to change the shape and curvature of the field lines where
the thread is allocated (based on Dahlburg et al., 1998). The configurations
are denoted by $A$, $B$ and $C$. The three configurations have been
represented with different total lengths for visualisation purposes only. The
parameters of each curve following the definition in Dahlburg et al. (1998)
are A (bottom curve): $s_{1}=0.05L$, $s_{2}=0.5L$, $d=0.1L$, B (middle curve):
$s_{1}=0.1L$, $s_{2}=0.5L$, $d=0.05L$, and C (top curve): $s_{1}=0.2L$,
$s_{2}=0.5L$, $d=0.01L$. $L$ is half the length of the magnetic field line.

##  3 Fully ionised plasma

For a static situation two equations for the force balance and thermal
equilibrium must be satisfied simultaneously. Gas pressure under static
equilibrium must satisfy that

$\displaystyle\frac{dp}{ds}(s)=\rho(s)g(s).$ (1)

When gravity is neglected gas pressure remains constant along $s$. This
equation is completed with the ideal gas law

$\displaystyle p(s)=\frac{1}{\tilde{\mu}}\frac{k_{\rm B}}{m_{\rm
p}}\rho(s)T(s),$ (2)

where $k_{\rm B}/m_{\rm p}$ is the ideal gas constant and $\tilde{\mu}$ is the
mean atomic weight. For a fully ionised hydrogen plasma $\tilde{\mu}=1/2$. The
modifications in the previous equation because of partial ionisation will be
introduced later in Sect. 4.

Using the ideal gas law we can eliminate the variable $p(s)$ to obtain,

$\displaystyle\frac{d\rho}{ds}(s)T(s)+\rho(s)\frac{dT}{ds}(s)=\tilde{\mu}\frac{m_{\rm
p}}{k_{\rm B}}\rho(s)g(s).$ (3)

For thermal equilibrium between conduction, radiative losses and heating the
next nonlinear ordinary differential equation has to be satisfied,

$\displaystyle\frac{d}{ds}\left(\kappa_{\parallel}(s)\frac{dT}{ds}(s)\right)-\rho(s)^{2}\Lambda(T(s))+E_{0}=0,$
(4)

where $\kappa_{\parallel}(s)=\kappa_{0}\,T(s)^{5/2}$, with $\kappa_{0}$ a
constant coefficient. $\Lambda(T)$ is the radiative loss function and $E_{0}$
is the background heating. It is instructive to write the previous energy
equation in integral form using the Gauss theorem in one dimension,

$\displaystyle\left[\kappa_{\parallel}(s)\frac{dT}{ds}(s)\right]_{s=\pm
L}+\int_{-L}^{L}\left(-\rho(s)^{2}\Lambda(T(s))+E_{0}\right)ds=0,$ (5)

where the first term of this equation represents the (minus) heat flux through
the footpoints of the field lines located at $\pm L$ and the second term is
the spatially integrated contribution of radiation and heating. If $E_{0}$ is
zero the heat flux through the boundaries because of thermal conduction, is
equal to the total energy lost by radiation in the whole domain.

Figure 2: Radiative loss functions used in the present work, Athay (1986)
(continuous line) and Hildner (1974) (dotted line). Athay’s function has
reduced losses in comparison with Hildner for typical temperatures in the
range 4,000-12,000 K.

For optically thin radiative losses the function $\Lambda(T)$ has several
parametrisations in the literature that are essentially valid for the corona
surrounding prominences. Nevertheless, threads are optically thick and
therefore radiative losses are expected to be greatly reduced. This effect can
be taken into account by artificially decreasing the values in the optically
thin radiative losses for low temperatures ($T<10^{4}K$), see for example
Milne et al. (1979), Schmitt & Degenhardt (1995). Here we consider two main
radiative loss functions, described in Athay (1986) and Hildner (1974) and
represented in Fig. 2 for comparison purposes. Athay’s radiation function has
the advantage of being an analytical function with continuous derivatives and
attains quite low values under prominence/thread conditions. As we show in the
next sections this has important consequences regarding the lengths of the
calculated threads. We have explored other radiative losses such as Klimchuk-
Raymond’s fits (Klimchuk & Cargill, 2001) and the parametrisation computed
from CHIANTI atomic database (Dere et al., 1997; Landi et al., 2012) but the
results are in general quite similar to those of Hildner’s function.

We derive the conditions to obtain a dense and cold thread at the centre of
the dip surrounded by a hot plasma. We denote the temperature and density at
this point as $T_{0}=T(0)$ and $\rho_{0}=\rho(0)$. By symmetry around $s=0$ we
impose that $dT/ds=0$, meaning that Eq. (4) at $s=0$ reduces to

$\displaystyle\frac{d^{2}T}{ds^{2}}(0)=\frac{\rho_{0}^{2}\,\Lambda(T_{0})-E_{0}}{\kappa_{\parallel}(0)}.$
(6)

Since we are seeking for solutions that represent a cold thread connecting
with the hot corona, the temperature at $s=0$ (where there is an extrema
because we have imposed that $dT/ds=0$) must have a minimum, i.e.,
$d^{2}T/ds^{2}>0$. According to Eq. (6) this condition is satisfied only if
$E_{0}<\rho_{0}^{2}\,\Lambda(T_{0})$. Conditions for the existence of
prominence solutions were already studied in some detail in the early work of
Milne et al. (1979). Note that even in the unlikely situation with no
background heating in the solar atmosphere, i.e., $E_{0}=0$, it is still
possible to obtain physically meaningful solutions. However, for
$E_{0}>\rho_{0}^{2}\,\Lambda(T_{0})$ no solutions representative of threads
(i.e., cold material surrounded by coronal plasma) are found. Instead, coronal
loop solutions of the type studied by, for example, Klimchuk et al. (2010) or
Mikić et al. (2013) are obtained. In fact, the stationary solutions calculated
by these last authors by solving the time dependent problem have been used as
a check of our numerical method based on standard routines. We have obtained
static equilibrium solutions that match well with the stationary solutions of
Mikić et al. (2013).

Typical numbers for background heating found in the literature are in the
range $10^{-4}-10^{-5}$ W m-3 (e.g., Karpen et al., 2001; Mikić et al., 2013)
however, for the specific radiative losses considered in this work, these
values are in general too high to satisfy the condition
$E_{0}<\rho_{0}^{2}\,\Lambda(T_{0})$ necessary to have a typical
thread/prominence with density $\rho_{0}=10^{-11}$ kg m-3 and temperature
$T_{0}=10,000$ K. For this reason, we have decided not to focus on a single
value of the background heating but changing the value of $E_{0}$ in the range
between 0 and slightly below $\rho_{0}^{2}\,\Lambda(T_{0})$. It is clear that
the value of the radiative loss function at $T_{0}$ is a relevant number that
determines the range of allowed values for $E_{0}$. If the radiative losses
are low for temperatures around $10,000$ K the background heating needs to be
decreased accordingly. For this reason we consider background heatings as low
as $E_{0}=2.75\times 10^{-9}$ W m-3 and $E_{0}=1.25\times 10^{-6}$ W m-3 for
Athay and Hildner’s radiation functions but in some cases the values are
slightly above these values (e.g., Fig. 3).

The system of Eqs. (3) and (4) are solved numerically for $T$ and $\rho$ under
given boundary conditions using standard numerical techniques based on a
variable-order, variable-step Adams method. At the centre of the thread
($s=0$) the temperature ($T_{0}$) and density ($\rho_{0}$) values are selected
and $dT/ds=0$ is imposed. The two coupled ordinary differential equations are
numerically integrated from $s=0$ to $s=L$ using numerical methods. When
necessary, adaptive mesh methods are used, this is especially important at the
PCTR. We do not need, in general, to implement a shooting method since the
three boundary conditions are imposed at $s=0$. The shooting method is only
required when density at $s=0$ has to be determined to match the temperature
at the corona ($s=L$), see end of Sect.3.3. In this last case we use a
Runge–Kutta–Merson method and a Newton iteration in the shooting and matching
technique.

Figure 3: Hydrostatic and thermal equilibrium along the field line with zero
gravity. The continuous line corresponds to Athays’s radiation function, the
dashed line represents the results for Hildner’s function, while the
background heatings are $E_{0}=2.75\times 10^{-9}$ W m-3 and $E_{0}=1.75\times
10^{-6}$ W m-3. In this particular example $T_{0}=10^{4}\,\rm K$,
$\rho_{0}=10^{-11}\,\rm kg\,m^{-3}$ and the total length of the field line is
$2L=200\,\rm Mm$. Model $A$ is used in this plot.

###  3.1 Constant gas pressure

For a better comprehension of the results we start with the case with zero
gravity. Gas pressure is constant along the field under such conditions, and
the solution of the equations is represented in Fig. 3 for the two radiation
functions. The solutions show a rather different behaviour although the
reference temperatures and densities are the same at $s=0$. For Athay’s
radiation function a single thread around $s=0$ is found and matches smoothly
through a PCTR a plasma that is close to coronal conditions (densities of the
order of $10^{-13}\,\rm kg\,m^{-3}$, and temperatures around $10^{6}\,\rm K$).
On the contrary, in the case of Hildner’s function, several (up to three) cold
and dense regions are found in the configuration. The densities and
temperatures of these thread-like solutions are chosen to be exactly the same
at the centre of the thread ($s=0$). The system displays periodic cold and
dense threads, and this behaviour also applies to Athay’s radiation function
but the corresponding spatial periodicity is much longer than the length of
the system.

The distance between successive threads depends on the reference values for
temperature and density and also on the value of the heating constant. It is
clear that for constant gas pressure the system has a characteristic spatial
scale and it is worth investigating the origin of this length. It turns out
that this scale is closely related to the so-called Field’s wavenumber after
Field (1965). This author found that the thermal mode in a uniform plasma with
constant temperature of the order of MK is always unstable if thermal
conduction is absent. But under the presence of thermal conduction there is a
critical length scale, that we denote hereafter by $L_{\rm C}$, for the
stabilisation of the thermal mode. In the model used by Field (1965), with
uniform temperature and density, the equilibrium satisfies that

$\displaystyle\rho_{0}^{2}\,\Lambda(T_{0})=E_{0},$ (7)

and for a given constant pressure value and using the gas law, density and
temperature are obtained from the previous equation, which must be solved
numerically in general. Since the equilibrium is isothermal there is no
conduction term in Eq. (7) but it is present in the perturbations.

An analysis of the perturbations in this configuration to understand the
features of the thermal mode has been done in the past by, for example, Field
(1965); van der Linden & Goossens (1991); Soler et al. (2011, 2012). The
obtained dispersion relation shows the existence of a critical length and only
those perturbations with wavelenghts below $L_{\rm C}$ are stable, the
expression for $L_{\rm C}$ (see Eq. (26a) in Field, 1965) is

$\displaystyle L_{\rm C}$ $\displaystyle=$ $\displaystyle
2\pi\sqrt{\frac{{\kappa_{\parallel}}}{{\frac{\rho}{T}}\left(\frac{\partial\rho^{2}\Lambda}{\partial\rho}\right)_{T}-\left(\frac{\partial\rho^{2}\Lambda}{\partial
T}\right)_{\rho}}}$ (8) $\displaystyle=$ $\displaystyle
2\pi\sqrt{\frac{{\kappa_{0}\,T_{0}^{5/2}}}{2\frac{\rho_{0}^{2}}{T_{0}}\Lambda(T_{0})-\rho_{0}^{2}\,\Lambda^{\prime}(T_{0})}}.$

Derivatives involving density and radiative losses at constant temperature and
constant density are present in the denominator of the previous expression and
are evaluated in the second part of the equation where
$\Lambda^{\prime}(T_{0})$ is the temperature derivative of the radiative
losses evaluated at $T_{0}$. It is worth to mention that in Eq. (8) the
explicit dependence on $E_{0}$ is absent, the reason being that in our model
the heating only affects the equilibrium but not the perturbations. Note that
Begelman & McKee (1990) define a Field’s length which is not the same as the
critical length used here (see also Koyama & Inutsuka, 2004; Sharma et al.,
2010). Begelman & McKee (1990) provide a relationship between the critical and
the Field’s length in their Eq. (4.16).

We apply Eq. (8) to the situation in Fig. 3 but keeping in mind that we are
going to compare the calculated inhomogeneous equilibrium with a model that
has constant density and temperature. In order to perform a reliable
comparison we chose the same gas pressure value (which is constant along the
tube) and the same background heating in the two models. Then we calculate the
corresponding density and temperature that satisfies Eq. (7) for the
homogeneous model. For example, for Hildner radiation function we find a
temperature of $342536\,\rm K$ and a density of $2.9\times 10^{-13}\,\rm
kg\,m^{-3}$ in the homogeneous model. If we compute the mean values for the
corresponding inhomogeneous model (Fig. 3) we obtain a mean temperature of
$418897\,\rm K$ and a mean density of $2.8\times 10^{-13}\,\rm kg\,m^{-3}$,
these numbers are similar to those calculated for the homogeneous case. The
obtained values of temperature and density for the homogeneous case are
introduced in Eq. (8). We find that the corresponding scale length for Hildner
radiation function is 36 Mm. This value agrees reasonable well with the
periodicity found in Fig. 3 for Hildner’s function, which is around 43 Mm.
Repeating the same procedure for Athay’s radiation function (calculating again
the temperature and density values) we find that in this case $L_{\rm
C}=22881$ Mm. This large value explains in Fig. 3 the lack of periodicity in a
length of 100 Mm. Therefore, we conclude that the critical length provides a
reasonable estimation of the expected periodicity that makes the system stable
regarding the thermal instability. Furthermore, we have computed different
numerical solutions by changing the value of $\kappa_{0}$ and have positively
checked that the obtained characteristic lengths are proportional to the
square root of the conduction coefficient, as expected from Eq. (8). We have
varied other parameters such as temperature and density and the results
confirm that Eq. (8) is an adequate approximation. Note that the role of sound
waves, and therefore, the effect of pressure variations is absent in the
definition of the critical length while in the full numerical solutions of
Fig. 3 the effect of gas pressure is included.

Now we focus on another characteristic scale in the system but of local
nature, the length of the individual threads. As we have mentioned in Section
1 the calculated lengths by Milne et al. (1979) and Degenhardt & Deinzer
(1993) are short in comparison with the measured thread lengths. To
investigate this question and since the dense plasma representing a thread is
located around the origin of our coordinate system we assume that the
temperature around $s=0$ can be written as a second-order series expansion of
the form

$\displaystyle
T(s)=T_{0}\left(1+b_{1}\frac{s}{L}+b_{2}\frac{s^{2}}{L^{2}}\right),$ (9)

where the dimensionless coefficients $b_{1}$ and $b_{2}$ need to be
determined. Additional terms in Eq. (9) are neglected since $s/L\ll 1$ and we
perform a local analysis. The boundary condition $dT/ds=0$ at $s=0$ yields to
$b_{1}=0$. For the plasma density we have that using the gas law for constant
pressure and in the situation $s/L\ll 1$ we can write

$\displaystyle\rho(s)\approx\rho_{0}\left(1-b_{2}\frac{s^{2}}{L^{2}}\right).$
(10)

We substitute these expansions in $s$ for temperature and pressure in the full
energy equation given by Eq. (4). Approximating the powers by the
corresponding Taylor series and using again the fact that $s/L\ll 1$ we find,
after some algebra and to zeroth order in $s$ that

$\displaystyle
b_{2}=\frac{1}{2}\frac{L^{2}}{\kappa_{0}\,T^{7/2}_{0}}\left(\rho_{0}^{2}\,\Lambda(T_{0})-E_{0}\right).$
(11)

This coefficient is in essence the ratio between the radiative minus the
heating term over the conduction term evaluated at $s=0$ and assuming that the
temperature changes on a spatial scale given by $L$. Equation (11) provides
the value of the term in front of the parabolic dependence on distance, $s$,
in Eq. (9) and it is used here as a proxy to obtain information about the
thread length. The smaller the value of $b_{2}$ the flatter the parabolic
curve and therefore the longer the length of the thread. The question that we
have to address is how to define a spatial scale, $l_{\rm th}$, associated
with a parabola. For this reason if we rewrite Eq. (9) as

$\displaystyle T(s)=T_{0}\left(1+\frac{s^{2}}{l_{\rm th}^{2}}\right),$ (12)

then

$\displaystyle l_{\rm
th}=\sqrt{\frac{L^{2}}{b_{2}}}=\sqrt{\frac{2\,\kappa_{0}\,T^{7/2}_{0}}{\rho_{0}^{2}\,\Lambda(T_{0})-E_{0}}},$
(13)

where we have used the expression for $b_{2}$ given in Eq. (11). The spatial
scale $l_{\rm th}$, can be understood as an approximation for the thread
length and provides valuable information about the dependence of this
parameter on: central density and temperature, radiative losses at this
temperature, conduction and background heating. The radiative losses are
dependent on the specific cooling table chosen as this will define the
behaviour and therefore shape of the PCTR. Since we are under the assumption
$E_{0}<\rho_{0}^{2}\,\Lambda(T_{0})$, the parameter $l_{\rm th}$ is always a
real positive number. According to Eq. (13) an increase in the conduction will
lead to longer threads while an increase in the radiative losses will shorten
them (if $E_{0}$ is constant). The dependence of $l_{\rm th}$ on the
conduction coefficient, $\kappa_{0}$, is exactly the same as for $L_{\rm C}$
(see Eq. (8)) although these two spatial scales have different physical
meanings. We will return to Eq. (13) in the forthcoming sections.

Figure 4: Hydrostatic and thermal equilibrium along the field line with
gravity. The continuous line corresponds to Athays’s radiation function, the
dashed line represents the results for Hildner’s function, while the
background heatings are $E_{0}=2.75\times 10^{-9}$ W m-3 and $E_{0}=1.25\times
10^{-6}$ W m-3. In this particular example $T_{0}=10^{4}\,\rm K$,
$\rho_{0}=10^{-11}\,\rm kg\,m^{-3}$ and the total length of the field line is
$2L=200\,\rm Mm$. Model $A$ is used in this plot.

###  3.2 Non-constant gas pressure, gravity included

When gravity is introduced in the system gas pressure is no longer constant
along the magnetic field line. We concentrate hereafter in Model $A$. For
large spatial scales we expect that gravity can significantly alter the
results in comparison to the constant pressure case, and indeed this is the
case. The numerical solution indicates that the second thread solution (the
one that was at $s\approx 45$ Mm in Fig. 3) for Hildner’s function tends to
have temperatures much below the minimum at the thread centre ($s=0$) and the
numerical integration of the equations fails to converge. The system does not
allow the periodicity found for the case with zero gravity for the same value
of the background heating. In this case no static solution is allowed in the
system for the selected parameters and a dynamical behaviour is expected
because of the thermal nonequilibrium.

However, by changing $E_{0}$ it is still possible to find, under some choice
of parameters, a situation with two threads under mechanical and thermal
balance. An example is shown in Fig. 4. Thus, it is possible to have some cold
material in equilibrium balance but not located at the dips. This cold and
dense material does not reach the temperature and density values of the
central thread though.

In Fig. 4, density and temperature as a function of position along the
magnetic field under gravity is also represented for Athay’s function.
Interestingly, in this case the differences with respect to the zero gravity
situation are not large, essentially the density tends to increase near the
footpoint as a result of the presence of gravity since around the footpoint
the gravity force is purely vertical and has the largest contribution (compare
with the continuous line in the top panel of Fig. 3). The thread length
obtained in this case from the simulations, hereafter denoted by $a$ (do not
confuse with $l_{\rm th}$, the analytical approximation), is $1.7\,\rm Mm$,
and it is calculated using the position where the density derivative with $s$
has a maximum ($a$ is twice this value). We find an extended prominence corona
transition region that eventually matches a plasma that is close to coronal
conditions (densities of the order of $10^{-13}\,\rm kg\,m^{-3}$ and
temperatures around $10^{6}\,\rm K$). For these reasons we conclude that the
obtained model is a good representation of a thread.

In Fig. 5 the contribution of the different terms in the energy equation is
represented for the same parameters of Fig. 4. The conduction term is always
positive and is essentially balanced by the radiative losses while the
background heating is very small in this example. Both conduction and
radiation have a strong peak at the PCTR where density and temperature change
abruptly.

Figure 5: Terms in the energy equation, Eq. (4), as a function of position at
the thread body. The continuous line corresponds to the conduction term, the
dashed to the radiation losses while the dotted line represents the constant
background heating. The same parameters as in Fig. 4 for Athay’s function have
been used. Model $A$ is used in this plot. Figure 6: Thread length, $a$, as a
function of the heating constant, $E_{0}$, for different reference densities
$\rho_{0}$. In this plot the reference temperature is $T_{0}=8\times 10^{3}$ K
and Athay’s radiative loss has been used. The analytical approximation given
by Eq. (13) is plotted with red dashed lines. Model $A$ is used in this plot.

###  3.3 A parametric survey

Once we know the main features of the solutions we study the influence of the
different parameters on the computed equilibrium for a wide range of values.
For this reason we have carried out a parametric study starting with the
dependence on $E_{0}$. Although from the results in Fig. 5 we may conclude
that the role of the background heating is small, it turns out that this
parameter has a strong influence on the obtained thread lengths, and as we
will show in Paper II, this affects significantly the damping times. This
effect of $E_{0}$ on the equilibrium was already noted by Schmitt & Degenhardt
(1995) and is further investigated here.

In Fig. 6 we have represented the numerically obtained thread lengths (twice
the position where the density derivative with $s$ has a maximum) as a
function of $E_{0}$ for a specific choice of temperature ($8\times 10^{3}$ K)
at the core of the thread. The thread length decreases when the heating is
reduced and the minimum thread length is achieved for zero background heating.
This agrees with the dependence of $l_{\rm th}$ on
$\rho_{0}^{2}\,\Lambda(T_{0})-E_{0}$ in the denominator of Eq. (13).
Analytical approximations for the thread length are more difficult to obtain
in this case but note that the projection of gravity at $s=0$ is precisely
zero, therefore as a first approximation we can still use the definition of
$l_{\rm th}$. The numerical results indicate that the longest thread lengths
are obtained for values of the heating tending to the radiative losses at the
thread centre, as expected from Eq. (13) too. Nevertheless, the analytical
approximation given by Eq. (13) and plotted in Fig. 6 with red dashed lines
overestimates the thread length in this last situation. It is worth to mention
that there are numerical problems regarding convergence in the calculation of
the solutions for values of the heating constant near to
$\rho_{0}^{2}\,\Lambda(T_{0})$, but the reason is clear from the denominator
of $l_{\rm th}$.

We focus now on the dependence of the thread length on the length of the field
lines. The three magnetic configurations, denoted by A, B and C, have been
analysed and the total length of the field lines has been changed in the range
$80-200$ Mm. We have found that for Athay’s radiative function the total
length of the threads, i.e. the parameter $a$, does not change much with the
length of the magnetic field lines and the values are typically of the order
of several Mm. Interestingly, the thread lengths are in the range of the
reported in the observations (e.g., Okamoto et al., 2007; Arregui et al.,
2018). The variation of the thread length with models A, B, and C is at most
around 1%. We conclude that the geometry has some influence on the length of
the threads but is not too significant. Note that this effect is not included
in Eq. (13) since no gravity was assumed in the derivation of this equation
(in our model the geometry of the field is only included through the
projection of gravity along the magnetic field lines).

On the other hand, when Hildner’s radiative function is considered, the thread
lengths are short, typically of the order of $150\,\rm km$ only. Again the
magnetic field geometry has a weak influence on the results in this case. The
short lengths of the threads for Hildner’s function were already reported by
Milne et al. (1979) and Schmitt & Degenhardt (1995) and the cause is that the
radiative losses are too high for typical thread temperatures. On the
contrary, Athay’s radiative losses are significantly lower for typical thread
temperatures (see the comparison in Fig. 2) and this eventually leads to much
longer threads.

Figure 7: Thread length, $a$, as a function of the central temperature for
different reference densities. The radiative function is based on Athay
(1986). The analytical approximation given by Eq. (13) is plotted with red
dashed lines. Model $A$ is used in this plot.

In Fig. 7 the obtained thread length is represented as a function of the
central temperature for different reference densities. In these calculations
we have used Athay’s function and have imposed that $E_{0}=0$, therefore, as
we have demonstrated in Fig. 6 we concentrate on the shortest threads. Figure
7 indicates that cool threads have longer lengths than hot threads.
Furthermore, light threads have larger sizes than heavy threads which is in
agreement again with Eq. (13) ($l_{\rm th}$ is inversely proportional to
$\rho^{2}_{0}$), see the red dashed lines in Fig. 7.

Figure 8: Thread density at $s=0$, as a function of the central temperature
for models $A$, $B$, and $C$, matching a coronal temperature of 1 MK. A
shooting method has been used. The radiative function is the same as in Athay
(1986).

So far we have integrated the ordinary differential equations from $s=0$ to
$s=L$ because we have imposed three boundary conditions at $s=0$ (temperature,
density and zero derivative of temperature) and the problem is well defined.
The calculation of the maximum heating is straight forward since it involves
the density and temperature at the centre of the thread. Although in some
cases we have obtained coronal temperatures slightly below 1 MK in the corona
this approach provides reasonable results. However, it is also interesting to
investigate the situation where, apart from the thread temperature, the
coronal temperature is imposed and the integration provides the rest of the
parameters. We are mainly interested in the calculation of the density at the
centre of the thread that produces a perfect match for a coronal temperature
at $s=L$. In this case a shooting method is required to solve the differential
equations. Now the maximum heating constant is not known since it involves the
thread density that has to be determined from the solution of the equations.
An iterative approach is required to calculate the maximum $E_{0}$. To avoid
this complication we have concentrated on the situation $E_{0}=0$ and studied
the dependence of the results on temperature for the three different
geometrical configurations. The results are shown in Fig. 8, where the
obtained thread densities are plotted as a function of the thread
temperatures. The range of thread temperatures for the three configurations
agree well with the densities estimated from the observations (e.g., Tandberg-
Hanssen, 1995; Patsourakos & Vial, 2002). The cooler the thread the higher the
density is the behaviour found in the configuration and it is what one would
expect (using the gas law) from a set of solutions that have essentially the
same gas pressure at the centre of the thread. We have investigated the
situation when the heating constant is different from zero, and again it leads
to longer threads than the situation for $E_{0}=0$, as expected from Eq. (13)
and demonstrated in Fig. 6.

###  3.4 Matching a chromospheric layer

Up to now we have described a system composed of a cold and dense material
representing a thread connecting with a plasma under coronal conditions. A
more realistic model should include the connection with the chromosphere near
the footpoints of the magnetic field lines that are supposed to be anchored in
the photosphere. The physics of the chromosphere is complex and it is out of
the scope of this work to include a detailed modelling of this layer. However,
it is worth investigating the physical conditions required to have a layer
akin to the chromosphere in our simplified 1D system. We have already seen
that it is possible to obtain cold and dense plasma regions along the field
lines depending on the periodicity of the condensations that can have similar
conditions to chromospheric plasmas, especially under the presence of gravity.
Nevertheless, these chromospheric regions can occur anywhere along the field
line, while we are mostly interested in a chromosphere localised near the
footpoints, i.e., near $s=\pm L$ in our model. The way to force the coronal
part of the previously obtained solutions to match with a chromosphere is to
include a localised heating near the footpoints. This approach has been used
in the past by several authors and the most common form for the heating
function used in the literature is,

$\displaystyle E(s)=E_{0}\left(1+h_{\rm ch}\,e^{(s-L)/\lambda}\right),$ (14)

where $\lambda$ is a spatial scale typically of the order of 10 Mm and $h_{\rm
ch}$ is a factor that is in the range 20-100 (see for example, Dahlburg et
al., 1998; Karpen et al., 2001). Using these values in the spatially dependent
heating function and solving the ordinary differential equations we obtain a
solution with a fast decrease in temperature and a quick increase in density
near the footpoint. However, the rapidly changing nature of the chromospheric
part of the solution requires, in general, a special treatment of this layer
(see Lionello et al., 2009; Johnston & Bradshaw, 2019; Zhou et al., 2021, for
modifications in the conduction coefficient to avoid such difficulties).
Nevertheless, the adaptive methods that are implicitly used in our numerical
treatment do not fail to resolve the steep gradients in the solution. It is
important to remark that the condition $E(s=0)<\rho_{0}^{2}\,\Lambda(T_{0})$
still needs to be satisfied to have a thread-like solution.

The approach used here to incorporate a chromospheric layer in the model is to
assume that at some point the temperature remains constant and takes typical
chromospheric values (see for example Mikić et al., 2013; Karpen et al.,
2001). In this situation no energy equation is required, gas pressure must be
continuous at this point and the presence of gravity leads to an exponentially
growing density and pressure as we move downwards in the chromospheric part.
The integration of the equations is performed from the thread centre to the
chromosphere. An example of such computed equilibrium is shown in Fig. 9 for
the case $h_{\rm ch}=35$. In this equilibrium a temperature threshold of
$10^{4}$ K has been imposed in the isothermal chromospheric layer which starts
at $s=95$ Mm and its thickness is therefore 5 Mm (larger than the real
thickness which is typically around 2 Mm). The exponential increase in density
(straight line in the logarithmic scale of the plot) is clear in the top panel
of the figure and leads to mean chromospheric density values.

Figure 9: Hydrostatic and thermal equilibrium along the field line with
gravity for Hildner’s function. The background heating is $E_{0}=5\times
10^{-7}$ W m-3 and $h_{\rm ch}=35,15,0$ for the continuous, dashed, and dotted
lines. In this particular example $T_{0}=10^{4}\,\rm K$,
$\rho_{0}=10^{-11}\,\rm kg\,m^{-3}$ and the temperature threshold for the
imposed chromosphere is $10^{4}$ K. Model $A$ has been used in this example.

Using Hildner’s function and for comparison purposes we have included in Fig.
9 two solutions with reduced localised heating to visualise the nature of the
solutions. For the situations with $h_{\rm ch}=15$ and $h_{\rm ch}=0$ we find
how temperature decreases near the footpoint but is still above the
temperature threshold ($10^{4}\,\rm K$) to have a chromosphere. An interesting
conclusion from Fig. 9 is that the presence of footpoint heating does not
affect the thread properties and the solution is essentially the same around
$s=0$ for the three cases. This is also true when Athay’s radiative function
is considered. It is worth mentioning that multiple condensations, discussed
in Sect. 3.2, are not present in Fig. 9 because the heating constant is always
lower than the one used in Fig. 4.

The model developed in this section including a dense and cold chromosphere
near the footpoint will be useful for understanding the effect of this layer
on the attenuation of the waves (Paper II) produced by the possible mechanism
of wave reflection at the chromosphere and wave leakage if the boundary is
open.

##  4 Partially ionised plasma

Let us assume that the plasma is not fully ionised. Based on 1D non-LTE
radiative transfer models, Heinzel et al. (2014) calculated among other
parameters, the ionisation degree in several prominence slabs. In particular,
Heinzel et al. (2015) provide tables for the ionisation degree for different
temperatures and pressures at the prominence. The idea here is to use these
values in our calculations of the equilibrium and study how our models are
modified by the presence of neutrals. We are again in the situation where the
description using an optically thin plasma is not fully correct but it can be
considered as a starting point.

The plasma is assumed to be composed of hydrogen and helium. The abundance of
helium is $10\%$ and is not ionised. The ionisation degree is defined here as
$i=n_{\rm e}/n_{\rm H}$ where $n_{\rm e}$ is the electron density and $n_{\rm
H}$ the total hydrogen density ($n_{\rm H_{I}}+n_{\rm p}$). The total particle
number is defined as $N=n_{\rm H}+n_{\rm He}+n_{\rm e}$ and using the previous
definitions and the helium abundance ($n_{\rm He}=0.1\,n_{\rm H})$ it is
written as

$\displaystyle N=n_{\rm e}\left(1+\frac{1.1}{i}\right).$ (15)

Using the ideal gas law we find that the electron density in terms of
pressure, temperature and ionisation degree is (now we write the explicit
dependence of the variables)

$\displaystyle n_{\rm e}(s)=\frac{p(s)}{k_{\rm
B}T(s)}\left(\frac{i(p(s),T(s))}{i(p(s),T(s))+1.1}\right).$ (16)

The ionisation degree, $i$, depends on $p$ and $T$ and it is calculated from
Table 1 of Heinzel et al. (2015). To simplify things hereafter we write
$i(p(s),T(s))$ as $i(s)$.

For the total density we have that $\rho=n_{\rm H}\,m_{\rm p}+n_{\rm
He}\,4\,m_{\rm p}+n_{\rm e}\,m_{\rm e}$ (being $m_{\rm p}$ and $m_{\rm e}$ the
proton and electron masses). It reduces to the following expression when the
electron mass is neglected in front of the proton mass,

$\displaystyle\rho(s)=\frac{m_{\rm p}}{k_{\rm
B}}\frac{p(s)}{T(s)}\left(\frac{1.4}{i(s)+1.1}\right).$ (17)

Now the term in the brackets plays the role of $\tilde{\mu}$ in Eq. (2). In
the partially ionised situation it is more convenient to solve the equilibrium
equations for gas pressure and temperature because the ionisation degree
calculated in Heinzel et al. (2015) depends on these two magnitudes.

The equation for hydrostatic equilibrium, Eq. (1), is written in terms of
pressure and temperature because density has been eliminated using Eq. (17).
Partial ionisation changes the equation for thermal equilibrium in the
conduction and radiation terms. The factor in front of the radiative term is
now

$\displaystyle n_{\rm e}(s)\,n_{\rm H}(s)=\frac{n_{\rm e}^{2}(s)}{i(s)}.$ (18)

In the conductivities the contribution of neutrals has to be added to the
electron contribution. We have that according to Soler et al. (2010, 2012),

$\displaystyle\kappa_{\rm e}(s)$ $\displaystyle=$
$\displaystyle\kappa_{0}\,\xi_{\rm p}(s)\,T^{5/2}(s),$ (19)
$\displaystyle\kappa_{\rm n}(s)$ $\displaystyle=$
$\displaystyle\left(\kappa_{1}\,\xi_{\rm H_{I}}(s)+\kappa_{2}\,\xi_{\rm
He_{I}}\right)T^{1/2}(s),$ (20)

where the relative density of species are in our case, $\xi_{\rm
p}(s)=i(s)/1.4$, $\xi_{\rm H_{I}}(s)=(1-i(s))/1.4$ and $\xi_{\rm
He_{I}}=0.4/1.4$. For the conductivities we have that $\kappa_{0}=1.1\times
10^{-11}$ W m-1 K-7/2, $\kappa_{1}=2.24\times 10^{-2}$ W m-1 K-3/2, and
$\kappa_{2}=3.18\times 10^{-2}$ W m-1 K-3/2. The effective conduction
coefficient is now the sum of electron and neutral conductivities

$\displaystyle\kappa_{\parallel}(s)=\kappa_{\rm e}(s)+\kappa_{\rm n}(s).$ (21)

The conduction term involves spatial derivatives of this coefficient. The
coefficient $\kappa_{\parallel}$ depends on temperature and ionisation degree,
and this last parameter also depends on temperature and pressure. Regarding
the ionisation, from the practical point of view we have adjusted a second
order polynomial to the function $i(p,T)$ given in Heinzel et al. (2015)
(table for an altitude of 20 Mm). The calculation of the partial derivatives
of $i$ with $p$ and $T$ are straight forward using the polynomial fit.

Figure 10: Hydrostatic and thermal equilibrium along a thread under partial
ionisation (continuous lines) and full ionisation (dashed lines). In this
particular example, temperature ($T_{0}=8\times 10^{3}$ K) and pressure at
($s=0$) are the same for the partially and fully ionised situation. In these
solution we have imposed that $E_{0}=0$. The radiative function is the same as
in Athay (1986). Model $A$ has been used in this plot.

We have repeated some of the previous calculations including the effect of
partial ionisation through the change in the radiation and conduction terms
and the modification of the ideal gas law. Since we have shown that the
inclusion of a chromosphere has little effect on the conditions within the
threads we do not include this pseudo-layer in the present calculations.

The main result is that partial ionisation significantly increases the size of
the thread. An example is shown in Fig. 10, where the length of the thread is
four times longer than in the fully ionised case. We have to bear in mind that
now the initial values of the thread parameters at $s=0$, required to perform
the integration of the differential equations, are different to those in the
fully ionised case because of the change introduced in the gas law by partial
ionisation (see Eq. (17)). In the present example temperature and gas pressure
are the same at the thread centre but the reference density is different. In
Fig. 10 (bottom panel) we can see how the ionisation degree changes inside the
thread because of the implemented model of Heinzel et al. (2015). The
ionisation degree raises smoothly from 0.72 at the thread centre to 1 at the
edge of the thread where high (coronal) temperatures are achieved.

The increment in the thread length produced by partial ionisation is a
consequence of the changes in thermal conduction, which is larger owing to the
contribution of neutrals on the conductivities. For a better understanding of
this feature we proceed as in the fully ionised situation when gravity is
absent. Using the same approximation as in Eq. (9) we find that under the
presence of partial ionisation the characteristic spatial scale of the thread
is

$\displaystyle l_{\rm
th}=\sqrt{2\frac{\kappa_{0}\frac{i_{0}}{1.4}T^{7/2}_{0}+\left(\kappa_{1}\frac{1-i_{0}}{1.4}+\kappa_{2}\frac{0.4}{1.4}\right)T^{3/2}_{0}}{\rho_{0}^{2}\Lambda(T_{0})\frac{i_{0}}{1.4^{2}}-E_{0}}},$
(22)

where $i_{0}=i(p_{0},T_{0})$ is the ionisation degree at the centre of the
thread. The value of $i_{0}$ is typically around 0.72 but the reference
density, $\rho_{0}$, also depends on $i_{0}$ (see Eq. (17)). Indeed, it is not
difficult to calculate the ratio of the radiation term for partial ionisation
relative to full ionisation for constant reference pressure, $p_{0}$ and
temperature, $T_{0}$. For partial ionisation the reference density is

$\displaystyle\rho_{\rm 0p}=\frac{m_{\rm p}}{k_{\rm
B}}\frac{p_{0}}{T_{0}}\left(\frac{1.4}{i_{0}+1.1}\right),$ (23)

while for full ionisation

$\displaystyle\rho_{\rm 0f}=\frac{m_{\rm p}}{k_{\rm
B}}\frac{p_{0}}{T_{0}}\left(\frac{1}{2}\right).$ (24)

The ratio of the two corresponding radiation terms is

$\displaystyle\frac{\rho_{\rm
0p}^{2}\,\Lambda(T_{0})\frac{i_{0}}{1.4^{2}}}{\rho_{\rm
0f}^{2}\,\Lambda(T_{0})}=\frac{4\,i_{0}}{\left(i_{0}+1.1\right)^{2}}.$ (25)

This ratio is always smaller than one, meaning that radiative losses under
partial ionisation are reduced in comparison with the fully ionised case.
Since this term appears in the denominator of Eq. (22) it means that it
produces an increase in the thread length.

Figure 11: Conduction terms as a function of position (in the range 0 to 30 Mm
for visualisation purposes). Conduction by neutral H (dashed line) goes to
zero at the edge of the thread where the plasma is fully ionised. Conduction
by neutral He (dotted line) dominates inside the thread and becomes zero in
the corona where there is conduction by electrons only (continuous line). Same
parameters as in Fig. 10 have been used. Figure 12: Thread length as a
function of temperature under partial ionisation and full ionisation. In these
solution we have imposed that $E_{0}=0$. The radiative function is the same as
in Athay (1986). The analytical approximations given by Eqs. (22) (partial
ionisation) and (13) (full ionisation) are plotted with orange and red dashed
lines. Same parameters as in Fig. 10 have been used.

For thermal conduction we need to evaluate the three terms that appear in the
numerator of Eq. (22). The sum of the three terms is typically around 4.5
times the value for the fully ionised situation. This increment leads to
longer threads under partial ionisation. It turns out that the increased
conduction under partial ionisation is mostly produced by the presence of
neutral helium. This is visualised in Fig. 11 where the different conduction
terms are plotted as a function of position for a typical case. Inside the
thread, conduction by neutrals dominates conduction by electrons. The
contribution of neutral helium is the largest, although its abundance is just
$10\%$. Conduction by neutral H goes to zero as we approach the edge of the
thread where plasma is fully ionised. In the coronal medium the contribution
of neutral He is artificially forced to go to zero since at coronal
temperatures is fully ionised. Neutral He is therefore responsible, together
with the reduced radiation, for the rise of the obtained thread lengths.

Finally, Fig. 12 shows how the thread length changes with the central
temperature under partial ionisation. The behaviour is the same as in Fig. 7
for the fully ionised case, and for comparison purposes we have also
represented the results for full ionisation for the same reference pressure at
$s=0$. Interestingly, as the central temperature increases the two curves
approach each other since the larger the temperature the larger the ionisation
degree in the prominence. We find again longer threads under the presence of
partial ionisation, being the difference relative to full ionisation a factor
that varies from 5.3 to 2.6 in the range of temperatures of the plot. The
analytical approximation given by Eq. (22) is also plotted in the figure and
the agreement with the numerical result is remarkable.

##  5 Summary and conclusions

In the present work we have analysed the features of one-dimensional
equilibrium thread models under hydrostatic and thermal balance. We have
started with the situation without gravity and have progressively increased
the complexity of the model. This has allowed us a better comprehension of the
results. Gravitational stratification, the presence of a chromosphere and
finally partial ionisation effects have been incorporated to our model. The
main outcome of our study is summarised in the following:

  1. 1.

The value of the background heating in comparison with the radiative losses at
the centre of the thread is crucial to lead to thread-like solutions
surrounded by coronal plasma. Only if $E_{0}<\rho_{0}^{2}\,\Lambda(T_{0})$
static, cold and dense plasma threads under thermal equilibrium exist. In the
opposite situation when the background heating is higher than the radiative
losses at the centre of the structure only loop-like solutions are achieved.

  2. 2.

Several thread-like condensations are in principle possible along the field
lines and not necessarily located at the dips of the magnetic field. We have
found that the physical origin of the secondary condensations is related to
the critical length introduced by Field (1965) and related to the
characteristic spatial scale for thermal stability.

  3. 3.

The presence of gravity in the model can produce that the secondary
condensations collapse and no equilibrium configuration on a given length of
the magnetic field line is achieved. Again the value of the background heating
plays a major role in the behaviour of secondary thread-like solutions,
located in general near the footpoints. The gravitational stability of the
obtained solutions needs to be investigated and it will be addressed in Paper
II.

  4. 4.

A parametric survey has been carried out to understand the dependence of the
thread length, on the different values of the parameters. The geometry of the
field lines is not especially important but the radiative losses for low
temperatures are crucial to obtain realistic thread lengths. Athay’s radiative
function, with reduced losses under typical thread conditions (see Fig. 2), is
the most suitable choice from the ensemble of radiative losses analysed in the
present work. In any case since the plasma is optically thick under typical
prominence/thread conditions the incorporation of more realistic physics
requires to properly solve the radiative transport problem.

  5. 5.

We have derived a simple analytical expression for the characteristic spatial
scale of the thread, $l_{\rm th}$, under static equilibrium that explains the
dependence of the computed thread lengths, $a$, with the different parameters
of the model, including partial ionisation. This could be used in a novel way,
for example, to infer the ratio between radiation and heating if the length of
the thread and the central temperature and density are known from
observations.

  6. 6.

Significantly long threads are obtained when partial ionisation is present.
This is a consequence of the reduced radiation and increased conduction
produced by neutral helium in comparison to the fully ionised case. It is
assumed that helium is totally neutral in the thread, but this is not
completely true under real conditions (see Soler et al., 2010) and in reality
this may decrease the length of the threads. This needs to be addressed in the
future.

  7. 7.

The connection of the thread-corona system with a chromosphere is obtained
when heating around a certain threshold is localised near the footpoints, in
agreement with Dahlburg et al. (1998). The chromospheric model used here (a
stratified and isothermal layer under full ionisation) is a first
approximation and more physics needs to be included in future studies.
However, the presence of a chromosphere is interesting for the applications
regarding the damped oscillations that will be discussed in Paper II.

  8. 8.

Related to the previous point, we have shown that localised footpoint heating
does not significantly alter the temperature and density of the cold plasma
representing the thread relative to the case without footpoint heating, at
least for the parameters considered in this work. In our model the localised
heating essentially leads to the existence of a chromospheric layer only.

The numerical solutions presented in this work will be used to compute the
corresponding eigenmodes to compare the damping rates of our calculations and
the reported in the observations (Paper II). First, this will allow us to
better understand the damping mechanism of the oscillations and second this
information will be used as a method to constrain or to infer some
modifications to the models (see Anzer & Heinzel, 2008) and most likely about
the dependence of the radiative losses on temperature. The comparison of
theory and observations will be used for seismological purposes.

###### Acknowledgements.

We acknowledge the support from grant AYA2017-85465-P (MINECO/AEI/FEDER, UE),
to the Conselleria d’Innovació, Recerca i Turisme del Govern Balear, and also
to IAC3. We are grateful to the ISSI Team led by Manuel Luna “Large-Amplitude
Oscillations as a Probe of Quiescent and Erupting Solar Prominences” for
inspiring this work in the fruitful meetings held in Bern in 2018 and 2019. We
also thank ISSI-Beijing for hosting the team ”The eruption of solar filaments
and the associated mass and energy transport” led by J. C. Vial and P. F. Chen
where part of the results of this paper were presented. We are grateful to
Rafel Prohens from the Departament de Ciències Matemàtiques i Informàtica,
Universitat de les Illes Balears (UIB) for his advise on phase analysis of
nonlinear differential equations. We thank the anonymous referee for their
useful comments that helped to improve the manuscript.

## References

  * Antiochos & Klimchuk (1991) Antiochos, S. K. & Klimchuk, J. A. 1991, ApJ, 378, 372
  * Anzer & Heinzel (2008) Anzer, U. & Heinzel, P. 2008, A&A, 480, 537
  * Arregui et al. (2018) Arregui, I., Oliver, R., & Ballester, J. L. 2018, Living Reviews in Solar Physics, 15, 3
  * Athay (1986) Athay, R. G. 1986, ApJ, 308, 975
  * Ballester & Priest (1989) Ballester, J. L. & Priest, E. R. 1989, A&A, 225, 213
  * Begelman & McKee (1990) Begelman, M. C. & McKee, C. F. 1990, ApJ, 358, 375
  * Dahlburg et al. (1998) Dahlburg, R. B., Antiochos, S. K., & Klimchuk, J. A. 1998, ApJ, 495, 485
  * Degenhardt & Deinzer (1993) Degenhardt, U. & Deinzer, W. 1993, A&A, 278, 288
  * Dere et al. (1997) Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149
  * Fiedler & Hood (1992) Fiedler, R. A. S. & Hood, A. W. 1992, Sol. Phys., 141, 75
  * Field (1965) Field, G. B. 1965, ApJ, 142, 531
  * Heinzel et al. (2015) Heinzel, P., Gunár, S., & Anzer, U. 2015, A&A, 579, A16
  * Heinzel et al. (2014) Heinzel, P., Vial, J. C., & Anzer, U. 2014, A&A, 564, A132
  * Hildner (1974) Hildner, E. 1974, Sol. Phys., 35, 123
  * Hillier & van Ballegooijen (2013) Hillier, A. & van Ballegooijen, A. 2013, ApJ, 766, 126
  * Jenkins et al. (2019) Jenkins, J. M., Hopwood, M., Démoulin, P., et al. 2019, ApJ, 873, 49
  * Johnston & Bradshaw (2019) Johnston, C. D. & Bradshaw, S. J. 2019, ApJ, 873, L22
  * Karpen et al. (2001) Karpen, J. T., Antiochos, S. K., Hohensee, M., Klimchuk, J. A., & MacNeice, P. J. 2001, ApJ, 553, L85
  * Klimchuk & Cargill (2001) Klimchuk, J. A. & Cargill, P. J. 2001, ApJ, 553, 440
  * Klimchuk et al. (2010) Klimchuk, J. A., Karpen, J. T., & Antiochos, S. K. 2010, ApJ, 714, 1239
  * Koyama & Inutsuka (2004) Koyama, H. & Inutsuka, S.-i. 2004, ApJ, 602, L25
  * Landi et al. (2012) Landi, E., Del Zanna, G., Young, P. R., Dere, K. P., & Mason, H. E. 2012, ApJ, 744, 99
  * Lionello et al. (2009) Lionello, R., Linker, J. A., & Mikić, Z. 2009, ApJ, 690, 902
  * Luna & Karpen (2012) Luna, M. & Karpen, J. 2012, ApJ, 750, L1
  * Luna et al. (2018) Luna, M., Karpen, J., Ballester, J. L., et al. 2018, ApJS, 236, 35
  * Luna et al. (2012) Luna, M., Karpen, J. T., & DeVore, C. R. 2012, ApJ, 746, 30
  * Mikić et al. (2013) Mikić, Z., Lionello, R., Mok, Y., Linker, J. A., & Winebarger, A. R. 2013, ApJ, 773, 94
  * Milne et al. (1979) Milne, A. M., Priest, E. R., & Roberts, B. 1979, ApJ, 232, 304
  * Okamoto et al. (2007) Okamoto, T. J., Tsuneta, S., Berger, T. E., et al. 2007, Science, 318, 1577
  * Patsourakos & Vial (2002) Patsourakos, S. & Vial, J.-C. 2002, Sol. Phys., 208, 253
  * Rempel et al. (1999) Rempel, M., Schmitt, D., & Glatzel, W. 1999, A&A, 343, 615
  * Schmitt & Degenhardt (1995) Schmitt, D. & Degenhardt, U. 1995, Reviews in Modern Astronomy, 8, 61
  * Sharma et al. (2010) Sharma, P., Parrish, I. J., & Quataert, E. 2010, ApJ, 720, 652
  * Soler et al. (2011) Soler, R., Ballester, J. L., & Goossens, M. 2011, ApJ, 731, 39
  * Soler et al. (2012) Soler, R., Ballester, J. L., & Parenti, S. 2012, A&A, 540, A7
  * Soler et al. (2010) Soler, R., Oliver, R., & Ballester, J. L. 2010, A&A, 512, A28
  * Tandberg-Hanssen (1995) Tandberg-Hanssen, E. 1995, The nature of solar prominences, Vol. 199
  * van der Linden & Goossens (1991) van der Linden, R. A. M. & Goossens, M. 1991, Sol. Phys., 134, 247
  * Zhang et al. (2012) Zhang, Q. M., Chen, P. F., Xia, C., & Keppens, R. 2012, A&A, 542, A52
  * Zhang et al. (2013) Zhang, Q. M., Chen, P. F., Xia, C., Keppens, R., & Ji, H. S. 2013, A&A, 554, A124
  * Zhou et al. (2014) Zhou, Y.-H., Chen, P.-F., Zhang, Q.-M., & Fang, C. 2014, Research in Astronomy and Astrophysics, 14, 581
  * Zhou et al. (2021) Zhou, Y.-H., Ruan, W.-Z., Xia, C., & Keppens, R. 2021, A&A, 648, A29