Temperature inversion in a gravitationally bound plasma: Case of the solar corona

The temperature of the solar atmosphere increases from thousands to millions of degrees moving from the lower layer (chromosphere) to the outermost one (corona), while the density drops accordingly. The mechanism behind this phenomenon, known as a temperature inversion, is still unknown. In this work, we model a coronal loop as a collisionless plasma confined in a semicircular tube that is subject to the Sun's gravity and in thermal contact with a fully collisional chromosphere behaving as a thermostat at the loop's feet. By using kinetic $N$-particle simulations and analytical calculations, we show that rapid, intermittent, and short-lived heating events in the chromosphere drive the coronal plasma towards a non-equilibrium stationary state. The latter is characterized by suprathermal tails in the particles' velocity distribution functions, exhibiting temperature and density profiles strikingly similar to those observed in the atmosphere of the Sun. These results suggest that a million-Kelvin solar corona can be produced without the local deposition of heat in the upper layer of the atmosphere that is typically assumed by standard approaches. We find that suprathermal distribution functions in the corona are self-consistently produced instead of postulated a priori, in contrast to classical kinetic models based on a velocity filtration mechanism.


Introduction
Temperature inversion, namely, an anticorrelation between temperature and density profiles, occurs in astrophysical systems, such as the filaments in molecular clouds (Arzoumanian et al. 2011;Toci & Galli 2015), the Io plasma torus around Jupiter (Meyer-Vernet et al. 1995), the Earth's magnetosphere (Ma et al. 2020), the hot gas in some galaxy clusters (Baldi et al. 2007;Wise et al. 2004) and, notably, the solar atmosphere, which is with no doubt the most striking and thoroughly studied example (see e.g., Golub & Pasachoff 2009, and references therein).The Sun's outermost layer, the corona, reaches temperatures above 10 6 K.It lies on top of a denser and cooler (roughly 10 4 K) layer, the chromosphere, with the two layers connected by the so-called transition region, a thin interface only hundreds of kilometers wide that sees temperature variations up to a factor of 50 and density drops by the same factor (Yang et al. 2009).The mechanism behind this "coronal heating problem" is still largely unknown (Klimchuk 2006;Parnell & De Moortel 2012).
Assuming local thermodynamic equilibrium, a corona can form only upon the local deposition of heat in the upper layer.The energy coming from the Sun is enough to heat the whole corona (Withbroe 1981), bringing on the question of how energy is transported and dissipated at coronal heights.Some proposed mechanisms involve the release of magnetic energy stored in the corona (Parker 1972;Dmitruk & Gomez 1997;Gudiksen & Nordlund 2005;Rappazzo et al. 2008;Rappazzo & Parker 2013;Wilmot-Smith 2015), transmission and damping of waves generated by photospheric motions (Heyvaerts & Priest 1983;Ionson 1978;Howson et al. 2020), or a direct insertion of hot plasma in form of chromospheric spicules (Pontieu et al. 2011).However, there are theoretical indications (Shoub 1983;Esser & Edgar 2000) and observational evidence (Dudík et al. 2017) to support that local thermodynamic equilibrium may not be satisfied in the transition region and in the corona.This allows for another way of obtaining temperature inversion, as first recognized by Scudder (1992a,b).If the velocity distribution functions (VDFs) of electrons and ions have suprathermal tails 1 in the chromosphere, then temperature must increase with height: faster particles are able to climb higher in the Sun's gravity well and temperature increases thanks to "velocity filtration."Notably, this implies no local deposition of heat in the upper layer.Unfortunately, the model does not produce a transition region and, more importantly, non-thermal distributions in the strongly collisional chromospheric plasma are difficult to justify.
Interestingly, recent numerical studies have shown that when an isolated system governed by long-range interactions is impulsively perturbed it undergoes collisionless relaxation, ultimately reaching a non-thermal stationary state with temperature inversion, whose origins can be traced back to velocity filtration (Casetti & Gupta 2014;Teles et al. 2015;Gupta & Casetti 2016).This might explain the inverted temperature-density profiles of filaments in molecular clouds (Di Cintio et al. 2017), but it cannot be applied as such to the coronal plasma, which is not isolated, but rather engaged in steady thermal contact with the chromosphere.However, 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 (Cauzzi et al. 2009;Ermolli et al. 2022), while its temperature is expected to fluctuate in space and time (Peter et al. 2014;Hansteen et al. 2014).
In this Letter, using both numerical simulations and analytical modelling, we show that rapid temperature fluctuations in the chromosphere are able to build up a hot corona.For simplicity, we modeled a coronal loop as a semicircular tube of collisionless plasma in thermal contact with a thermostat, the latter playing the role of the dynamic and collisional chromospheric plasma.The tunable part of the model is the statistics of the temperature fluctuations of the chromosphere.Temperature inversion appears as a very robust phenomenon and does not require any finetuning.However, matching the density and temperature values of the solar corona, while obtaining a (thick) transition region, requires a specific range of temperature fluctuations.Here, we briefly discuss whether the latter may be observationally tested.

Model and numerical simulations
We consider a model of the plasma in a coronal loop made up of N e = N i = 2N protons and electrons having masses m i and m e and charge e i = −e e = e.Particles are confined in a semicircular tube of length, 2L, and cross-section, S , and they are subject to a constant downward gravity and to an electric force that ensures charge neutrality, whose combined effect is proportional to (m e + m i )/2 (Pannekoek 1922;Rosseland 1924;Belmont et al. 2013).We assume that all quantities depend only on the coordinate along the loop axis, x ∈ [−L, L] and are symmetric with respect to the loop top located at x = 0. We expand the electrostatic interactions between the particles in a Fourier series, with wave vectors k n = πnx/L, only up to n = 1.This is a meanfield approximation, where each particle only experiences the average effect of all the others, analogous to that yielding the Hamiltonian mean field model for single-component long-range-interacting particle systems (Antoni & Ruffo 1995;Chavanis et al. 2005;Giachetti & Casetti 2019).With these assumptions, the equations of motion are: where j = 1, . . ., 2N numbers the particles, α = e, i denotes the species, g = GM /R 2 is the Sun's gravity, and is the self-consistent electric field 2 , with In Eq. ( 3), q α is the "stratification parameter" for each species, obtained by averaging the cosine of the particle angular position along the loop.As a result, q α = 0 corresponds to a uniform distribution of particles, q α = −1 (resp.+1) to a distribution concentrated at the base of the loop, in x = ±L (resp.at the top, in 2 We use electrostatic cgs units. x = 0).When q α ∈ (−1, 0) (resp.(0, 1)) the density decreases (resp.increases) from the bottom to the top of the loop.The difference q i − q e in Eq. ( 2) measures the charge imbalance giving rise to the electric field.The anchoring of the loop feet to the high chromosphere is modeled by coupling the coronal plasma to a thermostat at x = ±L: when a particle reaches the boundary, it is reinjected into the loop with a random velocity consistent with the flux from a thermal distribution at the temperature of the thermostat (Landi & Pantellini 2001;Lepri et al. 2021;Tehver et al. 1998).Once written using dimensionless variables, the equations of motion depend only on three parameters: the mass ratio, M = m i /m e , and the strengths of the interactions, C and g, expressed as the electrostatic and gravitational energy, respectively, normalized to the thermal energy, where n is the average number density of each species, k B is the Boltzmann constant, and T 0 = 10 4 K is the reference temperature of the thermostat and our unit of temperature.
Starting with the plasma in thermal equilibrium at the initial thermostat's temperature, T = 1 (in our units), the latter is allowed to vary so as to mimic the dynamic nature of the chromosphere.After a waiting time, t w , we increase the thermostat temperature by a quantity ∆T , we keep the thermostat at T = 1 + ∆T for a time, τ, and then we switch back to T = 1.These steps are iterated for the whole duration of the simulation run, drawing the values of t w and ∆T from two exponential distributions, The physical idea is that the chromospheric temperature fluctuates due to random heating events, among which those producing a larger ∆T are less likely to occur.We observe that when both t w and τ are sufficiently small and, in particular, when these values are smaller than the relaxation times to equilibrium, the plasma in the loop is kept away from the initial thermal state, never relaxes to a thermal distribution with T > 1; rather, it attains a non-equilibrium stationary state, with non-thermal distribution functions and inverted temperature-density profiles.We solved Eqs. ( 1)-( 3), coupled to the fluctuating thermostat, using a fourth-order symplectic algorithm (Candy & Rozmus 1991) with a time step of δt = 10 −4 .Due to the meanfield nature of the particles' equations of motion, the computational cost scales as N (instead of N 2 as in direct N-body codes) and in the limit of large N, the dynamics is collisionless and equivalent to a Vlasov one (Campa et al. 2014).

Simulations results
We now focus on the case of a coronal loop with average number density, n = 2.5 × 10 9 cm −3 , and half-length, L = π × 10 4 km, so that g = 16.64.For the mass ratio we chose the realistic value M = 1836.A realistic value of the parameter C ≈ 10 22 implies very fast plasma oscillations and prohibitively small integration time steps (≈1/ √ C) to reliably simulate the dynamics.However, we checked that in the stationary state the temperature and density profiles as well as the VDFs do not depend on the value of C, which only affects the amplitude and frequency of the oscillations around the means.Therefore, we arbitrarily chose C = 400 so that we may have oscillations with a periodicity greater than L5, page 2 of 6 the thermostat timescale (see below), meaning that our simulated particles have a much smaller charge than in the real world.The simulations were performed using N ≈ 2.1 × 10 6 , τ = 8 × 10 −4 , t w = 4 × 10 −2 , and T p = 90, i.e., 9 × 10 5 K in physical units.
Together with the choice of T 0 = 10 4 K, these values imply an average temperature T b ≈ 1.2 × 10 4 K at the base of the loop (see below).This sets the base at a height z b ≈ 2 × 10 3 km above the chromosphere, where the transition region begins, and the top of the loop at z t ≈ 2.2 × 10 4 km, well into the corona.
In the top panel of Fig. 1 we report the time evolution of the kinetic energies per particle where p j,α is the dimensionless momentum of the jth particle of species α, with M e = 1 and M i = M; in the bottom panel the stratification parameters q e and q i defined in Eq. ( 3) are reported, again as a function of time.After a transient, a stationary state is reached where all these quantities fluctuate around constant values, the transient being longer and the amplitude of the fluctuations smaller for the ions, as expected due to their larger mass.The initial value of q i = q e ≈ −0.99 corresponds to the stratification of an isothermal atmosphere at T = 1.In the stationary state, which is not a thermal equilibrium state (as we show in the following), a smaller value is obtained, q i ≈ q e ≈ −0.94, indicating a slightly milder stratification.In Fig. 2, we report the temperature and density profiles of electrons along the loop, obtained by time-averaging in the stationary state 3 .To make the comparison with the known profiles of the solar atmosphere easier, we report the results in physical units as a function of the height, z, above the chromosphere (temperatures in K, densities in cm −3 , height in km).A temperature inversion is clearly visible, as well as the presence of a transition region: indeed, there is a sharp increase (resp.decrease) of the temperature (resp.density) roughly between z = 2 × 10 3 and z = 5 × 10 3 km, with a jump from 1.2 × 10 4 to 6 × 10 5 K while density drops by two orders of magnitude, followed by a gentler change of both quantities.Both the values of T and n and the shape of the curves are strikingly similar to those of the Sun atmosphere (see e.g., Golub & Pasachoff 2009;Yang et al. 2009), aside from the fact that the transition region in our model is larger than in the real case, that is, roughly 3 × 10 3 km against less than 2 × 10 2 km.

Theory
The numerical results can be reproduced and explained computing the distribution functions (DFs) f α (ϑ, p) of ions and electrons in the non-equilibrium stationary state (ϑ is the dimensionless position).To this end, we first consider the time averaged incoming flux of energy at the feet of the loop, i.e., at ϑ = ±π, which can be written for each species as J IN,α = J T (t),α t , where we have  Fig. 2. Electron density (red) and temperature (blue) as a function of the height above the photosphere, obtained in a numerical simulation with parameters corresponding to the solar atmosphere (see text).The thick grey lines are the theoretical profiles computed using Eq. ( 11).Ion profiles (not shown) are the same.
over the probability distributions of waiting times and of thermostat temperatures, obtaining where A = τ/ τ + t w η and • γ and • η stand for the averages over the distributions γ and η that are given in Eq. ( 5).In our model the plasma in the loop is collisionless, so that f α obeys the Vlasov equation and in the stationary state it depends on ϑ and p only via the single-particle Hamiltonian, as implied by the Jeans theorem.We also note that an electrostatic interaction term proportional to the effective electric field should be present in Eq. ( 9), but it can be proven that it vanishes in the stationary state Barbieri et al., in prep.The DF at ϑ = ±π L5, page 3 of 6 is fixed imposing equality of the incoming and outgoing fluxes at the feet of the loop, that is, J IN,α = J OUT,α , with The DFs f α (ϑ, p) in the stationary state can be obtained using the distribution at the boundary combined with the Jeans theorem Barbieri et al. (in prep.) because the level sets of the singleparticle Hamiltonian ( 9) are always connected to the bottom boundary (i.e., the thermostat), expressed as where The interpretation of Eq. ( 11) is as follows: the stationary DF is given by a thermal distribution at temperature T = 1, plus a nonthermal contribution arising from the average of thermal distributions at T = 1 + ζ over the probability distribution γ(ζ) of the temperature fluctuations.The weight of the non-thermal contribution is proportional to A, the fraction of time in which the thermostat is not at temperature T = 1.The thermal population dominates at small heights z, and is depressed by the gravity term in H α when increasing z; conversely, the non-thermal contribution becomes more and more relevant at larger z due to velocity filtration, because faster particles can climb higher in the potential well, showing up as suprathermal tails in the distribution.This is shown in Fig. 3, where the VDFs of electrons normalized by the density are plotted at three increasing heights, z = 2.3, 3.9, 11×10 3 km from bottom to top, corresponding to the base of the transition region, the middle transition region, and the corona, respectively.Red, blue, and green curves are obtained from the simulation, while grey curves are the theoretical predictions of Eq. ( 11).We note that VDFs are plotted as functions of the signed kinetic energy sign(p) p 2 /2 in a semilogarithmic scale, so that a thermal distribution (a Gaussian) appears as a triangle symmetric about zero.The VDFs are always composed of a thermal core at small velocities plus a suprathermal tail at larger velocities.As the height increases, the thermal core progressively shrinks and almost disappears in the corona, where the VDF is basically suprathermal.These VDFs explain the shape of the temperature and density profiles reported in Fig. 2, where the grey curves are the theoretical prediction obtained with Eq. ( 11).The sharp temperature rise in the transition region (and density decrease) is caused by the dramatic change in the shape of the VDF that passes from being nearly thermal to almost completely suprathermal.Once coronal heights are reached, totally non-thermal VDFs produce a gentler variation of temperature and density, similar to the case described in Scudder (1992a,b).For completeness, the stationary values of kinetic energies and stratification parameters computed using Eq. ( 11) are drawn with black horizontal lines in Fig. 1.

Discussion and conclusions
We present a kinetic model of the solar atmosphere, where the collisonless coronal plasma is in steady contact with a thermostat mimicking a completely collisional chromosphere.The analytical and numerical results consistently show that in response to intermittent rapid and short-lived increments of the chromospheric temperature, suprathermal tails in the VDFs naturally form and gravitational filtering causes a sharp temperature rise and density decrease in the above atmosphere, consisting of a (thicker than observed) transition region, followed by an extended corona at roughly 10 6 K. Suprathermal electrons, along with a thermal population, are measured in situ in the solar wind (Pilipp et al. 1987;Halekas et al. 2020;Maksimovic et al. 2020).Their presence in the transition region (Dudík et al. 2017) or in flaring regions (Polito et al. 2018) is also compatible with remote sensing observations of non-thermal line widths, but their direct detection is still challenging.Our model supports the formation of a nonthermal population already at the base of the corona for closed loop geometry.We expect our mechanism to be valid also in the case of open field lines, thus favouring the formation of the solar wind by gravitational filtering as in exospheric models (Jockers 1970;Lemaire & Scherer 1971;Lamy et al. 2003;Maksimovic et al. 1997;Zouganelis et al. 2004).However, in an open spherical geometry, there are also escaping particles, forming the wind, and "trapped particles" (whose trajectories never reach the base); thus, a more complex formalism would be necessary to tackle this case (Zouganelis et al. 2004).
One aspect that is fundamental in achieving temperature inversion is that the coronal plasma is in a non-equilibrium stationary state.This requires the chromosphere to maintain any given temperature for intervals much shorter than the relaxation time of electrons in the corona, t R , which in our model, is expected to be the shortest between the sound travel time and the free fall time, that is, roughly 10 s for the parameters we considered above.To reach temperatures around 10 6 K in the corona, while keeping temperatures around 10 4 K at the base of the transition region, the mean of the chromospheric L5, page 4 of 6 temperature increments must be as large as the coronal temperature, namely, T p ≈ 10 6 K, but the ratio between their duration, τ, and the typical waiting time, t w , between them must be small, so that τ t w t R .Short-lived, intense, and small-scale brightenings are routinely observed on the Sun (Dere et al. 1989;Teriaca et al. 2004;Peter et al. 2014;Tiwari et al. 2019;Berghmans et al. 2021).Among the latter, the so-called campfires, recently observed in extreme UV imaging by Solar Orbiter, have temperatures of ≈10 6 K, while explosive events appearing in Hα line widths have smaller temperatures, around 2 × 10 5 K (Teriaca et al. 2004), but are ten times more frequent.This trend is compatible with the exponential distribution of temperature increments we assumed.However, current measurements have a temporal resolution which is at best of a few seconds, so that temperature increments able to form a corona according to our model remain unresolved.Anyhow, even in the case of unresolved events, a given temperature measurement brings a footprint of the underlying event distribution.For instance, assuming t w ≈ 1 s, with τ ≈ 0.1 s our model would predict an average temperature T ≈ 2 × 10 4 K at the top of the chromosphere, which is not observed; however, if τ ≈ 1/50 s we get T ≈ 1.2 × 10 4 K, as in the case we have shown in Fig. 2.This is consistent with ALMA observations of brightness temperatures as large as 1.2 × 10 4 K (Molnar et al. 2019), with measurements at a 2 s cadence.A possible physical mechanism at work is magnetic reconnection in the low atmosphere, as suggested by recent observations (Lee et al. 2020;Raouafi et al. 2023).
The main result in this Letter is that the average temperature profile of the solar atmosphere can be obtained keeping the coronal plasma in thermal contact with the chromosphere, without any local deposition of heat in the corona.At variance with previous results based on velocity filtration, a transition region appears, and non-thermal chromospheric VDFs are not needed: a corona can be formed out of a collisional, albeit dynamic, chromosphere.Clearly, this does not exclude that also other mechanisms contribute to coronal heating.
Potentially important physics is neglected in the model presented here, such as radiation losses, collisions in the corona, and magnetic fields.Although a detailed discussion of these effects is beyond the scope of this work, we address them in brief as follows.In our model, energy is transported in the corona by the particles' stream, and is exchanged with the thermostat on a timescale of the order of the crossing time of the loop, which is less than a minute and much smaller than the radiation timescale (around 30 m, see e.g., Serio et al. 1991;Shibata & Magara 2011;Antolin & Froment 2022), so that energy losses due to radiation can be safely neglected in our model.Moreover, radiation losses would not change the picture, because in the non-equilibrium state, a heat flux towards the radiating region would set in and restore the temperature profile (Shoub 1983;Landi & Pantellini 2001;Dorelli & Scudder 1999).Collisions between particles will surely occur in the corona.In any case, thanks to the fact that the mean free path of a particle in a plasma scales with velocity v as v 4 , we expect that only "cold" particles would be heavily affected by collisions; specifically, the VDFs would become closer to thermal at small energies, while non-thermal features related to the "hot" particles, which are the ones selected by gravitational filtering and are then able to reach coronal heights, would not be erased by collisions (Shoub 1983;Landi & Pantellini 2001).Velocity filtration is then expected to become more efficient in presence of some degree of collisionality; the latter effect, together with particle confinement by magnetic lines, may result in a transi-tion region sharper than the one obtained here and closer to the observed one.
The mechanism that ends up producing a temperature inversion in the solar corona in our model is general and it could prove relevant to other systems where hot coronae are thought to be present: stars other than the Sun, as well as active galactic nuclei and so on.

Fig. 1 .
Fig. 1.Kinetic energies (top) and stratification parameters (bottom) of ions and electrons in green and orange, respectively, as a function of time.The black horizontal lines are the theoretical stationary values predicted using Eq.(11).

Fig. 3 .
Fig.3.Electron VDFs (color curves) normalized by the electron densities as a function of the signed kinetic energy, at different heights (see labels).Theoretical VDFs computed from Eq. (11) are plotted as thick grey lines.Ion VDFs (not shown) are the same if plotted against sign(p)p 2 /2M.In the bottom panel, a magnification of the central region of the same VDFs is shown to highlight the disappearance of the Gaussian profile.