© N. Brughmans et al. 2022

1. Introduction

Prominences are cool, dense plasma structures in the solar corona, suspended by the magnetic field. They are usually found above polarity inversion lines (PILs) – features found in magnetogram observations where the normal component of the photospheric magnetic field changes sign (Mackay et al. 2010; Vial & Engvold 2015; Gibson 2018). In recent years, many high resolution simulations have been performed, both in 2D and 3D, that have investigated prominence formation through various mechanisms. These mechanisms usually require the magnetic field topology local to the condensations to be concave upwards with respect to gravity; levitation against gravity is facilitated by means of the magnetic Lorentz force (pressure and tension), as first shown by the self-similar model of Kippenhahn & Schlüter (1957). Generally speaking, two main formation processes are discerned: (i) the evaporation-condensation model, where localised footpoint heating of a magnetic arcade produces upflows which collect in any pre-existing dipped fields (Xia et al. 2012, 2014; Xia & Keppens 2016) – a theory shared with the phenomenon of coronal rain (Li et al. 2022); and (ii) the levitation-condensation model (Kaneko & Yokoyama 2015, 2018; Jenkins & Keppens 2021) or its variant reconnection-condensation (Kaneko & Yokoyama 2017), where a flux rope is dynamically formed through shearing and converging motions at lower altitudes – a sheared magnetic arcade, whose footpoints are subject to those motions, eventually changes internal connectivity via magnetic reconnection and builds a twisted, helical field (van Ballegooijen & Martens 1989). In doing so, the reconnection may redistribute or scoop up material from the chromosphere or lower corona, which remains suspended by the flux rope (levitation). The condensation aspect of all three formation mechanisms relies on a runaway cooling effect triggered by perturbations to the energy loss rate, in other words density and temperature variations, known as the thermal instability (Parker 1953; Field 1965).

The thermal instability is of paramount importance in many fields of astrophysics. The process leads to the formation of filamentary structures, providing an excellent alternative to the gravitational instability in the interstellar medium (Jennings & Li 2021) or galaxy clusters (Qiu et al. 2020). Numerical implementations commonly assume a purely optically thin ambient medium, such as is the case for the solar corona, where energy transported as radiation is able to ‘free-stream’ outwards. Since perturbations to density and temperature can lead to local increases in these radiative losses, variations feed back into the energy balance so as to yet further increase the density and decrease the temperature. The criteria that govern the (in)stability of this delicate equilibrium were first formulated by Parker (1953) and fully established through the seminal work of Field (1965). The latter author obtained an analytic dispersion relation in an idealised setting that governs the linear behaviour of the mode. Recently, Claes & Keppens (2019) and Claes et al. (2020) revisited these results and numerically studied a local coronal volume pervaded by interacting slow magneto-acoustic waves in 2D and 3D settings, respectively. The authors showed that the thermal mode eventually takes over, with flows along the magnetic field that drive material towards a thermally unstable location, leading to the formation of a filamentary structure initially oriented almost perpendicular to the magnetic field. Subsequent small differences in ram pressure from both sides force the redistribution of these structures along the magnetic field; in 3D, an intriguing misalignment was found between the condensations and the magnetic field. How different cooling curves impact the above process was investigated by Hermans & Keppens (2021), who introduced a bootstrap measure such that extremely thin filaments could be handled at unprecedented numerical resolutions. Further evidence for the importance of thermal instability within the solar atmosphere was recently provided by Claes & Keppens (2021), who examined the full magnetohydrodynamic (MHD) spectrum of a realistic stratified and magnetised atmosphere. The study highlighted that much of the solar chromosphere, transition region and corona can be liable to unstable or overstable magneto-thermal modes.

The recent study of Jenkins & Keppens (2021) considered the formation and evolution of 2.5D flux rope – prominence systems via the thermal instability as a component of the common levitation-condensation mechanism (Kaneko & Yokoyama 2015). Their solar coronal model adopted an exponential background heating profile meant to balance the local radiative and conductive losses of the system governed mechanically by hydrostatic equilibrium (Fang et al. 2013). Once the flux rope – prominence system formed, the simulation yielded high-resolution demonstrations of the dynamic formation of prominence condensations via the thermal instability. The authors clarified how the fully non-linear evolution can be related to sequential linear MHD mode stability considerations, with the thermal mode driving the condensation process. In their 2.5D simulations, the field-aligned thermal conduction term within their energy equation in essence thermally isolated the flux rope plasma from the background corona, but the time-independent, exponential background heating component still supplied a small but non-negligible contribution to this plasma at all times. This directly influences the thermal balance and in turn artificially lengthens the formation timescales of thermal instability (Hermans & Keppens 2021). This aspect is of critical importance when considering the dynamic formation, evolution, and decay processes in prominences, for instance the timescales of prominence mass recycling, as recently highlighted by the pioneering work of Kaneko & Yokoyama (2018). It should be stated that thermal instability in the context of solar prominence formation also depends on the presence of the chromosphere and transition region (Klimchuk 2019). Through thermal conduction, much of the energy of the solar corona is transported downwards, which can provide a more efficient mechanism for energy loss than through radiative losses alone. Moreover, the inclusion of the extra mass reservoir opens the aforementioned evaporation-condensation evolutionary pathways to modifying the eventual properties of the prominence material.

The thermal instability relies on the intricate balance between energy loss, gain, and transport mechanisms. Within the solar corona this balance is generally regulated by ambient heating, downwards thermal conduction and optically thin radiative cooling, with additional considerations for compression and Joule heating, for instance. However, as the cause of the 1 MK solar corona remains a long-standing open problem, often ad hoc background heating models are artificially imposed. Two classes of artificial heating models are generally discerned (Mandrini et al. 2000). The first class describes the heating rate as a function of height in the corona based on observations or theoretical models, typically formulated assuming an exponential decrease with height as in Jenkins & Keppens (2021) (see, e.g., Sturrock et al. 1996; for the observational basis, but also Fang et al. 2013; Xia & Keppens 2016; Zhao et al. 2017 and Fan 2017 for similar numerical implementations). The second class considers physically or observationally derived scaling laws inside coronal loops (e.g., Dahlburg et al. 2018), but the models can be applied more generally to the solar atmosphere. These scaling laws provide a parameterised heating rate based on local plasma conditions,

$$\begin{array}{c}\hfill \mathcal{H}\sim B^{\alpha }{\rho }^{\beta }{L}^{\gamma },\end{array}$$(1)

where B is the magnetic field strength, ρ the density, and L the length of the field line through the local volume. This heating model is characterised by the combination of powers (α, β, γ). A comprehensive summary of several known possible mechanisms behind these scaling laws is given in Table 5 of Mandrini et al. (2000), ranging from reconnection (DC models) to wave dissipation (AC models), with dependencies that sometimes include additional parameters beyond the three given above.

Each of these two classes (exponential or parameterised through Eq. (1)) of heating prescriptions has been successfully employed in simulations of both coronal loops (Mok et al. 2008, 2016) and prominences (Kaneko & Yokoyama 2015, 2017, 2018; Xia et al. 2012, 2014). Incorporating field line length dependence is computationally demanding, as it requires tracing field lines for every cell in the simulation domain during runtime. This is non-trivial to realise in a grid-adaptive, parallelised framework like MPI-AMRVAC (Keppens et al. 2021), although recent work addressing solar flares with electron beam ingredients along evolving reconnecting field lines used this in Ruan et al. (2021). Mok et al. (2016) provide an alternative in the context of coronal loop simulations through approximating the field line length L by the radius of curvature of the magnetic field R, which is the inverse of the magnetic curvature κ. This first order approximation is valid when the loops are almost circular. In the context of prominences, however, this approach proves unable to capture the complex twisted field that comprises a flux rope.

In this work, we implement a wide array of heating prescriptions and compare their influence on prominence formation through the 2.5D levitation-condensation models of Jenkins & Keppens (2021). For Eq. (1), a comparison between two such heating models, (α, β) = (2, 0) and (α, β) = (0, 1), with γ = 0, was previously presented by Kaneko & Yokoyama (2015) as they introduced the levitation-condensation model. They found that for those heating models that depend on magnetic field strength, condensations appear only when the flux rope is formed through anti-shearing motions, where the arcade footpoint motions are directed so as to decrease the magnetic shear of the initial condition. We extend these results here by looking at the morphology and physical properties of the formed prominences, doing so for each of the aforementioned heating prescriptions and senses of shear, in addition to a new ‘reduced heating’ model that furthermore approximates the influence of the L parameter for flux ropes. Moreover, a phase space distribution of the condensation process reveals a non-isobaric evolution of thermal instability, where constant pressure is eventually reached along the flux surfaces that contain the prominence.

In Sect. 2, we detail the simulation setup used throughout this work, together with the heating models that we considered and our method for flux rope tracking. In Sect. 3, we list the results from our simulations, which we further discuss in Sect. 4.

2. Numerical setup and equations

As already indicated, solar prominences are phenomena found embedded within the solar corona. Following the original levitation-condensation setup of Kaneko & Yokoyama (2015), we restrict our region of interest to the coronal volume and neglect the chromosphere below. We consider a 2D simulation domain in the x − y-plane (horizontal–vertical), with vector components in the invariant z-direction (horizontal) to complete the 2.5D description. To carry out the simulation, we make use of the fully open-source, adaptive-grid, parallelised MPI-AMRVAC toolkit¹ (Keppens et al. 2021). We simulate a cross-section of 24 × 25 Mm with x ∈ [ − 12, 12] Mm and y ∈ [0, 25] Mm, taking a base resolution of 96 × 96 cells with three additional levels of Adaptive Mesh Refinement (AMR). This yields a maximum resolution of approximately 24 Mm /(96 × 2³) = 31.25 km, which corresponds to the lower-resolution simulations of Jenkins & Keppens (2021). The AMR criteria are based on sharp gradients in the local density and the magnetic field components following the Löhner prescription (Löhner 1987).

2.1. Equations and physics

The full set of MHD equations solved by MPI-AMRVAC, including all non-ideal or non-adiabatic effects and source terms, are equivalent to (the code actually solves the related system using conservative variables),

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathit{v})=0,\hfill \end{array}$$(2)

$$\begin{array}{cc}& \rho \frac{\partial \mathit{v}}{\partial t}+\rho \mathit{v}·\mathrm{\nabla }\mathit{v}+\mathrm{\nabla }p-\mathit{j}\times \mathbf{B}-\rho \mathbf{g}=0,\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathbf{B}}{\partial t}-\mathrm{\nabla }\times (\mathit{v}\times \mathbf{B})+\mathrm{\nabla }\times \eta \mathit{j}=0,\hfill \end{array}$$(4)

$$\begin{array}{cc}& \rho \frac{\partial T}{\partial t}+\rho \mathit{v}·\mathrm{\nabla }T+(\gamma -1)[p\mathrm{\nabla }·\mathit{v}+\rho \mathcal{L}-{\eta |\mathit{j}|}^2-\mathrm{\nabla }·(\stackrel{\rightleftarrows }{\kappa }·\mathrm{\nabla }T)]=0,\hfill \end{array}$$(5)

as given in Goedbloed et al. (2019)². The thermodynamic variables are related through an equation of state, for which we use the ideal gas law,

$$\begin{array}{c}\hfill p\mu =\mathcal{R}\rho T.\end{array}$$(6)

Here, μ is the mean molecular mass in proton masses m_p and ℛ ≈ 3.81 J K⁻¹ is the ideal gas constant. We assume a monatomic gas so that the ratio of specific heats γ = 5/3. Governing the magnetic field are Gauss’ law for magnetism and Ampère’s law,

$$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{B}=0,\mathrm{\nabla }\times \mathbf{B}={\mu }_0\mathit{j},\end{array}$$(7)

where μ₀ = 4π is the permeability of vacuum. The (dependent) variables appearing in this closed system of equations, are the density ρ, velocity v, gas pressure p, current density j, magnetic field B and temperature T.

The temperature Eq. (5) contains energy sinks and sources within the square brackets. Firstly, the heat-loss function ρℒ = n_Hn_HeΛ(T)−ℋ contains the net volumetric energy loss from optically thin radiative cooling n_Hn_HeΛ(T) and heating ℋ, with units ergs cm⁻³ s⁻¹ (Field 1965). We note that the cooling rate depends on the hydrogen and helium number densities n_H, He and temperature T through a cooling curve Λ, constructed using tabulated energy losses. The volumetric heating rate ℋ – although currently unknown – must necessarily depend on a number of quantities, such as density ρ or magnetic field strength B, and we consider a few possibilities throughout this work (Mandrini et al. 2000). Secondly, Ohmic heating ∼η|j|² is included in line with the consideration of magnetic resistivity. The coefficient η has an artificially high value of 2.33 × 10⁸ m² s⁻¹(=0.002 in code units), ensuring that we dominate any possible numerical diffusion, and allowing the reconnection necessary during the early formation phase of the flux rope. Thirdly, anisotropic thermal conduction $\mathrm{\nabla }·(\stackrel{\rightleftarrows }{\kappa }·\mathrm{\nabla }T)$ is included, where for the parallel component we take the Spitzer conductivity κ_∥ = 8 × 10⁻⁷T^5/2 ergs cm⁻¹ s⁻¹ K⁻¹ (Spitzer 2006). Due to being eight orders of magnitude smaller than the parallel component under coronal conditions, perpendicular thermal conduction is neglected throughout this study.

The plasma is assumed to be completely ionised (an assumption that may be invalid inside cool and dense prominences Hillier 2019; Braileanu et al. 2021a,b), and is partly composed of helium ions with assumed abundance n_He = 0.1n_H. The mean molecular mass is hence,

$$\begin{array}{c}\hfill \mu =\frac{\rho /m_{\mathrm{p}}}{n}=\frac{n_{\text{H}}+4n_{\text{He}}}{2n_{\text{H}}+3n_{\text{He}}}\approx 0.6.\end{array}$$(8)

For consistency with earlier work (Jenkins & Keppens 2021), we adopt a non-constant gravitational acceleration,

$$\begin{array}{c}\hfill \mathit{g}(y)=-g_{\text{cor}}\frac{R_{\odot }^2}{{(R_{\odot }+y)}^2}\widehat{\mathit{y}},\end{array}$$(9)

with y the height above the solar surface and g_cor = 274 m s⁻² the gravitational acceleration at the solar surface, located at radius R_⊙ = 695 700 km. In our simulation domain, the gravitational acceleration falls from g_cor at the bottom to 0.93g_cor at the top (y = 25 Mm).

The MHD Eqs. (2)–(5) are solved using a threestep, third-order Runge–Kutta method with the HLL flux scheme (Harten et al. 1983). To limit spurious oscillations, we use the third-order asymmetric CADA3 slope limiter proposed by Čada & Torrilhon (2009). As a divergence fix for the magnetic field, we use the linde method to perform parabolic cleaning (Keppens et al. 2003) and reduce monopole errors through the bottom boundary driving where a second order central differences approximation for ∇ ⋅ B = 0 is implemented. MPI-AMRVAC allows for a split treatment of the magnetic field, fixing a steady background field B₀ while performing calculations for the perturbed field B₁ as described by Tanaka (1994), which is also possible for force-free fields (Xia et al. 2018). Magnetic field splitting ensures accuracy and efficiency by better capturing behaviour at low plasma-beta. We also solve an auxiliary energy equation for the internal energy density so as to handle typical coronal low-beta regions and replace faulty pressure values when necessary. For the radiative cooling physics, we use the Colgan_DM cooling curve by Colgan et al. (2008), modified for lower temperatures by Dalgarno & McCray (1972), which is suitable for this study (Hermans & Keppens 2021). The cooling curve is interpolated at 12 000 points and calculated at runtime using the exact integration method (Townsend 2009). A minimal allowed temperature of 10³ K is enforced throughout the simulations, well below the typical minimum temperatures found in our simulations (∼7000 K). MPI-AMRVAC uses normalisation of the plasma variables, which is fixed by the following choice: number density ${\stackrel{\sim }{n}}_{\text{H}}={10}^8$ cm⁻³, length $\stackrel{\sim }{L}={10}^8$ cm and temperature $\stackrel{\sim }{T}={10}^6$ K.

2.2. Simulation setup

As initial conditions, we set up a stratified solar coronal atmosphere in hydrostatic and thermal near-equilibrium. Assuming an isothermal corona at a temperature of T₀ = 1 MK and assuming a plane-parallel atmospheric stratification, the momentum Eq. (3) reduces to a separable ODE if the magnetic field is force-free. The resulting equation is solved numerically with the trapezoid rule in our initial conditions, which is numerically convenient and equivalent to setting the analytic hydrostatic solution. In hydrostatic equilibrium, it follows from the ideal gas law Eq. (6) that pressure and density feature the same exponential variation with height,

$$\begin{array}{c}\hfill p_0(y)=p_{\text{cor}}exp(-\frac{y}{H(y)}),{\rho }_0(y)={\rho }_{\text{cor}}exp(-\frac{y}{H(y)}),\end{array}$$(10)

where,

$$\begin{array}{c}\hfill H(y)=H_0\frac{R_{\odot }+y}{R_{\odot }},\text{with}{H}_0=\frac{\mathcal{R}{T}_0}{\mu g_{\text{cor}}},\end{array}$$(11)

is a local scale height based on the hydrostatic pressure scale height H₀ ≈ 50 Mm and p_cor = 0.434 dyn cm⁻² and ρ_cor = 3.2 × 10⁻¹⁵ g cm⁻³ are the pressure and density values at the base of the corona, respectively.

The magnetic configuration within the simulation starts with the linear force-free sheared magnetic arcade B₀ given by Kaneko & Yokoyama (2015, 2017, 2018) and Jenkins & Keppens (2021). Then the magnetic field strength varies as,

$$\begin{array}{c}\hfill B_0(y)=B_{a}exp\left(\frac{-y}{H_0}\right).\end{array}$$(12)

Here, B_a = 10 G is the field strength at the bottom of the domain. The magnetic arcade has exactly one half period in the simulation domain with the PIL at x = 0 Mm, setting an initial shear angle of ≈9° with the out-of-plane normal. Our choice of B_a corresponds to the strongest magnetic field considered by Jenkins & Keppens (2021). The field strength, gas pressure and density share approximately the same height profile such that magnetic pressure (∼B²) decreases twice as fast. The initial sheared arcade prescription is taken to be the steady background field B₀ used in the aforementioned magnetic field splitting. An image of this initial setup can be found in Jenkins & Keppens (2021, Fig. 1).

2.3. Boundary conditions

We take the same boundary conditions described in Jenkins & Keppens (2021), where there is no outflow out of the domain and the magnetic field is vertical (in line with the periodicity) at the left and right edges. At the lower and upper boundary, the magnetic field is extrapolated. To form the flux rope, we impose the shearing and converging motions at the bottom boundary as originally described in Kaneko & Yokoyama (2015). Following these authors, there is a distinction to be made between shearing and anti-shearing motions, which increase and decrease the magnitude of the magnetic field component in the invariant direction, respectively. The motion is prescribed as,

$$\begin{array}{c}\hfill V_z(x,t)=\{\begin{array}{cc}-V_x(x,t),\hfill & \text{shearing,}\hfill \\ V_x(x,t),\hfill & \text{anti-shearing},\hfill \end{array}\end{array}$$(13)

where V_x(x, t) is prescribed following Eqs. (13)–(17) of Jenkins & Keppens (2021); the associated 12 km s⁻¹ magnitude for driving and shearing motions is then far below the characteristic coronal Alfvén speed of ≈442 km s⁻¹. We adopt the driving prescription in time and space from this work but do not redefine the moment they are initiated to t = 0. The driving hence begins at t = 1300 s and lasts until t = 2800 s.

2.4. Heating models

We now include an overview of the heating models adopted in this work, and describe how detailed balance can be achieved against the initial radiative losses. Since these initial losses scale as ρ₀(y)²Λ(T₀), they are explicitly known at t = 0 s using Eq. (10) and T₀ = 1 MK. For the required initial thermal equilibrium, any heating model has to fulfill the condition that there are no initial net energy losses or gains:

$$\begin{array}{c}\hfill 0={\rho }_0{\mathcal{L}}_0={\rho }_0^2\mathrm{\Lambda }(T_0)-{\mathcal{H}}_0.\end{array}$$(14)

2.4.1. Exponential and mixed heating models

The first heating model is time independent and is hence completely determined by condition (14), as in for instance Claes & Keppens (2019) and Hermans & Keppens (2021). The heating rate then falls exponentially with height, since it is the absolute value of the initial cooling rate. We henceforth refer to this model as the ‘exponential heating model’, or model E. The heating prescription, which has been successfully implemented by many previous authors (Fang et al. 2013; Xia & Keppens 2016; Zhao et al. 2017; Jenkins & Keppens 2021) then becomes:

$$\begin{array}{c}\hfill \mathcal{H}={\rho }_{\mathrm{cor}}^2\mathrm{\Lambda }(T_0)exp(-\frac{2y}{H(y)}).\end{array}$$(15)

A second class of heating models is based on those heating scaling laws deduced for coronal loops. The ‘mixed heating model’, or model M, is of the general form,

$$\begin{array}{c}\hfill \mathcal{H}=c{B}^{\alpha }{\rho }^{\beta },\end{array}$$(16)

where (α, β)∈ℝ² and c(y) is a function determined by the equilibrium condition (14). Instead of using Eq. (12) for the initial magnetic field strength, we approximate it with the exponential profile with varying scale height for convenience in the expression below. This sets up an approximate initial equilibrium which relaxes during the first 1300 s, with the initial background heating given by,

$$\begin{array}{c}\hfill {\mathcal{H}}_0=c(y)B_a^{\alpha }exp(-\frac{\alpha y}{H(y)}){\rho }_{\mathrm{cor}}^{\beta }exp(-\frac{\beta y}{H(y)}).\end{array}$$(17)

The function c(y) can then be found as,

$$\begin{array}{c}\hfill c(y)=\frac{{\rho }_{\mathrm{cor}}^{2-\beta }}{B_a^{\alpha }}\mathrm{\Lambda }(T_0)exp(-\frac{(2-\alpha -\beta )y}{H(y)}),\end{array}$$(18)

where the adopted approximate magnetic scale height allows the exponentials to be combined. We note that choosing α + β = 2 makes c a constant. This is preferable both from a theoretical and practical point of view: if α + β ≠ 2 then there is still a steady, exponential component to the heating, losing the truly local nature of the heating model and maintaining the potentially unphysical residual heating to the flux rope. Hence, all the mixed heating models considered in this work are chosen to satisfy α + β = 2, as in the model (α, β) = (2.5, −0.5) used by Mok et al. (2008).

Indeed, heating prescriptions depending on local parameters have been used before, most notably by Mok et al. (2008) and Kaneko & Yokoyama (2017) in its two-parameter form above. Several other combinations of parameters (α, β) have been used throughout the literature, most often with β = 0 so as to obtain a heating rate that depends purely on the local magnetic energy density, hereafter model M0 (magnetic heating) (Kaneko & Yokoyama 2015, 2018; Xia et al. 2012, 2014). We investigate the influence of models using a small, but increasingly, negative β on the prominence formation, which we term models M1 (quasi-magnetic heating, β = −0.1), M2 (β = −0.2), and M3 (Mok et al. 2008, β = −0.5) – an overview of the models used in this work and their prescription is given in Table 1. With the inverse density dependence, the heating prescription is meant to be magnetic energy-based, but with an additional ‘penalty’ for regions of increased density. When performing prominence simulations, this is a desirable property since then the heating rate will decrease in any locations subject to a density enhancement. We note that the previous exponential heating model E is equivalent to the mixed heating model with choice (α, β) = (0, 0).

2.4.2. Reduced heating model

Finally, we propose a new heating model pertaining to 2.5D prominence formation simulations that effectively approximates the true 3D nature of a flux rope. We want to account for the ignored variation in z, and note that a flux rope within the solar atmosphere is typically characterised by extended field lines of some hundred Mm. Kaneko & Yokoyama (2017) demonstrated that inside a flux rope, longer field lines preferentially cool down and host condensations notwithstanding a heating prescription that did not explicitly trace field lines. By assuming the action of heating instead scales exponentially along field lines rather than simply with height (many 1D hydro models of prominence formation exist which inherently support the relevance of the field line length to the condensation process; see, e.g., Karpen & Antiochos 2008; Klimchuk & Luna 2019; Pelouze et al. 2022), we may argue that the flux rope should experience a reduced heating rate in its central portion compared to the footpoint regions. In the 2.5D MHD scenario studied here, we thereby want to incorporate how field lines in 3D have varying lengths throughout our flux rope cross-section.

To this aim, we apply a time-dependent heating reduction that follows the flux rope cross-sectional shape throughout the simulations. This requires a dynamic detection of the flux rope during runtime, for which we propose a method that relies on the 2.5D local magnetic field curvature (i.e., curvature defined using the 3D magnetic field vector but omitting gradients in the z-direction). Such a local prescription removes the need for tracing the global field associated with the flux rope. This resulting ‘reduced heating model’ is applied to both the exponential heating model (now termed RE) and the mixed model (RM) so as to approximate the thus far ignored dependence of the general heating model (1) on field line length L inside the flux rope.

The flux rope cross-section at time t is modelled as an ellipse centred at (x_O(t),y_O(t)), coinciding with the central flux rope axis, which appears as a magnetic null point (O-point) in the 2D projection. The ellipse’s major (a(t)) and minor (b(t)) axes are given as the distance from the centre to the vertical and horizontal edges of the flux rope, respectively. After these parameters have been determined in the current timestep (see Appendix A), the volumetric background heating rate ℋ_b is transformed following,

$$\begin{array}{c}\hfill \mathcal{H}=f(\mathit{u})\times {\mathcal{H}}_b,f(\mathit{u})=1-\delta {\mathrm{e}}^{-{|\mathit{u}|}^8+{|\mathit{u}|}^4}/1.284,\end{array}$$(19)

where,

$$\begin{array}{c}\hfill \mathit{u}=(\frac{x-x_{\mathrm{O}}}{0.85b},\frac{y-y_{\mathrm{O}}}{0.85a}),\end{array}$$(20)

is a coordinate transformation centred on the flux rope axis and scaled according to its dimensions. The reduction function f is shown in Fig. 1 and takes values between [1 − δ, 1] since the exponential is normalised by its maximum of 1.284. To ensure that f does not influence the regions external to the flux rope, we scale a and b by 0.85 in the coordinate transformation. The above form for the reduction function was chosen so as to reflect the structure of a theoretical twisted flux rope, that is, the field lines at the centre of the flux rope are less twisted and hence shorter than those field lines that form the ‘boundary’ with the corona (cf. Titov & Démoulin 1999). For a heating rate that varies exponentially along field lines, this leads to modest ≈60% reduction along the central axis of the flux rope that increases radially away from the central axis up to the flux rope boundary where f returns smoothly to a value of 1 (see Fig. 1). Figure 2 features several snapshots of the reduced heating model RM1 applied during runtime, where the flux rope appears as the region of closed poloidal field lines. The volumetric heating rate is seen to be of the order of 10⁻⁴ − 10⁻³ ergs cm⁻³ s⁻¹, which corresponds to values typically taken in literature (e.g., Antiochos & Klimchuk 1991; Dahlburg et al. 1998). At the PIL however, the heating rate greatly increases due to the increased magnetic pressure from bringing the loop footpoints together during the shearing and converging period.

Fig. 1. Cross-section f(x, yO) of f along the horizontal flux rope axis. Heating is reduced most at the flux rope edges since the field lines are longer. Outside the flux rope, there is no reduction as f = 1.

Fig. 2. Demonstration of parameterised heating combined with the reduced background heating (RM1) inside the flux rope, which is tracked during the simulation. Panels a and b feature zoom-ins of the domain.

The maximal reduction in heating rate inside the flux rope is controlled by the parameter δ, which is loosely connected to the 3D length of the flux rope. The relation between δ and L can be approximated assuming that the heating rate ℋ falls exponentially along field lines: let s be the arc-length coordinate along the central axis of the flux rope, which has length L, where our simulation domain consists of a vertical cross-section at s₀ = 0. Then the loop footpoints are at $s=\pm \frac{L}{2}$. Assuming the decay length of ℋ equals the pressure and magnetic scale height H₀ ≈ 50 Mm, we can derive L from δ:

$$\begin{array}{cc}& (1-\frac{\delta }{1.284})=\mathcal{H}(s_0)/\mathcal{H}(\pm \frac{L}{2})=exp(-\frac{|s_0\pm \frac{L}{2}|}{H_0})\hfill \\ \hfill & \iff L=-2H_{0}log(1-\delta /1.284).\hfill \end{array}$$(21)

For the reduced exponential heating model, RE (δ = 0.8), this amounts to a flux rope length of L = 96.6 Mm, while for the reduced quasi-magnetic heating model, RM1 (δ = 0.9), we have L = 119.5 Mm. These numbers fit well within the observed range of prominence lengths summarised by Parenti (2014).

3. Results

In our simulations, a flux rope forms through multiple reconnection events, creating a 2D structure of poloidal nested flux surfaces, alternately hotter and cooler as a result of mergers of multiple subordinate flux ropes – a detailed description can be found in Jenkins & Keppens (2021). With the flux rope fully formed after 2800 s, the energy injection into the plasma as a consequence of the formation process ceases. From then on, the thermodynamic properties of the dense material suspended by the concave upwards field are thus governed primarily by the ratio between the heating and cooling terms of the energy equation. Under the conditions that the local cooling term dominates over the combined influence of the heating or conduction terms, the associated plasma may become unstable to the thermal instability. However, and as we shall henceforth show, when, how, and whether this happens at all depends directly on the adopted heating model in combination with the (anti-)shearing motions involved in the initial formation of the flux rope.

3.1. Condensation phase

Table 1 summarises the different combinations of heating models and footpoint motions under consideration and whether thermal instability occurs, in analogy with (Kaneko & Yokoyama 2015, Table 1). In Sects. 3.1.1–3.1.4, we look at how the flux rope material reacts to these different prescriptions.

3.1.1. Exponential heating model

Model E contains a steady background heating rate that decreases exponentially with height. We consider a flux rope formed through positive shearing motions, as in Jenkins & Keppens (2021), and find similar overall evolution: a dense cloud suspended by the magnetic field collapses following the condensation process described in Claes & Keppens (2019). That is, radiative losses increase, resulting in a decrease of temperature and increase of density, which again enhances the radiative losses and leads to a runaway evolution; matter is accreted along each flux surface, producing a condensation perpendicular to the magnetic field, typical of the thermal instability in 2D (Claes & Keppens 2019; Hermans & Keppens 2021). After the two main condensations form at around y = 3 Mm and 4.5 Mm, some smaller condensations appear at various locations along the outer flux surfaces of the flux rope. These smaller blobs then fall down along the magnetic field under the effect of gravity, reaching speeds up to 190 km s⁻¹ and eventually join the larger condensation. Imbalances in the velocities of inflowing material deform the smaller condensations as they move downwards. The momentum of the plasma gained through falling results in sloshing motions of the prominence body, with distinct oscillation periods along each flux surface. Towards the end of the simulation, the material slows down and the condensed material reaches relative equilibrium, during which the prominence sinks as a whole due to enhanced (η = 0.002) mass slippage over the field lines (Low et al. 2012; Jenkins & Keppens 2021).

3.1.2. Mixed heating models

The similar study of Kaneko & Yokoyama (2015) previously concluded that in situ prominence formation is inhibited when positive shearing motions are coupled with a heating model ℋ ∼ B, while coupling the motions with a model ℋ ∼ ρ does produce condensations. In agreement with their results, we find no condensations in our positive shearing simulations with mixed −0.2 ≤ β ≤ 0 models M0, M1, M2, but rather hot material at a temperature of about 6.3 MK. Since the temperature of the plasma within the flux rope constructed using the positive shearing M2 model was so high, we did not test a positive shearing M3 configuration with α = 2.5. Forming the flux rope through anti-shearing motions, however, does produce condensations. These have similar characteristics for all models, yet are quite distinctive from those obtained with the exponential heating model. In particular, the results for the quasi-magnetic heating model (α, β) = (2.1, −0.1) closely match those for the magnetic heating model (α, β) = (2, 0). For M0–M3, thermal instability initiates at around the same time (t = 3250 s). Figure 3 shows the density and temperature evolution for the quasi-magnetic M1 model with anti-shearing motions. Here, one large monolithic condensation forms at the bottom of the flux rope. Post-formation, the prominence is seen to expand and descend (panels c and d of Fig. 3), while its top moves due to asymmetric inflows. Unlike in the other simulations, parts of the bottom atmosphere become thermally unstable due to the background heating decrease by the anti-shearing motions. As the resulting coronal rain forms outside of the flux rope, it falls down along the open field lines. However, the bottom boundary in our setup does not allow the material to leave the coronal domain as a chromosphere has not been included, and hence the condensations remain unphysically suspended above the bottom boundary. As the prominence body grows, its lower end makes contact with the flux surface on which these condensations reside. To allow for an analysis of the prominence characteristics, we applied an additional consistency check (for ∇ ⋅ B = 0) to the bottom boundary, which slightly increases the temperature inside the subordinate flux ropes due to increased dissipation of localised currents. This subsequently disconnects the prominence material from these ‘coronal rain’ blobs, and enables an analysis of the prominence material up to up to t = 6741 s. After this time, the blobs do connect to the flux surface of the main prominence body and we end the analysis.

Fig. 3. Density and temperature evolution for the simulation performed with quasi-magnetic heating (M1) and anti-shearing motions. One monolithic condensation forms, growing and descending at a later stage in its evolution. A movie of this figure, where the coronal rain blobs outside of the flux rope are also featured, is available with the online version of this manuscript.

3.1.3. Reduced exponential heating model

Model RE combines the exponential heating model with the reduction function f that dynamically modifies the heating rate inside the flux rope, as explained in Sect. 2.4.2. One can thus expect that the simulations will be identical up until the flux rope formation is initiated (t = 1300 s). We therefore initialise the RE simulation with the one for model E at that time. Since we choose δ = 0.8, the heating rate is reduced in the core of the flux rope to approximately 40% of its original value, and to 20% near the edges.

Figure 4 shows the first large condensation to form around 3000 s at a height of around 8 Mm. Generally speaking, condensations occur in every other flux surface, as a result of the temperature alternating between hotter-cooler due to the preceding mergers of several secondary flux ropes. Some condensations consist of only one blob, some consist of a blob on either side of the central vertical axis which fall down and collide. After the collision, a shock-like feature is seen to propagate away from the impact location which resembles the shocks seen at the collapse of condensations within previous simulations of the thermal instability (cf. Claes et al. 2020). Thermal instability remains relevant up to t = 4700 s, after which no additional condensations form as most of the dense matter bound by the flux rope has already cooled down.

Fig. 4. Density and temperature for flux rope formed by shearing motions with reduced exponential heating model. All over the flux rope, condensations form and fall down along the magnetic field, collecting in one large prominence body. A movie of this figure is available with the online version of this manuscript.

Most condensations fall down rapidly along field lines under the effect of gravity, overshoot the concave-up portion of the domain due to their finite kinetic energy, before subsequently performing damped oscillations about x = 0 with different periods along each field line. Some smaller blobs move upwards as a result of upflows and changing pressure conditions as a result of material evolving elsewhere within the flux rope. This is, however, a consequence of the 2.5D description of the flux rope, where the poloidal field line of a flux surface connects plasma elements that would not necessarily be connected in a 3D description of the flux rope. Strikingly, both of these condensations eventually evaporate while falling down at around 5000 s before reaching the main prominence body. After these condensations disappear, much of the material inside the flux rope has collected in the extended prominence, whereas the density of the hot material in the flux rope drops to about 20% of the ambient coronal density. As highlighted well in the online movie that accompanies Fig. 4, the large prominence appears to have a significant effect on the hosting flux rope, clearly deforming the field line topology despite the strong 10 G field. As in the simulations with the mixed heating model, the bottom of the prominence sinks and compresses field lines, leading to an increased perpendicular current density. The outer flux rope shell is heated upon compression by the material at the bottom of the domain, as is evident from panel d of Fig. 4. At the same time, the temperature at the flux rope centre greatly increases through resistive heating and mass slippage of the prominence top.

3.1.4. Reduced quasi-magnetic heating model

Finally, we perform a simulation with the reduced quasi-magnetic heating model (RM1) and a flux rope formed through shearing motions. In Sect. 3.1.2, it was shown that models M combined with shearing motions produce too hot a flux rope for thermal instability to occur, analogous to Kaneko & Yokoyama (2015). However, observational studies have long asserted that shearing motions are ubiquitous to PILs and hence transfer complimentary shear to any associated magnetic arcades that straddle this divide (e.g., Athay et al. 1986). Remarkably, through the dynamic cross-sectional heating reduction, this contradiction is overcome, although outside of the flux rope the heating model does still lead to high temperatures. As can be seen in Fig. 2, the heating rate at the vertical edges at |x| = 12 Mm increases up to order 10⁻³ ergs s⁻¹ cm⁻³ due to low density and high magnetic pressure, leading to downflows that further decrease the local density. One filamentary condensation appears at the bottom of the flux rope (Fig. 5) around 3850 s.

Fig. 5. Density and temperature for flux rope formed by shearing motions with reduced quasi-magnetic heating model. One monolithic condensation forms at the bottom of the flux rope, visibly containing more mass than for the unreduced mixed heating model in Fig. 3. A movie of this figure is available with the online version of this manuscript.

The condensation closely resembles those obtained with models M and anti-shearing motions, but is larger. Again, we observe shock-like features emanating from the initial cloud out along the magnetic field when it collapses. Even though the heating rate is reduced, thermal instability remains absent from the central region of the flux rope, which consequently only contains hot material. The prominence is seen to grow from its lower end, probably as a consequence of the high-density region at the bottom, which reduces the local heating rate in model M1, in combination with the mass-slippage mechanism previously remarked upon within these models (Low et al. 2012; Jenkins & Keppens 2021).

3.2. Comparison of the heating models

Table 1 summarises the outcomes for the combinations of heating models and footpoint motions considered, analogous to Table 1 in Kaneko & Yokoyama (2015). We find those heating models based on local parameters to require anti-shearing motions in order to produce conditions suitable for prominence formation. Applying the ad hoc, masked heating reduction inside the flux rope re-enables the possibility of condensation formation with shearing motions.

The locations of condensation formation vary between the heating prescriptions but remain governed by the imbalance between heating and cooling and the stabilising effect of thermal conduction. Figure 6 features the net heating rate along a vertical cut over the y-axis, immediately after the flux rope has fully formed. The regions of negative net heating rate at this instant indeed correspond well to those flux surfaces that ultimately experience the in situ condensation shown in the first panels of Figs. 3–5. In particular, model RE is seen to lead to a net cooling rate over the entire flux rope, which enables most of the levitated material to cool and rain down.

Fig. 6. Cut along the y-axis showing the net heating rate, −ρℒ, when the driving motions end. Regions where the energy transfer rate is negative are more prone to thermal instability.

The evolution of the mean magnetic energy density over the domain for the four selected simulations that produced condensations, detailed in Fig. 7, shows clearly the distinction between shearing and anti-shearing footpoint motions: the former store energy in the magnetic field while the latter decrease the magnetic shear and hence also the value of B_z. In the initial relaxation phase the magnetic energy evolution is comparable for all simulations. The abrupt change at 1300 s is due to the activation of the driven boundary motions, bringing the footpoints of the initial arcade together and hence driving an energetic evolution of the system. On a related note, there appears to be a relationship between the direction of shearing motions and the number of additional subordinate flux ropes appearing with the current value of η: shearing motions consistently lead to the formation of four large subordinate flux ropes that merge with the main flux rope, while anti-shearing motions produce only two.

Fig. 7. Mean magnetic energy density evolution over the domain for models E, M1, RE and RM1. Shearing motions increase the magnetic energy, anti-shearing motions lead to it decreasing.

An analysis of the maximum temperature over the simulation domain reveals the influence of the heating models on the flux rope and the ambient atmosphere. The mixed models can be ranked according to the power α for the magnetic field strength: a higher α produces a hotter atmosphere and flux rope regardless of shearing direction. Ordered by increasing maximum temperature, we have M0–M1–M2–M3. When the flux rope has fully formed after 2800 s, the maximal temperatures within the mixed models are of the order of 8 MK. For the models with reduced heating, the temperatures are then much lower at around 4 − 5 MK. This is clearly the influence of the diminished heating rate inside the flux rope. For model RE, the maximum temperature also increases at a slower pace due to the lower heating rate inside the flux rope: the centre of the flux rope therefore heats mostly due to resistive Ohmic dissipation (Joule heating). Another effect is that the maximum temperature of model RE levels out soon after the shearing motions end, while it continues to increase for the other heating models. However, at around 5000 s, resistive heating in the flux rope centre becomes more significant and increases the temperature out of equilibrium for model RE as well.

We present the minimum temperature evolution for models E, M1, RE and RM1 within Fig. 8. This temperature is exclusively reached inside the condensations and prominence body, which we limited to |x|≤3 Mm for the run with model M1 as to exclude the coronal rain from the analysis. The late phase, post-condensation temperatures are found in the range 6000 − 10 000 K, which are values well above the artificial minimum of 10³ K adopted in our simulations. For both the original and the reduced heating models, the minimum temperature decreases non-exponentially, although a satisfactory exponential fit can be obtained for models E and M1 nonetheless. This is not unexpected, as the reduced heating models impose conditions out of thermal equilibrium, which could explain the sharp decrease in temperature as a non-linear perturbation. It is hence non-trivial to map the initial growth rate of thermal instability to the evolution into the non-linear regime, as was previously possible in the thermal instability simulations of Claes & Keppens (2019). Nevertheless, we find approximate growth rates ω_i = 0.002 s⁻¹ (model M1), ω_i = 0.004 s⁻¹ (model E) and ω_i ∼ 0.007 (R-models), which agrees well with an estimate from linear theory based on the local parameters, and neglecting thermal conduction and resistivity, of ω_i ∼ 0.001 s⁻¹.

Fig. 8. Evolution of the minimum temperature for models E, M1, RE and RM1. The minimum temperature decreases non-exponentially to less than 104 K, with apparently higher growth rates for the reduced models.

The evolution of the average prominence, and maximum domain densities for models E, M1, RE, and RM1 is presented in Fig. 9, wherein the occurrence of the first condensation is marked with a vertical black line. We see that the reduced heating model RE leads to the development of condensations at an earlier time than model E. For models M1 and RM1, this behaviour is flipped as a result of the accompanying shearing motions in the second case, which lead to less favourable conditions to in situ condensation since the field strength and hence heating rate increases.

Fig. 9. Evolution of maximum (solid) and average (dashed) density for models E, M1, RE and RM1. Black line indicates the appearance of the first condensation.

To calculate the average prominence density, we define the prominence material as having: T < 25 000 K and ρ > 1.171 × 10⁻¹⁴ g cm⁻³ (=50 code units), and in the case of model M1 limit the analysis to |x|≤3 Mm. In all cases, the maximum density is reached at the bottom of the prominence, as a result of the larger-area flux surface encapsulating more material. The small peak towards the end for model RE is, however, reached in the prominence top when it compresses due to field slippage. We find both the maximum density of condensations produced by model RM to be consistently higher than those obtained with the other models. This is most likely related to the nature of the mixed model, which reduces heating in regions of increased density and hence the material is more free to cool and contract in turn. A similar phenomenon may occur with model M if the simulation were extended without interference by the coronal rain. In terms of average prominence density, models E and M are found to outclass those for the reduced models by a factor two. We find this to be a consequence of the larger condensations created by the reduced heating models, where the additional volume is provided by material of lower density. Hence, the distribution of prominence densities with these modified models is concentrated around lower values (see also Sect. 4.3).

We present the evolution of both prominence mass and area for models E, M1, RE, and RM1 in Fig. 10. We note immediately that the two reduced heating models produce larger and more massive prominences compared to the ‘original’ models, with a mass increase of about 30% for model E and 160% for model M1, and increase in area by a factor three. A similar observation applies to models (R)E compared to models (R)M: the mixed models produce significantly less massive condensations, with even the reduced heating model RM1 being barely able to outpace model E. The resulting condensations are however comparable in terms of area. Mass increase with model RE stabilises around 6000 s when the last condensations have either rained down or evaporated. For models (R)E, the prominence mass stabilises only briefly before decreasing due to evaporation of material within the prominence ‘monolith’, at rates of −0.2 and −0.15 g cm⁻¹ s⁻¹, respectively. In the later stages of evolution, both simulations with the exponential models show a decreasing area as a consequence of the concentration of plasma at the concave upwards portions of the magnetic topology over time. Under these conditions the extent of the plasma along a flux surface is increasingly governed solely by the local pressure-scale height.

Fig. 10. Evolution of the total prominence mass per unit length and its 2D surface area, for the exponential and quasi-magnetic heating models (dashed) and their reduced counterparts (solid). Both in terms of area and mass, the reduced heating models lead to larger prominences.

In addition to the size of condensations, the total area evolution also contains information on the oscillations of condensed material in the dipped magnetic field as mentioned above and visible in Figs. 4 and 5. Extracting the period of the largest oscillatory motions from the periodicity in the area evolution, we find dominant periods of 335 s for model E, 904 s for model RE and 840 s for model RM1. We explore this in more detail in Sect. 4.4. The quasi-magnetic heating model M1 does not lead to any oscillations since the single condensation forms in-place at the locations of concave-up magnetic topology.

3.3. Phase space evolution

Criteria for the initiation of the thermal instability have been derived under both isochoric and isobaric criteria (Parker 1953; Field 1965). To ascertain which of these approximations is more consistent with the evolution found within our prominence condensations, we now look at the ρ⁻¹ − T phase space, shown in Fig. 11, at several snapshots for the simulation with model RE. The coloured diagonal lines indicate the pressure isocontours governed by the normalised ideal gas law p = ρT. The distribution is coloured by total volume of cells in the simulation domain that correspond to a certain bin in phase space, and hence contains information on the physical extent of regions in phase space. A similar analysis was performed by Waters & Proga (2019) for unmagnetised astrophysical fluids. In contrast to their analysis, we do not use the concept of an ‘S-curve’ of thermal equilibrium conditions, but instead the qualitative description of the phase evolution of thermal instability.

Fig. 11. Evolution in ρ−1 − T state space for model RE, with pressure isocontours diagonally. Distribution given in total cell volume. A movie of this figure is available with the online version of this manuscript.

Initially, the atmosphere is isothermal with a varying density profile due to hydromagnetostatic equilibrium, which translates to a horizontal line in the ρ⁻¹ − T diagram (panel a). As time progresses, the atmosphere begins to heat and the distribution in phase space takes on complex shapes when the coronal loop footpoints are driven. After some time, the thermal instability is initiated and a branch towards the cold, dense lower left corner of the state space domain splits off (panels c and d). The abrupt cooling of the material seems to follow the pressure isocontours fairly well in the initial stage at t = 2903 s, but the branch is seen to cross several isocontours soon afterwards and follows a more vertical and hence more isochoric path. In panel e, the evolution of the condensation at its coolest point happens under almost isothermal conditions as a consequence of the significant drop in radiative losses at ≈10⁴ K. In the final snapshot on panel f, the hot, tenuous flux rope and relatively colder, massive condensation appear as two distinct populations in the lower left and upper right corners of the distribution, respectively. In fact, two hot and tenuous populations can be distinguished as two lines of increased volume aligned with the pressure isocontours and hence at approximately uniform pressure. The left, upper population corresponds to the atmosphere surrounding the flux rope, while the right, lower population corresponds to the flux surfaces on which the prominence resides. The latter are tenuous because much of the material has condensed. The hot flux rope centre appears as a vertical extension in the right upper corner of the phase diagram.

4. Discussion

We performed 2.5D simulations of prominence formation through the levitation-condensations mechanism, where the basic simulation setup is taken as in Jenkins & Keppens (2021). The original model of Kaneko & Yokoyama (2015) has been further modified to include a density stratification profile that is fully consistent with the varying gravitational acceleration Eq. (10) over the domain. A value of B_a = 10 G was taken for the background magnetic field at the bottom of our simulation domain while the simulation of Jenkins & Keppens (2021) at the same resolution had a 3 G background field. This has implications for the magnitude of Ohmic heating at the flux rope centre and hence the amount of material that is able to condense.

4.1. Influence of heating models on prominence formation

4.1.1. Exponential heating models

The simulation performed with the exponential heating model Eq. (15) closely matches the results from Jenkins & Keppens (2021). However, in our case the material collects in one main prominence body while it remained disconnected in different flux surfaces in the 10 G field simulations of Jenkins & Keppens (2021). This is most likely a consequence of the lower resolution taken in our simulations, which leads to a slightly lower heating input from reconnection of subordinate flux ropes. Nevertheless, the large central part of the flux rope remains as a dense and hot region insusceptible to thermal instability, as evident from the heating-cooling balance in Fig. 6. Indeed, following the formation of the flux rope, the material bound by the flux rope should be isolated from the strongest contributions of the background heating as thermal conduction can no longer transport this energy along field lines. Therefore, we introduced a method to implement an ad hoc modification to the heating profile so as to approximate the influence of field line connectivity to the heating experienced within the 2.5D flux rope core.

This reduced exponential heating model leads to a significantly cooler flux rope than its unmodified counterpart, enabling radiative cooling to easily dominate over background heating. Strikingly, nearly all of the material in the flux rope is found to cool and rain down into the prominence, turning the flux rope into a cavity with densities of the order of 10⁻¹⁶ g cm⁻³. Although similar results were found before, this is the first time it has been demonstrated with a field strength as high as 10 G. The flux rope centre is much cooler than in previous simulations without reduced heating and even hosts in situ condensations. In the final stages of the simulation, the flux rope temperature increases to order 10 MK which is rather high compared to observed cavity temperatures of 1.5 − 2.2 MK (Bąk-Stȩślicka et al., 2019). Such an evolution is compounded by the high resistivity value in our simulations leading to increased Joule heating in the cavity combined with increased mass slippage at the prominence top.

4.1.2. Mixed heating models

The original study of Kaneko & Yokoyama (2015) investigated the effect of background heating ℋ ∼ B² and ℋ ∼ ρ in a similar setup. They used a background magnetic field of 3 G and an artificial cut-off value on the minimum temperature to ensure numerical stability. In this work, we performed simulations with the mixed heating model, depending on both B and ρ at the same time (model M), as in Mok et al. (2008). Kaneko & Yokoyama (2017) also briefly discuss a 3D simulation with anti-shearing motions and (α, β) = (1, 1), although their main result focuses on an (α, β) = (2, 0) configuration. In full agreement with the results of Kaneko & Yokoyama (2015), model M only leads to condensations in flux ropes formed through anti-shearing motions. This can be explained by the increased magnetic energy stored in the field by shearing motions, which act to amplify B_z (Fig. 7) and hence increase the local heating rate within the flux rope. Moreover, regions throughout the simulation domain where the magnetic pressure is high, like the converging loop footpoints and the flux rope edges, also experience a high heating rate that leads to unfavourable conditions for the condensation process. Hence, for those simulations involving anti-shearing motions, only the bottom of the flux rope and lower atmosphere with increased density could provide sufficient radiative losses to initiate the thermal instability and lead to the formation of one large monolithic condensation in the flux rope and a string of coronal rain blobs in the surrounding atmosphere. Then, as already mentioned, our adoption of a 10 G background field strength, although bringing us closer to observational values (Casini et al. 2003; Wang et al. 2020), increases the core heating rate even further.

As mentioned above, attempts to form condensations with the combined use of shearing motions and B-dependent heating models had thus far proven unsuccessful. Here, we introduced the physically motivated dynamic masked heating reduction throughout the flux rope interior, and this is shown to lead to successful prominence formation. Indeed, as anti-shearing motions are not commonly observed along PILs (Gibb et al. 2014; Mackay 2015, where differential rotation usually increases the shear of magnetic arcades), it was previously highly puzzling as to how a prominence could form under ‘typical’ solar conditions since 2.5D models indicated such behaviour would categorically prohibit prominence formation when adopting B-dependent background heating. Moreover, the spontaneous coronal rain phenomenon reported here indicates that anti-shearing motions may decrease magnetic energy too substantially to consistently explain prominence formation.

4.2. (Additional) influence of Joule heating

For all of the setups presented here, we adopted a constant η value for the resistivity throughout the simulation space and time. This permits a range of resistive evolutions within the simulations, such as magnetic reconnection (but see also Sect. 4.4), with perhaps the most ubiquitous of them being the dissipation of strong currents as Joule heating η|j|². In Table 2, we quantify the magnitude of this Joule heating at two typical locations, and compare it against the contribution from the assumed background heating model.

In the centre of the flux rope at the end of the simulation, we find the Joule heating contribution to be of the same order of magnitude, but in all cases smaller than the contribution of the background heating model. During the preceding formation process, however, the picture is quite the opposite with the Joule heating contributing a considerable amount of energy to the local plasma (cf. the alternating hot shells in Jenkins & Keppens 2021), most significantly at the footpoint locations during their migration towards the PIL.

It is clear that the evolutions within the simulations lead to the generation of strong currents. In all cases, the influence that these currents have on the thermodynamic properties are directly proportional to the magnitude of the assumed resistivity η since the Joule heating scales linearly with it. In the current simulation, we deliberately overestimated the magnitude of this resistivity so as to assist in the flux rope formation process – a magnitude that was necessarily above the numerical resistivity limit set by the resolution. In future simulations we may relax this overestimation, as an increase in spatial resolution will similarly decrease the restriction on the resolvable resistivity. Subsequent current enhancements may peak to higher values with smaller extents, depositing the equivalent energy, but we nevertheless suggest this leading to an overall reduction in the amplitude of the associated Joule heating across the simulation domain. The non-linear nature of these simulations requires, however, dedicated studies so as to accurately conclude the influence of these modifications. For more dynamic scenarios, one may also explore the influence of the spatially varying ‘anomalous’ prescriptions, as in Zhao et al. (2017) and Ruan et al. (2020).

4.3. Towards increased realism in prominence simulations

Despite the ever-increasing body of work on the coronal heating problem, no single heating mechanism or combination thereof has yet been deemed the definitive solution (remaining key questions are highlighted by Klimchuk 2015). How this heating affects prominence formation is here shown to differ between various parametric models. Specifically, the timescale on which instability occurs, decreases for parametric models, as seen on Fig. 9, which also indicates that imposing a time-independent exponential heating background without differentiating between the flux rope and the surrounding atmosphere imposes a non-negligible residual heating rate inside the flux rope, affecting the onset of thermal instability: it was shown that applying a heating reduction inside the flux rope decreases this time significantly. Then, for models (R)M1, we found the directional choice of shearing motions to have a large impact, explaining why condensations appear at a slightly later time for the simulation that was paired with positive shearing despite employing the reduced heating model. It is thus expected that the same reduced heating setup with anti-shearing motions would instead lead to a considerably earlier development of condensations compared to the standard M1 setup.

In Table 1, we list the n_max and T_min characteristics of the prominences towards the end of the simulations, as formed by the different heating models. Despite the different models employed, the maximum number density and minimum temperature are comparable for most heating prescriptions, indicating that the physics governing the prominence material, once formed, is independent of the heating prescriptions explored here. Furthermore, prominence observations summarised by Parenti (2014) list number densities ranging from 10⁹ to 10¹¹ cm⁻³, which shows that the values obtained in all our simulations are in good agreement. The same may be said for the core temperature of the prominences of our simulations, which is reported by the same author to have a range of 7500 − 9000 K. The spatial dimensions of our simulated prominences are, however, on the smaller side of observed widths between 1 and 10 Mm resulting from the relatively small computational domain. Pressures inside the prominence are found in the range 0.08 − 0.41 dyn cm⁻² for models E and M1 and in the range 0.05 − 0.64 dyn cm⁻² for their reduced counterparts, in perfect correspondence with observations (Parenti 2014).

Panel a of Fig. 10 clearly shows evaporation of prominence material for models (R)E, which is found to occur at the prominence top where the cool material resides next to the hot flux rope core. This is possibly a result of Joule heating through the large perpendicular currents, which could provide a net heating effect in the neighbouring parts of the prominence.

One distinct difference arises between the average prominence densities of models with and without heating reduction in Fig. 9: the former have lower average densities. We calculated the distributions of prominence density and pressure to seek an explanation. Models RE, RM1 produce very peaked distributions around the (lower) average values, indicating that the resulting prominences are largely uniform in pressure and density. The uniform region takes up most of the prominence, which does show stratification at top and bottom ends. For model M1, the distributions have a larger standard deviation but are still rather peaked, although around multiple higher density or pressure values. For model E, the distributions are much broader, pointing to prominences with a strongly pronounced internal density and pressure stratification, which clearly appears on a vertical profile taken within the prominence. This stratification can be solely attributed to the steady exponential heating background, which imposes a vertical stratification inside the prominence body instead of a stratification along field lines as would be expected for plasma bound to the magnetic field (Blokland & Keppens 2011a,b). Indeed, a similar stratification is observed in the bottom part of the prominence that escapes the modified heating mask of model RE – a comparison between Figs. 2 and 5 indicates how the lowest-most extent of the prominence lies outside the ellipse. In particular, the part of the prominence for model RE inside the reduced heating mask – which has almost uniform density and pressure – is seen to extend further horizontally in the magnetic dips compared to the bottom end due to the different thermodynamic properties of the plasma, more specifically due to a different pressure scale height. A vertical cut along the prominence, featured on Fig. 12, clearly shows an exponential density profile in the prominence material outside the reduced heating mask, and almost uniform conditions above the location where the heating mask is applied. Such a feature, in particular, provides a particularly strong argument as to why an exponential background should not be applied indiscriminately to 2.5D prominence formation in flux ropes: not only does the exponential influence the timescales, but also the stratification of plasma properties within the flux rope – a key component in synthesising any resulting simulation as an observation (e.g., Jenkins & Keppens 2022).

Fig. 12. Vertical profile of density, temperature and heating rate through the prominence for model RE. The density varies exponentially in the lower part of the prominence and is more uniform in the upper part inside the heating mask.

To compare our obtained prominence masses against observations, we multiply the mass per unit length by a typical flux rope or prominence length of 100 Mm – this assumption agrees with the δ parameter used in the reduced heating prescription (21), although the prominence plasma rarely extends down towards the chromospheric footpoints of a flux rope and this estimation hence gives an upper limit. The estimated prominence masses then lie between 2.3 and 8.3 × 10¹³ g, a similar order of magnitude to previous simulation results (Jenkins & Keppens 2021), but still below the typical inferred values from observations of 10¹⁴ − 2 × 10¹⁵ g (Parenti 2014). Hence, the inclusion of more-realistic considerations for the heating, however valid, do not address the outstanding issue of creating more ‘realistic’ prominences in terms of order-of-magnitude total mass content. Still, the reduced heating approach does increase the prominence masses compared to the other heating models. The larger masses can generally be attributed to the increased net energy loss rate due to the decrease in background heating, which in turn makes a larger portion of the flux rope more prone to thermal instability. The correspondence between observational and simulation number densities and temperatures, yet the persistent underestimation for total mass content indicates that the simulation domain and subsequent dimensions of the ab initio prominences are smaller than present within the actual solar atmosphere. The consideration of both a larger domain, and one that includes a chromosphere, would provide the coronal flux rope with an abundance of additional material, as already shown with the realistic prominence mass values found in the work of Zhao et al. (2017). Moreover, the justification for the artificial heating reduction employed in our 2.5D simulations is automatically accounted for in 3D simulations, since field-aligned thermal conduction is unable to efficiently transport energy towards the middle of the flux rope. We thus envision the combination of a 3D setup as in Kaneko & Yokoyama (2017) that implements localised heating considerations and drives the formation of a flux rope via positive shearing motions with the presence of a chromosphere, to be the natural next step. Moreover, the reduced heating model presented in this work based on common assumptions in 1D prominence formation models can benefit from recent progress, for example by Huang et al. (2021), who combined evaporation and injection into one model. Extending our setup as described above could additionally incorporate these effects in the current levitation-condensation model by supplying localised heating in the chromosphere. It is anticipated that the above aspects will be crucial in both further advancing these models towards self-consistency, and facilitating effective comparisons against observations. The results obtained in this paper hence will, and should, one day be superseded by more realistic models of a complete solar atmosphere.

4.4. Additional results: slippage and oscillations

As reported in earlier work with the same setup, the high resistivity value in our simulations enhances the effect of mass slippage over the field lines (Jenkins & Keppens 2021; Low et al. 2012). Due to this effect, the prominence body as a whole performs motions perpendicular to the poloidal magnetic field, where through resistive dissipation it sinks towards the bottom of the simulation domain. We observe this phenomenon in all our simulations, with its typical signature of increased current density at the slippage sites. In a benchmark simulation, where resistivity was switched off after the formation of the flux rope, this slippage effect was all but eliminated.

The observed mass slippage also affects the prominence oscillations reported in this work. Condensations form at different heights throughout the flux rope on distinct poloidal field lines. After falling down, they perform oscillatory motions about x = 0 with different periods, as is clearly visible in Fig. 4. The periodicity is governed by the pendulum model, with dependence on the local radius of curvature of the magnetic field (Arregui et al. 2018). The oscillations found in this work occur spontaneously in response to evolving condensation formation and do not require an imposed perturbation. Longitudinal oscillations have been invoked to explain counter-streaming motions in prominences (Zirker et al. 1998), as they have been found in other simulations of prominence formation (e.g., Xia et al. 2011). An analysis of the oscillations along selected field lines – using the method from Liakh et al. (2020) – for model RE indeed reveals a period dependence on height, where the period increases from top of the prominence (region closest to the centre of the host flux rope) to bottom with decreasing the radius of field curvature (the accompanying analysis is included in a recent conference proceeding, Brughmans & Keppens 2022). In particular, oscillation damping times also show a height dependence. Oscillations are even temporarily amplified inside the top regions of the prominence as well as in its bottom part, as was reported for the first time by Liakh et al. (2020, 2021). Again from the benchmark simulation without resistivity, we find that in general, mass slippage decreases the apparent periods and increases the apparent damping time. We note that the oscillations observed here are plane-projected manifestations of longitudinal oscillations. Observations of longitudinal prominence oscillations yield typical periods of 50 − 60 min (Arregui et al. 2018), somewhat longer than the periods of 5 − 10 min reported in our work, which again underlines the need to consider larger simulation domains and full 3D models with line-tied flux rope ends in future works.

4.5. Condensation velocities and phase space

An analysis of the maximum velocity over the domain provides further insight in the formation and evolution of condensations. During the formation phase, inflow velocities of 55 − 100 km s⁻¹ are found. These velocities decrease after the formation of the initial condensation for models M1, RM1, but for the exponential models, subsequently, velocities up to 190 km s⁻¹ are observed within those falling condensations that originate near the apex of the flux rope. These high velocities, in excess of those anticipated from free fall acceleration, are a consequence of pressure evolution associated with the formation of condensations at other locations along the same flux surface. This surplus of kinetic energy driven by dramatic pressure evolution within the exponential models is the reason why the material performs damped oscillations with large initial amplitudes after formation.

The phase space evolution in Fig. 11 featured a non-isobaric evolution of the thermal instability for model RE. Models E, M1, and RM1 also lead to a very similar evolution, although in the very onset of thermal instability, the distribution follows the pressure isocontours a bit longer than for model RE. The non-isobaric evolution of the condensation process does not come as a surprise: the inclusion of gravity introduces non-isobaric dynamics – Jenkins & Keppens (2021) already found baroclinicity localised to the condensations. In short, the onset of thermal instability visually matches isobaric conditions, but not far into the linear evolution, the isochoric criterion becomes more appropriate as the temperature decreases faster than the density increases (Xia et al. 2012). Moschou et al. (2015) also arrived at this conclusion by identifying the dominant instability criterion in the thermally unstable regions associated with coronal rain.

5. Conclusions

We performed 2.5D simulations of prominence formation through the levitation-condensation mechanism with two classes of heating models, and an additional ad hoc modification for 2.5D simulations that approximates the 3D flux rope structure, and tracks the flux rope during runtime based on magnetic curvature. Exponential heating models lead to larger prominences but have a high residual heating rate inside the flux rope and prominence, while the mixed models have the more realistic assumptions but do not produce condensations when combined with positive shearing motions. Both of these inconsistencies are overcome by the reduced heating mask: the reduced exponential model leads to large condensations, almost at the lower edge of observed masses for whole prominences, and simultaneously eliminates the residual flux rope heating, thus creating a much cooler flux rope and prominences with almost uniform density and pressure. Crucially, those models with both magnetic-field-strength-dependent heating and the ad hoc reduction mask, combined with shearing motions, now lead to the successful formation of condensations. As both classes of heating models are shown to lead to prominences with different evolutions and morphologies but very similar physical properties, taking these models to either a larger 2D or even fully 3D domain will likely overcome many of the problems and inconsistencies that continue to reside within our simulations, and in turn further increase the simulated prominence masses. Finally, a phase-space visualisation of the condensation process describes neither isobaric nor isochoric behaviour, but rather a combination of the two, with a state of constant pressure along flux surfaces recovered once force-balance is achieved within the flux rope.

A natural progression of this work is to extend the setup to 3D, using a larger simulation domain with the inclusion of a chromosphere. This will enable us to assess the effect of the mixed heating model on a flux rope formed through positive shearing motions and whether the resulting prominences can reach realistic masses. Indeed, the inclusion of a chromosphere and transition region will drastically modify the energy balance of the model atmosphere, which will have a large effect on stability provided by thermal conduction of energy towards the transition region, where it is radiated away. The results of this work, which only considers a coronal domain, could then be modified significantly in view of this added realism, where perhaps thermal instability is less likely to occur except when evaporation from the chromosphere is taken into account. In fact, such a setup has been successfully implemented by Xia & Keppens (2016). It remains unclear how the cooling of the longest field lines obtained by Kaneko & Yokoyama (2017) would be affected by the presence of a chromosphere and transition region.

Movies

Movie 1 associated with Fig. 3 (aa44071-fig3_movie) Access here

Movie 2 associated with Fig. 4 (aa44071-fig4_movie) Access here

Movie 3 associated with Fig. 5 (aa44071-fig5_movie) Access here

Movie 4 associated with Fig. 11 (aa44071-fig11_movie) Access here

Acknowledgments

N.B. acknowledges support by his Fonds Wetenschappelijk Onderzoek (FWO) fellowship (grant number 11J2622N) and would like to thank Valeriia Liakh for her advice on analysing prominence oscillations. J.M.J. and R.K. received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 833251 PROMINENT ERC-ADG 2018), and from Internal Funds KU Leuven, project C14/19/089 TRACESpace. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government – department EWI. Simulation visualisations were created using yt (https://yt-project.org/).

References

Appendix A: Flux rope tracking

We here set forth to describe the algorithm used to track the flux rope during runtime. First, the 3D field curvature is calculated over the domain using,

$$\begin{array}{c}\hfill \mathit{\kappa }=\mathbf{b}·\mathrm{\nabla }\mathbf{b},\end{array}$$(A.1)

where $\mathbf{b}=\frac{\mathbf{B}}{B}$ is the unit vector along the magnetic field. Since there is no variation in the z-direction in our 2.5D setup, the gradient in (A.1) has only non-trivial components along x and y. The centre of the flux rope is then detected as a point of minimal curvature as it corresponds to a straight field line that locally is oriented directly out-of-plane in the 2.5D representation. As we do not deliberately displace the central axis of the flux rope in our simulations (cf. Jenkins & Keppens 2021), we can assume that the flux rope centre lies directly above the PIL and search for a point (0, y_O). This detection method is in accordance with Yang et al. (2019), who found that weak curvature correlates with a strong normal force (i.e. straight field lines in 2.5D) and strong curvature occurs in the neighbourhood of points with a weak magnetic field (i.e. 3D O-points). The flux rope edges are detected instead as local maxima in the curvature, resulting from the dynamical evolution of the flux rope in the surrounding atmosphere.

A.1. Detailed algorithm for tracking

Suppose an N × N grid with resolution Δx × Δy is given by a collection of points {(x_i, y_i) | i = 1, …, N}. The characteristics of the model ellipse are fixed by three values: let y_O be the height of the flux rope centre and x_b = x_O + b,  y_a = y_O + a the horizontal and vertical coordinate of the corresponding edges, ascertained at the previous simulation time. The new locations in the current timestep are then determined by solving optimisation problems of the locations and curvature, which are solved in a discrete manner by looping over the AMR grids.

To find the flux rope centre, we solve the following optimisation problem,

$$\begin{array}{cc}& \text{Maximise}{y}_i\hfill \\ \hfill & \begin{array}{cc}\text{such that}& |x_i|\le \mathrm{\Delta }x,\hfill \\ & \kappa (x_i,y_i)\text{is minimal and}\le 0.01,\hfill \\ & |y_i-y_O|<d,\hfill \\ & 0<i<N.\hfill \end{array}\hfill \end{array}$$

If no new location of minimal curvature satisfying these conditions is found, the maximal allowed jump d is increased by 10% and the problem solved again, up to a maximum of d = 0.2. Here, d varies throughout the simulation to accommodate the more dynamic periods of flux rope formation and prevent ‘overshooting’. Its value in code units depends on the variable timestep Δt and is given by

$$\begin{array}{c}\hfill d=\{\begin{array}{cc}1\hfill & \text{if}{y}_O,x_b,y_a\text{have not all been found,}\hfill \\ 10\mathrm{\Delta }t\hfill & \text{if}{y}_O,x_b,y_a\text{have all been found.}\hfill \end{array}\end{array}$$(A.2)

After the coordinate y_O has been found, b is found by looking for maxima of κ over a horizontal strip at height y_O while a is found in the same way over a vertical strip at x_O = 0. To find b, the following problem is solved,

$$\begin{array}{cc}& \text{Maximise}{x}_i\hfill \\ \hfill & \begin{array}{cc}\text{such that}& |y_i-y_O|\le \mathrm{\Delta }y,\hfill \\ & x_i\ge 0,\hfill \\ & \kappa (x_i,y_i)\text{is maximal},\hfill \\ & |x_i-x_b|<5d,\hfill \\ & 0<i<N.\hfill \end{array}\hfill \end{array}$$

If the above problem has no solution, d is again increased in steps of 10% until it reaches d = 0.1. While approximately d < 45Δt, we first look for solutions x_i ≥ x_b that increase the dimensions of the flux rope since it is expanding for most of its dynamic evolution, but this condition is relaxed when no solution is found. The optimisation problem for finding the location of the vertical edge y_a of the flux rope is similar, with x_i, b interchanged with y_i, a.

A.2. Flux rope tracking results

The tracked location of the flux rope centre is shown to oscillate as a result of mergers between the primary and subsequent subordinate flux ropes (cf. Jenkins & Keppens 2021). We find an average initial period of 280 − 300 s.

The magnetic flux through the tracked flux rope is approximated as Φ = B_zπab, where B_z is taken to be the field strength reached in the flux rope centre. The obtained values are hence rather an order-of-magnitude estimate, more so since the horizontal length b is underestimated after the end of the flux rope formation. We find the estimated magnetic flux bound by the flux rope, across models, to be of order 10¹⁹ Mx, consistent with earlier simulations (e.g. Zhao et al. 2017) yet several orders of magnitude below observational results (e.g. 10²¹ Mx in Dissauer et al. 2018) which leads to the conclusion that we are working with a small flux rope in these simulations.

The tracked flux rope is not always in a 1-to-1 correspondence with the actual flux rope. First, the tracking algorithm can suffer from mis-identifications whenever the field topology undergoes a sudden evolution during, for instance, flux rope mergers, but the shape is recovered after a few time units. Second, the assumption of a perfectly elliptic poloidal field will by definition exclude some portions of the flux rope immediately above the X-point as the topology is deformed by a combination of magnetic pressure and tension associated with the ongoing reconnection (cf. Jenkins & Keppens 2021). Third, once the footpoint driving motions are switched off the flux rope is assumed to be fully formed. However, due to the large value of η maintained here on out, a slower diffusion of field lines leads to the expansion of the flux rope partially out of the upper domain boundary. After this time, the previous methods and justifications for defining the edge of the flux rope, through curvature features alone, are no longer applicable. To avoid inconsistencies and mis-identifications, we opt to fix the horizontal half-length b from the moment the driving velocities cease, and update b if and only if a higher value is detected. Hence, as time progresses, the detected dimensions tend to underestimate the actual size of the flux rope. An example of this underestimation can be observed in Fig. 2, panel d.

All Tables

All Figures

	Fig. 1. Cross-section f(x, yO) of f along the horizontal flux rope axis. Heating is reduced most at the flux rope edges since the field lines are longer. Outside the flux rope, there is no reduction as f = 1.
In the text

	Fig. 2. Demonstration of parameterised heating combined with the reduced background heating (RM1) inside the flux rope, which is tracked during the simulation. Panels a and b feature zoom-ins of the domain.
In the text

	Fig. 3. Density and temperature evolution for the simulation performed with quasi-magnetic heating (M1) and anti-shearing motions. One monolithic condensation forms, growing and descending at a later stage in its evolution. A movie of this figure, where the coronal rain blobs outside of the flux rope are also featured, is available with the online version of this manuscript.
In the text

	Fig. 4. Density and temperature for flux rope formed by shearing motions with reduced exponential heating model. All over the flux rope, condensations form and fall down along the magnetic field, collecting in one large prominence body. A movie of this figure is available with the online version of this manuscript.
In the text

	Fig. 5. Density and temperature for flux rope formed by shearing motions with reduced quasi-magnetic heating model. One monolithic condensation forms at the bottom of the flux rope, visibly containing more mass than for the unreduced mixed heating model in Fig. 3. A movie of this figure is available with the online version of this manuscript.
In the text

	Fig. 6. Cut along the y-axis showing the net heating rate, −ρℒ, when the driving motions end. Regions where the energy transfer rate is negative are more prone to thermal instability.
In the text

	Fig. 7. Mean magnetic energy density evolution over the domain for models E, M1, RE and RM1. Shearing motions increase the magnetic energy, anti-shearing motions lead to it decreasing.
In the text

	Fig. 8. Evolution of the minimum temperature for models E, M1, RE and RM1. The minimum temperature decreases non-exponentially to less than 104 K, with apparently higher growth rates for the reduced models.
In the text

	Fig. 9. Evolution of maximum (solid) and average (dashed) density for models E, M1, RE and RM1. Black line indicates the appearance of the first condensation.
In the text

	Fig. 10. Evolution of the total prominence mass per unit length and its 2D surface area, for the exponential and quasi-magnetic heating models (dashed) and their reduced counterparts (solid). Both in terms of area and mass, the reduced heating models lead to larger prominences.
In the text

	Fig. 11. Evolution in ρ−1 − T state space for model RE, with pressure isocontours diagonally. Distribution given in total cell volume. A movie of this figure is available with the online version of this manuscript.
In the text

	Fig. 12. Vertical profile of density, temperature and heating rate through the prominence for model RE. The density varies exponentially in the lower part of the prominence and is more uniform in the upper part inside the heating mask.
In the text