1. Introduction
The lower solar atmosphere consists of a dense, collisional plasma in thermal equilibrium at temperatures around 10 000 K. In contrast, the outer atmospheric layer, the solar corona, is composed of a much more rarefied but significantly hotter plasma, with temperatures reaching 1–2 million K. Understanding the physical mechanisms responsible for this dramatic temperature increase remains one of the central unresolved challenges in solar and plasma physics, commonly referred to as the coronal heating problem (Ionson Reference Ionson1978; Heyvaerts & Priest Reference Heyvaerts and Priest1983; Scudder Reference Scudder1992a, Reference Scudderb; Gudiksen & Nordlund Reference Gudiksen and Nordlund2005; Klimchuk Reference Klimchuk2006; franco Rappazzo et al. Reference Rappazzo, franco, Velli, Einaudi and Dahlburg2008; Pontieu et al. Reference De Pontieu, McIntosh, Carlsson, Hansteen, Tarbell, Boerner, Martinez-Sykora, Schrijver and Title2011; Parnell & De Moortel Reference Parnell and De Moortel2012; Rappazzo & Parker Reference Rappazzo and Parker2013; Wilmot-Smith Reference Wilmot-Smith2015; Howson, De Moortel & Reid Reference Howson, De Moortel and Reid2020; Hau et al. Reference Hau, Chang, Lazar and Poedts2025).
A recent study introduced a novel kinetic
$N$
particle model of coronal loops in the solar atmosphere (Barbieri et al. Reference Barbieri, Casetti, Verdini and Landi2024a), showing through numerical simulations that rapid, stochastic temperature increments in the high chromosphere, faster than the electron crossing time of the loop, can drive the plasma into a non-equilibrium stationary state characterised by anti-correlated temperature and density profiles, a phenomenon referred to as temperature inversion. In a subsequent work (Barbieri et al. Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b), the model was analysed within a kinetic framework, leading to the development of a temporal coarse-graining formalism capable of describing the long-time behaviour of the system. It was proven that the plasma dynamics can be captured by a pair of coupled Vlasov equations for coarse-grained distribution functions, with stationary solutions exhibiting suprathermal velocity distributions. These suprathermal tails are self-consistently produced by the temperature increments, and the temperature inversion naturally arises as a result of gravitational filtering (i.e. more energetic particles in the suprathermal tails can rise higher than the others in the gravity well, Scudder Reference Scudder1992a, Reference Scudderb). Remarkably no fine-tuning of parameters is required to produce temperature inversion, as long as the characteristic time scale of thermal fluctuations remains shorter than the electron crossing time. Moreover, these results are independent of the self-consistent electrostatic field. The formalism was subsequently applied to solar-type main-sequence stars, predicting the presence of a hot corona for all such stars (Barbieri et al. Reference Barbieri, Casetti, Verdini and Landi2025).
In the present study, we extend this analysis by relaxing the restriction on fluctuation time scales and generalising the temporal coarse-graining approach accordingly. We show that the dynamics can still be described by two coupled kinetic equations, but with an additional term arising from the coarse-graining procedure itself. This term becomes relevant when temperature fluctuations occur on time scales comparable to the electron crossing time, resulting in species-dependent density and temperature profiles and a suppression of gravitationally driven temperature inversion. If the fluctuation time scales become much larger than the proton crossing time, the system exhibits oscillatory behaviour between multiple thermal configurations. In this regime, the stationary coarse-grained distribution corresponds to a superposition of thermal states, and temperature inversion can still be recovered under specific definitions of the coarse-grained temperature.
The paper is organised as follows. In § 2, we review the kinetic
$N$
-particle model introduced by Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024a). In § 3, we develop the extended temporal coarse-graining formalism and show that it recovers the earlier formulation as a limiting case under appropriate assumptions. § 4 contains a detailed numerical analysis of plasma dynamics across different fluctuation time scales, interpreted within the new theoretical framework. Finally, § 5 summarises the main results and outlines potential future developments.

Figure 1. Schematic diagram of the two-component plasma loop model. The vertical axis, designated as
$z$
, represents the altitude within the atmosphere, while the curvilinear abscissa is denoted as
$x$
. The symbols
$\sigma _{m,\alpha }= m_{\alpha } n_S$
and
$\sigma _{\alpha }=e_{\alpha } n_S$
, with
$\alpha =\{e,i\}$
are the surface mass density and the surface charge density of the species
$\alpha$
, while
$n_S$
is the surface number density.
2. Two-component gravitationally bound plasma model
Let us now briefly describe the model for geometrically confined plasma structures, specifically the coronal loops that are ubiquitous in the Sun’s atmosphere (see, e.g. Aschwanden Reference Aschwanden2005) introduced by Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024a, Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casettib). The loop is modelled as a semicircular tube of length
$2L$
and cross-sectional area
$S$
, with the charge distribution discretised into
$n_S$
density sections. These consist of
$2N$
sections, each containing electrons of charge
$-e$
and mass
$m_e$
, and protons of charge
$e$
and mass
$m_i$
. All particles are subjected to an external field comprising the gravitational field and the Pannekoek–Rosseland electrostatic field (see Pannekoek Reference Pannekoek1922; Rosseland Reference Rosseland1924). For an in-depth treatment, see the work of Belmont et al. (Reference Belmont, Grappin, Mottez, Pantellini and Pelletier2013). It is assumed that all quantities are symmetric with respect to the apex of the loop. Electrostatic interactions are treated using a multimode approximation known as the Hamiltonian mean field (HMF) model (Antoni & Ruffo Reference Antoni and Ruffo1995; Chavanis, Vatteville & Bouchet Reference Chavanis, Vatteville and Bouchet2005; Elskens & Escande 2019), in which the electrostatic potential is expanded in Fourier series and truncated in the first mode, as discussed by Barbieri et al. (Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b). A scheme of the two-component loop plasma model is shown in figure 1. Under these assumptions, the equations of motion for each particle
$j$
are given by

where
$e_{\alpha }$
is equal to
$e$
for protons and
$-e$
for electrons, while the electric field
$E$
reads as

the parameter
$Q$
is given by

the quantities
$q_{n,\alpha }$
are given by

$\alpha \in \{e,i\}$
denotes the species (here electrons or protons),
$g={G M_{\odot }}/{R_{\odot }^2}$
is the gravitational field at the surface of the Sun where
$M_{\odot }$
is the solar mass and
$R_{\odot }$
is the solar radius, and
$x$
is the spatial coordinate (i.e. the curvilinear abscissa of the loop). Henceforth, the parameter
$Q$
is referred to as the charge imbalance parameter and the
$q_{\alpha }$
as the stratification parameters. If the majority of particles of species
$\alpha$
are symmetrically concentrated at the base of the loop, then
$q_{\alpha } \approx -1$
; a uniform distribution yields
$q_{\alpha } \approx 0$
and a concentration at the apex results in
$q_{\alpha } \approx 1$
. A non-zero
$Q$
therefore implies a net charge imbalance in the system.
The system is assumed to be in ideal thermal contact with a thermal boundary that mimics a fully collisional chromosphere.
2.1. System of units
The following units are respectively defined for velocity, mass and length:

Using these units, the equations of motion given in (2.1) can be rewritten in dimensionless form as

where the external and electrostatic forces are

with the parameters
$Q$
and
$q_{\alpha }$
given by

In these equations,
$M_{\alpha }$
equals the mass ratio
$M=m_i/m_e$
for protons and
$M=1$
for electrons. The coordinate
$\theta$
represents the dimensionless spatial position. The dynamics of the system is thus fully characterised by the three dimensionless parameters

where
$n_0=n_S N/L$
is the average density of a given species. The parameters
$C$
and
$\tilde {g}$
quantify the strengths of the electrostatic interaction and the external field, respectively, in units of thermal energy. Unless otherwise specified, all equations and plotted quantities hereafter are expressed in these dimensionless units.
2.2. Vlasov dynamics and thermal equilibrium solution
In the mean-field limit, the phase-space dynamics of the system is governed by two coupled Vlasov equations, one for each species

where
$f_\alpha$
are the single-particle distribution functionsFootnote 1 of species
$\alpha$
, and
$F_{\alpha }$
is the mean-field force derived from the Hamiltonian

where the charge imbalance parameter
$Q$
is a functional of the distribution functions

The thermal equilibrium solution corresponds to the isothermal atmosphere, where the distribution function takes the form

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

3. Extended temporal coarse-graining
In this section, we develop a general theory of temporal coarse-graining and subsequently show how it reduces to the specific case presented by Barbieri et al. (Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b).
3.1. Extended temporal coarse-graining: kinetic formalism
If the boundary temperature
$T_0$
is held constant, the loop evolves to thermal equilibrium at the same temperature. However, as discussed by Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024a), the chromosphere is highly dynamic. We thus assume its temperature fluctuates due to heating events of amplitude
$\Delta T$
, drawn from a distribution
$\gamma (\Delta T)$
, with duration
$\tau$
, separated by waiting times
$t_w$
drawn from a distribution
$\eta (t_w)$
, during which the boundary returns to temperature
$T_0$
. A schematic representation of this time-dependent thermal boundary is shown in figure 2. In what follows, we present a general formulation of the temporal coarse-graining theory. Specifically, we consider the time average of the Vlasov equations over a coarse-graining interval
$\tilde {t}$
, such that


Figure 2. Scheme of the time series of the temperature of the thermal boundary. During the time intervals of duration
$\tau$
, the temperature increases by an amount
$\Delta T$
, and during the waiting times
$t_w$
, it returns to the value
$T_0$
.
We also focus on the stationary condition

Averaging the Vlasov equation over the time interval
$\tilde {t}$
yields

where the coarse-grained distribution function
$\tilde {f}_{\alpha }$
is defined via time-averaging as

We then decompose the distribution function as
$f_\alpha = \tilde {f}_\alpha + \delta f_\alpha$
, where
$\tilde {f}_\alpha$
is the coarse-grained component and
$\delta f_\alpha$
represents the fluctuations around it. Substituting this decomposition into the averaged equation, we obtain

Here,
$\delta E$
is the fluctuating part of the electrostatic field, given by

Time-averaging the fluctuating boundary conditions leads to a coarse-grained energy reservoir described by

with

We can now split the left-hand side of (3.5) into two physically distinct contributions:
-
(i) mean-field term (MF), describing the standard Vlasov evolution:
(3.9)\begin{equation} \biggl (\frac {\partial \tilde {f}_{\alpha }}{\partial \tilde {t}}\biggl )_{MF} = -\frac {p}{M_\alpha }\frac {\partial \tilde {f}_{\alpha }}{\partial \theta }-F_{\alpha }[\tilde {f}_{\alpha }]\frac {\partial \tilde {f}_{\alpha }}{\partial p}\,; \end{equation}
-
(ii) coarse-graining term (CG), which originates from the time-averaging of the fluctuations
(3.10)\begin{equation} \biggl (\frac {\partial \tilde {f}_{\alpha }}{\partial \tilde {t}}\biggl )_{CG}=-\mathrm{sign}(e_{\alpha }) C \biggl \langle \delta {E} \frac {\partial \delta f_{\alpha }}{\partial p} \biggl \rangle _{\tilde {t}}\,. \end{equation}
This second term is absent in the restricted formulation described by Barbieri et al. (Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b), and introduces corrections that depend explicitly on the electrostatic coupling and on the magnitude of the fluctuations.
We now aim to estimate under which conditions the coarse-graining correction becomes relevant. To this end, we introduce the characteristic time scale over which the energy of each species is redistributed along the loop, referred to as the relaxation time
$t_{R,\alpha }$
. This quantity can be approximated as the minimum between the thermal and gravitational crossing times:

If the heating duration
$\tau$
and the waiting time
$t_w$
are both much smaller than
$t_{R,\alpha }$
, we do not expect significant energy fluctuations within the coarse-graining interval
$\tilde {t}$
, since the system has not sufficient time to redistribute energy along the loop. In this regime, we expect that the time-averaged distribution
$\tilde {f}_{\alpha }$
closely approximates the instantaneous distribution
$f_{\alpha }$
. Conversely, if
$\tau$
and/or
$t_w$
are comparable to or larger than
$t_{R,\alpha }$
, the system cannot fully equilibrate between heating events, leading to significant energy fluctuations induced by the time-dependent thermal boundary. In this case,
$f_{\alpha } \neq \tilde {f}_{\alpha }$
and the deviation can be estimated as

Given this estimate for the amplitude of fluctuations, we can also assess the magnitude of the coarse-graining correction terms in the kinetic equation

Equation (3.13) thus provides a criterion for estimating the relevance of fluctuations, analogous in spirit to the Ginzburg criterion used to assess fluctuation effects near phase transitions.
3.2. Restricted formulation of the temporal coarse-graining
Using (3.13), we can estimate the magnitude of the coarse-graining term in the regime

in which the heating duration and waiting times are both much shorter than the electron relaxation time. In this limit, the coarse-graining correction is negligible, effectively of order zero. As a result, the plasma dynamics can be accurately described by a set of Vlasov equations for the coarse-grained distribution functions,

In conclusion, under condition (3.14), we recover the specific stationary-state regime previously discussed by Barbieri et al. (Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b). Analytical stationary solutions compatible with the energy reservoir described by (3.7) are

The constants
$\mathcal{A}_{\alpha }$
are normalisation constants so that
$\tilde {f}_{\alpha }$
are normalised to
$1$
. The stationary-state distribution functions can be interpreted as the sum of two components: a thermal distribution at the base temperature
$T_0=1$
(i.e. the reference temperature set by the thermostat) and a non-thermal contribution arising from the average of thermal distributions at higher temperatures
$T\gt T_0$
, weighted by the probability distribution
$\gamma (T)$
of temperature fluctuations. The amplitude of this non-thermal contribution is proportional to
$A$
, which represents the fraction of time the thermostat deviates from
$T_0$
, that is, the fraction of time during which the chromosphere is actively heated. At low altitudes
$z$
, the thermal component dominates; however, its influence decreases with height due to the suppression by the gravitational potential embedded in
$H_{\alpha }$
. In contrast, the non-thermal component becomes increasingly significant at larger
$z$
, where velocity filtration allows higher-energy particles to populate greater heights, manifesting as suprathermal tails in the distribution functions. From these expressions, one can compute the number density and kinetic temperature as moments of
$\tilde {f}_{\alpha }$
. For instance, the number density is given by

and the kinetic temperature is

It is important to note that in this regime, both the density and temperature profiles are the same for electrons and protons, and are independent of the electrostatic coupling parameter
$C$
. However, this symmetry breaks in the extended formulation as discussed previously. This implies that, when
$C \neq 0$
, the stationary density and temperature profiles of electrons and protons may differ, and the mechanism of gravitational filtering may no longer guarantee temperature inversion.
In the next section, we will develop a two-fluid formalism based on this extended kinetic description to assess the role of electrostatic fluctuations more precisely.
3.3. Extended temporal coarse-graining: two-fluid formalism
We can derive a set of two-fluid equations directly from the coarse-grained kinetic equations established in § 3.1. The method relies on computing the moments of the distribution functions. Specifically, for any function of momentum
$\psi (p)$
, we define the moment

Multiplying the kinetic equation by
$\psi (p)$
and integrating over
$p$
, we obtain a hierarchy of fluid-like equations

For
$\psi (p) = 1$
, we obtain the continuity equation

where the coarse-grained density
$n_\alpha$
and mean velocity
$u_\alpha$
are defined as

For
$\psi (p) = p / M_\alpha$
, we obtain the momentum balance equation

where
$P_\alpha$
is the coarse-grained kinetic pressure of species
$\alpha$
, defined as

and
$R_\alpha$
is the contribution from the coarse-graining correction term

This additional term
$R_\alpha$
quantifies the strength of the fluctuations due to the electrostatic field, and provides a way to measure the influence of coarse-graining on the macroscopic momentum balance. In regimes where
$R_\alpha$
is non-zero, the coarse-graining term modifies the dynamics, leading to deviations from the standard Vlasov-based mean-field predictions. In § 4.4, we will use this two-fluid formalism in conjunction with numerical simulations to evaluate
$R_\alpha$
, and examine how the presence of fluctuations changes the stationary profiles of density and temperature across different regimes.
3.4. Superposition of many thermal solutions
We now consider the opposite regime with respect to the restricted coarse-graining scenario presented in § 3.2. Here, we assume that the system has sufficient time to relax to thermal equilibrium during each heating or cooling phase. That is, the characteristic relaxation times
$t_{R,\alpha }$
are much shorter than both the duration of heating events
$\tau$
and the average waiting time between them
$\langle t_w \rangle$
, that is,

This corresponds to a regime where the plasma continuously oscillates between well-defined thermal equilibria. In this regime, the long-time dynamics of the system, on time scales much larger than
$\tau$
and
$\langle t_w \rangle$
, can be described by a coarse-grained distribution that results from averaging over the sequence of these thermal states. The effective distribution is then given by

where
$f_{\alpha ,T}(\theta ,p)$
is the thermal solution at the temperature
$T$
computed via (2.13). These distribution functions exhibit suprathermal tails, arising from the superposition of multiple thermal equilibrium statesFootnote 2 contained in the first term of (3.27). As a result, the density and temperature profiles derived from such distributions necessarily exhibit temperature inversion, a signature feature of gravitational filtering of higher-energy tails. The inverted profiles can be explicitly calculated analytically and we get for the number density

where
$n_{T}(\theta )$
is given by

while for the kinetic temperature, we get

It is important to emphasise that the temperature obtained from (3.30) is the kinetic temperature of the time-averaged distribution function: since the system evolves through a sequence of distinct thermal states, this expression does not represent the time-average of the individual equilibrium temperatures. The latter, which characterises the mean energy content of the sequence of thermal states, is given by

We note that, just as in the opposite regime described in § 3.2, the coarse-grained temperatures computed via (3.30) and (3.31), as well as the number densities obtained from (3.28), are independent of the self-consistent electrostatic interaction, and are thus identical for both species.
4. Dynamics of the plasma model
We now apply the coarse-graining formalism to the model of a gravitationally confined plasma described in § 2. Barbieri et al. (Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b) demonstrated that the plasma develops a hot corona, i.e. relaxes towards a non-thermal stationary configuration exhibiting temperature inversion, if two key conditions are met:
-
(i) the time scales of thermal fluctuations (
$\tau$ and
$t_w$ ) must be much shorter than the electron relaxation time;
-
(ii) the heating events must be rare and intense, i.e.
$\tau \ll \langle t_w \rangle$ and
$\Delta T \gg T_0$ .
In what follows, we explore other dynamical regimes, i.e. we relax the condition
$\tau , t_w \ll t_{R,e}$
, while keeping the constraints on
$\Delta T$
and
$\langle t_w \rangle$
. For clarity, we present a simplified case in which the boundary temperature alternates between two fixed values,
$T = 1$
and
$T = 1 + \Delta T$
, with fixed waiting time
$t_w$
. The associated distributions are

Although this case is specific, the physical conclusions that follow are general.
4.1. Computation of the relevant quantities
To characterise the dynamical evolution of the plasma system, we compute both the kinetic energy and the stratification parameter for each species. These quantities are given by

To compute the moments of the distribution functions of a physical quantity
$\psi (p)$
given by (3.4), we use a two-step time-coarse-graining procedure; furthermore, the moments of the distribution functions of a physical quantity
$\psi (p)$
have been computed at a temporal coarse-graining level as follows.
-
(i) Instantaneous moments of the distribution functions are computed at each time step.
-
(ii) These moments are then time-averaged over a coarse-graining window
$\tilde {t} \gg \tau ,t_w$ , ensuring smoothing over thermal boundary fluctuations.
Using the definition of the coarse-grained distribution functions in (3.4), this procedure recovers the kinetic moments defined in (3.19), ensuring that all coarse-grained quantities depend solely on the time-averaged distributions.
As for temperature, since it is defined through the ratio of two moments, there is no single way to define it. In particular, we can define temperature at the level of temporal coarse-graining in two different ways.
-
(i) The first defines temperature by time-averaging the instantaneous kinetic temperature
(4.3)where\begin{equation} T_{ST,\alpha }(\theta ,\tilde {t})=\frac {1}{\tilde {t}}\int _{\tilde {t}} T_{\alpha }(\theta ,t)\,\text{d}t,\quad T_{\alpha }(\theta ,t)=\frac {P_{\alpha }(\theta ,t)}{n_{\alpha }(\theta ,t)} , \end{equation}
$P_{\alpha }$ is kinetic pressure and
$n_{\alpha }$ is the number density at time
$t$ . Here, we assume that the measurements of the moments can be done with a resolution much shorter than the coarse-graining time-scale allowing to recover the microscopic dynamic of the system.
-
(ii) If the dynamics of the system below the coarse-graining time-scale is not accessible by the measurements, it is still possible to define an averaged temperature as the ratio of the time-averaged pressure and density,
(4.4)From a physical point of view, this definition is meaningful when the system dynamics below the coarse-graining time scale\begin{equation} \tilde {T}_{\alpha }(\theta ,\tilde {t}) = \frac {\tilde {P}_{\alpha }(\theta ,\tilde {t})}{\tilde {n}_{\alpha }(\theta ,\tilde {t})} . \end{equation}
$\tilde {t}$ is not accessible. In this regime, the behaviour of the two species is described by the coarse-grained distribution functions
$\tilde {f}_{\alpha }$ , and the kinetic temperature is defined using the kinetic formulation applied to
$\tilde {f}_{\alpha }$ , as specified by the present definition.
Both definitions are meaningful in the coarse-graining context, but they may yield different results depending on the magnitude of fluctuations
$\delta f_{\alpha }$
. Indeed, by decomposing the distribution function as
$f_{\alpha } = \tilde {f}_{\alpha } + \delta f_{\alpha }$
, it is possible to establish a link between the two definitions of temperature given by (4.3) and (4.4). After some algebraic manipulation, we obtain

where
$\delta n_{\alpha }$
is given by

$\delta P_{\alpha }$
is given by

and
$\delta u_{\alpha }$
is given by

The functional
$\mathcal{K}_{\alpha }[\delta f_{\alpha }]$
measures the ratio between the two temperature definitions. Furthermore, we have that
$\mathcal{K}_{\alpha }[\delta f_{\alpha }]\rightarrow 1$
when
$\delta f_{\alpha } \rightarrow 0$
. This implies that the two temperature definitions are equal only if the fluctuations are negligible. This difference becomes important in regimes where fluctuations are non-negligible and will be discussed in detail in the subsequent sections.
4.2. Numerical simulation
The
$N$
-particle dynamics described by (2.6) is solved numerically using a fourth-order symplectic integrator (Candy & Rozmus Reference Candy and Rozmus1991). The thermal flux from the fluctuating boundary is modelled using a standard approach (see Tehver et al. Reference Tehver, Toigo, Koplik and Banavar1998; Landi & Pantellini Reference Landi and Pantellini2001). Specifically, when a particle interacts with the thermal boundary, it is re-injected with a momentum
$p$
sampled from a flux-weighted distribution:

where
$r \in [0,1]$
is a uniformly distributed random variable,
$T$
denotes the instantaneous temperature of the thermal boundary and
$p$
denotes the momentum of the particle. We note that naively reintroducing the particle with a new velocity drawn from a half-Gaussian distribution at temperature
$T$
would violate the stationary thermal equilibrium, as re-injected particles would, by construction, have a higher probability of possessing velocities close to zero.
To impose central symmetry, we restrict the simulation domain to
$\theta \in [-\pi ,0]$
. When a particle reaches the upper boundary at
$\theta = 0$
, central symmetry is enforced by applying the transformation


Figure 3. Left panels: time evolution of the kinetic energies
$K_{\alpha }$
(top) and stratification parameters
$q_{\alpha }$
(bottom) for protons (green) and electrons (orange), computed from simulations using (4.2). Right panels: time-averaged temperature and density profiles for electrons (red dashed lines) and protons (blue lines). The upper panel shows the coarse-grained temperatures calculated using (4.4), while the lower panel uses the kinetic temperature definition from (4.3). Theoretical predictions from (3.17) (for density) and (3.18) (for temperature) are shown as grey lines.
4.3. Stationary-state regime
As previously discussed, when the time scales associated with temperature fluctuations, namely,
$\tau$
and
$t_w$
, satisfy the condition

the plasma undergoes a relaxation process towards a non-thermal stationary configuration characterised by inverted density and temperature profiles. Results from a simulation in this regime are shown in figure 3. To generate these results, simulation parameters were fixed as follows:
$\tau =0.05, t_w = 0.5, g=1, M=100, C=400, \Delta T=90, N=2^{16}$
, and time step
$\text{d}t=4\times10^{-3}$
. The left panels in figure 3 illustrate the temporal evolution of the kinetic energies
$K_{\alpha }$
and the stratification parameters
$q_{\alpha }$
for both electrons and protons. As shown, both quantities tend towards a stationary regime. The coarse-grained temperature and density profiles, obtained by time averaging during the steady phase, as described in § 4.1, are shown in the right panels of figure 3. Theoretical predictions were computed using (3.18) for the temperature profiles and (3.17) for the number densities. The top right panel shows temperatures computed according to the definition in (4.4), while the bottom right panel uses the definition in (4.3). In both cases, the simulations produce identical results for both species, in excellent agreement with the theoretical predictions. This agreement stems from the fact that fluctuations are negligible in this regime. As a result, according to (4.5), the two definitions of coarse-grained temperature must coincide. In this regime, the general temporal coarse-graining formalism described in § 3.1 simplifies substantially. Specifically, the fluctuation term
$\delta f_{\alpha }$
is negligible and the additional coarse graining term in the kinetic equation (3.10) vanishes. Consequently, the system reduces to the stationary form of the Vlasov equation, and the stationary solutions are those given by (3.16). Remarkably, these solutions are independent of the self-consistent electrostatic interaction parameter
$C$
, as confirmed by the numerical simulations.
4.4. Hybrid state regime
By increasing the time scales
$\tau$
and
$t_w$
, the system reaches the following regime:

In this regime, electrons do not have sufficient time to reach thermal equilibrium at temperature
$T = 1 + \Delta T$
during a heating event of duration
$\tau$
, but they do relax to the base temperature
$T_0 = 1$
during the longer waiting time
$t_w$
. Protons, by contrast, remain in the regime described in the previous section, since both
$\tau$
and
$t_w$
are still much smaller than their relaxation time
$t_{R,i}$
. We refer to this scenario as the ‘hybrid state regime’. Numerical simulations in this regime reveal that the plasma dynamics depends sensitively on the presence or absence of self-consistent electrostatic interactions. Accordingly, we first analyse the case
$C = 0$
before introducing the effects of a finite
$C$
.
4.4.1. Independent atmospheres
We begin by examining the dynamics in the absence of self-consistent electrostatic interactions, i.e. with
$C = 0$
. For the simulations shown in figure 4, the model parameters were set as follows:
$\tau = 0.5$
,
$t_w = 25$
,
$\tilde {g} = 1$
,
$M = 1836$
,
$C = 0$
,
$\Delta T = 90$
,
$N = 2^{16}$
and time step
$\text{d}t = 4 \times10^{-3}$
. The left panels of figure 4 show that both the electron kinetic energy and the stratification parameter exhibit picks associated with the heating events of duration
$\tau$
. During the waiting time
$t_w \gg t_{R,e}$
, the electrons relax back towards a thermal configuration at temperature
$T_0 = 1$
. As a result, the electron dynamics is no longer truly stationary, in contrast to the protons, which remain effectively stationary due to
$\tau , t_w \ll t_{R,i}$
. Nevertheless, when observed on a coarse-grained time scale,
$\tilde {t}$
, as defined in (4.4), are identical for electrons and protons, and also match the analytical prediction given by (3.18). This can be understood as follows. First of all,
$C=0$
and, in addition, the coarse-grained dynamics is stationary so that (4.11) is satisfied. These conditions imply that the general kinetic equation given in (3.5) reduces to a stationary system of Vlasov equations for the two species, as in the steady-state regime described in § 4.3. Consequently, the temperature profiles computed via (4.4) must be identical for both species and match the theoretical profile given in (3.18). Moreover, this also justifies the equality of the coarse-grained number densities for electrons and protons, which coincide with the analytical profile from (3.17). These results are fully consistent with the temporal coarse-graining two-fluid formalism presented in § 3.3. In this specific case, the two-fluid equations of motion, (3.23), reduce to the hydrostatic balance:

Given the condition
$\tilde {n}_e = \tilde {n}_p$
, we can subtract the equations for the two species to obtain

This relation admits the equality of coarse-grained pressures as a valid solution. Combining the equalities of pressure and density, we conclude that the coarse-grained temperatures must also be equal, i.e.
$\tilde {T}_e = \tilde {T}_p$
. However, since the electron dynamics is no longer locked in the non-equilibrium stationary state, the condition
$\tilde {f}_e = f_e$
no longer holds, and the fluctuation component
$\delta f_e$
becomes non-zero. According to (4.5), this implies a mismatch between the two definitions of temperature, so that
$T_{ST,e} \neq \tilde {T}_e$
. This difference is clearly visible in the right panels of figure 4.
4.4.2. Role of the self-electrostatic interaction
Figure 5 displays the temperature profiles (solid lines) and number density profiles (dashed lines) obtained from numerical simulations. The blue curves correspond to electrons, while the red curves represent protons. In the left panel, the temperatures are calculated using (4.4), while in the right panel, they are calculated using (4.3). The theoretical temperature and density profiles are plotted as thick grey lines: the decreasing density profile is given by (3.17), and the increasing temperature profile by (3.18). The simulation parameters were fixed as follows:
$\tau = 0.5$
,
$t_w = 25$
,
$\tilde {g} = 1$
,
$M = 1836$
,
$C = 400$
,
$\Delta T = 90$
,
$N = 2^{16}$
and the integration time step was set to
$\text{d}t = 4 \times 10^{-3}$
. As evident from figure 5, the coarse-grained temperature and density profiles no longer agree with the analytical predictions of (3.17) and (3.18), for both species. Moreover, the electron profiles, computed using both (4.4) and (4.3), no longer exhibit a monotonic increase in temperature with respect to the spatial coordinate
$\theta$
. These results suggest that, in this regime, the electrostatic interaction plays a non-negligible role in shaping the temperature and density profiles of the plasma. The breakdown of the agreement with the analytical profiles is a direct consequence of the presence of non-negligible fluctuations, whose effects will be analysed in detail later in this section.

Figure 5. Left panel: coarse-grained density (dotted lines) and temperature (solid lines) profiles for protons (red) and electrons (blue), calculated using (4.4). Right panel: same quantities, but temperature profiles are computed using (4.3). In both panels, grey curves show theoretical predictions from (3.17) (density) and (3.18) (temperature). Green curves represent proton-only theoretical profiles from (4.24) (density) and (4.25) (temperature), which account for the electrostatic contribution.
Since
$C \neq 0$
, the ‘coarse-graining’ correction terms described by (3.10) are non-zero. As a result, the system of kinetic equations can no longer be described solely by the mean-field Vlasov dynamics, represented by (3.9). Consequently, the distribution functions given by (3.16) are no longer valid solutions of the full kinetic equations. To assess the presence and magnitude of fluctuations, we make use of the two-fluid formalism introduced in § 3.3. In particular, the momentum-balance equations (3.23) are central to this analysis. Due to the assumption of central symmetry, all odd moments of the coarse-grained distribution functions vanish. This property has been confirmed via numerical simulations. As a result, only the pressure gradient terms contribute to the left-hand sides of (3.23). To quantify the intensity of the fluctuations
$\delta f_{\alpha }$
, we compute the residual terms
$R_{\alpha }$
appearing in (3.23). This can be done by evaluating, from simulations, the pressure term as

and the total force term as

Therefore, the difference of (4.15) and (4.16) gives the ‘coarse-graining’ force


Figure 6. Top left: electron pressure gradient
$\mathcal{F}_{e,1}$
(red), total force
$\mathcal{F}_{e,2}$
(blue) and coarse-graining correction
$\mathcal{F}_{e,3}$
(green), as functions of the spatial coordinate
$\theta$
. The black curve shows the total force balance, which vanishes as expected. Top right: same quantities for protons. Bottom left: electrostatic force (green), gravitational force (blue), and total force (red) acting on electrons, calculated respectively from (4.20) to (4.22). Bottom right: same as left, but for protons.
In figure 6, we plot the three force components
$\mathcal{F}_{\alpha ,j}$
for
$j = 1,2, 3$
as functions of
$\theta$
: the top left panel corresponds to electrons, and the top right to protons. For electrons, it is evident that the pressure gradient force
$\mathcal{F}_{\alpha ,1}$
(shown in red) differs from the mean-field force
$\mathcal{F}_{\alpha ,2}$
(blue curve). As a result, the coarse-graining correction
$\mathcal{F}_{\alpha ,3}$
(green) is non-zero and of comparable magnitude to the other contributions, indicating the presence of significant fluctuations in this regime. In contrast, for protons, since
$\tau , \langle t_w \rangle \ll t_{R,p}$
, the system remains effectively in the stationary-state regime. Accordingly, we observe that
$\mathcal{F}_{\alpha ,1} = \mathcal{F}_{\alpha ,2}$
, implying that the coarse-graining term
$\mathcal{F}_{\alpha ,3}$
vanishes. The sum of all three force components, shown in black in both panels, satisfies the force balance condition imposed by (3.23), and thus must be identically zero.
Now that the intensity of the fluctuations has been quantified, we shall return to the kinetic approach. Using (3.5) and imposing stationarity, we obtain the following kinetic equation:

for the electrons and

for the protons. In (4.18) and (4.19), the coarse-graining correction term due to fluctuations has been retained for electrons but omitted for protons. As a result, the electron temperature and density profiles cannot coincide with the analytical mean-field predictions given by (3.18) for temperature and (3.17) for density. This discrepancy arises because the electron dynamics is governed by the full kinetic equation, including the fluctuation term in (4.18). For protons, by contrast, the kinetic equation (4.19) includes only the Vlasov term. However, unlike in the purely gravitational case described by (3.16), the electrostatic field is now non-zero. This is clearly visible in the two bottom panels of figure 6, where we plot the individual force components acting on each species as functions of the spatial coordinate
$\theta$
. Specifically, in these panels, we show:
-
(i) in green,
(4.20)\begin{equation} F_{\alpha ,1} = \mathrm{sign}(e_{\alpha }) , C \boldsymbol{\cdot }Q \boldsymbol{\cdot }\sin {\theta }; \end{equation}
-
(ii) the external gravitational force in blue,
(4.21).\begin{equation} F_{\alpha ,2} = \tilde {g} , \sin \left (\frac {\theta }{2}\right ); \end{equation}
-
(iii) the total force acting on a particle of species
$\alpha$ in red,
(4.22)\begin{equation} F_{\alpha ,3} = F_{\alpha ,1} + F_{\alpha ,2}. \end{equation}
As can be seen, the electrostatic force is of the same order of magnitude as the gravitational force for both electrons and protons, and thus plays a significant role in the system’s dynamics. Consequently, the density and temperature profiles obtained from simulations cannot match the analytical predictions of (3.17) and (3.18). However, an improved analytical expression can still be derived in this regime by incorporating the electrostatic potential into the single-particle Hamiltonians, as in (2.11). Replacing
$\tilde {H}_\alpha$
with the full Hamiltonian
$H_i$
in (3.16), we obtain the corrected coarse-grained distribution function

From this distribution, we compute the number density

and the kinetic temperature

Here, the total potential
$V$
acting on protons is

The density (a decreasing function) and temperature (an increasing function) profiles given by (4.24) and (4.25) are plotted in green in figure 5. These show good agreement with the numerical results, demonstrating that the electrostatic field significantly influences proton dynamics. Furthermore, since the distribution functions in (4.23) remain suprathermal and are governed by a mean-field Vlasov dynamics, the observed temperature inversion continues to be explained by the velocity filtration mechanism. Finally, even when the temperature is computed using the alternative definition in (4.3), a clear distinction between electron and proton behaviour remains, as shown in figure 5. Nevertheless, the proton temperature profile computed via (4.3) coincides with that obtained using (4.4). This agreement confirms that proton fluctuations are negligible, and according to (4.5), the two temperature definitions must converge. For electrons, by contrast, fluctuations are non-negligible, and the two definitions yield distinct temperature profiles. Moreover, since the dynamics for electrons are no longer governed by a purely Vlasov-type evolution, temperature inversion via gravitational filtering is no longer guaranteed across the entire spatial domain. This breakdown is clearly visible in the simulation results presented in the right-hand panels of figure 5.

Figure 7. Top panel: time series of the kinetic energies
$K_{\alpha }$
for protons (orange) and electrons (blue), computed numerically using (4.2). Bottom panels: spatial profiles of the coarse-grained density and temperature for protons (dashed red) and electrons (solid blue), evaluated following the procedure outlined in § 4.1. In the left panel, temperature profiles are computed using (4.4); in the right panel, using (4.3). Theoretical predictions from (3.17) (for density) and (3.18) (for temperature) are shown in grey. Dashed yellow curves correspond to the analytical profiles obtained by superposing multiple thermal configurations, as described in § 3.4. Specifically, densities are computed using (3.28) (bottom panels), and temperatures using (3.30) (left) and (3.31) (right).
4.5. Superposition regimes
By increasing the time scales
$\tau$
and
$t_w$
of the temperature fluctuations, the electrons will first reach a configuration of oscillation in time between various thermal equilibria, and then the protons will do so until the following regime is reached,

In this case, both the kinetic energies of protons and electrons oscillate around different thermal solutions, as illustrated in figure 7. In this regime, both species are not in a stationary state, but since we are interested in observing the dynamics on a coarse-grained time scale given by (3.1), then the dynamics becomes stationary at such a coarse graining time scale. All parameters have been fixed at the following numerical values:
$\tau =400, t_w = 4000, \tilde {g}=1, M=100, C=400, \Delta T=90, N=2^{16}$
and integration time step
$\text{d}t=4\times 10^{-3}$
. As illustrated in figure 7, the coarse-grained numerical temperature profiles of both protons (in red) and electrons (in blue) calculated with both (4.4) and (4.3) exhibit significant discrepancies from the analytical solution (in grey) provided by (3.18) and they yield the profiles obtained from the superposition of disparate thermal solutions (in yellow), specifically calculated with (4.4) for the left panel and with (4.3) for the right panel. The same considerations made in (4.3) apply to the densities. In this instance, and in accordance with the theory outlined in § 3.4, the results are independent on the value of electrostatic coupling constant
$C$
. We also note that in this regime, a temperature inversion does or does not occur depending on how the coarse-grained temperature is calculated. In particular, if the temperature is calculated using the definition given by (4.3), then during the dynamics, the various temperatures are added together, resulting in a total temperature that is the superposition of the temperatures of many thermal states, as given by (3.31). Conversely, if the temperature is calculated using the definition given by (4.3), it corresponds to the temperature associated with the temporal coarse graining distribution function, which in this case is the superposition of many thermal distributions and is described by (3.27). As previously outlined in § 3.4, the distribution functions in question exhibit suprathermal tails. Given that it is determined at each point by the conservation of energy of a single particle,
$\tilde {H}_{\alpha }$
, gravitational filtering leads to the emergence of a monotonically increasing temperature profile that is anti-correlated to the density profile.
5. Summary and perspectives
In this work, we have extended the analysis of the dynamics of a stratified, collisionless plasma gravitationally confined in a loop and in thermal contact with a fluctuating boundary beyond the previously studied regime where the fluctuation time scales are much smaller than the relaxation times of electrons and protons. Building upon previous investigations by Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024a, Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casettib), we have developed a generalised temporal coarse-graining framework that remains valid when the characteristic time scales of the temperature fluctuations are no longer much shorter than the electron crossing time. By deriving a set of coarse-grained kinetic equations, we have identified a novel contribution, arising directly from the coarse-graining procedure that affects the stationary profiles of the plasma. Our numerical results show that when the time scale of temperature fluctuations becomes comparable with the electron crossing time, this additional term leads to a species-dependent separation in the temperature and density profiles. In this regime, the stationary state of the system is no longer described by the effective distribution function introduced by Barbieri et al. (Reference Barbieri, Papini, Di, Pierfrancesco, Simone and Casetti2024b) and gravitational filtering is no longer sufficient to guarantee temperature inversion. This effect is interpreted through the emergence of electrostatic fluctuations that generate a correction to the standard Vlasov dynamics. We have also investigated the opposite limit, in which the time scales of temperature fluctuations are much longer than the crossing times of both species. In this regime, the system evolves through a sequence of thermal states, and its long-term behaviour can be described by a coarse-grained distribution that is a superposition of thermal solutions corresponding to different boundary temperatures.
We have also introduced two definitions of temperature at the temporal coarse-graining level. The first is defined as the ratio between the coarse-grained pressure and the coarse-grained density. The second is obtained by computing the instantaneous kinetic temperature at each time step and then averaging over time. The first definition assumes that measurements of the moments can be performed with a resolution much shorter than the coarse-graining time scale, thereby allowing access to the system’s microscopic dynamics. In contrast, the second definition is appropriate when the system dynamics below the coarse-graining scale
$\tilde {t}$
are not accessible. In this regime, the behaviour of the two species is described by the coarse-grained distribution functions
$\tilde {f}_{\alpha }$
, and the kinetic temperature is computed from these via the standard kinetic definition. The two definitions yield the same result only when fluctuations are negligible, that is, in the stationary regime. In other regimes, particularly in the superposition regime, a temperature inversion persists if the temperature is computed as the ratio of coarse-grained pressure to coarse-grained density. However, if the coarse-grained temperature is instead defined as the average of the instantaneous kinetic temperatures, the resulting profile becomes isothermal. Such an inversion is expected to be observed when the integration time of the measurement significantly exceeds the relaxation time.
As an initial avenue for future investigation, it would be worthwhile to extend the HMF approximation for the electrostatic field by incorporating multiple Fourier modes in the expansion. Similarly, the temporal coarse-graining procedure could be generalised to account for a multi-mode electrostatic interaction. This extension would allow us to address several key questions. Does the overall dynamics depend on the specific form of the electrostatic interaction? And, more importantly, does the plasma always reach the stationary regime described in § 4.3? This latter question is intriguing not only from a theoretical perspective but also for its implications in solar physics. As highlighted earlier, our model applied to the solar atmosphere (see, e.g. Barbieri et al. Reference Barbieri, Casetti, Verdini and Landi2024a) demonstrates that rapid, intense, intermittent and short-lived heating events can drive the plasma towards a stationary configuration characterised by anti-correlated density and temperature profiles similar to those observed in the Sun. An affirmative answer to this question would provide strong evidence of the robustness of our results. Another promising extension is the inclusion of a magnetic field along the loop axis in the plasma model. Such an addition is of interest not only for modelling solar coronal loops, but also in the context of fusion plasmas confined in Tokamak devices (see, e.g. Goedbloed et al. Reference Goedbloed, Keppens and Poedts2010; Ciraolo et al. Reference Ciraolo2018). Introducing an axial magnetic field, which would naturally decrease in strength from the base to the top of the loop due to expansion with altitude, would impose the additional constraint of magnetic moment conservation on the particle motion. This conservation would tend to increase the parallel temperature while decreasing the perpendicular temperature, thereby enhancing gravitational filtering along the magnetic field and diminishing it in the orthogonal direction. With this extra ingredient, several new questions arise. How does the temperature inversion change in the steady-state configuration? For what values of the loop expansion rate does the perpendicular temperature still increase monotonically with height? And how does the overall plasma dynamic picture change?
Finally, it is important to note that throughout this discussion, the plasma has been treated as collisionless. The framework presented in this paper may be affected by the introduction of collisions in the first two regimes, but not in the superposition regime described in § 4.5, since in that case, the system oscillates between various thermal configurations that are also solutions of the Boltzmann equation. In the context of solar coronal applications, the introduction of Coulomb collisions becomes particularly relevant. The main questions to be addressed are the following. Beyond the superposition regime, how does plasma dynamics change with the inclusion of collisions? Is it always possible to identify a stationary state regime comparable to that described in § 4.3? If so, how do the properties of this regime evolve under collisional effects? Under what conditions do temperature inversion and non-thermal distributions persist in the presence of Coulomb collisions?
Acknowledgements
The authors wish to thank the anonymous reviewers for their insightful comments, which improved the presentation of the work. L.B. wishes to thank Arnaud Zaslavsky for valuable discussions.
Editor Francesco Califano thanks the referees for their advice in evaluating this article.
Funding
We acknowledge the Fondazione CR Firenze under the projects HIPERCRHEL. This research was partially funded by the European Union – Next Generation EU – National Recovery and Resilience Plan (NRRP) – M4C2 Investment 1.4 – Research Programme CN00000013 ‘National Centre for HPC, Big Data and Quantum Computing’ – CUP B83C22002830001 and by the European Union – Next Generation EU – National Recovery and Resilience Plan (NRRP) – M4C2 Investment 1.1 – PRIN 2022 (D.D. 104 del 2/2/2022) – Project ‘Modeling Interplanetary Coronal Mass Ejections’, MUR code 31. 2022M5TKR2, CUP B53D23004860006. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. L.B. wants to thank the Sorbonne Université in the framework of the Initiative Physique des Infinis for financial support.
Declaration of interests
The authors report no conflict of interest.