1. Introduction
Stationary states out of thermal equilibrium occur in nature and often exhibit unexpected properties. An example is given by the so-called inverted temperature–density profile, or temperature inversion, that corresponds to an increase in temperature while the density is decreasing. Numerical simulations have shown how temperature inversion can be achieved in non-equilibrium stationary states, usually referred to as quasi-stationary states, attained after processes of violent relaxation (see e.g. Lynden-Bell (Reference Lynden-Bell1967), Kandrup (Reference Kandrup1998), Levin, Pakter & Teles (Reference Levin, Pakter and Teles2008) and Ewart et al. (Reference Ewart, Brown, Adkins and Schekochihin2022), see also Barbieri et al. (Reference Barbieri, Di Cintio, Giachetti, Simon-Petit and Casetti2022) and references therein) of many-body systems with long-range inter-particle potentials (see Di Cintio, Ciotti & Nipoti (Reference Di Cintio, Ciotti and Nipoti2013) and Campa et al. (Reference Campa, Dauxois, Fanelli and Ruffo2014) and references therein) brought out of thermal equilibrium by energy injection (Casetti & Gupta Reference Casetti and Gupta2014; Teles et al. Reference Teles, Gupta, Di Cintio and Casetti2015; Gupta & Casetti Reference Gupta and Casetti2016).
Temperature inversion occurs in several astrophysical systems such as, for example, filaments in molecular clouds (Arzoumanian et al. Reference Arzoumanian, André, Didelon, Konyves, Schneider, Menoshchikov, Sousbie, Zavagno, Bontemps and Di Francesco2011; Toci & Galli Reference Toci and Galli2015; Di Cintio, Gupta & Casetti Reference Di Cintio, Gupta and Casetti2018) the Io plasma torus around Jupiter (Meyer-Vernet, Moncuquet & Hoang Reference Meyer-Vernet, Moncuquet and Hoang1995) and the hot gas in galaxy clusters (Wise, McNamara & Murray Reference Wise, McNamara and Murray2004; Baldi et al. Reference Baldi, Ettori, Mazzotta, Tozzi and Borgani2007). The most relevant and widely studied case is represented by the atmosphere of the Sun (Golub & Pasachoff Reference Golub and Pasachoff2009). The problem of temperature inversion in the atmosphere of the Sun has been approached in several ways (see Klimchuk (Reference Klimchuk2006) and Parnell & De Moortel (Reference Parnell and De Moortel2012) and references therein) but remains far from being completely solved.
Among all the existing theoretical models, the one that is relevant for the present work is the kinetic approach introduced by Scudder (Reference Scudder1992a,Reference Scudderb) and dubbed by the author ‘velocity filtration’. In this approach the solar atmosphere is treated as a system of non-interacting particles under the action of the Sun's gravity. Thanks to the conservation of energy, only particles with a sufficient amount of kinetic energy can climb the gravity well and reach a given altitude, thus forming the atmosphere (for this reason the mechanism was dubbed ‘velocity filtration’). If the particle's velocity distribution functions (VDFs) at the base of the atmosphere are thermal, then the stationary configuration is the standard isothermal stratified atmosphere (Aschwanden Reference Aschwanden2005), with a density profile exponentially decreasing with altitude and a constant temperature all over the atmosphere. If, as in the original Scudder model, the velocity distribution functions are suprathermal (e.g. they are given by a kappa distribution, see Lazar & Fichtner Reference Lazar and Fichtner2021), then the temperature profile is a growing function of the atmosphere's height, thus resulting in stationary anti-correlated density and temperature profiles. In the Scudder formalism, however, the chromospheric VDFs are assumed fixed and suprathermal, even if the chromosphere is highly collisional, posing the problem of how such distributions can be sustained in such an environment. Moreover, modelling the chromosphere with a fixed kappa function (as in Scudder Reference Scudder1992a,Reference Scudderb) produces a linear increase of temperature rather than a transition region and corona, as observed in the solar atmosphere. To overcome the difficulty of sustaining suprathermal populations in the strongly collisional chromosphere, mechanisms related to turbulence and wave–particle interactions (e.g. Parker & Tidman Reference Parker and Tidman1958) or statistical processes based on Levy's flights (Collier Reference Collier1993) have been often evoked. In the present work, as well as in Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024), we want to overcome this difficulty by modelling the chromosphere as fully collisional so that it can be schematised at every instant of time as a thermal boundary. We then ask if rapid heating events, modelled as sudden temperature increments of the thermal boundary, are able to produce suprathermal distribution functions. More precisely, in Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024) it has been shown that rapid, intermittent and short-lived heating events are able to drive the coronal collisionless plasma towards a non-equilibrium stationary state characterised by inverted temperature–density profiles, with a transition region and a corona similar to those observed in the solar atmosphere. As pointed out in Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024), the idea behind studying the role of temperature fluctuations in the chromosphere lies in the fact that, while the VDFs of the chromospheric plasma are likely to be thermal due to collisionality, the chromosphere is a very dynamic environment, showing fine-scale structures down to instrumental resolution (see Molnar et al. Reference Molnar, Reardon, Chai, Gary, Uitenbroek, Cauzzi and Cranmer2019; Ermolli et al. Reference Ermolli, Giorgi, Murabito, Stangalini, Guido, Molinaro, Romano, Guglielmino, Viavattene and Cauzzi2022), hence its temperature is expected to fluctuate in space and time (Hansteen et al. Reference Hansteen, De Pontieu, Carlsson, Lemen, Title, Boerner, Hurlburt, Tarbell, Wuelser and Pereira2014; Peter et al. Reference Peter, Tian, Curdt, Schmit, Innes, De Pontieu, Lemen, Title, Boerner and Hurlburt2014).
In this paper we investigate this matter further, by analysing the problem in detail from the point of view of kinetic plasma physics. We introduce a novel temporal coarse-graining technique to model the physics of the non-equilibrium plasma dynamics described above and in Barbieri et al. (Reference Barbieri, Casetti, Verdini and Landi2024). Here, we show how the dynamics can be described by a set of two coupled Vlasov equations for the temporal coarse-graining version of the distribution functions of the two species (here electrons and protons). This system of equations is coupled to effective coarse-grained distributions at its boundary describing the temperature fluctuations of the chromosphere. We stress the fact that, at variance with Scudder (Reference Scudder1992a,Reference Scudderb), we do not consider non-interacting particles, but rather we explicitly take into account their mutual electrostatic interaction, although only at the mean-field level. We derive analytical expressions for the coarse-grained distribution functions of the two species in the non-equilibrium stationary configuration. These distribution functions exhibit suprathermal high-energy tails because of the dynamics induced by the temperature fluctuations of the thermal boundary, that is, they are self-consistently produced by our modelling of the chromosphere. Of course, extra parameters are introduced to specify the distribution of temperature fluctuations, but no fine tuning of these parameters is necessary to produce the temperature inversion and to allow the use of a coarse-graining approach. The only requirement is a fluctuation time scale that is smaller than the relaxation time scale of the system. As a consequence, our coarse-graining formalism lets us interpret our results in terms of Scudder's velocity filtration mechanism, but without imposing a priori suprathermal tails. We note that suprathermal tails in the distribution functions can be obtained even in stationary states of an isolated collisionless plasma as a consequence of the dynamical phase mixing in the single-particle phase space (see Ewart, Nastac & Schekochihin (Reference Ewart, Nastac and Schekochihin2023) and references therein).
The paper is structured as follows. In § 2 we present the model used to describe the collisionless plasma atmosphere in contact with the fluctuating thermal boundary. In § 3 we present the temporal coarse-graining approach describing the plasma dynamics and we obtain the coarse-grained version of the distribution functions of the two species in the non-equilibrium stationary configuration. In § 4 we test the theory against kinetic numerical simulations and we discuss its limits of applicability. In§ 5 we apply our method to the specific case of the atmosphere of the Sun and discuss the results, also comparing our findings with the observed density and temperature profiles of Yang et al. (Reference Yang, Zhang, Jin, Li and Duan2009). Moreover, we discuss some observational evidence to support our result and the physical limits of modelling the coronal plasma neglecting Coulomb collisions. Finally, in § 6 we summarise and discuss the future perspectives of this study.
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:
(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).
(ii) The symmetry is cylindrical along the loop-following coordinate $x$.
(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).
(iv) Every observable is symmetric with respect to the top of the loopFootnote 1. 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 2.
The first two assumptions imply that the total energy of the system is
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
where the charge density
is related to the potential $\phi$ by the Poisson equation
Here, $U_g$ is the total external potential per section and reads
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:
We then define the charge imbalances as
where the vectors $\boldsymbol {q}_{n,\alpha }$ are given by
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
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.
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
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
where the electric field $E$ decomposed in Fourier modes is
We now implement the HMF assumption by truncating the Fourier expansion at the first mode. The equations above become
and the electric field becomes
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
where
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:
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
where the expressions for the external forces and the electrostatic field are
and the charge imbalance and stratification parameters are given by
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
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:
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
where the (2.19a,b) become
and (2.19a,b) become
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
where $f_\alpha$ are the single-particle distribution functions for both species, and $H_{\alpha }$ are the mean-field Hamiltonians
where
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:
where $f_{\alpha }$ is normalised to 1 for both species. By combining the above equation with (2.28a,b), we obtain
and integration with respect to $p$ gives
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
where $\tilde {H}_{\alpha }$ are the mean-field Hamiltonians (2.27) with $\phi =0$, that is
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
and for the kinetic temperatures
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 3 , 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
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
The temporal coarse-grained phase-space distribution function, defined as the time average of $f_{\alpha }$ over $\tilde {t}$, reads as
and the time-averaged Vlasov dynamics over $\tilde {t}$ becomes
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
The thermal boundary at the bottom boundary of the system induces an incoming flux defined at a given time by
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
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
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
so that the coarse-grained flux becomes
where $\tilde {f}_{\alpha }(p)$ is defined by
and $A$ is
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
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
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
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
while for the kinetic temperatures
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
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
together with their theoretical value (indicated by the straight solid lines) in the stationary configuration, given by
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$.
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.
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.
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
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.
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
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)$.
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
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$.
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).
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:
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).
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).
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.
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.
5.1. Observational evidence of temperature increments and physical limits of our modelling
We now discuss the limits of our modelling approach. In our model, we only considered one-sided temperature fluctuations, based on the assumption that a physical process heats up the chromosphere for a finite amount of time. Such an assumption justifies the use of positive temperature increments only. Observational evidence of physical processes that can produce such increments will be presented later in this section. We stress the fact that, even if we decrease the temperature below the background value, we still have a temperature inversion. In this case, the ‘hot’ particles are the ones produced by the background temperature and the ‘cold’ ones are produced by the temperature decrements. Therefore, velocity filtration is still active, and the temperature starts from a value that is smaller compared with the background due to the colder particles but rises to the background value.
In conclusion, the temperature fluctuations drive the above collisionless plasma environments towards a non-equilibrium stationary configuration with temperature inversion that can be described by our temporal coarse-graining method.
This is true only if the time scales of the temperature fluctuations are much smaller than the relaxation time $t_R$ (computed via (3.2), roughly 10 s for the parameters that we have considered above).
Moreover, to reach coronal temperatures of $10^6$ K and maintain, at the same time, the temperature at the base of the transition region at around $10^4$ K, the fraction of time in which the chromospheric temperature increments must be as large as one million degrees must be around $A \sim 0.02$. More precisely, with this value of $A$, the temperature at the base is around $T_b \sim 1.2 \times 10^4$ this is compatible with ALMA observations of temperature $1.2 \times 10^4$ K (Molnar et al. Reference Molnar, Reardon, Chai, Gary, Uitenbroek, Cauzzi and Cranmer2019).
Indeed, short-lived and intense brightenings are observed on the Sun's surface (Dere, Bartoe & Brueckner Reference Dere, Bartoe and Brueckner1989; Teriaca et al. Reference Teriaca, Banerjee, Falchi, Doyle and Madjarska2004; Peter et al. Reference Peter, Tian, Curdt, Schmit, Innes, De Pontieu, Lemen, Title, Boerner and Hurlburt2014; Tiwari et al. Reference Tiwari, Panesar, Moore, De Pontieu, Winebarger, Golub, Savage, Rachmeler, Kobayashi and Testa2019; Berghmans et al. Reference Berghmans, Auchère, Long, Soubrié, Mierla, Zhukov, Schühle, Antolin, Harra and Parenti2021). Among all of these the so-called campfires, observed in extreme UV imaging, have temperatures of $\approx 10^6$ K. Other explosive events appearing in $\mathrm {H}\alpha$ line widths have smaller temperatures, around $2 \times 10^5$ K, but are ten times more frequent (Teriaca et al. Reference Teriaca, Banerjee, Falchi, Doyle and Madjarska2004). This trend justifies the use of a decreasing exponential distribution for the temperature increments. However, current measurements have a time resolution of the order of a few seconds (i.e. of the order of the relaxation time $t_R$) so that, even if the intensity of such increments could be compatible with recently observed explosive events, they remain unresolved because of their rapid time duration. Interestingly, recent EUV solar observations (Raouafi et al. Reference Raouafi, Stenborg, Seaton, Wang, Wang, DeForest, Bale, Drake, Uritsky and Karpen2023) show that a possible physical mechanism that is at work in the lower atmosphere is magnetic reconnection. In particular, they show how jetting activity at the base of the solar corona driven by discrete small-scale magnetic reconnection events (i.e. continuous sources of temperature increments) can produce a flow of matter that propagates from the corona up to the heliosphere so that they could be at the origin of the coronal heating and of the solar wind acceleration.
The most questionable aspect of our treatment is that the chromosphere is modelled as fully collisional and concurrently the overlying plasma (the transition region and the corona) is modelled as collisionless. Due to the low values of the Knudsen number $K_n$ ($K_n\sim 10^{-2}\unicode{x2013}10^{-3}$ in the transition region and $K_n\sim 10^{-1}$ in the corona) the collisions are certainly present in both the transition region and the corona. Clearly, a fully self-consistent treatment of Coulomb collisions in the whole solar atmosphere is beyond the scope of this work. However, it is worth briefly addressing this problem. Thanks to the $1/r$ (where $r$ is the inter-particle distance) nature of electrostatic potential, the electron–electron mean free path in a plasma strongly depends on the particles velocity $v$, more specifically it is proportional to $v^4$. In our procedure we have ‘cold’ particles generated by the background temperature $T_0$ and ‘hot’ particles generated by the temperature increments. We expect that only ‘cold’ particles would be strongly affected by collisions so that the VDFs would be closer to thermal at low energies, while non-thermal features associated with the ‘hot’ particles, which are the ones selected by gravitational filtering, would be able to reach the corona (Landi & Pantellini Reference Landi and Pantellini2001). Moreover, velocity filtration should become more efficient in presence of Coulomb collisions, thus producing a transition region even steeper than the one produced by our model and closer to the observed one.
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).
Acknowledgements
We wish to thank G. Cauzzi for very useful discussions.
Editor Francesco Califano thanks the referees for their advice in evaluating this article.
Funding
We acknowledge partial financial support from the MIUR-PRIN2017 project Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST) n. 201798CZL, the National Recovery and Resilience Plan, Mission 4 Component 2 – Investment 1.4 - NATIONAL CENTER FOR HPC, BIG DATA AND QUANTUM COMPUTING – funded by the European Union – NextGenerationEU – CUP B83C22002830001, the Solar Orbiter/Metis program supported by the Italian Space Agency (ASI) under the contracts to the National Institute of Astrophysics (INAF), Agreement ASI-INAF N.2018-30-HH.0, and the Fondazione CR Firenze under the projects HIPERCRHEL and THE SWITCH. 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.
Declaration of interests
The authors report no conflict of interest.