2. Charge neutrality in a gravitationally stratified stellar atmosphere

We consider a plane-parallel atmosphere composed of a plasma of electrons and protons stratified by the gravitational field. For each species, the hydrostatic equilibrium condition is expressed as

(2.1) \begin{equation} \frac {\text{d} P_{\alpha }}{\text{d}z} = n_{\alpha }F_{\alpha } +R_{\alpha }, \quad F_{\alpha }=\mathrm{sign}(e_{\alpha }) e E -g m_{\alpha } , \end{equation}

where $\alpha \in \{e,p\}$ denotes the species (electrons or protons), $m_{\alpha }$ is the mass of a particle of species $\alpha$ , $e$ is the elementary charge, $g$ is the gravitational acceleration, $z$ is the height, $n_{\alpha }$ is the number density, $P_{\alpha }$ is the pressure, $F_{\alpha }$ is the total force per unit volume, $E$ is the electric field and $R_{\alpha }$ is the friction force due to collisions. Using the ideal gas law $P_{\alpha } = n_{\alpha } k_B T$ , and assuming quasi-neutrality, $n_{e}(z)=n_{p}(z)=n(z) \, \forall z$ , and the equality of temperatures, $T_{e}(z)=T_{p}(z)=T(z) \, \forall z$ , we subtract the equilibrium equations for each species to isolate the electric field:

(2.2) \begin{equation} E= E_{\textrm{PR}} - \frac {1}{2ne} \sum _{\alpha \in \{e,p\}}\mathrm{sign}\left (e_{\alpha }\right )R_{\alpha } + \frac {1}{2ne}\sum _{\alpha \in \{e,p\}}\mathrm{sign}\left (e_{\alpha }\right )\frac {\text{d} \left (T_{\alpha } n\right )}{\text{d}z} , \end{equation}

where $E_{\textrm{PR}}$ is the Pannekoek–Rosseland field defined as

(2.3) \begin{equation} E_{\textrm{PR}} = -\frac {\left (m_p-m_e\right )}{2e} g . \end{equation}

Even if the electron and proton temperatures are equal, the system may not be in thermal equilibrium, and thus $R_{\alpha } \neq 0$ due to collisional effects. Therefore, the total electric field must include an additional ambipolar contribution to restore charge neutrality. Taking the divergence of the electric field yields

(2.4) \begin{equation} \frac {\text{d}E}{\text{d}z} = \frac {1}{2e}\frac {\text{d}}{\text{d}z}\left [\frac {1}{n}\sum _{\alpha \in \{e,p\}} \mathrm{sign}\left (e_{\alpha }\right ) \left (\frac {\text{d} (T_{\alpha } n)}{\text{d}z}-R_{\alpha }\right )\right ] , \end{equation}

indicating that a self-consistent, spatially varying electric field must arise from the plasma to sustain quasi-neutrality. This ambipolar electric field becomes essential when collisional terms are included, as the classical Pannekoek–Rosseland field alone is insufficient. Given that Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024a, Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casettib) utilise the HMF model to approximate electrostatic interactions, a key question arises: can the HMF model, which retains only the first Fourier mode, reproduce the ambipolar effect and maintain charge neutrality in the presence of gravity? If not, can the model be suitably modified to include such effects by extending the Fourier mode expansion?

Figure 1.

Schematics of the loop model. The coronal plasma in the loop is treated as collisionless and in thermal contact with a fully collisional chromosphere (modelled as a thermal boundary).

3. The two-component gravitationally bound plasma model

We consider geometrically confined plasma structures, specifically coronal loops, which are prevalent throughout the solar atmosphere (e.g. Aschwanden Reference Aschwanden2005). Each loop is modelled as a semicircular tube of length $2L$ and cross-sectional area $S$ (figure 1). The plasma inside the loop is treated as a one-dimensional, collisionless, two-species electrostatic system subject to solar gravity. The electrostatic interaction is modelled via a multi-mode Fourier expansion of the potential, truncated at a chosen mode $N_n$ , while enforcing central symmetry at the loop apex. The equations of motion for each particle $j$ are given by

(3.1) \begin{equation} m_{\alpha } \ddot {x}_{j,\alpha }=eE \left (x_{j,\alpha }\right )+g m_{\alpha }\sin {\left (\frac {\pi x_{j,\alpha }}{2L}\right )}, \end{equation}

where $x \in [-L,L]$ denotes the curvilinear coordinate along the loop. The self-consistent electric fieldFootnote ¹ is expressed as

(3.2) \begin{equation} E(x)=8\,\mathrm{sign}{\left (e_{\alpha }\right )}e \boldsymbol{\cdot }n_SN\sum _{n=1}^{N_n} \frac {Q_{n}}{n}\sin {\left (\frac {\pi n x}{L}\right )} , \end{equation}

where $2N$ is the total number of particles of the species $\alpha$ , $n_S$ is the surface number density and $Q_n$ are the charge imbalance parameters defined as

(3.3) \begin{equation} Q_n=\sum _{\alpha \in \{e,i\}}\mathrm{sign}{\left (e_{\alpha }\right )}q_{n,\alpha }, \end{equation}

with the stratification parameters $q_{n,\alpha }$ given by

(3.4) \begin{equation} q_{n,\alpha }=\frac {1}{N}\sum _{j=1}^{N} \cos {\left (\frac {\pi nx_{j,\alpha }}{L}\right )}. \end{equation}

Physically, $q_{n,\alpha } \approx -1$ corresponds to particles concentrated near the base of the loop, $q_{n,\alpha } \approx 0$ to a uniform distribution and $q_{n,\alpha } \approx 1$ to particles clustered near the loop top. A non-zero $Q_n$ implies a charge imbalance at spatial scale $L/n$ . The Fourier components $\rho _n$ of the charge density $\rho (x)$ are related to $Q_n$ via the Poisson equation:

(3.5) \begin{equation} \frac {\partial ^2 \phi }{\partial x^2}=-4\pi e\rho \left (x\right ) \quad \rightarrow\quad \rho _n = -e n_0 Q_n, \quad n_0=\frac {N n_S}{L} , \end{equation}

indicating that $Q_n$ directly determines the spatial structure of the charge density at the spatial scale $L/n$ .

The loop is assumed to be in ideal thermal contact with a lower boundary representing the fully collisional chromosphere. However, as discussed in § 1, the upper chromosphere and the base of the transition region are highly dynamic environments, characterised by intense, short-lived and small-scale brightenings (Dere et al. Reference Dere, Bartoe and Brueckner1989; Teriaca et al. Reference Teriaca, Banerjee, Falchi, Doyle and Madjarska2004; Peter et al. Reference Peter2014; Tiwari et al. Reference Tiwari2019; Berghmans et al. Reference Berghmans2021).

We model these brightenings as temperature pulses as illustrated in figure 2. When the boundary temperature is held constant, the loop plasma reaches thermal equilibrium (isothermal loop). However, the dynamic behaviour of the chromosphere introduces stochastic heating events, which are analysed in § 5.2.

Figure 2.

Time evolution of the temperature at the thermal boundary. During intervals of duration $\tau$ , the temperature is increased by an amount $\Delta T$ , while during the waiting times $t_w$ , it is at the baseline value $T_0$ .

3.1. The system of units

We now proceed to express the equations of motion in dimensionless form. We introduce characteristic units

(3.6) \begin{equation} \begin{gathered} v_0=\sqrt {\frac {k_B T_{0}}{m_e}}, \quad m_0=m_e, \quad L_0=\frac {L}{\pi } , \end{gathered} \end{equation}

and define the dimensionless variables

(3.7) \begin{equation} M=\frac {m_p}{m_e}, \quad C=\frac {8 e^2 L^2 n_0}{\pi k_B T_0}, \quad \tilde {g}=\frac {g L m_p}{\pi k_B T_{0}}. \end{equation}

The quantities $C$ and $\tilde {g}$ quantify the strength of the self-electrostatic energy and of the gravitational energy of the electrons in units of thermal energy. In these units, denoting with $\theta$ the dimensionless spatial coordinate, the equations of motion become

(3.8) \begin{equation} M_{\alpha } \ddot {\theta }_{j,\alpha }= \mathrm{sign}{\left (e_{\alpha }\right )} E \left (\theta _{j,\alpha }\right )+\tilde {F}\left (\theta _{j,\alpha }\right ) , \end{equation}

the expressions for the external forces and the electrostatic field are

(3.9) \begin{equation} \tilde {F} \left (\theta _{j,\alpha }\right )= \tilde {g}_{\alpha } \sin {\left (\frac {\theta _{j,\alpha }}{2}\right )}, \quad E(\theta )= C\sum _{n=1}^{N_n}\frac {Q_n}{n}\sin {\left ( n\theta \right )} \end{equation}

and

(3.10) \begin{equation} Q_n=\sum _{\alpha \in \{p,e\}} \mathrm{sign}{\left (e_{\alpha }\right )}q_{n,\alpha },\quad q_{n,\alpha }=\frac {1}{N}\sum _{j=1}^{N}\cos {\left (n\theta _{j,\alpha }\right )} . \end{equation}

In the aforementioned equations, $M_{\alpha }$ is equal to the mass ratio $M=m_p/m_e$ for protons and $1$ for electrons, while $\tilde {g}_{\alpha }$ is equal to $\tilde {g}$ for protons and $\tilde {g}/M$ for electrons. Henceforth, all quantities are reported in dimensionless units unless specified otherwise.

3.2. Vlasov dynamics

In the mean-field limit, the dynamics of the phase-space distribution functions is governed by a system of two Vlasov equations:

(3.11) \begin{equation} \frac {\partial f_\alpha }{\partial t}+\frac {p}{M_{\alpha }}\frac {\partial f_{\alpha }}{\partial \theta }+F_{\alpha }[f_{\alpha }]\frac {\partial f_{\alpha }}{\partial p}=0, \quad F_\alpha =-\frac {\partial H_{\alpha }}{\partial \theta } , \end{equation}

where $f_\alpha$ are the distribution functions for both species and $H_{\alpha }$ are the mean-field Hamiltonians:

(3.12) \begin{equation} H_\alpha =\frac {p^2}{2M_\alpha }+ V_{\alpha }, \qquad V_{\alpha } \left ( \theta \right )=\mathrm{sign}{\left (e_\alpha \right )} \phi \left (\theta \right )+2\tilde {g}_{\alpha }\cos {\left (\frac {\theta }{2}\right )} . \end{equation}

Here $\phi$ is given by

(3.13) \begin{equation} \phi \left (\theta \right )= C \sum _{n=1}^{+\infty } \frac {Q_n}{n^2} \left (\cos {\left (n\theta \right )}+\left (-1\right )^{n+1}\right ) , \end{equation}

and $Q_n$ are given by

(3.14) \begin{equation} Q_n[f_{\alpha }]=\sum _{\alpha \in \{e,p\}} \mathrm{sign}{\left (e_{\alpha }\right )} q_{n,\alpha }[f_{\alpha }], \quad q_{n,\alpha }[f_{\alpha }]=\int _{-\pi }^{\pi } \text{d}\theta \int _{-\infty }^{\infty }\text{d}p \cos {\left (n\theta \right )} f_{\alpha }\left (\theta ,p\right ) . \end{equation}

4. Reproducing the Pannekoek–Rosseland field and charge neutrality with the multi-mode model

In this section, we demonstrate that the multi-mode electrostatic interaction model is capable of reproducing both the Pannekoek–Rosseland potential and the charge neutrality condition in the limit of a large number of Fourier modes. According to Jeans’ theorem, all stationary solutions of the Vlasov equation must depend solely on the mean-field Hamiltonians. Consequently, the distribution functions for each species can be expressed as

(4.1) \begin{equation} f_{\alpha } \left (\theta ,p \right )=f_{\alpha }\left (H_{\alpha }\right ) . \end{equation}

By combining (4.1) with (3.14) we obtain

(4.2) \begin{equation} Q_n=\sum _{\alpha \in \{e,p\}} \mathrm{sign}{\left (e_{\alpha }\right )} \int _{-\infty }^{+\infty } \text{d}p \int _{-\pi }^{+\pi } \text{d}\theta \cos {\left (n \theta \right )} f_{\alpha }\left (H_{\alpha }\right ) . \end{equation}

Using the kinetic definition of number density, (4.2) can be rewritten as

(4.3) \begin{equation} Q_n=\sum _{\alpha \in \{e,p\}} \mathrm{sign}{\left (e_{\alpha }\right )} \int _{-\pi }^{+\pi } \text{d}\theta \cos {\left (n\theta \right )} n_{\alpha }\left (V_{\alpha }\right ) , \end{equation}

where $V_{\alpha }$ is the total potential for each species. Since $V_{\alpha }$ depends on all $Q_n$ , this results in a system of coupled nonlinear equations. In general, this system cannot be solved analytically for the number densities $n_{\alpha }$ , but an analytical solution can be obtained under broad assumptions. First, we assume that the functional form of the number density is identical for both species:

(4.4) \begin{equation} n_{\alpha }\left (V_{\alpha }\right )=n \left (V_{\alpha }\right ) . \end{equation}

This assumption is satisfied by several well-known distribution functions, including thermal distributions, kappa distributions (Scudder Reference Scudder1992a, Reference Scudderb) and superstatistics (Beck & Cohen Reference Beck and Cohen2003), all widely used in theoretical modelling and observations of stellar atmospheres (Dudík et al. Reference Dudík, Polito, Dzifčáková, Zanna and Testa2017; Lazar & Fichtner Reference Lazar and Fichtner2021; Barbieri et al. Reference Barbieri, Casetti, Verdini and Landi2024a).

We now parametrise the Pannekoek–Rosseland potential along the loop’s curvilinear coordinate $\theta$ as

(4.5) \begin{equation} \phi _{\textrm{PR}}\left (\theta \right )= \sum _{\alpha \in \{e,p\}} \mathrm{sign}\left (e_{\alpha }\right ) \frac {\tilde {g}_{\alpha }}{2} \mathrm{cos}\left (\frac {\theta }{2}\right ) , \end{equation}

which can be expanded in a Fourier series as

(4.6) \begin{equation} \phi _{\textrm{PR}}\left (\theta \right )=\sum _{k=0}^{+\infty }c_k \cos {\left (k\theta \right )}, \quad c_k= \sum _{\alpha \in \{e,p\}} \mathrm{sign}\left (e_{\alpha }\right ) \tilde {g}_{\alpha }\frac {\left (-1\right )^k}{\left (1-4k^2\right )} . \end{equation}

The Pannekoek–Rosseland potential can be reproduced by the electrostatic interactions if the following condition holds:

(4.7) \begin{equation} C\sum _{n=1}^{N_n}\frac {Q_n}{n^2}\left (\cos {\left (n\theta \right )}+\left (-1\right )^{n+1}\right )=\sum _{k=0}^{+\infty }c_k \cos \left (k\theta \right ) . \end{equation}

This can be rewritten as

(4.8) \begin{equation} C\left (\sum _{n=1}^{N_n}\frac {Q_n}{n^2}\cos {\left (n\theta \right )}+Q_0\right ) = \sum _{k=0}^{+\infty }c_k \cos \left (\frac { \pi k x}{L}\right ), \qquad Q_0=\frac {\left (-1\right )^n Q_n}{n^2} . \end{equation}

This condition is satisfied only if the number of Fourier modes $N_n \rightarrow +\infty$ , and the $Q_n$ coefficients satisfy

(4.9) \begin{equation} Q_n = C_q \frac {n^2 \left (-1\right )^n}{\left (1-4n^2\right )},\quad C_q= \frac {1}{C} \sum _{\alpha \in \{e,p\}} \mathrm{sign}\left (e_{\alpha }\right ) \tilde {g}_{\alpha } . \end{equation}

This ensures that the electrostatic potential matches the Pannekoek–Rosseland potential exactly. Furthermore, if (4.4) is satisfied, then the right-hand side of (4.3) vanishes, implying $Q_n=0\ \forall n$ . However, (4.9) shows that $Q_n \neq 0$ , unless $C_q$ is negligible.

For typical coronal plasma conditions (e.g. $2L \sim 60 \,\mathrm{Mm}$ , $n_0 \sim 10^{9}\, \mathrm{cm}^{-3}$ , $T_0 = 10^6 \,\text{K}$ ), we find $C_q \sim 10^{-20}$ , which justifies the approximation $Q_n \sim 0$ , and hence the quasi-neutrality condition is satisfied. The coefficients $Q_n$ are related to the Fourier modes of the charge density $\rho (x)$ through (3.5). Since $|Q_n| \rightarrow C_q/4$ for large $n$ , the Fourier series of $\rho (x)$ diverges unless it is truncated at a finite mode $\hat {n}$ , satisfying

(4.10) \begin{equation} \hat {n} Q_{\hat {n}} \sim 1 \quad \mathrm{and} \quad \hat {n} \gg 1 , \end{equation}

leading to

(4.11) \begin{equation} \hat {n} \sim 4/C_{q} \quad \forall n\lt \hat {n} . \end{equation}

Now, since $C_q \sim 10^{-20}$ , for (4.11) we have that $\hat {n} \sim 4 \times 10^{20}$ . This implies that the solution of the system of (4.9) is given by (4.3) up to $n\lt \hat {n}$ . Now, since $Q_n \sim 0$ for all $n\lt \hat {n}$ , and given that the charge density components $\rho _n$ are proportional to $Q_n$ via (3.5), it follows that $\rho _n \sim 0$ for all $n\lt \hat {n}$ as well.

In conclusion, we have shown that electrostatic interactions can reproduce the Pannekoek–Rosseland potential and enforce quasi-neutrality down to a very small spatial scale proportional to $L/\hat {n}$ , in the realistic limit where electrostatic forces dominate over gravity, i.e. $C_q \ll 1$ . This behaviour appears to be general and not limited to loop geometry, as further demonstrated in Appendix A for a plane-parallel atmosphere. Now a natural question arises: where does this electric field come from?

To investigate where any residual charge imbalance occurs, we consider the case $C_q \sim 400$ , implying $\hat {n}=100$ according to (4.11). The electrostatic potential is then approximated by

(4.12) \begin{equation} \phi \left (\theta \right ) \cong \sum _{k=0}^{\hat {n}} c_{k}\cos {\left ( k\theta \right )} , \end{equation}

and the total charge density becomes

(4.13) \begin{equation} \rho \left (\theta \right )\cong -\sum _{n=1}^{\hat {n}}Q_n \cos {\left (n\theta \right )} . \end{equation}

Figure 3 illustrates that the self-consistent electrostatic potential (blue) closely approximates the Pannekoek–Rosseland potential (grey) throughout the loop. The middle panel shows that the discrepancies are more pronounced near the base of the loop. These differences are associated with a localised charge imbalance, as depicted in the right-hand panel. This analysis indicates that the Pannekoek–Rosseland potential primarily arises from a thin charge layer at the base of the loop, while the upper regions remain approximately quasi-neutral. This outcome is anticipated, as the gravitational force attains its maximum at the base, whereas the electrostatic field vanishes at that point.

Figure 3.

Left: the self-consistent electrostatic potential $\phi$ (blue) and the Pannekoek–Rosseland potential $\phi _{\textrm{PR}}$ (grey), plotted as functions of the curvilinear coordinate $\theta$ , both normalised by the factor $\sum _{\alpha \in \{e,p\}} \mathrm{sign} (e_{\alpha }) \tilde {g}_{\alpha }$ . The resulting rescaled potentials are denoted as $\tilde {\phi }$ and $\tilde {\phi }_{\textrm{PR}}$ , respectively. The self-consistent potential is computed from (4.12), while the Pannekoek–Rosseland potential is obtained from (4.5). Centre: the absolute difference $|\tilde {\phi }-\tilde {\phi }_{\textrm{PR}}|$ as a function of $\theta$ , highlighting the spatial deviation between the two potentials. Right: the charge density $\rho (\theta )$ plotted as a function of $\theta$ , calculated using (4.13).

Figure 4.

Top row, left-hand panel: electron (red) and proton (blue) number densities as functions of the curvilinear coordinate along the loop. The densities are computed using (5.2) with a single Fourier mode ( $N_n=1$ ), as indicated in the subplot title. The grey curve shows the reference density profile corresponding to the Pannekoek–Rosseland potential. Top row, right-hand panel: the Pannekoek–Rosseland potential (grey), calculated using (2.3), and the self-consistent electrostatic potential (blue), computed using (3.13) with one Fourier mode. Middle and bottom rows: same as the top row, but for $N_n=2$ and $N_n=9$ , respectively.

5. Validation of the multi-mode electrostatic model

In the previous section, we demonstrated that under the following assumptions:

1. (i) the number density satisfies $n_{\alpha } (V_{\alpha })=n (V_{\alpha })$ ,

2. (ii) the parameter $C_q \ll 1$ , that is, electrostatic interactions are much stronger than gravitational forces,

the multi-mode electrostatic model increasingly approximates charge neutrality and the Pannekoek–Rosseland potential as the number of Fourier modes increases. We now validate this conclusion for different types of distribution functions.

5.1. Thermal solution for multi-mode models

Among all the stationary solutions of the Vlasov equation, the thermal solution is given by

(5.1) \begin{equation} f_{\alpha }\left (H_{\alpha }\right )=\frac {1}{Z_{\alpha }}e^{- {H_{\alpha }}/{T}} \quad \mathrm{and} \quad Z_{\alpha }=\int _{-\infty }^{+\infty }\text{d}p \int _{-\pi }^{\pi }\text{d}\theta e^{-{H_{\alpha }}/{T}} , \end{equation}

where $Z_{\alpha }$ is the partition function.

Using the standard kinetic definitions, the number density for each species becomes

(5.2) \begin{equation} n_{\alpha }\left (\theta \right )= n \left (V_{\alpha }\right )=\frac {e^{- {V_{\alpha }(\theta )}/{T}}}{\int _{-\pi }^{\pi }\text{d}\theta e^{-{V_{\alpha } (\theta )}/{T}}} . \end{equation}

From this expression, it is evident that thermal distributions satisfy the density condition in (4.4). Substituting (5.2) into (4.3) and performing some algebra yields

(5.3) \begin{equation} Q_n= \sum _{\alpha \in \{e,p\}}\mathrm{sign}(e_{\alpha }) \frac {\int _{-\pi }^{0} \text{d}\theta \cos {\left (n\theta \right )}e^{-{V_{\alpha }(\theta )}/{T}}}{\int _{-\pi }^{0}\text{d}\theta e^{- {V_{\alpha } (\theta )}{/T}}}. \end{equation}

Solving this nonlinear system provides the values of $Q_n \ \forall n$ . To ensure the recovery of charge neutrality and the Pannekoek–Rosseland potential, the condition $C_q \ll 1$ must be satisfied.

Figure 4 presents the numerical results obtained with dimensionless parameters satisfying $C_q \ll 1$ , specifically: $\tilde {g}=32,\ M=m_p/m_e,\ C=10^4,\ T=90$ . The left-hand panels show the electron (red) and proton (blue) number density profiles computed via (5.2), where the coefficients $Q_n$ are derived by solving (5.3) using $N_n=1$ (top row), $N_n=2$ (middle row) and $N_n=9$ (bottom row) Fourier modes. The grey curves represent the density profiles computed using the Pannekoek–Rosseland potential from (4.5). As the number of Fourier modes increases, the species densities converge towards the Pannekoek–Rosseland profile, indicating improved charge neutrality. The right-hand panels of figure 4 show the corresponding electrostatic potential computed from (3.13) (blue), alongside the Pannekoek–Rosseland potential (grey). The convergence of the self-consistent potential towards the Pannekoek–Rosseland solution with increasing $N_n$ confirms the analytical predictions of § 4.

5.1.1. Validation via numerical simulations

We now demonstrate that the dynamic model introduced in § 3 relaxes towards the thermal equilibrium configuration described above, provided that the thermal boundary is held at a constant temperature $T$ . We initialise the system with isothermal distributions for both electrons and protons, accounting for their respective gravitational stratifications:

(5.4) \begin{equation} f_{\alpha }(\theta ,p)= \frac {e^{- ({1}/{T})( {p^2}/{(2M_{\alpha })}+2 \tilde {g}_{\alpha } \cos {(\theta /2 )} )}}{\sqrt {2\pi M_{\alpha }T}\int _{-\pi }^{\pi }\text{d}\theta e^{- ({2 \tilde {g}_{\alpha }}/{T}) \cos {({\theta }/{2} )}}} . \end{equation}

This configuration is not dynamically stable, as it yields distinct stratifications for the two species. As a result, an ambipolar electric field emerges and progressively acts to equalise the species densities, enforcing charge neutrality. This effect can only be reproduced if a sufficient number of Fourier modes is included in the electrostatic model.

We integrate numerically (3.8)Footnote ² using $N_n=9$ modes and the same parameter values as in figure 4, for $N=2^{17}$ particles. The top-left panel of figure 5 shows the time evolution of the kinetic energies of electrons (orange) and protons (green), defined as

(5.5) \begin{equation} K_{\alpha } = \frac {1}{N}\sum _{j=1}^{N} \frac {p_{j,\alpha }^2}{2M_{\alpha }} . \end{equation}

Both kinetic energies relax to an asymptotic value, consistent with thermal equilibrium. In this state, the electron and proton densities coincide and match the Pannekoek–Rosseland density profile (grey), shown in the top-right panel. The bottom-left panel illustrates the time evolution of the total electrostatic energy:

(5.6) \begin{equation} E_{el}= \frac {C}{2} \sum _{n=1}^{N_n} \frac {Q_n^2}{n^2} . \end{equation}

This quantity also relaxes to an asymptotic value. The electrostatic potential corresponding to this asymptotic state (blue), shown in the bottom-right panel, closely matches the Pannekoek–Rosseland potential (grey), confirming the model’s ability to self-consistently reproduce the equilibrium solution.

Figure 5.

Top left: time evolution of the kinetic energies $K_{\alpha }$ of protons (green curve) and electrons (orange curve). Top right: number densities of electrons (blue) and protons (red). The grey curves correspond to the theoretical density profile obtained from (5.2), where the electrostatic potential $\phi$ is replaced by the Pannekoek potential given in (4.5). Bottom left: time evolution of the total electrostatic energy $E_{el}$ , evaluated numerically via (5.6). Bottom right: self-consistent electrostatic potential $\phi$ , computed from simulations using (3.13). The grey curve represents the Pannekoek potential given by (4.5).

5.2. Multi-temperature solution for multi-mode models

The thermal solution presented above is valid under the assumption of a static thermal boundary at fixed temperature $T$ . However, as discussed in § 3 as well as in § 1, our primary interest lies in the case where the thermal boundary exhibits time-dependent temperature fluctuations. As shown in Barbieri et al. (Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b), if the time scales of these fluctuations are shorter than the relaxation time of the loop at each corresponding temperature, then the system relaxes towards a non-equilibrium stationary state. In this regime, the distribution functions of the plasma particles are described by the following analytical expression:

(5.7) \begin{equation} f_{\alpha }\left (\theta ,p\right )=\mathcal{N}_{\alpha }\left (A\int _1^{+\infty }\text{d}T\frac {\gamma \left (T\right )}{T}e^{- {H_{\alpha }}/{T}}+\left (1-A\right )e^{-H_{\alpha }}\right ), \quad A=\frac {\tau }{\tau + \langle t_w \rangle _{\eta }} , \end{equation}

where $\mathcal{N}_{\alpha }$ is the normalisation constant ensuring that $f_{\alpha }$ integrates to unity, $\gamma (T)$ is the probability distribution of the temperature increments, $\tau$ is the duration of the heating events and $\langle t_w \rangle$ is the average waiting time between events over the probability distribution $\eta (t_w)$ . Using the standard kinetic definitions, we compute the number density and kinetic temperature for each species:

(5.8) \begin{equation} n_{\alpha }\left (\theta \right )=n \left (V_{\alpha } \right )=\frac {A\int _{1}^{+\infty }\text{d}T{\gamma (T )}/{\sqrt {T}}e^{- {V_{\alpha }}/{T}}+\left (1-A\right )e^{-V_{\alpha }}}{A\int _{1}^{+\infty }\text{d}T{\gamma (T )}/{\sqrt {T}}\int _{-\pi }^{\pi } \text{d}\theta e^{- {V_{\alpha }}/{T}}+\left (1-A\right )\int _{-\pi }^{\pi } \text{d}\theta e^{-V_{\alpha }}} , \end{equation}

(5.9) \begin{equation} T_{\alpha }\left (\theta \right )=T \left (V_{\alpha } \right )=\frac {A\int _{1}^{+\infty }\text{d}T{\gamma (T )}/{\sqrt {T}}e^{- {V_{\alpha }}/{T}}+\left (1-A\right )e^{-V_{\alpha }}}{A\int _{1}^{+\infty }\text{d}T{\gamma (T )}/{\sqrt {T}}e^{-{V_{\alpha }}/{T}}+\left (1-A\right )e^{-V_{\alpha }}} . \end{equation}

These multi-temperature distributions also satisfy the condition given by (4.4), required for the system to approach charge neutrality. By inserting (5.8) into the system of (4.3), and after some manipulation, we obtain

(5.10) \begin{align} Q_n&=\sum _{\alpha \in \{e,p\}} \mathrm{sign}\left (e_{\alpha } \right )\nonumber \\&\quad \times \frac {A\int _{1}^{+\infty }\text{d}T({\gamma \left (T \right )}/{\sqrt {T}})\int _{-\pi }^{0} \text{d}\theta \cos {\left (n\theta \right )}e^{- {V_{\alpha }}/{T}}+\left (1-A\right )\int _{-\pi }^{0} \text{d}\theta \cos {\left (n\theta \right )}e^{-V_{\alpha }}}{A\int _{1}^{+\infty }\text{d}T{\gamma (T )}/{\sqrt {T}}\int _{-\pi }^{0} \text{d}\theta e^{- {V_{\alpha }}/{T}}+\left (1-A\right )\int _{-\pi }^{0} \text{d}\theta e^{-V_{\alpha }}} . \end{align}

Solving this system yields the values of $Q_n$ for all modes. As with the thermal case, we demonstrate this behaviour using dimensionless parameters that satisfy the condition $C_q \ll 1$ . The results are shown in figure 6, where we assume a delta-distribution for the temperature increments,

(5.11) \begin{equation} \gamma \left (T\right ) = \delta \left (T-\left (1+\Delta T\right )\right ) , \end{equation}

corresponding to a boundary temperature that alternates between two discrete values. The chosen parameters are: $A=0.1,\ \tilde {g}=2,\ M=m_p/m_e,\ C=10^4,\ \Delta T=90$ . In figure 6, the left-hand panels display the temperature and number density profiles of electrons (blue) and protons (red dashed), computed from (5.8) and (5.9), respectively, using $N_n=1,\ 2,\ 9$ . The grey lines represent the profiles obtained using the Pannekoek–Rosseland potential. As in the thermal case, increasing the number of Fourier modes leads to convergence of the electron and proton profiles towards the Pannekoek–Rosseland solutions, indicating that charge neutrality and the correct electrostatic potential structure are recovered. The right-hand panels show the corresponding electrostatic potentials computed using (3.13) (blue) compared with the Pannekoek–Rosseland potential (grey). Once again, convergence is observed as the number of modes increases. In Appendix B we show a case in which $\gamma (T)$ is not a simple delta-distribution proving the robustness of the present result.

Figure 6.

Top row, left-hand panel: electron (red) and proton (blue) number densities (globally decreasing functions) and temperatures (globally increasing functions) as functions of the curvilinear coordinate along the loop. The densities and temperature are computed using (5.2) with a single Fourier mode ( $N_n=1$ ), as indicated in the subplot title. The grey curve shows the reference temperature and density profiles corresponding to the Pannekoek–Rosseland potential. Top row, right-hand panel: the Pannekoek–Rosseland potential (grey), calculated using (2.3), and the self-consistent electrostatic potential (blue), computed using (3.13) with one Fourier mode. Middle and bottom rows: same as the top row, but for $N_n=2$ and $N_n=9$ , respectively.

5.2.1. Validation against numerical simulations

The dynamics governed by the equations of motion (3.1), coupled with the fluctuating thermal boundary, has been numerically simulated following the methodology described in Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024a, Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casettib), with the notable difference that the electric field is computed using (3.2). In this section, we illustrate how the fluctuating thermal boundary drives the overlying collisionless plasma towards a stationary state characterised by a temperature inversion. To this end, the system is initialised in a thermal equilibrium configuration, in which the distribution functions for both species are given by

(5.12) \begin{equation} f_{\alpha }\left (\theta ,p\right )= \frac {e^{- {p^2}/{(2M_{\alpha })}-2 \tilde {g} \cos {({\theta }/{2} )}}}{\sqrt {2\pi M_{\alpha }}\int _{-\pi }^{\pi }\text{d}\theta e^{-2 \tilde {g} \cos { ({\theta }/{2} )}}} \quad \forall \alpha . \end{equation}

This configuration is not stationary under the influence of the dynamically fluctuating thermal boundary at the base. However, as discussed in § 3, if a sufficiently large number of Fourier modes is retained in the modelling of the self-consistent electrostatic interactions, we expect the system to recover both charge neutrality and the Pannekoek–Rosseland potential in the stationary state. We report here the results of a simulation conducted with $N_n=9$ modes and using the same set of numerical parameters employed for the production of figure 6. The total number of particles is set to $N=2^{17}$ . In analogy with figure figure 5, the upper-left panel of figure 7 shows the time evolution of the total kinetic energy of each species $\alpha$ , computed according to (5.5). In this case as well, the kinetic energies of both species relax towards an asymptotic value. The upper-right panel displays the final stationary profiles of temperature and density computed via the numerical simulation together with their theoretical counterparts, which clearly exhibit a temperature inversion. The theoretical profiles, shown in grey, are computed using (5.8) and (5.9), where the electrostatic potential $\phi (\theta )$ is replaced by the Pannekoek–Rosseland expression given in (4.5). The lower-left panel illustrates the time evolution of the total electrostatic energy, evaluated using (5.6). As observed in the case discussed in § 5.1, this quantity also relaxes to an asymptotic value. Finally, the lower-right panel confirms that the self-consistent electrostatic potential generated by the plasma coincides with the Pannekoek–Rosseland potential at equilibrium.

Figure 7.

Same quantities and colour scheme as in figure 5, but in the multi-temperature case regime described in § 5.2.

6. Summary and perspectives

The main plasma physics result of the present work is the identification of a self-consistent mechanism to reproduce the ambipolar effect that gives rise to the Pannekoek–Rosseland electric field and ensures charge neutrality in a gravitationally stratified, collisionless plasma. This is achieved by extending the electrostatic interaction model from a single Fourier mode (as in the standard HMF approximation) to multiple modes. Under very general assumptions on the form of the distribution functions, assumptions that are satisfied by thermal, kappa and superstatistical distributions, we demonstrated that in the limit of many modes, the self-consistent electrostatic potential converges to the classical Pannekoek–Rosseland potential. At the same time, the system exhibits quasi-neutrality apart from a thin layer of charge at the base of the atmosphere which originates the Pannekeok–Rosseland potential. These theoretical predictions were confirmed numerically. We have shown that, under suitable conditions, the distribution functions of both species relax towards a stationary configuration in which charge neutrality and the Pannekoek–Rosseland potential are simultaneously recovered. This was observed both in the thermal equilibrium case (where the boundary temperature is fixed) and in the more physically relevant multi-temperature case, where the thermal boundary undergoes random energy fluctuations.

While the results of this work do not aim to provide a general theory of ambipolar electric fields nor a complete description of coronal heating, they represent a foundational step in that direction. Now that we have established how to incorporate the ambipolar effect in a self-consistent way into the model originally introduced in Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024a, Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casettib), we are in a position to extend it. In particular, we can include additional ingredients such as collisions, ionisation processes and asymmetric temperature fluctuations at the loop footpoints. These effects will allow us to determine how the ambipolar electric field predicted by the multi-mode Fourier approach is modified, and how the picture of coronal heating evolves in this more complete framework.

As discussed in § 2, we expect collisions to modify the self-generated electric field with respect to the classical Pannekoek–Rosseland expression. For example, we aim to investigate the following questions: How many Fourier modes are required to reproduce the modified ambipolar field in the collisional case? In which regions of the atmosphere are the deviations from the Pannekoek–Rosseland field most significant? How does this modified electric field act on the temperature inversion produced by gravitational filtering in the collisionless limit? Finally, by including collisions, we can study the interplay between the collisionless gravitational filtering, which produces the temperature inversion effect, and the collisional thermal conduction, which allows heat to flow from the corona into the chromosphere.

Another important extension of the present model concerns the role of partial ionisation, especially in the lower part of the transition region, where the ionisation process is known to occur. In this regime, deviations from the classical Pannekoek–Rosseland field may become more significant, as the plasma is not yet fully ionised. As altitude increases and ionisation progresses, the assumptions of a fully ionised and collisionless plasma, underlying the present treatment, become more accurate. Understanding how the ambipolar electric field evolves across this ionisation gradient represents an important development of our model.

We note that in our model we have central symmetry with respect to the loop apex. This will ensure a zero net flux in the system. By introducing asymmetric temperature fluctuations at the two footpoints, a net flux across the loop appears.

However, the resulting asymmetry does not alter the physical mechanism responsible for the temperature inversion. The inversion still arises because ‘hot’ particles, generated by temperature increments, can climb higher in the gravitational potential well, while ‘cold’ particles remain concentrated near the base and the temperature profile along the loop is no longer symmetric. Although this adds complexity to the model, it does not change the physical origin of the temperature inversion. It represents an interesting extension that we plan to explore in a future work.

We note that generating different temperature increments at the two footpoints of the loop produces a flow into the system, altering the ambipolar electric field, which would no longer be the classical Pannekoek–Rosseland field. Another way to obtain a different ambipolar field could be through the introduction of distinct temperature fluctuations for the two species. This leads to different temperature profiles for electrons and protons, and consequently, as evident from (2.2), an ambipolar electric field that deviates from the classical Pannekoek–Rosseland expression.

Finally, as an additional possible follow-up, it would be interesting to include alpha particles in the modelling, as they contribute significantly to the mass density.