2. The two-component gravitationally bound plasma model

We model a gravitationally bound plasma, focusing on geometrically confined plasma structures, with specific reference to the coronal loops that are common in the Sun's atmosphere (see e.g. Aschwanden Reference Aschwanden2005). We model the loop as a semicircular tube of length $2L$ and section $S$, where the charge distribution is discretised onto $N_e$ sections with surface number density $n_S$, surface charge density $\sigma _e=-e n_S$ and surface mass density $\sigma _{m,e}=m_e n_s$; and $N_p$ sections with surface charge density $\sigma _p=e n_S$ and surface mass density $\sigma _{m,p}=m_p n_S$, representing the electron and proton components, respectively. Moreover, the charge of the electron is represented by the constant e, while the masses of an electron and a proton are denoted by $m_e$ and $m_p$, respectively. We assume that $N_e=N_p=2N$ in order to ensure global charge neutrality. A scheme of a two-component loop plasma model is sketched in figure 1. We now introduce the following assumptions:

1. (i) Particles are subject to an external force field pointing towards the loop's end, representing the Sun's (uniform) gravity plus the Pannekoek–Rosseland electrostatic field (Pannekoek (Reference Pannekoek1922) and Rosseland (Reference Rosseland1924), see also Belmont et al. (Reference Belmont, Grappin, Mottez, Pantellini and Pelletier2013) for an extended review).

2. (ii) The symmetry is cylindrical along the loop-following coordinate $x$.

3. (iii) Electrostatic self-interactions are modelled by truncating to the first-order Fourier expansion as in the so-called Hamiltonian mean-field (HMF) models (see Antoni & Ruffo Reference Antoni and Ruffo1995; Chavanis, Vatteville & Bouchet Reference Chavanis, Vatteville and Bouchet2005; Giachetti & Casetti Reference Giachetti and Casetti2019).

4. (iv) Every observable is symmetric with respect to the top of the loopFootnote ¹. As a consequence, if the top of the loop coincides with $x = 0$, all the functions of the canonical coordinates are centrally symmetric with respect to the origin of the two-dimensional single-particle phase spaceFootnote ².

The first two assumptions imply that the total energy of the system is

(2.1)\begin{equation} H=S\left(\sum_{j=1}^{2N}\sum_{\alpha \in \{e,p\}} \frac{1}{2}\sigma_{m,\alpha} \dot{x}^2_{j,\alpha} + E_{el} + U_{g}\right), \end{equation}

where $\alpha \in \{e,p\}$ denotes the species (here electrons or protons), $g={G M_{\odot }}/{R_{\odot }^2}$ is the gravitational field at the surface of the Sun, $M_{\odot}$ is the mass of the Sun, G is Cavendish gravitational constant, $R_{\odot}$ is the radius of the Sun and $x$ is the spatial coordinate (i.e. the curvilinear abscissa of the loop). The total electrostatic energy per section $E_{el}$ reads

(2.2)\begin{equation} E_{el}=\frac{1}{2}\int_{{-}L}^{L} {{\rm d}x} \rho(x)\phi(x), \end{equation}

where the charge density

(2.3)\begin{equation} \rho(x)= \sum_{j=1}^{2N}\sum_{\alpha \in \{e,p\}} \sigma_{\alpha} \delta (x-x_{j,\alpha}), \end{equation}

is related to the potential $\phi$ by the Poisson equation

(2.4)\begin{equation} \frac{\partial^2 \phi}{\partial x^2}={-}4{\rm \pi}\rho. \end{equation}

Here, $U_g$ is the total external potential per section and reads

(2.5)\begin{equation} U_g=\frac{gL\sum_{\alpha \in \{e,p\}} \sigma_{m,\alpha}}{\rm \pi} \sum_{j=1}^{2N}\sum_{\alpha \in \{e,p\}} \cos{\left(\frac{{\rm \pi} x_{j,\alpha}}{2L}\right)}. \end{equation}

In order to use the HMF modellisation we Fourier analyse the electrostatic energy $E_{el}$. Hereafter, we will use the following definition for the Fourier transform:

(2.6a,b)\begin{equation} f(x)=\sum_{n\in \mathbb{Z}} f_n \exp\left({{\rm i}\frac{{\rm \pi} n x}{L}}\right),\quad f_n=\frac{1}{2L}\int_{{-}L}^L {{\rm d}x} f(x) \exp\left({-{\rm i}\frac{{\rm \pi} n x}{L}}\right). \end{equation}

We then define the charge imbalances as

(2.7)\begin{equation} \boldsymbol{Q}_{n}=(Q_{n,x},Q_{n,y})=\sum_{\alpha \in \{e,p\}} \mathrm{sign}{(\sigma_{\alpha})}\boldsymbol{q}_{n,\alpha}, \end{equation}

where the vectors $\boldsymbol {q}_{n,\alpha }$ are given by

(2.8)\begin{equation} \boldsymbol{q}_{n,\alpha}=\frac{1}{2N}\sum_{j=1}^{2N} \left[\cos{\left(\frac{{\rm \pi} n x_{j,\alpha}}{L}\right)}, \sin{\left(\frac{{\rm \pi} n x_{j,\alpha}}{L}\right)}\right]. \end{equation}

Performing the Fourier transform (2.6a,b) of the Poisson equation (2.4) and applying (2.7)–(2.8) we obtain the modes of the potential–density pair as

(2.9)\begin{equation} \phi_n=\frac{4L^2}{{\rm \pi} n^2}\rho_n,\quad \rho_n={-}\frac{\sigma_e N}{L}(Q_{n,x}-i Q_{n,y})\quad \forall\ n \neq 0,\enspace \rho_0=0. \end{equation}

Since each density mode $\rho _n$ defined above is related to a specific $\boldsymbol {Q}_n$, any vanishing density mode $\rho _n$ corresponds to a zero charge imbalance. This can be also seen directly from (2.7) and (2.8). Indeed, if most of the $2N$ particles of a given species are symmetrically concentrated at the bottom of the loop, then $\boldsymbol {q}_{n,\alpha } \approx -1$. Indeed, if they are uniformly distributed then $\boldsymbol {q}_{n,\alpha } \approx 0$ ($\boldsymbol {q}_{n,\alpha } \approx +1$ corresponds to a situation where most of the particles are concentrated at the top of the loop). As a consequence, if their difference $\boldsymbol {Q}_n$ does not vanish, there is a charge imbalance in the system.

Figure 1. Sketch of the two-component plasma loop model. The vertical axis $z$ is the altitude in the atmosphere; $x$ is the curvilinear abscissa of the loop. Here, $\sigma _{m,\alpha }$,$\sigma _{\alpha }$ with $\alpha =\{e,p\}$ are respectively the surface mass density and the surface charge density of the species $\alpha$.

From now on we will refer to the quantities $\boldsymbol {q}_{n,\alpha }$ as the stratification parameters. By performing the Fourier expansion (2.6a,b) of the charge density (2.3) and of the electrostatic potential $\phi (x)$ and by using (2.7)–(2.9), after some algebra we get the decomposition in Fourier modes of the electrostatic energy per section (2.2) as

(2.10)\begin{equation} E_{el}=\frac{8L (N\sigma_e)^2}{\rm \pi}\sum_{n=1}^{+\infty} \frac{||\boldsymbol{Q}_n||^2}{n^2}. \end{equation}

The proportionality to $(N\sigma _e)^2$ arises from the fact that each particle of the system feels the electrostatic forces due to all other particles. Its proportionality to $L$ means that, the larger the size of the system is, the greater is $E_{el}$. The electrostatic energy $E_{el}$ depends on all the charge imbalances $||\boldsymbol {Q}_n||^2$ scaled by $n^2$, that is, large-scale modes make a larger contribution to the electrostatic energy compared with the small-scale ones. Using (2.10) combined with the total energy (2.1) we derive the equations of motion in the form

(2.11)\begin{equation} m_{\alpha} \ddot{x}_{j,\alpha}=eE(x_{j,\alpha})+g \frac{\displaystyle\sum_{\alpha \in \{e,p\}}m_{\alpha}}{2} \sin{\left(\frac{{\rm \pi} x_{j,\alpha}}{2L}\right)}, \end{equation}

where the electric field $E$ decomposed in Fourier modes is

(2.12)\begin{equation} E(x)=8\,\mathrm{sign}{(e_{\alpha})}e_{\alpha}n_SN\sum_{n=1}^{+\infty} \left[\frac{Q_{n,x}}{n}\sin{\left(\frac{{\rm \pi} n x}{L}\right)}- \frac{Q_{n,y}}{n}\cos{\left(\frac{{\rm \pi} n x}{L}\right)}\right]. \end{equation}

We now implement the HMF assumption by truncating the Fourier expansion at the first mode. The equations above become

(2.13)\begin{equation} m_{\alpha} \ddot{x}_{j,\alpha}=e_{\alpha}E(x_{j,\alpha})+ g \frac{\displaystyle\sum_{\alpha \in \{e,p\}}m_{\alpha}}{2}\sin{\left(\frac{{\rm \pi} x_{j,\alpha}}{2L}\right)}, \end{equation}

and the electric field becomes

(2.14)\begin{equation} E(x)=8\,\mathrm{sign}(e_{\alpha})e_{\alpha}n_SN \left[Q_{x}\sin{\left(\frac{{\rm \pi} x}{L}\right)}-Q_{y}\cos{\left(\frac{{\rm \pi} x}{L}\right)}\right]. \end{equation}

In the equations above, $\boldsymbol {q}_{\alpha }=\boldsymbol {q}_{1,\alpha }$ are the first mode (large-scale) stratification parameters and $Q_x=Q_{1,x}$, $Q_y=Q_{1,y}$ are the first mode charge imbalance components. From now on, unless explicitly stated, we will refer to $\boldsymbol {q}_{\alpha }$ and $\boldsymbol {Q}$, respectively, as the stratification parameters and the charge imbalance. As one can see, the HMF assumption reduces the computational cost of the electrostatic interactions from $N^2$ to $N$.

Let us now comment further on the physics behind the expression of the self-electrostatic interactions. To do so, we use again the electrostatic energy in terms of its Fourier modes (2.10) truncating the expression at the first mode. After some elementary algebra we get

(2.15)\begin{equation} E_{el}=E_0\sum_{j,k=1}^{2N}\sum_{\alpha,\beta \in \{e,p\}} \tilde{E}_{\alpha,\beta}\boldsymbol{q}_{j,\alpha} \boldsymbol{\cdot}\boldsymbol{q}_{k,\beta}, \end{equation}

where

(2.16)\begin{equation} \boldsymbol{q}_{j,\alpha}=\left[\cos{\left(\frac{{\rm \pi} x_{j,\alpha}}{L}\right)}, \sin{\left(\frac{{\rm \pi} x_{j,\alpha}}{L}\right)}\right], \end{equation}

where $\tilde {E}_{\alpha,\beta }=\mathrm {sign}(e_{\alpha })\,\mathrm {sign}(e_{\beta })$, and $E_0={4{\rm \pi} ^{-1} L e^2n_S^2N}$. From these expressions it is apparent that the interactions among particles of the same species are described by an antiferromagnetic HMF term while the interactions among particles of different species are, conversely, described by a ferromagnetic HMF term. The antiferromagnetic HMF tends to increase the angle between $\boldsymbol {q}_{j,\alpha }$ and $\boldsymbol {q}_{k,\beta }$ in order to minimise the energy. In terms of physical quantities in the loop, this interaction tends to increase the physical distance $x_{j,\alpha }-x_{k,\alpha }$ between two particles of the same species $\alpha$ up to the value $L$ and thus such a term mimics the electrostatic repulsion between charges of the same sign. The ferromagnetic HMF tends to minimise the energy by decreasing the distance $x_{j,\alpha }-x_{k,\beta }$ between two particles of different species $\alpha,\beta$ up to zero and this mimics the electrostatic attraction between charges of the different sign.

In the following, we consider this plasma model in thermal contact with a thermal boundary of temperature $T_0$ (the base of the plasma atmosphere).

2.1. Normalisation and central symmetry assumption

To study the dynamics defined by (2.13), we define the following units for velocity, mass and length:

(2.17a–c)\begin{equation} v_0=\sqrt{\frac{k_B T_{0}}{m_e}},\quad m_0=m_e,\quad L_0=\frac{L}{\rm \pi}. \end{equation}

The choice of units given above implies that the unit of energy is $E_0 = k_B T_0$, where $T_0$ is the unit of temperature and $k_B$ is the Boltzmann constant. For the specific case of the Sun the unit of temperature $T_0$ is set equal to the reference temperature of the thermal boundary (chromosphere), that is $T_0=10^4$ K. The equations of motion are expressed in dimensionless form as

(2.18)\begin{equation} M_{\alpha} \ddot{\theta}_{j,\alpha}= \mathrm{sign}{(e_{\alpha})} C E(\theta_{j,\alpha})+\tilde{F}(\theta_{j,\alpha}), \end{equation}

where the expressions for the external forces and the electrostatic field are

(2.19a,b)\begin{equation} \tilde{F}(\theta_{j,\alpha})= \tilde{g} \sin{\left(\frac{\theta_{j,\alpha}}{2}\right)},\quad E(\theta)=Q_{x}\sin{\theta}-Q_{y}\cos{\theta}, \end{equation}

and the charge imbalance and stratification parameters are given by

(2.20a,b)\begin{equation} \boldsymbol{Q}=(Q_x,Q_y)=\sum_{\alpha \in \{p,e\}} \mathrm{sign}{(e_{\alpha})}\boldsymbol{q}_{\alpha},\quad \boldsymbol{q}_{\alpha}=\frac{1}{2N}\sum_{j=1}^{2N}(\cos{(\theta_{j,\alpha})},\sin{(\theta_{j,\alpha})}). \end{equation}

In the equations above, $M_{\alpha }$ equals $M_p ={m_p}/{m_e}$ for the protons and $M_e = 1$ for the electrons, while $\theta$ is the dimensionless spatial coordinate. The dynamics is therefore fully identified by the three parameters

(2.21a–c)\begin{equation} M=\frac{m_p}{m_e},\quad C=\frac{8 e^2 L^2 n_0}{{\rm \pi} k_B T_0},\quad \tilde{g}=\frac{g L (m_p+m_e)}{2 {\rm \pi}k_B T_{0}}, \end{equation}

where $n_0=Nn_S/L$ is the average density of a given species. The quantities $C$ and $\tilde {g}$ measure the strength of the electrostatic interaction and of the external field in units of thermal energy, respectively. Hereafter, unless explicitly stated, we will use dimensionless quantities in all the equations and for all the quantities to be plotted in the figures.

Assuming that all quantities are symmetric with respect to the top of the loop, in terms of the particles’ phase-space coordinates, corresponds to imposing the following symmetry rules:

(2.22a,b)\begin{equation} \theta_{j+N,\alpha}={-}\theta_{j,\alpha},\quad \dot{\theta}_{j+N,\alpha}={-}\dot{\theta}_{j,\alpha};\quad \forall j=1 \cdots N, \end{equation}

where the first $j=1 \cdots N$ particles populate one half of the loop $\theta _{j,\alpha } \in [-{\rm \pi},0]$ and the remaining particles populate the other half. With such an assumption the system of (2.18) can be reduced to a system of equations for only the first half of the loop as

(2.23)\begin{equation} M_{\alpha} \ddot{\theta}_{j,\alpha}= \mathrm{sign}{(e_{\alpha})} C E(\theta_{j,\alpha})+\tilde{F}(\theta_{j,\alpha}), \end{equation}

where the (2.19a,b) become

(2.24a,b)\begin{align} \tilde{F}(\theta_{j,\alpha})= \tilde{g} \sin{\left(\frac{\theta_{j,\alpha}}{2}\right)},\quad E(\theta)=Q\sin{\theta}, \end{align}

and (2.19a,b) become

(2.25a,b)\begin{equation} Q=\sum_{\alpha \in \{p,e\}} \mathrm{sign}{(e_{\alpha})}q_{\alpha},\quad q_{\alpha}=\frac{1}{N}\sum_{j=1}^{N}\cos{(\theta_{j,\alpha})}. \end{equation}

The positions and velocities of particles in the other half of the loop are then determined using (2.19a,b).

2.2. Vlasov dynamics and thermal equilibrium solution

In the continuum limit, in terms of phase-space distribution functions the mean-field dynamics of the normalised plasma model given by (2.23) is governed by a system of two Vlasov equations

(2.26a,b)\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 single-particle distribution functions for both species, and $H_{\alpha }$ are the mean-field Hamiltonians

(2.27)\begin{equation} H_\alpha=\frac{p^2}{2M_\alpha}+\mathrm{sign}{(e_\alpha)}C \phi(\theta)+2\tilde{g} \cos{\left(\frac{\theta}{2}\right)};\quad \phi(\theta)=Q\cos{\theta}, \end{equation}

where

(2.28a,b)\begin{equation} Q=\sum_{\alpha \in \{e,p\}} \mathrm{sign}{(e_{\alpha})} q_{\alpha},\quad q_{\alpha}=\int_{-{\rm \pi}}^{\rm \pi} {\rm d}\theta \int_{-\infty}^{\infty}{\rm d} p \cos{\theta} f_{\alpha}(\theta,p). \end{equation}

In virtue of the Jeans theorem (see e.g. Nicholson Reference Nicholson1983), all stationary solutions of the system (2.24a,b) are functions the sole mean-field Hamiltonians (2.27). Technically speaking, in general, the stationary solution of a Vlasov equation is a function of all the constants of motion. Here, since our model is one-dimensional the only independent constants of motion are the mean-field Hamiltonians $H_\alpha$. Among all the solutions of the system above, the thermal equilibrium one has the following expression:

(2.29a,b)\begin{equation} f_{\alpha}(\theta,p)=\frac{\exp\left({-\dfrac{H_{\alpha}}{T}}\right)}{Z_{\alpha}},\quad Z_{\alpha}=\int_{-{\rm \pi}}^{\rm \pi}{\rm d} \theta \int_{-\infty}^{+\infty}{\rm d} p \exp\left({-\frac{H_{\alpha}}{T}}\right), \end{equation}

where $f_{\alpha }$ is normalised to 1 for both species. By combining the above equation with (2.28a,b), we obtain

(2.30)\begin{equation} Q=\sum_{\alpha \in \{e,p\}} \mathrm{sign}{(e_{\alpha})} \int_{-\infty}^{+\infty} {\rm d} p \int_{-{\rm \pi}}^{\rm \pi} {\rm d}\theta \cos{\theta} \frac{\exp\left({-\dfrac{H_{\alpha}[Q]}{T}}\right)}{ \displaystyle\int_{-\infty}^{+\infty} {\rm d} p\int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta \exp\left({-\dfrac{H_{\alpha}[Q]}{T}}\right)}, \end{equation}

and integration with respect to $p$ gives

(2.31)\begin{align} Q & = \int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta \cos{\theta} \frac{\exp\left({-\dfrac{CQ\cos{\theta}+\tilde{g} \cos{\left(\dfrac{\theta}{2}\right)}}{T}}\right)}{\displaystyle\int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta \exp\left({-\dfrac{CQ\cos{\theta}+\tilde{g}\cos{\left(\dfrac{\theta}{2}\right)}}{T}}\right)} \nonumber\\ & \quad -\int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta\cos{\theta}\frac{\exp\left({-\dfrac{-CQ\cos{\theta}+\tilde{g} \cos{\left(\dfrac{\theta}{2}\right)}}{T}}\right)}{\displaystyle\int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta \exp\left({-\dfrac{-CQ\cos{\theta}+\tilde{g}\cos{\left(\dfrac{\theta}{2}\right)}}{T}}\right)} , \end{align}

whose solution is $Q=0$ (and so $\phi =0$). Therefore, the equilibrium solution of a two-component plasma model, where the two species are subject to the same external potential, requires a vanishing self-electrostatic potential, (i.e. $\phi =0$). The thermal solution for the plasma atmosphere model is the so-called isothermal atmosphere given by

(2.32a,b)\begin{equation} f_{\alpha}(\theta,p)=\frac{\exp\left({-\dfrac{\tilde{H}_{\alpha}}{T}}\right)}{Z_{\alpha}},\quad Z_{\alpha}=\int_{-{\rm \pi}}^{\rm \pi} {\rm d}\theta \int_{-\infty}^{+\infty} {\rm d} p\exp\left({-\dfrac{\tilde{H}_{\alpha}}{T}}\right), \end{equation}

where $\tilde {H}_{\alpha }$ are the mean-field Hamiltonians (2.27) with $\phi =0$, that is

(2.33)\begin{equation} \tilde{H}_\alpha=\frac{p^2}{2M_\alpha}+2\tilde{g}\cos{\left(\frac{\theta}{2}\right)}. \end{equation}

The knowledge of the distribution functions allows us to compute the number densities $n_\alpha$ and the temperatures profiles $T_\alpha$ using the standard kinetic definitions (see e.g. Nicholson Reference Nicholson1983). For the number densities we get

(2.34)\begin{equation} n_{\alpha}(\theta)=\frac{\exp\left({-\dfrac{2\tilde{g}}{T} \cos{\left(\dfrac{\theta}{2}\right)}}\right)}{\displaystyle\int_{-{\rm \pi}}^{\rm \pi}{\rm d} \theta \exp\left({-\dfrac{2\tilde{g}}{T}\cos{\left(\dfrac{\theta}{2}\right)}}\right)}, \end{equation}

and for the kinetic temperatures

(2.35)\begin{equation} T_{\alpha}(\theta)=T. \end{equation}

3. Temporal coarse graining

Starting from a thermal equilibrium state (2.32a,b) with temperature $T_0=1$, the out-of-equilibrium state is induced by introducing fluctuations in the value of the temperature of the thermal boundary. In practice, we first increase the temperature from the initial equilibrium value $T_0$ up to $T>T_0$ (sorted from a probability distribution $\gamma (T)$) for a finite fixed time interval $\tau$Footnote ³ , of length smaller than the relaxation time $t_R$ to equilibrium at temperature $T$. After a time $\tau$ the temperature of the thermal boundary is reverted to $T_0$ for a finite time interval $t_w$ sorted from a probability distribution $\eta (t_w)$, again much smaller than the relaxation time $t_R^\prime$ back to $T_0$. Iterating this procedure prevents the system from relaxing to a thermal equilibrium at either $T_0$ or $T$. Under the conditions prescribed above we now define the coarse-graining time scale $\tilde {t}$ such that

(3.1)\begin{equation} \tau,\langle t_w \rangle \ll \tilde{t} \ll \tilde{t}_{R}, \end{equation}

where $\tilde {t}_R\equiv \min (t_R;t_R^\prime )$. The value of $\tilde {t}_R$ is evaluated as the minimum between the thermal crossing time and the gravity crossing time of electrons

(3.2)\begin{equation} \tilde{t}_R\equiv \min(t_T;t_{\tilde{g}})\quad t_T=\frac{2{\rm \pi}}{\sqrt{T}} \quad t_{\tilde{g}}=\sqrt{\frac{\rm \pi}{\tilde{g}}}. \end{equation}

The temporal coarse-grained phase-space distribution function, defined as the time average of $f_{\alpha }$ over $\tilde {t}$, reads as

(3.3)\begin{equation} \tilde{f}_{\alpha}(\theta,p,\tilde{t})=\langle f_{\alpha} \rangle_{\tilde{t}}= \frac{1}{\tilde{t}}\int_{\tilde{t}} {\rm d} t f_{\alpha}(\theta,p,t), \end{equation}

and the time-averaged Vlasov dynamics over $\tilde {t}$ becomes

(3.4)\begin{equation} \frac{f_{\alpha}(\theta,p,\tilde{t})-f_{\alpha}(\theta,p,0)}{\tilde{t}}+ \frac{p}{M_{\alpha}}\frac{\partial \tilde{f}_{\alpha}}{\partial \theta}+ \left\langle \tilde{F}_{\alpha}\frac{\partial f_{\alpha}}{\partial p}\right\rangle_{\tilde{t}}=0. \end{equation}

Given the condition $\tilde {t} \ll \tilde {t}_{R}$, during the time interval $\tilde {t}$ the system energy cannot be redistributed along the loop. We can therefore approximate $f_{\alpha }$ with its time average $f_{\alpha }=\tilde {f}_{\alpha }$, while the finite difference on the left-hand side of (3.4) becomes a time derivative. We then get

(3.5)\begin{equation} \frac{\partial \tilde{f}_{\alpha}}{\partial \tilde{t}}+\frac{p}{M_\alpha} \frac{\partial \tilde{f}_{\alpha}}{\partial \theta}+\tilde{F}_{\alpha}[\tilde{f}_{\alpha}] \frac{\partial \tilde{f}_{\alpha}}{\partial p} =0. \end{equation}

The thermal boundary at the bottom boundary of the system induces an incoming flux defined at a given time by

(3.6)\begin{equation} J_{T,\alpha}=\int_0^{+\infty}{\rm d} p\frac{p}{TM_{\alpha}}\exp\left({-\frac{p^2}{2TM_{\alpha}}}\right), \end{equation}

where the value of $T$ is drawn from the distribution $\gamma (T)$ during the time intervals of length $\tau$; while $T = 1$ during time intervals of length $t_w$ drawn from $\eta (t_w)$.

We now compute the temporal coarse-grained flux by taking the time average of (3.6) obtaining

(3.7)\begin{equation} J_{\alpha} = \frac{t_{t_w}}{\tilde{t}}\langle J_{T,\alpha} \rangle_{t_{t_w}} + \frac{t_{\tau}}{\tilde{t}}\langle J_{1,\alpha} \rangle_{t_{\tau}}, \end{equation}

where $\tilde {t}=t_{t_w}+t_{\tau }\ t_{t_w}=\sum _{i=1}^{n_p} t_{w_i}\ t_{\tau }=n_p \tau$ and $n_p$ is the total number of temperature increments within $\tilde {t}$. Assuming $\tau,t_{w_i}\ll \tilde {t}$ implies that the number of temperature fluctuations in such an interval is large enough to justify the use of an ergodic assumption. We can therefore replace the time average over $t_{t_w}$ with the average with respect to its probability distribution, as

(3.8)\begin{equation} \langle J_{T,\alpha} \rangle_{t_{t_w}}=\int_1^{+\infty}{\rm d} T\gamma(T)J_{T,\alpha}. \end{equation}

Furthermore, the assumption $\tau,t_w\ll \tilde {t}$ also implies $t_{t_w}=\sum _{i=1}^{n_p} t_{w_i} = n_p \langle t_w \rangle _{\eta }$, where $\langle t_w \rangle _{\eta }$ is given by

(3.9)\begin{equation} \langle t_w \rangle_{\eta} = \int {\rm d} t_w t_w \eta(t_w), \end{equation}

so that the coarse-grained flux becomes

(3.10)\begin{align} J_{\alpha}=\int_0^{+\infty}{\rm d} p p \tilde{f}_{\alpha}(p), \end{align}

where $\tilde {f}_{\alpha }(p)$ is defined by

(3.11)\begin{equation} \tilde{f}_{\alpha}(p)=\frac{A}{M_{\alpha}}\int_1^{+\infty}{\rm d} T \frac{\gamma(T)}{T}\exp\left({-\frac{p^2}{2TM_{\alpha}}}\right)+ \frac{1-A}{M_{\alpha}}\exp\left({-\frac{p^2}{2M_{\alpha}}}\right), \end{equation}

and $A$ is

(3.12)\begin{equation} A=\frac{\tau}{\tau +\langle t_w \rangle_{\eta}}. \end{equation}

In practice, the dynamics of a two-component plasma coupled to a time fluctuating thermal boundary at a given boundary is described, at the temporal coarse-grained level, by a set of Vlasov equations (3.5) for the two species coupled to two incoming effective fluxes (3.10) at the boundary. The latter are generated by the non-thermal distributions $\tilde {f}_\alpha$ given in (3.11).

In the formalism introduced above, it is possible to evaluate an analytic expression for the coarse-grained phase-space distribution functions in the stationary configuration. To do so, one has first to determine such phase-space distributions at the boundary in the stationary configuration. This can be carried out by imposing the stationarity and the continuity conditions at the boundary, together with the symmetry condition in the momentum space

(3.13a–c)\begin{equation} \frac{\partial \tilde{f}_{\alpha}}{\partial \tilde{t}}=0,\quad \mathcal{N}_{\alpha}\tilde{f}_{\alpha}(p)=\tilde{f}_{\alpha}(-{\rm \pi},p),\quad \tilde{f}_{\alpha}(-{\rm \pi},p)=\tilde{f}_{\alpha}(-{\rm \pi},-p),\quad p>0. \end{equation}

In the expressions above $\mathcal {N}_{\alpha }$ are normalisation constants. We can now build up the solution for the entire phase space using the Jeans theorem. We then get

(3.14)\begin{equation} \tilde{f}_{\alpha}(\theta,p)=\mathcal{N}_{\alpha}\left(A\int_1^{+\infty}{\rm d} T \frac{\gamma(T)}{T}\exp\left({-\frac{H_{\alpha}}{T}}\right)+(1-A)\exp\left({-H_{\alpha}}\right)\right), \end{equation}

where $H_{\alpha }$ are the mean-field Hamiltonians defined in (2.27). The constants $\mathcal {N}_{\alpha }$ are fixed by the normalisation condition for the two species

(3.15)\begin{equation} \int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta\int_{-\infty}^{+\infty}{\rm d} p \tilde{f}_{\alpha}(\theta,p)=1. \end{equation}

Since $H_{\alpha }$ depends on $\tilde {f}_{\alpha }(\theta,p)$ only through the electrostatic potential $\phi$, this is a self-consistent problem that can be solved with respect to $\phi$ as done for the equilibrium solution (2.30), yielding $\phi =0$. As a consequence, the mean-field-Hamiltonians $H_{\alpha }$ reduce to $\tilde {H}_{\alpha }$ as in (2.33). The first term on the right-hand side of (3.14) is induced by the temperature fluctuations of the thermal boundary and is a superposition of Boltzmann exponential distributions, each of which has a temperature $T>1$ weighted with the probability $\gamma (T)$. Such a term introduces the suprathermal contribution and from now on will be referred to as the multitemperature population.

The second term on the right-hand side of (3.14) is proportional to a Boltzmann exponential at $T=1$, due to the fact that, along each interval of length $t_w$, the thermal boundary is brought back to the temperature $T=1$. From now on we will refer to this term as the thermal population. The relative contribution of the multitemperature and thermal populations is weighted by the factor $A$ given by (3.12), corresponding to the fraction of time in which the thermal boundary is set at one of the temperatures of the multitemperature population.

The shape of the distribution functions depends only on the mean-field Hamiltonians $\tilde {H}_{\alpha }$, since the system is ruled by a Vlasov-type dynamics.

Our formalism allows one to discuss the physics in the stationary configuration in terms of velocity filtration. Given the presence of suprathermal tails in the distribution functions (as induced by the multitemperature population) and an external field, the temperature inversion process is expected to take place. This can be easily verified by computing the number densities $n_{\alpha }$ and the kinetic temperatures $T_{\alpha }$ from (3.14), using the standard kinetic definitions (see e.g. Nicholson Reference Nicholson1983). For the number densities we get

(3.16)\begin{equation} \tilde{n}_{\alpha}(\theta)=\frac{A \displaystyle\int_{1}^{+\infty}{\rm d} T\dfrac{\gamma(T)}{\sqrt{T}} \exp\left({-\dfrac{2\tilde{g}}{T}\cos{\left(\dfrac{\theta}{2}\right)}}\right)+ (1-A)\exp\left({-2\tilde{g}\cos{\left(\dfrac{\theta}{2}\right)}}\right)}{A\displaystyle \int_{1}^{+\infty}{\rm d} T\dfrac{\gamma(T)}{\sqrt{T}}\int_{-{\rm \pi}}^{\rm \pi} {\rm d}\theta \exp\left({-\frac{2\tilde{g}}{T}\cos{\left(\frac{\theta}{2}\right)}}\right)+(1-A) \int_{-{\rm \pi}}^{\rm \pi} {\rm d}\theta \exp\left({-2\tilde{g}\cos{\left(\frac{\theta}{2}\right)}}\right)}, \end{equation}

while for the kinetic temperatures

(3.17)\begin{equation} \tilde{T}_{\alpha}(\theta)=\frac{\displaystyle A\int_{1}^{+\infty}{\rm d} T\gamma(T) \sqrt{T}\exp\left({-\frac{2\tilde{g}}{T}\cos{\left(\frac{\theta}{2}\right)}}\right)+ (1-A)\exp\left({-2\tilde{g}\cos{\left(\frac{\theta}{2}\right)}}\right)}{\displaystyle A\int_{1}^{+\infty}{\rm d} T\frac{\gamma(T)}{\sqrt{T}} \exp\left({-\frac{2\tilde{g}}{T}\cos{\left(\frac{\theta}{2}\right)}}\right) +(1-A)\exp\left({-2\tilde{g}\cos{\left(\frac{\theta}{2}\right)}}\right)}. \end{equation}

In the next section the above formulas will be evaluated for different choices of model parameters and compared with the output of kinetic numerical simulations.

4. Comparison with numerical simulations

We validated our model with $N-$particle numerical simulations where we integrated the equations of motion (2.23) of the system. In the numerical experiments discussed here, as a rule, we fixed $N=2^{21}$ and used a fourth-order symplectic algorithm (see e.g. Candy & Rozmus Reference Candy and Rozmus1991) with fixed time step $\delta t=10^{-4}$. We modelled the incoming energy flux from the fluctuating thermal boundary with the standard technique (see Tehver et al. Reference Tehver, Toigo, Koplik and Banavar1998; Landi & Pantellini Reference Landi and Pantellini2001) such that, when a particle of a given species crosses the boundary, it is re-injected in the system (at the bottom) with a positive velocity sampled from the flux density argument of (3.6). We note that naively re-introducing the particle with a new velocity extracted from a (half) Gaussian at temperature $T$ would break the stationary thermal equilibrium solution, as re-injected particles would have, by construction, a higher chance of having a near zero velocities (see e.g. Lepri et al. (Reference Lepri, Ciraolo, Di Cintio, Gunn and Livi2021) and references therein).

4.1. Stationary state and model parameters

Motivated by the modellisation of the Sun's atmosphere, we set the distributions of the temperature fluctuations $\gamma (T)$ and waiting times $\eta (t_w)$ as

(4.1a,b)\begin{equation} \gamma(T)=\frac{1}{T_p}\exp\left({-\frac{(T-1)}{T_p}}\right),\quad T>1,\quad \eta(t_w)=\frac{1}{\langle t_w \rangle_\eta}\exp\left({-\frac{t_w}{\langle t_w \rangle_\eta}}\right). \end{equation}

The idea is that the chromospheric (i.e. the thermal boundary) temperature fluctuates randomly in time, and the strong temperature increments (i.e. those producing a large $T$), are rather rare (see Barbieri et al. Reference Barbieri, Casetti, Verdini and Landi2024). This established, we tune the remaining parameters ($A, \tilde {g},C,T_p$) accordingly. In figures 2–6 (left panels) we show the time evolution of the kinetic energies $K_{\alpha }$ (top panels) and of the stratification parameters $q_{\alpha }$ (bottom panels) evaluated in the numerical simulations for different choices of the parameter set as

(4.2a,b)\begin{equation} K_{\alpha}=\frac{1}{N}\sum_{i=1}^N \frac{p_{j,\alpha}^2}{2M_{\alpha}},\quad q_{\alpha}=\frac{1}{N}\sum_{j=1}^N \cos{(\theta_{j,\alpha})}, \end{equation}

together with their theoretical value (indicated by the straight solid lines) in the stationary configuration, given by

(4.3a,b)\begin{equation} q_{\alpha,SS}=\int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta\int_{-\infty}^{+\infty}{\rm d} p \cos{\theta}\tilde{f}_{\alpha}(\theta,p), \quad K_{\alpha,SS}=\int_{-{\rm \pi}}^{\rm \pi}{\rm d}\theta\int_{-\infty}^{+\infty}{\rm d} p \frac{p^2}{2M_{\alpha}}\tilde{f}_{\alpha}(\theta,p). \end{equation}

In all cases, both $K_{\alpha }$ and $q_{\alpha }$ relax towards their predicted stationary value. In the right panels of the same figures we show the kinetic temperature $T_{e}$ together with the number density $n_{e}$ for the electrons – the profiles for the protons having the same shape – as a function of the spatial coordinate $\theta$. The latter were computed numerically using the standard kinetic definitions (see e.g. Nicholson Reference Nicholson1983) for a single snapshot in the stationary state and then time averaging over subsequent snapshots. In addition, using (3.16) and (3.17) we have also computed their theoretical expressions. In all cases presented here, we observe a very good match between the numerical and analytical curves.

4.1.1. The fraction of time $A$

As apparent from (3.12), the parameter $A$ controls the fraction of time in which the thermal boundary is at a given value $T$ sorted from $\gamma (T)$. In figure 2 we present two cases: $A=0.5 (\tau =0.01,\langle t_w \rangle _\eta =0.01)$ and $A=0.25(\tau =0.01,\langle t_w \rangle _\eta =0.03)$. All the other parameters are fixed to $\tilde {g}=1, T_p=4,C=400,M=100$.

Figure 2. (a,b) Evolution of proton (red for $A=0.25$ and green for $A=0.5$) and electron (blue for $A=0.25$ and orange for $A=0.5$) kinetic energies $K_\alpha$ (a) and stratification parameters $q_\alpha$ (b) as numerically computed from simulations via (4.1a,b), together with theoretical prediction for their mean value (black horizontal lines, see (4.1a,b)). (c) Numerical density (in blue for $A=0.25$ and in red for $A=0.5$) and temperature (same rule of colours) of electrons vs curvilinear abscissa of the loop $\theta$. Grey curves denote the corresponding theoretical predictions from analytical formulas (3.16)–(3.17).

Increasing $A$ corresponds to an increment of the weight of the multitemperature population in (3.14) and therefore more energy is injected into the plasma. As a consequence, the plasma becomes less stratified in density (see variation of $\tilde {n}_{\alpha }$ and of the stratification parameters $q_{\alpha }$ in figure 2). For the same reason, the mean value of the kinetic energies $K_{\alpha }$ in the non-equilibrium stationary configuration increases, as the profile of the kinetic temperature $\tilde {T}_{\alpha }$ does at each in point of the space.

4.1.2. The intensity of temperature increments $T_p$

The parameter $T_p$ (see (4.1a,b)) controls the mean value of the (exponential) distribution of the temperature fluctuations. In figure 3 we present two cases: $T_p=1$ and $T_p=4$. As above, all the other parameters are fixed and are $A=0.5 (\tau =0.01,\langle t_w \rangle _{\eta }=0.01), \tilde {g}=1,C=400,M=100$. Also, in this case, increasing the value of $T_p$ makes the energy injected by the multitemperature population increase. We thus obtain the same behaviour as in the previous case (figure 2) for all quantities.

Figure 3. (a,b) Evolution of proton (red for $T_p=4$ and green for $T_p=1$) and electron (blue for $T_p=4$ and orange for $T_p=1$) kinetic energies $K_\alpha$ (a) and stratification parameters $q_\alpha$ (b) as numerically computed from simulations via (4.1a,b), together with theoretical prediction for their mean value (black horizontal lines, see (4.1a,b)). (c) Numerical density (in blue for $T_p=4$ and in red for $T_p=1$) and temperature (same rule of colours) of electrons vs curvilinear abscissa of the loop $\theta$. Grey curves denote the corresponding theoretical predictions from analytical formulas (3.16)–(3.17).

4.1.3. The strength of the electrostatic interactions $C$

The parameter $C$ controls the strength of the electrostatic interactions between particles. In figure 4 we report the two cases $C=100$ and $C=400$. All the other parameters are fixed such that $A=0.5 (\tau =0.01,\langle t_w \rangle _{\eta }=0.01), T_p=4,\tilde {g}=1,M=100$. This figure clearly shows that the electrostatic interaction $C$ does not play any significant role in shaping the inverted density–temperature profiles, as obvious from the theory since $\phi =0$, cf. § 3.

Figure 4. (a,b) Evolution of proton (red for $C=400$ and green for $C=100$) and electron (blue for $C=400$ and orange for $C=100$) kinetic energies $K_\alpha$ (a) and stratification parameters $q_\alpha$ (b) as numerically computed from simulations via (4.1a,b), together with theoretical prediction for their mean value (black horizontal lines, see (4.1a,b)). (c) Numerical density (in yellow for $C=400$ and in black for $C=100$) and temperature (same rule of colours) of electrons vs curvilinear abscissa of the loop $\theta$. Red curves denote the corresponding theoretical predictions from analytical formulas (3.16)–(3.17).

4.1.4. The stratification parameter $\tilde {g}$

This parameter determines the strength of the external stratification field in units of the thermal energy $k_B T_0$ of the thermal boundary (i.e. the chromosphere). In figure 5 we present two cases: $\tilde {g}=0.5$ and $\tilde {g}=1$. All of the other parameters are as before. The expressions of the distribution functions of the theory reported in (3.14) depend only on $\tilde {H}_\alpha$. At the bottom (i.e. $\theta =-{\rm \pi}$), the mean-field Hamiltonians $\tilde {H}_{\alpha }$ are the same for both species and read

(4.4)\begin{equation} \tilde{H}_{\alpha}(p,\theta={-}{\rm \pi})=\frac{p^2}{2M_{\alpha}}. \end{equation}

Therefore, the two distribution functions have a fixed width in momentum space that is independent of the specific value of $\tilde {g}$. Hence, also the kinetic temperature at $\theta =-{\rm \pi}$ is independent of $\tilde {g}$, as apparent in figure 5.

Figure 5. (a,b) Evolution of proton (red for $\tilde {g}=1$ and green for $\tilde {g}=0.5$) and electron (blue for $\tilde {g}=1$ and orange for $\tilde {g}=0.5$) kinetic energies $K_\alpha$ (a) and stratification parameters $q_\alpha$ (b) as numerically computed from simulations via (4.1a,b), together with theoretical prediction for their mean value (black horizontal lines, see (4.1a,b)). (c) Numerical density (in blue for $\tilde {g}=1$ and in red for $\tilde {g}=0.5$) and temperature (same rule of colours) of electrons vs curvilinear abscissa of the loop $\theta$. Grey curves denote the corresponding theoretical predictions from analytical formulas (3.16)–(3.17).

In the right plot of figure 5 we show that, by increasing the value of $\tilde {g}$, the kinetic temperature grows with a stronger gradient going from the bottom ($\theta =-{\rm \pi}$) to the top ($\theta =0$). In practice, as implied by (3.14), with higher values of $\tilde {g}$, only the Boltzmann factors with large values of the temperature in the multitemperature population can give a contribution at higher altitudes. In practice, from all the possible value extracted from $\gamma (T)$, only high $T$ temperature fluctuations give sufficient kinetic energy to the particles so that they can climb the gravitational well up to the top of the system $\theta =0$. We note that this is essentially Scudder's gravitational filtering mechanism.

4.1.5. Varying the temperature fluctuation distribution $\gamma (T)$

Finally, we change the distribution of temperature increments $\gamma (T)$. We consider here the case of a one-sided Gaussian $\gamma (T)$ distribution

(4.5)\begin{equation} \gamma(T)=\frac{2}{\sqrt{2{\rm \pi} T_p^2}}\exp\left({-\frac{(T-1)^2}{2T_p^2}}\right),\quad T>1. \end{equation}

In figure 6 we show the case for the following values of the system's parameters: $A=0.5 (\tau =0.01,\langle t_w \rangle _\eta =0.01),\tilde {g}=1,T_p=4,C=400,M=100$. As expected, temperature increments still give origin to temperature inversion independently of the specific shape of $\gamma (T)$.

Figure 6. (a,b) Evolution of proton (green) and electron (red) kinetic energies $K_\alpha$ (a) and stratification parameters $q_\alpha$ (b) as numerically computed from simulations via (4.1a,b), together with theoretical prediction for their mean value (black horizontal lines, see (4.1a,b)). (c) Numerical density and temperature (red) of electrons vs curvilinear abscissa of the loop $\theta$. Grey curves denote the corresponding theoretical predictions from analytical formulas (3.16)���(3.17).

4.2. Limits of the theoretical approach

We now discuss the limits of the theoretical approach, using numerical simulations. For the sake of simplicity, we focus on the case in which the temperature oscillates between the two values $T=1$ and $T=T_p$ with a fixed waiting time $t_w$. In this set-up the distributions $\gamma (T)$ and $\eta (t_w)$ are given by

(4.6a,b)\begin{equation} \gamma(T)=\delta(T-T_p),\quad \eta(t_w)=\delta(t-t_w). \end{equation}

In figure 7, we report the results obtained by varying $\tau$, with $\tau =t_w$. The other parameters are fixed to $C=400,T_p=4,\tilde {g}=1,M=100$.

Figure 7. (a–d) Kinetic energies of the electrons for four different couples of values of $\tau =t_w={0.1,20,200,600}$ from left to right. (e–h) Electron temperature–density profiles numerically evaluated time averaging over the interval $t_1=300$ and $t_2=1000$ (red for temperature and blue for density) together with the theoretical predictions in grey and calculated via (3.16)–(3.17).

We observe how, on increasing both $\tau$ and $t_w$ and retaining the ratio $A={\tau }/{(\tau + t_w)}$ fixed to $0.5$, the theoretical equations computed using (3.14) start failing in reproducing the density–temperature profiles evaluated from the numerical simulations. In particular, in the left panels the time scales $\tau$ and $t_w$ are so small that the particles do not have enough time to relax to a thermal configuration, neither at $T=T_p$ or at $T=1$. The system is therefore locked in a non-equilibrium stationary configuration and the inverted density–temperature profiles computed numerically are well reproduced by the theoretical counterparts. Technically speaking, the inverted density–temperature profiles shown in figure 7 are obtained by time averaging over several snapshots from $t_1=300$ to the end time $t_2=1000$. Since the system is in a stationary state these profiles are independent of the length and of the location along the stationary state of the interval $t_2-t_1$.

For the other cases, the density–temperature profiles are numerically evaluated in the same way. As both $\tau$ and $t_w$ increase, the kinetic energy starts to develop large amplitude oscillations. This arises from the fact that, when moving towards a regime (shown in the two right columns of the figure 7) in which both $\tau$ and $t_w$ are large enough (i.e. comparable to the relaxation time), the system evolves in time passing from a thermal configuration at $T=T_p$ to a thermal configuration at $T=1$ (see figure 8).

Figure 8. (a) Electron kinetic energy as a function of time. (b) Electron densities and temperatures for the two thermal configurations $T=T_p=4$ and $T=1$. Blue and magenta, respectively, show the numerical temperature and the numerical density for the case $T=T_p=4$, together with their respective theoretical equilibrium counterparts in grey and calculated via (2.34)–(2.35). Red and black, respectively, show the numerical temperature and the numerical density for the case $T=1$, together with their respective theoretical equilibrium counterparts in grey and calculated via (2.34)–(2.35).

Of course, in all cases except when $\tau =t_w=0.1$, since the system is no longer locked in a stationary configuration, the theory cannot be applied and thus the numerically recovered density–temperature profiles evaluated via time averaging do not match the theoretical counterparts. Moreover, since the system is no longer locked in a stationary state, the numerical density–temperature profiles depend explicitly on the length and on the location of the time average interval $t_2-t_1$ for the specific run.

5. Application to the solar atmosphere

Here, we show the application of the temporal coarse graining to predicting the inverted density–temperature profiles of the Sun's atmosphere. To compute such profiles we have used (3.16)–(3.17), and (3.14) for the distribution functions. We put the base of the model at $z_b=2\times 10^3$ km (i.e. at the base of the transition region in the Sun's atmosphere) and the top at $z_t=2.2\times 10^4$ km (i.e. in the corona). We then fix the length $2L$ of half of the loop via the following equation:

(5.1)\begin{equation} 2L=(z_t-z_b){\rm \pi}=2{\rm \pi} \times 10^4\,{\rm km},\quad z=\frac{2L}{\rm \pi}\cos{\left(\frac{{\rm \pi} x}{2L}\right)}, \end{equation}

that corresponds to the typical length of a coronal loop in the solar atmosphere. With such a choice, the dimensionless parameter $\tilde {g}$ is fixed to $\tilde {g}=16.64$. For the distribution of the temperature fluctuations and the waiting times, we use (4.1a,b) discussed in § 4 above. Thus, we are now left with the two free parameters $T_p$ and $A$. We initially fix the value of $A$ to $1$ and vary $T_p$.

As shown in § 4, increasing the value of $T_p$ results in increasing the values of the temperature profile. In order to have a temperature of the order of $10^6$ K at the top of the loop, we must increase $T_p$ up to a value around $10^6$ K (as shown in figure 10b). Therefore, from now on, we fix $T_p=90$ (that correspond to a temperature of $9\times 10^5$ K). From figures 10(a) and 9(a) we observe that, for $A=1$, we have a temperature at the base of the order of $T_b=5\times 10^5$ K. Such a value is far greater than the observed one, the latter being smaller than $T_{b,obs}=1.2\times 10^4\,{\rm K}$, as the chromospheric temperature varies between $8\times 10^3$ and $1.2\times 10^4$ K, see Molnar et al. (Reference Molnar, Reardon, Chai, Gary, Uitenbroek, Cauzzi and Cranmer2019).

Figure 9. Kinetic temperature (in K, (a)) computed via (3.17) and number density computed via (3.16) (b) scaled by the mean number density $n_0=2.5\times 10^9\,\mathrm {cm}^{-3}$ as a function of the height (in km) for different values of $A$ passing from $A=1$ to $A=0.01$ as listed in the legend.

In the previous section we have shown how the temperature at the base of the loop $T_{b}$ decreases with decreasing $A$ (see figure 2). In figure 10 we plot $T_b$ as a function of $A$. We observe that, at around A=$0.022$, the temperature $T_b$ crosses the upper limit value $T_{b,obs}$. Combining these two properties that relate $T_p$ to $A$, in figure 9 we show how the inverted density–temperature profiles with $A=0.02$ and $0.01$ (i.e. very rare temperature increments with $T_p = 90$) are similar to the one observed in the Sun's atmosphere. These profiles start from $T_b \sim 1.1\times 10^4$ K and have an initial rapid rise in temperature (the transition region) followed by a slower increase in the above region (i.e. the solar corona).

Figure 10. (a) Temperature at base $T_b$ of the atmosphere $z=0$ computed via (3.17) against the parameter $A$ (3.12). The red horizontal line is $T_{b,{\rm obs}}=1.2\times 10^4$ K. (b) The temperature at the top $T_t$ computed via (3.17) as function of the mean value of the temperature fluctuations $T_p$ for $A=1$.

This change of the shape caused by varying the parameter $A$ can be understood in terms of velocity filtration, as is apparent from the structure of the electron VDFs at the base of the model ($z=z_1=2.3\times 10^3$ km) as shown in figure 11 where we plot in a semilogarithmic scale the distribution of the signed electron kinetic energy $\tilde {f}_e (\text {sign}(p)p^2/2)$ divided by the density. By doing this, it becomes clearer whether the distribution function in question is thermal or not, since for a thermal distribution (i.e. a Gaussian) one would get two straight lines symmetric with respect to zero. As can be seen from figure 11(a,b), the thermal (Gaussian) part of the electron VDFs increasingly dominates the core of the distribution, as the value of $A$ decreases. This is the reason why the temperature at the base decreases crossing $T_{b,{\rm obs}}=1.2\times 10^4$ K and tends to $T_0=10^4$ K. Moreover, by fixing $A=0.02$, from figure 12 we observe that, in virtue of the energy conservation (thanks to the velocity filtration mechanism), the cold central thermal core of the electron DFs at the base of the model $z=z_1=2.3\times 10^3$ km progressively disappears by passing through the transition region $z=z_2=2.3\times 10^3$ km and is totally filtered out in the corona $z=z_3=1.1\times 10^4$ km. As a consequence of this, one has a rapid increase of temperature through the transition region and a slow increase of temperature in the corona.

Figure 11. (a) Electron VDFs computed via (3.14) normalised by the electron number densities (computed via (3.16)) at the base of the atmosphere at $z=z_1=2.3\times 10^3\,\mathrm {km}$ (the base of the transition region) for different values of $A$. In (b), a magnification of the central region of the same VDFs is shown to emphasise the central thermal (Gaussian) core. In (c) are depicted the electron VDFs in the corona $z=z_3$ for the same values of $A$ as the other panels. In this panel all the electron VDFs are collapsed onto each other because of velocity filtration.

Figure 12. (a) Electron VDFs (computed via 3.14 normalised by the electron number densities (computed via 3.16) for $A=0.02$, computed at three increasing heights ($z=z_1,z_2,z_3$) listed in the legend as a function of the signed kinetic energy. In (b), a magnification of the central region of the same VDFs is shown to appreciate the disappearance of the Gaussian profile with height, as expected when velocity filtration is at work.

In figure 11(c) we can observe that all the electron VDFs rescaled by the density in the corona at $z=z_3=1.1\times 10^4$ km for all the values of $A$ collapse onto each other, due to the fact that the central thermal core has been filtered out in the corona itself. Since the rescaled distribution functions are the same for all values of $A$, then the temperatures collapse to the same curve in the corona, as shown in figure 9(a).

The number density profiles are inverted with respect to the temperature profiles with an initial strong negative gradient (the transition region) followed by a slower decrease within the corona. By decreasing the value of $A$, the number density at the bottom of the system increases and at the same time the number density at the top decreases. Again this is due to the fact that the cold population in (3.14) becomes more and more pronounced for decreasing values of $A$, that is, we now have more particles that are not able to climb higher in the gravity well. Since the total number of particles is fixed independently of the values of $A$ as dictated by the normalisation condition (3.15), then, if the density at the bottom increases, the density at the top (the corona) should naturally decrease.

6. Summary and perspectives

In this work, we have addressed the question of generating suprathermal distribution functions and temperature inversion in a confined plasma structure by means of temperature increments at its edges produced by some heating events. We tackled this problem by introducing a novel formalism to treat the plasma dynamics in contact with a thermal boundary whose temperature fluctuates in time recently studied by Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024). We have shown how the multi-component (in the case discussed here, electrons and protons) plasma dynamics can be efficiently modelled in terms of the temporal coarse-grained distribution functions $\tilde {f}_{\alpha }$ for the two species. The dynamics just mentioned is governed by the effective system of Vlasov equations (3.5) in contact with the non-thermal distributions $\tilde {f}_{\alpha }(p)$ at its boundary (3.11).

We have obtained analytical expressions for the distribution functions of the two species in the non-equilibrium stationary configuration (3.14), in terms of the mean-field Hamiltonians $\tilde {H}_{\alpha }$ (2.33). In this set-up we can interpret the anti-correlated density–temperature profiles in terms of velocity filtration, in analogy to the Scudder mechanism. However, here, the suprathermal tails are not imposed as in the Scudder approach but are self-consistently produced by the temporal variations of the temperature of the thermal boundary (the chromosphere) yielding particles with different temperatures that mix together giving rise to the multitemperature contribution to the coarse-grained distribution functions in (3.14).

We have tested the theoretical predictions for temperature inversion against numerical simulations, finding an excellent agreement between the theory and the numerical results. In addition, we have applied our theoretical formalism to the specific case of the Sun's atmosphere. We have shown how decreasing the percentage of time of duration of the temperature increments (i.e. decreasing $A$) a transition region naturally forms in the inverted density–temperature profile (see figure 9). Moreover, in the condition of very rare temperature fluctuations $A\sim 0.02$ at $T_p \sim 10^6$ K we were able to produce an inverted temperature–density profile very similar to the one observed in the Sun's atmosphere with a base temperature below $1.2\times 10^4$ K (as observed e.g. by Molnar et al. Reference Molnar, Reardon, Chai, Gary, Uitenbroek, Cauzzi and Cranmer2019), transition region (wider than the observed one) and then a million kelvin corona.

In the present work we have derived the coarse-grained Vlasov dynamics (3.5) coupled to the effective incoming flux (3.10) in the one-dimensional case and using the HMF modellisation. The extension of our procedure to systems in higher dimensionality is expected to be possible because the temporal coarse graining does not depend explicitly on the dimensionality of the system and on the HMF modellisation for the self-electrostatic interactions at hand. This, in principle, opens up the possibility of including the contribution of a magnetic field along the curvilinear abscissa for the simple loop plasma model discussed here. We note that such an additional ingredient is interesting not only for the coronal loops in the solar atmosphere but also in the context of fusion plasmas confined inside a tokamak machine (see Goedbloed, Keppens & Poedts Reference Goedbloed, Keppens and Poedts2010; Ciraolo et al. Reference Ciraolo, Bufferand, Di Cintio, Ghendrih, Lepri, Livi, Marandet, Serre, Tamain and Valentinuzzi2018).

The formalism of temporal coarse graining could be also useful to exospheric models of heliospheric plasmas (i.e. the plasma that populates the interplanetary space) that describe interplanetary plasma as a non-collisional medium in contact with a fixed distribution at its boundary (Chamberlain Reference Chamberlain1960; Jockers Reference Jockers1970; Lemaire & Scherer Reference Lemaire and Scherer1971; Maksimovic, Pierrard & Lemaire Reference Maksimovic, Pierrard and Lemaire1997; Lamy Reference Lamy2003; Zouganelis et al. Reference Zouganelis, Maksimovic, Meyer-Vernet, Lamy and Issautier2004). Such systems have the boundary at the base of the heliosphere (outer zone of the solar atmosphere), while the theoretical formalism discussed here allows one to implement an exospheric model that has the base in the high chromosphere and that is able to reproduce the plasma of the solar atmosphere together with the heliospheric plasma. This approach can potentially give a model and a mechanism to produce suprathermal electrons in the heliospheric plasma, whose presence is suggested by in situ measurements (see for example Pilipp et al. Reference Pilipp, Muehlhaeuser, Miggenrieder, Montgomery and Rosenbauer1987; Berčič et al. Reference Berčič, Larson, Whittlesey, Maksimović, Badman, Landi, Matteini, Bale, Bonnell and Case2020; Halekas et al. Reference Halekas, Whittlesey, Larson, McGinnis, Maksimovic, Berthomier, Kasper, Case, Korreck and Stevens2020; Maksimovic et al. Reference Maksimovic, Bale, Berčič, Bonnell, Case, Dudok de Wit, Goetz, Halekas, Harvey and Issautier2020).

Finally, our formalism could also be used in the context of laser–plasma interaction, where a region of the system is in thermal contact with another where the temperature fluctuates as a consequence of the external laser pumping (see e.g. Gibbon & Förster Reference Gibbon and Förster1996; Atzeni Reference Atzeni2001).