© ESO 2021

1. Introduction

The mechanisms suggested to be responsible for the formation of solar prominences within the million degree solar corona are historically split into three main processes: levitation, injection, and condensation (cf. Mackay et al. 2010). The ‘levitation’ mechanism thereby describes the lifting of already cool but very heavy material (i.e., plasma-β ⪆ 1) from the solar chromosphere up into coronal heights. In that process, the magnetic dips are already present and provide the upward balance against gravity (see e.g., Kippenhahn & Schlüter 1957; van Ballegooijen & Martens 1989; Zhao et al. 2017, 2019; Zhao & Keppens 2020). The ‘injection’ mechanism usually invokes some magnetic reconnection within the chromospheric volume so as to propel material into coronal heights. Sustained capture and suspension of this propelled plasma is then only achieved if the correct topology already exists within the corona and is located along field lines that have propelled material; otherwise, this material would simply fall back to the chromosphere (Wang 1999; Chae et al. 2001). The final method, ‘condensation’, takes advantage of the optically thin emission processes occuring within the coronal volume (as eluded to by Kiepenheuer 1953). Herein, coronal plasma cools into condensations that collect within, or themselves create, the topological dips required for their longer-term stability (Antiochos & Klimchuk 1991). In recent years, modern variants or mixtures of these mechanisms have come into focus, with the ‘evaporation-condensation’ process especially invoking chromospheric plasma evaporations due to localised, sustained footpoint heating followed by in situ condensations in arcades and flux ropes (e.g., Xia et al. 2014; Keppens & Xia 2014; Xia & Keppens 2016). In this work, we revisit the coronal ‘levitation-condensation’ mechanism proposed in Kaneko & Yokoyama (2015); in that work a coronal-only volume is considered, where a flux rope forms and coronal plasma is lifted and finally condenses in situ due to TI.

The condensation and TI mechanisms have been studied extensively over the last half century (Parker 1953; Field 1965). However, for prominences it was quickly realised that this in situ condensation process alone could not provide enough plasma to the prominence to explain their characteristic masses since the associated coronal volume simply does not contain sufficient mass (e.g., Saito & Tandberg-Hanssen 1973). Subsequently, authors explored pairing the condensation mechanism with localised heating at the footpoints of the associated magnetic fields. Such approaches are also commonly invoked in the thermal non-equilibrium (TNE) cycle models to explain the mass supply to the coronal rain phenomenon (cf. Klimchuk 2019; Antolin 2020). These TNE models have thus far been explored parametrically in 1D hydrodynamic settings and require a magnetic connectivity to the lower chromosphere along with parametrised heating near transition region heights. The resulting additional ‘evaporation’ introduces a method of providing more material to the coronal volume without the requirement for any reconnection, which is needed in the injection model (Antiochos et al. 1999; Klimchuk 2019). The continuous supply of material to the coronal volume in this way overcomes the limitations of the condensation model. Indeed, this evaporation-condensation mechanism has been adopted as the basis for many recent numerical simulations of coronal rain (e.g., Fang et al. 2013; Moschou et al. 2015; Xia et al. 2017; Klimchuk & Luna 2019) and prominence formation (e.g., Xia et al. 2014; Keppens & Xia 2014; Xia & Keppens 2016; Fan 2017, 2018, 2020). In contrast to these works, we will here focus on a levitation-condensation configuration, where TNE cycles are completely absent, since there will be no connectivity to lower layers, and everything is happening due to TI and other magnetohydrodynamic (MHD) processes. It should be noted that the pure TI mechanism has also been studied in more idealised settings by Claes & Keppens (2019) and Claes et al. (2020), who point out that the initial stages of the TI process do not lead to a perfect alignment between condensations and magnetic field lines, bringing further complications to interpretations of observed fine-structure.

Prominences, once formed, come in a wide variety of shapes, sizes, and with dynamics understood to be driven by their formation, internal, and global evolution (Labrosse et al. 2010; Mackay et al. 2010; Vial & Engvold 2015). Structurally, stable prominences often appear as ‘clouds’ within the solar atmosphere with internal striations aligned approximately to the radial direction (e.g., Berger et al. 2008). In comparison with their on-disk counterparts, or ‘filaments¹’, similar striations or ‘threads’ appear but are instead aligned predominantly with the solar surface (e.g., Lin et al. 2005). For both the prominence and filament projections, the internal dynamics of the associated material are oriented in line with their internal structures, which are themselves conflicting despite being identical phenomena (e.g., Zirker et al. 1998; Berger et al. 2008). Furthermore, the observed magnetic field within prominences favours the picture laid out by general filament observations in that their associated magnetic field, responsible for their suspension and therein their appearance, is oriented predominantly horizontal to the solar surface (as predicted by Kippenhahn & Schlüter 1957, see e.g., Casini et al. 2003; Merenda et al. 2007; Orozco Suárez et al. 2014; Wang et al. 2020).

For more than a decade now, this discrepancy between the horizonal field of a prominence and their vertical internal structuring and dynamics has been interpreted as a form of a magneto-convection process (Berger et al. 2008, 2010, 2011; Liu et al. 2012a). Analysis of prominence observatons provided by the Solar Optical Telescope (SOT; Tsuneta et al. 2008) instrument on board the Hinode (Hinode; Kosugi et al. 2007) spacecraft led Ryutova et al. (2010) to interpret the vertical structuring as combinations of the ‘plumes’ and ‘fingers’ associated with the Rayleigh-Taylor instability (RTI). Numerous authors have successfully explored this interpretation further with both observations, numerical models, and combinations therein (e.g., Hillier et al. 2011, 2012a,b; Keppens et al. 2015; Kaneko & Yokoyama 2018). However, Low et al. (2012a) presented an alternative theory, where this observed behaviour is instead a consequence of the quasi-periodic breakdown in the frozen-in condition, allowing prominence plasma to slip across field lines towards the surface (see also Low et al. 2012b; Low & Egan 2014). In this paper, we present the first numerical demonstration where this slippage is indeed occuring, and where the vertical motions found are not directly related to the RTI. By restricting our setup to 2.5D, true interchange or RTI activity is strongly inhibited by means of the in-plane magnetic field, while we can make direct contact with advanced linear MHD stability criteria (Goedbloed et al. 2019).

Finally, pioneering 3D numerical simulations of Xia & Keppens (2016) effectively demonstrated the ability of condensations to spontaneously form throughout a prominence-hosting flux rope. The subsequent and complex evolution of these plasma packets to lower heights both within and out-of the global flux rope, reproduced the complex absorption and emission evolution observed within prominences as seen by the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) instrument on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012) spacecraft (cf. Régnier et al. 2011). Here, the observed dynamics are a consequence of the condensation formation process. In this paper, we perform extreme resolution 2.5D simulations, where we are able to identify a variety of both ideal and non-ideal processes to explain similarly complex prominence-internal structure and dynamics.

2. Numerical setup

To explore the small-scale structure and dynamics involved in the formation of solar prominences, we performed a series of 2.5 dimensional (2D + invariance in the 3rd dimension) MHD simulations using the parallelised Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al. 2012, 2020; Porth et al. 2014; Xia et al. 2018) toolkit. The domain adopts a cartesian (x, y) grid, ignoring the curvature of the Sun, oriented such that the solar surface is positioned along the x axis at y = 0. The flux rope hosting the prominence was dynamically formed using the common technique of converging motions imposed within the bottom boundary so as to approximate the cumulative effect of magnetic diffusion driven by photospheric granulation (van Driel-Gesztelyi et al. 2003). In order to achieve the formation of cool condensations ab initio, we included the effects of optically thin radiative cooling, field-aligned thermal conduction, and gravity (see Xia et al. 2018, for the full description of how MPI-AMRVAC handles parallel and perpendicular thermal conduction).

This paper is split largely into two separate studies. We begin in Sect. 3.1 by expanding on the ab-initio prominence formation study of Kaneko & Yokoyama (2015). In Sect. 3.2, we make full use of the adaptive mesh refinement (AMR) capabilities of MPI-AMRVAC to explore the dynamics of prominence condensations at scales an order of magnitude smaller than ever before. Within this study, we also varied the strength of the background magnetic field so as to explore the resulting effect on prominence formation and evolution.

2.1. Governing equations

The MPI-AMRVAC toolkit is used to solve the conservative form (plasma density, momentum density, energy, and magnetic field) of the following ‘primitive’² system of non-linear, non-adiabatic MHD equations:

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

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

$$\begin{array}{cc}& \rho \frac{\partial T}{\partial t}+\rho \mathit{v}·\mathbf{\nabla }T+(\gamma -1)[p\mathbf{\nabla }·\mathit{v}+\rho \mathit{g}·\mathit{v}+\rho \mathcal{Q}\hfill \\ \hfill & -\eta {\mathit{J}}^2-\mathbf{\nabla }·(\mathit{\kappa }·\mathbf{\nabla }T)]=0,\hfill \end{array}$$(3)

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

where ρ, v, p, and B are the primitive variables for plasma density, velocity, pressure, and magnetic field, respectively, and bold font indicates a vector quantity of three components. The final primitive variable, T, is the plasma temperature defined with,

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

where ℛ and μ are the gas constant and mean molecular mass, respectively, thus closing Eqs. (1)–(4). γ is the ratio of specific heats, under the assumption of a monoatomic ideal gas this is taken as 5/3. To account for both hydrogen and helium populations, the simulated plasma is considered fully ionised with an assumed abundance ratio of 10:1 (H:He) leading to the relation ρ = 1.4m_pn_H, where m_p is the proton mass and n_H the number density of hydrogen. ρ𝒬 then represents the combined effect of optically thin radiative losses ℒ and background heating ℋ as,

$$\begin{array}{c}\hfill \rho \mathcal{Q}=\underset{\rho }{\underbrace{{\rho }^2\mathrm{\Lambda }(T)}}L-\underset{\mathcal{H}}{\underbrace{{\rho }_0^2\mathrm{\Lambda }(T_0)e^{(-y/2a)}}},\end{array}$$(6)

where Λ(T) is the radiative loss per unit mass as a function of temperature, a the vertical magnetic scale height of the simulation, and ρ₀, T₀ are, respectively, the initial density and temperature set at the bottom of the simulation to ensure the system is initiated in thermodynamic equilibrium. This formulation implicitly assumes that the contribution of coronal heating to the solar atmosphere is an exponential with the maximum taken at photospheric heights (e.g., Fang et al. 2013; Xia & Keppens 2016). κ is the thermal conductivity, typically taken as the Spitzer conductivity (Spitzer 2006) equal to κ_∥ = 8 × 10⁻⁷T^5/2 erg cm⁻¹ s⁻¹ K⁻¹ (we note that for solar coronal conditions, κ_⊥ may be considered comparatively negligible and therefore neglected from the computation). Finally, η represents the resistivity of the system and is assumed constant across the simulation domain.

2.2. Initial conditions

The 2.5D simulation domain is initialised at t = 0 with dimensions −12 ≤ x ≤ +12 Mm and 0 ≤ y ≤ 25 Mm. The density distribution across the domain is stratified according to gravity i.e., hydrostatic equilibrium, with

$$\begin{array}{c}\hfill g(y)=g_0{(r_{\odot }/(r_{\odot }+y))}^2\end{array}$$(7)

and the magnitude of gravitational acceleration at the solar surface g₀ set to 27 400 cm s⁻². The density therefore decreases exponentially with height, from a maximum value of 3.2 × 10⁻¹⁵ g cm⁻³ at −12 ≤ x ≤ +12, y = 0. The background temperature of the stratified atmosphere was set to a constant value of 1 MK. As such, the pressure stratification is fixed according to the equation of state (5). Initially, the (v_x, v_y, v_z) components of the velocity were set to 0.

Following Kaneko & Yokoyama (2015), we adopted an initial linear force-free background magnetic field distribution of,

$$\begin{array}{cc}\hfill B_x& =-\left(\frac{2L_a}{\pi a}\right)B_{a}cos\left(\frac{\pi }{2L_a}x\right)exp[-\frac{y}{a}],\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill B_y& =B_{a}sin\left(\frac{\pi }{2L_a}x\right)exp[-\frac{y}{a}],\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill B_z& =-\sqrt{1-{\left(\frac{2L_a}{\pi a}\right)}^2}{B}_{a}cos\left(\frac{\pi }{2L_a}x\right)exp[-\frac{y}{a}],\hfill \end{array}$$(10)

where L_a, a, and B_a are the lateral scale height, vertical scale height, and photospheric magnitude of the linear force free field, respectively. In all of the cases that follow, L_a = 12 Mm and a ≈ 50 Mm. Finally, in line with the simulations of Zhao et al. (2017), η was set in dimensionless units at a constant value of 0.002 across the whole domain. This choice of initial conditions ensures the simulation began in a state of near-equilibrium.

2.3. Boundary conditions

The boundary behaviour was set for the primitive (and converted to conservative) variables as described in Table 1.

Cases ‘symm’ and ‘asymm’ refer to symmetric and asymmetric boundary conditions, respectively, which use mirror reflections about the boundary to fill the ghost cells. The ⟨ρ/p⟩₀(y = 0) and (y = y_max) cases describe the enforcement of the initial values. The extrapolation indicates the use of a linear, second-order, zero-gradient extrapolation method,

$$\begin{array}{cc}\hfill \mathrm{Bottom}:& u\left(x_i\right)=\left(-u\left(x_{i+2}\right)+4u\left(x_{i+1}\right)\right)/3,\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill \mathrm{Top}:& u\left(x_i\right)=\left(-u\left(x_{i-2}\right)+4u\left(x_{i-1}\right)\right)/3.\hfill \end{array}$$(12)

Following Kaneko & Yokoyama (2015) and Zhao et al. (2017), the construction of the flux rope topology from the initial linear force-free field was achieved through a combination of converging and shearing footpoint motions imposed on the bottom boundary, and magnetic reconnection controlled using the numerical diffusion η parameter of Eq. (4). The former is prescribed according to observational theory (e.g., van Ballegooijen & Martens 1989) and practically, here, in a similar way as Kaneko & Yokoyama (2015),

$$\begin{array}{c}\hfill \begin{array}{ccc}& (0<x<L_a/4,y<0)& (L_a/4⩽x⩽L_a,y<0)\\ \left\{\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right\}=& \left\{\begin{array}{c}-v_0(t)\frac{x}{L_a/4}\\ 0\\ v_0(t)\frac{x}{L_a/4}\end{array}\right\}& \left\{\begin{array}{c}-v_0(t)\frac{L_a-x}{3L_a/4}\\ 0\\ v_0(t)\frac{L_a-x}{3L_a/4}\end{array}\right\}\end{array}\end{array}$$(13)

$$\begin{array}{c}\hfill \begin{array}{ccc}& (-L_a/4<x<0,y<0))& (-L_a⩽x⩽-L_a/4,y<0)\\ \left\{\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right\}=& \left\{\begin{array}{c}{v}_0(t)\frac{|x|}{L_a/4}\\ 0\\ -v_0(t)\frac{|x|}{L_a/4}\end{array}\right\}& \left\{\begin{array}{c}{v}_0(t)\frac{L_a-|x|}{3L_a/4}\\ 0\\ -v_0(t)\frac{L_a-|x|}{3L_a/4}\end{array}\right\}\end{array}\end{array}$$(14)

with the velocity v₀(t) a function of time,

$$\begin{array}{cc}\hfill v_0(t)=v_{00},& (0<t<t_1),\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill v_0(t)=v_{00}(t_2-t)/(t_2-t_1),& (t_1⩽t<t_2),\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill v_0(t)=0,& (t⩾t_2),\hfill \end{array}$$(17)

where t, t₁, and t₂ are the current time and the times the converging and shearing velocity begins to slow (t₁ = 1000 s) and stop (t₂ = 1500 s), respectively, with v₀₀ = 12 km s⁻¹. The near-equilibrium state achieved in the initial conditions is permitted a short relaxation period of ≈200 s to ensure that any initial relaxation effects are smoothed out before the driving motions at the boundary are imposed. Time t = 0 is defined at the moment these driving motions are initiated.

2.4. Numerical methods

The Eqs. (1)–(4) are solved using a three-step, third-order Runga-Kutta method, with a HLL approximate Riemann solver making use of the third-order slope limiter presented by Čada & Torrilhon (2009). In our explicit time-stepping procedure, we adopt a courant value of 0.8. To negate the difficulties in handling large magnetic field strengths, MPI-AMRVAC permits the treating of the magnetic field in a split way following Tanaka (1994), but extended to the linear force free case as explained in Xia et al. (2018), where the analytic equilibrium current also enters. As such, the magnetic field of Eqs. (8)–(10) was specified as a time-independent background field, and we evolved the deviation from this field in a fully non-linear MHD evolution. ∇ ⋅ B = 0 is ensured on an AMR grid by using the diffusive approach as explained in Keppens et al. (2003), wherein a source term was added in a split manner to the energy Eq. (3) to smoothly diffuse/advect any present monopole errors and reduce the signature of 0 ⪅ ∇ ⋅ B ≪ 1. We note that MPI-AMRVAC has many different means to control monopole errors discretely, such as by using source terms, or by solving an elliptic projection equation (see Teunissen & Keppens 2019), or by exploiting a constrained transport approach.

3. Analysis methods and results

3.1. Baseline 31 km resolution simulation

To establish a baseline for comparison between this work and that of previous authors, the first aim of this study is to compare and contrast the results of Kaneko & Yokoyama (2015) with the one obtained using MPI-AMRVAC under the numerical schemes described in the preceeding sections. To that end, we adopt a grid of 96 × 96 with three additional levels of refinement permitting a maximum effective grid size of ≈31 km (approximately equal to the highest-resolution case explored by Kaneko & Yokoyama 2015). We also adopt the same piecewise radiative cooling curve of Hildner (1974). We deviate from the Kaneko & Yokoyama (2015) work in two main aspects (1) the exponential heating function of Eq. (6) is chosen in line with our recent works for consistency (e.g., Fang et al. 2013; Xia & Keppens 2016; Zhao et al. 2017), (2) and we remove the minimum temperature approximation that was used in their study (likely for numerical stability reasons). The strength of the initial surface magnetic field is set to 3 G, as in Kaneko & Yokoyama (2015).

3.1.1. Global evolution

Figure 1 presents the evolution of the density, temperature, and magnetic field during the ‘build phase’ of the magnetic flux rope. Material is elevated from the bottom of the ‘atmosphere’ and collects within the flux rope alongside the reconnection process. Here, the current sheet coincides with a decrease in density and an increase in temperature at the neutral line (x = 0). The persistence of an inverted ‘Y’ configuration below the forming flux rope enables reconnection to continue to supply both plasma and flux. By the time the footpoint motions have stopped and the flux rope has formed, we see a clear difference in the distribution of temperature and density at t ≈  1500 s in comparison with t = 0.

Fig. 1. Distribution of plasma density (top row) and temperature (bottom row) within the forming flux rope. Magnetic field lines have been overplotted in white to outline the isolated flux region commonly referred to as a magnetic flux rope.

After the cessation of the converging boundary motions at t = 1300 s, the system enters a relaxation phase where the flux rope that was accelerated upwards through the atmosphere undergoes damped periodic oscillations. The system remains in this phase for approximately 6200 s. Figure 2 details how the optically thin radiative losses Λ(T) evolve during this time, indicating a strong increase in emission around 7000 s as expected when the material within the flux rope begins to cool. It should be noted that the cooling process is linked directly to the linear TI (cf. Field 1965; Xia et al. 2014; Xia & Keppens 2016; Claes & Keppens 2019; Claes et al. 2020), as in this scenario there is no TNE cycle in our 2.5D simulation.

Fig. 2. Comparison between both the distribution and the magnitude of radiative losses within the flux rope immediately after the formation of the flux rope (left) and just prior to the appearance of condensing material (right).

Shortly thereafter, and as shown in the top row of Fig. 3, a large condensation begins to form in the wake of the increased radiative losses, the associated material has an initial velocity of ≈30 km s⁻¹ at t = 7264.6 s. Around 7264.6 <  t <  7350.4 s, the shape and density of the condensation is similar to that described by Kaneko & Yokoyama (2015). In contrast to this work, we find a subsequent collapse of this dense cloud leading to smaller-scale structure. This difference can be ascribed to the artificial and fairly high temperature cut-off adopted by these authors. Indeed, in Fig. 3 much lower temperatures of order < 10⁴ K are recorded within the forming condensations. As time progresses, we see the condensation narrowing at a velocity of ≈80 km s⁻¹ to form a single, vertically oriented monolithic structure. This is reminiscent of the assumptions adopted in many radiative transfer equation (RTE) inversion methods (e.g., Beckers 1964; Heinzel & Anzer 2001; Asensio Ramos et al. 2008, 2011; Schwartz et al. 2019). The velocity of the material involved in the condensation process is up to an order of magnitude greater than that observed within post-formation prominence condensations (cf. Labrosse et al. 2010), but of the same order of magnitude as previous, fully 3D simulations of this process (e.g., Xia & Keppens 2016).

Fig. 3. Evolution of the condensing material within the formed flux rope of the 31 km resolution case. The first row details the evolution of the density, the second details the temperature, the third details the Brunt-Väisälä instability ($N_{\mathrm{m},\mathrm{p}}^2$), and the fourth details the baroclinicity (positive → clockwise rotation). Velocity quivers of relative magnitude are overplotted as white arrows to indicate the direction of, and increase in, condensation rate with time. A movie of this figure is available online.

3.1.2. The formation of a ‘monolithic’ prominence

In the presence of a gravitational field, plasma can fall in the direction of decreasing gravitational potential. A gravitationally stratified atmosphere in hydrostatic equilibrium will have a density stratification that decreases exponentially with a scale height equal to the pressure scale height H,

$$\begin{array}{c}\hfill H=\frac{k_{\mathrm{B}}T}{\mu g},\end{array}$$(18)

where k_B is the Boltzmann constant, T the gas temperature, μ the mean molecular mass, and g the acceleration due to gravity. A small packet of this gas displaced out of hydrostatic equilibrium is then said to oscillate at the Brunt-Väisälä (BV) frequency. The presence of damping will then return this displaced packet of gas to equilibrium.

In the presence of a magnetic field, there is the additional magnetic pressure (+ tension) that, depending on the magnetic topology, may levitate charged gas to altitudes far in excess of the gas pressure scale height (cf. Hillier & van Ballegooijen 2013). To be explicit, the magnetic force F_m supplied to ionised gas (plasma) may be decomposed as follows,

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{m}}=\mathit{J}\times \mathit{B}=(\mathbf{\nabla }\times \mathit{B})\times \mathit{B}=\underset{B_{\mathrm{pressure}}}{\underbrace{-\mathrm{\nabla }\left(\frac{B^2}{2}\right)}}+\underset{B_{\mathrm{tension}}}{\underbrace{(\mathit{B}·\mathbf{\nabla })\mathit{B}}},\end{array}$$(19)

wherein,

$$\begin{array}{c}\hfill B_{\mathrm{tension}}=B(\widehat{\mathit{b}}·\mathbf{\nabla })(B\widehat{\mathbf{b}})=\widehat{\mathit{b}}\widehat{\mathit{b}}·\mathrm{\nabla }\left(\frac{B^2}{2}\right)+B^2\frac{\widehat{\mathbf{n}}}{R_{\mathrm{c}}},\end{array}$$(20)

and hence magnetic pressure acts in all directions whereas magnetic tension acts to nullify the parallel force, and has a component only in the direction $\widehat{\mathit{n}}$ of the radius of curvature R_c. We note we have chosen units such that the value of the permeability of free space μ₀ = 1.

Blokland & Keppens (2011a) studied the application of this to prominences in hydrostatic equilibrium that represent the flux-rope embedded prominence stage that we reached in this simulation. Indeed, the 2.5D nature of the problem has nicely nested poloidal flux surfaces that can levitate plasma above the solar surface by magnetic pressure and tension. After the prominence forms, we reach a phase of much slower evolution, and this can be seen as a succession of 2.5D magnetostatic equilibria with prominence plasma nested along flux surfaces.

Blokland & Keppens (2011b) went on to analyse these perfectly force-balanced states, in terms of linear MHD stability. A crucial factor in this linear stability analysis is played by the flux-surface projected forms of the BV frequencies, given by,

$$\begin{array}{cc}\hfill N_{\mathrm{B},\mathrm{p}}^2& =\left[\frac{\mathit{B}·\mathbf{\nabla }p}{\rho B}\right][\frac{\mathit{B}}{\rho B}·(\mathbf{\nabla }\rho -\frac{\rho }{\gamma p}\mathbf{\nabla }p)],\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill N_{\mathrm{m},\mathrm{p}}^2& =\left[\frac{\mathit{B}·\mathbf{\nabla }p}{\rho B}\right][\frac{\mathit{B}}{\rho B}·(\mathbf{\nabla }\rho -\frac{\rho }{\gamma p+B^2}\mathbf{\nabla }p)],\hfill \end{array}$$(22)

where N_B, p and N_m, p are the traditional and magnetically modified (i.e., including magnetic pressure) BV frequencies, respectively. When (parts of) flux surfaces exist where this N² <  0, the force-balanced equilibrium state is found to be unstable, since the organised continuous parts of the MHD spectrum become overstable. The authors refer to this as driving a convective continuum instability (CCI), which is possible for those magnetostatic equilibria where the density variation in the flux rope behaves as a flux function. For our simulations here, the density is indeed virtually constant on every nested flux surface for a long time during the simulation (prior to about 7000 s), as seen in Figs. 1–3 and the online animations.

Moschou et al. (2015) studied the application of the CCI theory to the coronal rain phenomenon within a 3D numerical simulation. The authors note that the evolution of their condensations could not be directly related to the CCI as that setup was done in an arcade and hence had no nested flux surfaces, also due to lacking the required translational symmetry. Nevertheless, the spatial coincidence of N² <  0 regions with actual condensations remained consistent, concluding that the motions of the plasma were driven at all times by some variant of the CCI. Following these authors, the third row of Fig. 3 details the evolution of the magnetically modified BV frequency within a cutout of the simulation domain. At t = 7264.6 s we find that the $N_{\mathrm{m},\mathrm{p}}^2<0$ regions trace the majority of those flux surfaces on which material is condensing. The material present along these specific flux surfaces may therefore be considered CCI-unstable. After this time we find that the $N_{\mathrm{m},\mathrm{p}}^2<0$ signatures are highly structured and the plasma located in these $N_{\mathrm{m},\mathrm{p}}^2<0$ regions will slide along their own magnetic surface under the influence of gravity, until the pressure balance is restored.

During the very early stages of condensation, the plasma within a flux rope is being redistributed as a consequence of the unstable flux projected BV distribution. At any given moment, a short-lived localisation of density may occur that increases the radiative energy losses ρℒ in this location. As such, the associated plasma will cool. For plasma at a temperature of 1 MK, the cooling curve of Hildner (1974) indicates the radiative losses would increase as a result of this cooling, leading to a feedback process. However, for sufficiently small temperature perturbations, thermal conduction smooths out these localised variations. Hence, the (mostly temperature) gradient of the cooling curve must be steep enough in combination with a sufficiently large temperature perturbation. Furthermore, as before, the localised decrease in temperature is accompanied by a localised increase in density, where the ρ² dependence on the radiative losses then ensures this cooling process cannot be easily halted or reversed. Parker (1953) and Field (1965) derived conditions for this TI under isochoric and isobaric conditions, respectively. Xia et al. (2012) studied the formation of prominence condensations through the TI and noted that the temperatures and pressures of the condensations dropped both faster and earlier than their densities, as was alluded to above. Therein they find the isochoric assumption more relevant for prominence condensations than the isobaric (an identical conclusion was also reached by Moschou et al. 2015). The TI-unstable criteria C for the isochoric assumption is given by,

$$\begin{array}{c}\hfill C=k^2-\frac{1}{\kappa }(\frac{\partial \mathcal{H}}{\partial T}-\frac{\partial (\rho \mathcal{L})}{\partial T})<0,\end{array}$$(23)

where k is the wavenumber of the perturbation and κ is the coefficient of thermal conduction. These authors also noted the relative insensitivity of C to the value of k in the isochoric case. In line with previous studies, k takes the value of double the condensation size (∼2 Mm). The top panel of Fig. 4 demonstrates that the temperature decreases are responsible for driving the TI evolution (bottom panel) in the early stages of the prominence studied here. We find that the localisation of the largest negative values for C coincide with the locations of the condensations that form at a later time as shown in Fig. 3.

Fig. 4. Locations of TI within the 31 km gridstep prominence. Top: ratio of thermal to density influence on the varying energy loss 𝒬 term in Eq. (3). Bottom: spatial map of the isochoric linear stability metric C (Eq. (23)) for the thermal mode instability.

Prior to the coalescence of the final monolithic structure, local density enhancements appear to form within the general condensation (for example see the density panel at t = 7436.2 s). Considering the background, vertically stratified atmosphere at this stage, it is clear that denser material finds itself positioned above less dense material. In hydrodynamics, such density inversions may become unstable to the RTI. The behaviour of material related to the RTI is then typically quantified by,

$$\begin{array}{c}\hfill \mathcal{A}=\frac{{\rho }_{+}-{\rho }_{-}}{{\rho }_{+}+{\rho }_{-}},\end{array}$$(24)

where 𝒜 is the Atwood number and ρ₊ and ρ₋ are the more and less dense regions, respectively. For values of 𝒜 ⪆ 0, we find that the mixing is somewhat symmetrical and the less dense plasma penetrates into the denser plasma to an equal degree as the reverse. For 𝒜 ≈ 1, the more dense plasma penetrates further into the less dense region with relatively little of the reverse (Chandrasekhar 1955). Moschou et al. (2015) explored the distribution of 𝒜 within their 3D simulations of coronal rain, calculated as a search for a density gradient in the y-direction of their domain. However, the 𝒜 parameter quantifies the density stratification in only the y-direction with particular emphasis on the locations of density gradients. The motions of plasma subject to the RTI are also dependent on the pressure gradient that the evolving plasma encounters and, in fact, the misalignment between these pressure and density gradients leads to the vortices commonly observed in the wake of the main falling fingers of RTI-unstable material. However, the 2.5D configuration used here is oriented across the flux loops of the flux rope and hence these vortices are heavily stabilised by the magnetic tension within the topological dips. As such, we do not expect to observe similar large-scale structures to those recorded in previous simulations of prominence evolution involving B oriented across the 2.5D plane, or in full 3D (e.g., Hillier et al. 2012a,b; Terradas et al. 2015; Keppens et al. 2015; Kaneko & Yokoyama 2018; Popescu Braileanu et al. 2021). Nevertheless, we expect the misalignments between the pressure and density gradients to be present in 2.5D and therein induce motions inline with the RTI but on scales of the order of a flux tube. Evolution of the plasma in this way may also be characterised using baroclinicity, a component of the evolving vorticity ω that arises in a flow with anisotropic temperature,

$$\begin{array}{c}\hfill \frac{\partial \mathit{\omega }}{\partial t}+(\mathit{v}·\mathbf{\nabla })\mathit{\omega }=(\mathit{\omega }·\mathbf{\nabla })\mathit{v}-\mathit{\omega }(\mathbf{\nabla }·\mathit{v})+\underset{\mathrm{baroclinicity}}{\underbrace{\frac{1}{{\rho }^2}\mathbf{\nabla }\rho \times \mathbf{\nabla }p}}.\end{array}$$(25)

This baroclinicity considers that a misalignment between the pressure and density gradients gives rise to the rotation of plasma as it slides down a constant-pressure gradient. In general, plasma unstable to the RTI in a direction other than vertical may be identified using this approach. The fourth row of Fig. 3 details the evolution of the baroclinicity within the condensing material. At t = 7264.6 s, the signature of baroclinicity covers a broad region but one that is focussed at the periphery of the condensation region and forms a ‘U’ shape. As the evolution progresses, this ‘U’ shape is amplified in magnitude, but also becomes flat and increasingly localised to the bottom of the forming monolithic structure.

In the density panels, at t = 7436.2 s we see that the two early condensations develop a curvature as they move towards the concave-up portions of the field topology. This is consistent with the baroclinicity-induced rotation. It is important to emphasise here that vorticity maps quantify, generally, the rotation of the velocity field (induced by some physical process), whereas the baroclinicity maps show the source of vorticity and represent a physical process in itself. Hence, maps of the vorticity would likely also show this rotation but would not offer information on the driving process.

3.1.3. Evolution of the monolithic prominence

Figure 5 details the evolution of the now-formed prominence after the initial condensation period. Initially, between 7607.8 <  t <  8294.2 s, it can be seen that the condensation region significantly deforms the field topology of the host flux rope as the concentration of highest density falls from an altitude of ≈9–6 Mm with an average velocity of 5 km s⁻¹. The largest magnetic field strength is recorded within the condensation at this time, with a value of ≈4.2 G. The falling material drags the associated field lines down and those field lines that lie directly below the condensation are compressed. Shortly after t = 8294.2 s, reconnection is triggered just above the densest portion of the falling monolith, creating a magnetically isolated region of high density that continues to fall, now with a velocity of ≈10 km s⁻¹. Behind this falling ‘plasmoid’³, another high-density region forms in its wake, and at t ≈ 8500 s reconnection is again triggered behind this new region. Once the new density enhancement is isolated by the reconnection, it begins to fall much faster (v ≈ 7 km s⁻¹) and eventually approaches and coalesces with the first falling plasmoid, which has slowed to ≈5 km s⁻¹ by this time due to increasing magnetic tension see the t = 8809.0 s panel of Fig. 3).

Fig. 5. Evolution of the monolithic prominence once formed by catastrophic, optically thin radiative cooling. Top row: evolution of the density with magnetic field lines indicated using white streamlines. Middle row: value of plasma β throughout the flux rope and prominence. Bottom row: evolution of the baroclinicity.

The separations of the plasmoids from the main monolithic, prominence condensation are comparable to AIA observations reported by Liu et al. (2012a). The authors related the apparent evolution of dynamic prominence plasma observed in the cooler 304 Å passband to the RTI. Indeed, the simulation here and that of Kaneko & Yokoyama (2018) contain plasma stratifications that are remarkably similar to these observations, albeit differing from the small-scale dynamics produced in dedicated simulations of the RTI in prominences (e.g., Hillier et al. 2012b; Keppens et al. 2015). To explore this, the bottom row of panels in Fig. 5 shows the evolution (7607.8 <  t <  8809.0 s) in the distribution of baroclinicity within a cut-out field of view (FOV) covering only the region of the plasmoid motions.

As in Fig. 3, the strong signatures of baroclinicity are focussed at the edges of both the main condensation and the falling plasmoids. The predominantly negative-left and positive-right signatures on either side of the main condensation are caused by a relatively strong (weak) horizontal (vertical) pressure gradient across the condensation boundary (recognise that ∇ρ × ∇p = [∂_xρ∂_yp − ∂_yρ∂_xp]e_z, so we always show the out-of-plane component in figures). The same is true for the left- and right-sides of the falling plasmoids, the relative importance of this is discussed in more detail in Sect. 4.

According to Low et al. (2012a), the initial compression of the magnetic field by the heavy prominence should create a large but finite current between the flux surfaces that, through dissipation, allows material to pass from one field line to another in the direction of gravity (cf. their Fig. 8). Hence, material may slip across field lines, a theory invoked to offer a potential explanation to why quiescent prominences projected above the limb appear to have a more vertical structuring than their on-disk, seemingly horizontal counterparts (see also Low et al. 2012b). Low & Egan (2014) then constructed a 1D steady-state model of a Kippenhahn-Schlüter-type prominence slab and demonstrated that under the action of gravity, a current layer was formed accompanied by a local velocity enhancement indicating the motion of material across the field lines.

To be more explicit, one of the cases explored by Low et al. (2012b) considered a temperature variation across those flux surfaces being compressed by the heavy prominence material. These authors also briefly considered that the compression may be significant enough for two adjacent field lines to be brought very close together in space. The perpendicular thermal conduction may therein also become important. However, the authors also note that the ratio ϵ between the perpendicular thermal conductivity κ_⊥ and resistivity η has, of course, the limits of ϵ ≪ 1 and 1 ≪ ϵ, corresponding to resistive or conductive behaviour, respectively. In our simulation, the perpendicular thermal conductivity is assumed to be negligible and therein completely neglected, leading to ϵ = 0 everywhere. However, they also found that ϵ = 0.64 × (plasma − )β and thus provides a quick check for whether our assumption remains valid within the region between these compressed flux surfaces. The distribution of β within the simulation domain is presented in the middle row of Fig. 5. It is shown that the value of β is of the order of one or less in the entire simulation domain, including the condensation region, and ≪1 within the remainder of the evacuated flux rope. At the interface between the falling condensations and the lower, compressed field we find localised regions of β ≈ 0.1. Hence, ϵ may be considered (much) less than one at the positions of these interfaces and hence resistive behaviour dominates over conductive.

In the left panel of Fig. 6 we show a snapshot at t = 8294.2 s of the density and orientation of the magnetic field lines within the region of the falling material (cf. Low et al. 2012a). The panels on the right then detail the evolution (7522.0 <  t <  8294.2 s) of the density, current density, vertical velocity, and temperature along a cut taken at x = 0. The four cuts are presented as zoom-ins on the bow i.e., the leading edge, of the falling material, indicated by the increase/decrease in density/temperature by several orders of magnitude. At t = 7522.0 s, the current density is enhanced at ≈9 Mm as the local field lines begin to be compressed ahead of the growing condensation. By t = 7779.4 s we see that the vertical velocity of the falling monolithic structure has decreased and the current enhancement is now more localised to the bow of the falling material. As the simulation progresses further, a small and localised negative velocity enhancement is observed to grow just below the current enhancement.

Fig. 6. Magnetic field and density map, and the evolution of plasma density, current density, vertical velocity, and temperature along a cut taken vertically through x = 0. Left: zoom-in of the region containing the bottom of the prominence condensation. Magnetic field lines are plotted in red with density contours overlaid according to the colourbar. Right: four panels detailing the evolution of the aforementioned physical parameters across the leading edge of the falling condensation, times are indicated as titles to each plot.

Figure 7 then explores the same quantities as Fig. 6 but for a later time t = 8654.56 s. At this time, and as can be seen in the left panel of this figure, the density concentrations have formed three distinct fragments or isolated islands. The top panel then details the variation in the density, current density, vertical velocity, and temperature in a cut from 0 Mm to 8 Mm in height at x = 0. The bottom three panels of Fig. 7 focus on the leading edges of each distinct density concentration, the lowest referred to with a (1), the middle with a (2), and the highest with a (3).

Fig. 7. Magnetic field and density map, and stratification of plasma density, current density, vertical velocity, and temperature along a cut taken vertically through x = 0. Left: a zoom-in of the region containing the multiple prominence condensations. Magnetic field lines are plotted in red with density contours overlaid according to the colourbar. Top: stratification of the aforementioned plasma parameters between 2–8 Mm, isolating the three prominence consensations. Bottom: three panels focussing on the distribution of parameters across the cusp of each blob.

3.2. 5.7 km resolution and a continous cooling curve

3.2.1. Modified setup

In this section we revisit the formation and evolution of prominences by using a much higher numerical resolution, and for two different field strengths. To that end, we make full use of the MPI-AMRVAC toolkit and adopt a base grid of 528 × 528 with three additional levels of refinement permitting a maximum effective gridcell size of ≈5.7 km. We also improve on our assumption of optically thin radiative losses by exchanging the piecewise cooling curve of Hildner (1974) for the continuous solar metallicity cooling curve of Schure et al. (2009), wherein the lower temperatures are governed by the behaviour of H I (see Dalgarno & McCray 1972, and the MPI-AMRVAC website⁴). For the majority of the simulation time, the Eqs. (1)–(4) are solved using the same numerical methods as described in Sect. 2.4. However, we also adopt a more stable two-step, second-order Total Variation Diminishing Lax-Friedrich (TVDLF) method, with a monotonised central (or ‘woodward’) limiter. We also artificially enforce an overall minimum temperature of 1000 K, far below the anticipated temperature of the prominence condensations to ensure, at all times, positive pressure. Similarly, we solve an auxiliary internal energy equation in tandem with the total energy equation such that we can guarantee positive pressures; in very low-β plasmas the numerical error for total energy is of the order of the internal energy and hence the total energy may become negative. For β <  0.005 we trust the solution to the internal energy, apply a linear combination of internal and total energy for 0.005 <  β <  0.05, and simply continue with total energy for β >  0.05. The two surface magnetic field cases explored in this section adopt a B_a value of 3 G and 10 G, labelled ‘3G’ and ‘10G’, respectively. Finally, we choose to break the artificial left-right symmetry adopted in Sect. 3.1 with the imposition of an arbitrary perturbation to both the 3G and 10G cases at t = 1945 s of the form v_x = −v_maxsin(πy/25), where v_max = 12 km s⁻¹.

3.2.2. Prominence formation and evolution

In each case 3G and 10G, the formation of the flux rope follows the previous example wherein the linear, force-free coronal magnetic field initially relaxes before the bottom boundary drives their footpoints together. The initiation of reconnection at the point x = 0 then builds the twisted magnetic field that is lifted above the surface as reconnection supplies flux to the forming rope. The increased resolution also allows us to resolve further reconnection-related fine-structure, formed within the current sheet underneath the main flux rope, and in both 3G and 10G cases, additional subordinate flux ropes are formed in this sheet.

The top panels of Fig. 8 capture these ropes after formation and during their launch into the main flux rope, wherein they are later incorporated (such dynamic mergers have previously been recorded in high-resolution numerical simulations e.g., Zhao et al. 2017, 2019). The particularly large dimensions of these ropes in comparison with previous studies is influenced by the fairly large η value of 0.002 across all simulations within this work. As such, reconnection progresses much more easily than expected on the actual Sun. Nevertheless, this is analagous to the observations of ‘flux feeding’ suggested to play a role in the mass-supply of prominences or in CME eruptions (e.g., Liu et al. 2012b; Zhang et al. 2014, 2020). The lower panels of Fig. 8 detail the distribution of electric current density (into the plane i.e., J_z) and mass density along a cut taken vertically at x = 0.

Fig. 8. Production of subordinate flux ropes below the main rope and their involvement with the supply of material to the prominence. Top panels: distribution of Jz current density. Bottom panels: current and plasma density distribution in a vertical cut through the current sheet at x = 0. Left and right: construction of two additional flux ropes (t = 915.4, 1215.7 and t = 735.2, 1001.2 s) within the 3G and 10G cases, respectively.

The flux rope formation process provides a non-negligible upward acceleration to the flux rope and trapped plasma. Once the boundary driving and reconnection has ceased, the formed flux ropes begin oscillating vertically with a damped amplitude. For the 10G case, the oscillation has a period of ≈350 s and persists for only a short while whereas the 3G case has a longer period of ≈900 s and continues to oscillate for an extended period due to the relatively weaker magnetic field and restoring tension force therein. The symmetry is then broken in line with the aforementioned perturbation applied at t = 1945 s. Shortly thereafter, the temperature within the core of the 3G flux rope begins to decrease rapidly and condensations form quickly (t ≈ 2743 s). For the 10G case, the larger magnetic field strength/current and associated Joule heating ηJ² maintains the temperature of the plasma for an extended period of time until finally the temperature in one of the outer layers of the flux rope drops at t ≈ 5162 s, accompanied by the formation of condensations. In the following, we apply the same analysis as explored in Sect. 3.1.

Figure 9 presents the evolution of plasma density and temperature within the 3G case. The core of the flux rope is the first region to cool as it is the first region to be disconnected from the surface; thermal conduction is no longer transporting energy from the bottom boundary where the magnitude of the heating function ℋ is largest. Indeed, the flux addition to the flux rope appears to be mimicked in the ordering of runaway cooling from the core to the outer-edges of the flux rope. Nearly identical to Fig. 4, we find that maps of ∂𝒬/∂T ÷ ∂𝒬/∂ρ and C for the 3G case indicate the decreasing temperature of plasma is driving the TI from the core of the flux rope to its outer edges with time. At the time of initial condensation, the symmetry-breaking perturbation of the flux rope has not completely damped out and so the forming condensations carry some of this momentum. Nevertheless, the evolution of all cooled plasma is ubiquitously towards the concave-up and dipped portions of the flux rope topology. As in Sect. 3.1, we find velocities of material entering the condensations as high as ≈90 km s⁻¹. Some packets of condensed plasma inherit significant momentum from the host rope. In combination with gravity these packets are observed to quickly slide to the bottom and are thrown back up to higher heights within the flux rope with velocities of the order 35 km s⁻¹ as they pass the topological dips. The condensing material displays a variety of ‘sloshing’ and oscillatory motions as its kinetic energy is damped and eventually the material comes to near-rest at the bottom of the rope. Interestingly, throughout the evolution of the system we see material condensing at locations throughout the flux rope i.e., not just as a continuous process but rather as isolated condensations even at higher heights (see, for example, the t = 5377 s panel and the localised temperature decreases at h ≈ 10 and 21 Mm and compare with Kaneko & Yokoyama 2015).

Fig. 9. Evolution of the condensing material within the 3G formed flux rope. First row: evolution of the density, second: temperature. Velocity quivers of relative magnitude are overplotted as before. A movie of this figure, including $N_{\mathrm{m},\mathrm{p}}^2$ and baroclinicity maps, is available online.

At t = 3661 s we can see that the condensing material initially cools the rope, almost in its entirety, down to an average temperature of between 10 000–100 000 K. Once the material begins to collect in the concave-up dips of the magnetic flux rope (t = 4519 s), we see that the background temperature of the rope returns to ≈1 MK. At the later time of t = 5377 s, the core of the flux rope has begun to heat to ≈4 MK. Towards the end of the condensation process, unlike as observed with the configuration in Sect. 3.1, the whole prominence begins to fall with an average velocity of ≈2 km s⁻¹, likely due to the cross-field slippage and reconnection in its wake. It is this reconnection that is responsible for heating the core of the 3G flux rope seen in the snapshot at t = 5377 s. The peak magnetic field strength is recorded within the condensed plasma portion of the 3G flux rope with a value of ≈4.1 G.

Comparing the distribution of the $N_{\mathrm{m},\mathrm{p}}^2<0$ locations with regions of high baroclinicity, the $N_{\mathrm{m},\mathrm{p}}^2<0$ solutions appear to be distributed throughout both the condensing region and the surrounding flux rope whereas the baroclinitity signatures appear localised to the condensing region itself (see the online movie associated with Fig. 9). At later times, the $N_{\mathrm{m},\mathrm{p}}^2<0$ signatures remain largely spread out whereas the baroclinicity has become localised once more to the positions of the prominence condensations.

Towards the end of the simulation, the condensations of the 3G case coalesce at the bottom of the magnetic topology and form once more (cf. the lower-resolution case) a monolithic structure. Recall that we chose to break the symmetry by applying a perturbation. At the time of the formation of the monolithic structure (see the t = 5377.0 s snapshot of Fig. 9), residual yet small-scale oscillations persist, requiring delicate consideration when applying the same analysis as presented in Figs. 6 and 7. Figure 10 presents the cut of plasma density, current density, vertical velocity, and temperature along a vertical cut at x = 0 at a time when the bulk oscillations in both the vertical and horizontal directions are minimal. The complete monolithic structure is captured at this time, clearly indicated in the middle panel with the increase (decrease) in density (temperature) between 4–9 Mm. The right panel then focusses on the bottom 3.95–4.05 Mm of the condensation wherein the field lines are compressed as a result of the plasma weight. We note here that the compression appears to be to a lesser degree than the case explored in Sect. 3.1, nevertheless a clear current enhancement is also present at the cusp in addition to the local velocity enhancement in the negative vertical direction. The increase in the resolution of the simulation in the 3G case has enabled the local enhancements in both current and vertical velocity to be better resolved than in Figs. 6 or 7.

Fig. 10. Magnetic field and density map, and stratification of plasma density, current density, vertical velocity, and temperature within the 3G case along a cut taken vertically through x = 0 at t = 5136.76 s. Left: a zoom in of the bottom of the prominence condensation. Magnetic field lines are plotted in red with density contours overlaid at thresholds according to the colourbar. Middle: stratification of the aforementioned plasma parameters between 0–12.5 Mm, isolating the entire prominence condensation. Right: a zoom-in of the cusp of the condensation, showing the behaviour of the plasma parameters across the boundary.

Figure 11 presents the evolution of plasma density and temperature within the 10G case. The core of the 10G flux rope is significantly hotter than the 3G case (up to 5 MK), owed to the larger magnetic field driving enhanced Joule heating. The large currents involved in the formation process i.e., the subordinate flux ropes presented in Fig. 8, provide significant current to a shell around the core. Hence, both of these regions remain hotter than the corresponding regions in the 3G case. This has the primary consequence of preventing these regions from cooling and hence condensations do not begin in the core but instead in one of the non-Joule-heated outer layers. The structure and shape of the initial condensations within the 3G case appear to be more global before becoming localised, whereas the initial condensation in the 10G case begins far more localised. The formation of the condensation starts to collect material along a flux surface with a velocity of order 70 km s⁻¹, before gravity accelerates it to an average velocity of 25 km s⁻¹ towards lower heights. Most strikingly, the orientation of this initial filamentary condensation is oriented perpendicular to the host magnetic field, a zoom-in is presented in the top panel of Fig. 12. At t = 6878.5 s, we see that two additional condensations are in the process of falling to the concave-up portions of the magnetic topology after having formed at two very different heights. The first and lowest of these two additional condensations forms at a similar height to the first condensation at t = 5162 s, whereas the second occurs at a far more elevated height of ≈20 Mm, before they both fall to the bottom of their flux surface. Both of these condensations have a contraction velocity of order 90 km s⁻¹ (driven by varations in (thermal/ram) pressure cf. Claes et al. 2020). The lower of the two condensations falls with a maximum velocity of ≈20 km s⁻¹. Interestingly, the initial contraction and subsequent fall of this lower condensation leaves a very low pressure region in its wake. The higher condensation now falling under gravity is also accelerated by this low pressure void, reaching velocities as high as 70–80 km s⁻¹. As was found for the 3G case the peak magnetic field strength is recorded within the condensed plasma portion of the 10G flux rope, with a value of ≈10.7 G.

Fig. 11. Evolution of the condensing material within the 10G formed flux rope. First row: evolution of the density, second: temperature. Velocity quivers of relative magnitude are overplotted as before. A movie of this figure, including $N_{\mathrm{m},\mathrm{p}}^2$ and baroclinicity maps, is available online.

Figure 13 presents the maps of ∂𝒬/∂T ÷ ∂𝒬/∂ρ and C (wherein k = 0.1 Mm, cf. Fig. 12) for the 10G case at times just prior to the two episodes of condensation formation (t = 5134.34 and 6629.68 s). The maps of the TI, constructed as maps of the isochoric C metric, elegantly isolate the exact locations (both position and extent) of the forming condensations. Furthermore, and consistent with the other cases explored in the paper, we find that the early stage of the TI within the condensing plasma of the 10G case is driven by the temperature variations. In fact, the resulting filamentary condensations of order ρ = 10⁻¹³ g cm⁻³, T = 10 000 K, p = 0.2 dyn cm⁻² are directly comparable to the density structures formed by MHD slow-wave induced TI modes (this will be discussed in the Sect. 4 but see Claes & Keppens 2019; Claes et al. 2020).

Fig. 12. Density, temperature, and pressure cross-section of the initial filamentary condensation within the 10G prominence formation case. Top panel: distribution of density within the forming filamentary condensation. Velocity quivers of relative magnitude are overlaid in white. Magnetic field lines are overlaid as solid-red lines. Bottom panel: density, temperature, and pressure profiles across the condensation cross-section, taken along the dashed-black line in the top panel. Profiles computed using yt (yt-project: Turk et al. 2011).

Finally, those locations evolving with $N_{\mathrm{m},\mathrm{p}}^2<0$ are once again distributed throughout the flux rope as noted previously for the 3G case, including those nested flux surfaces that lack discrete condensations (cf. CCI). For the large baroclinicity values (and variations therein), one finds satisfactory agreement with the locations of the prominence condensations (see the online movie associated with Fig. 11).

3.2.3. Visualising small-scale evolutions

The increased resolution employed in the 3G and 10G cases lead to complex non-linear evolutions. Figures 9 and 11 are accompanied by movies in the online version of this paper. We include some Hydrogen-α (H-α) visualisations of the evolution following the methods of Heinzel et al. (2015) (see the explicit implementation in the appendix of Claes et al. 2020). For all visualisations, the prominence material is synthesised using the 10 Mm table of Table 1 in Heinzel et al. (2015).

Figure 14 presents the (H-α) appearance of the forming and evolving 3G and 10G cases from both a top-down (filament) and side-on (prominence) perspective. The flux rope formation period 0 <  t <  1516 s, associated oscillations, induced perturbation t = 1945 s, and their damping can be visually identified as intensity variations within each visualisation. During the formation phase of the 10G case, the corresponding filament displays longitudinal oscillations across the width of the flux rope indicated by the periodic broadening of its cross section between 1000 <  t <  1800 s. The transverse, horizontal, and persistent oscillations induced by the symmetry-breaking perturbation are well captured within the filament visualisation of the 3G case. The transverse, vertically oriented oscillations induced by this formation process are present within the prominence visualisations of either case. To emphasise these vertical oscillations and their relation to the forming condensations, the position of each flux rope’s axis (hereafter ‘O-point’) are overplotted on the prominence visualisations. In each case, the formation and rise of the flux rope is captured along with the resulting vertical oscillations that are weakly damped in the 3G case, or strongly so in the 10G case. The sudden jump in the height of the 10G O-point is a result of the injection of the larger of two subordinate ropes presented in Fig. 8. In addition, we have overplotted the time evolution of total mass per unit length within the forming prominences, defined as the mass with a temperature below T = 10⁵ K.

4. Discussion

4.1. What drives the evolution of prominences in the early stages of their formation

This paper focusses on the mechanisms for, and drivers of, the formation and evolution of prominences within the solar atmosphere through numerical simulations. The host magnetic flux rope is constructed in-line with the assumption of photospheric diffusion and magnetic cancellation at the intersection of two separate flux regions of opposite polarity (cf. Martin et al. 1985; van Ballegooijen & Martens 1989; van Driel-Gesztelyi et al. 2003). The footpoints are converged and sheared by imposing a velocity on the lower boundary following Kaneko & Yokoyama (2015) to investigate the levitation-condensation process. The result, as detailed in Fig. 1, is a flux rope with a left-handed twist that is progressively elevated through the solar atmosphere as a consequence of continued reconnection. In all cases explored here, the formation process leads to a vertical oscillation of the constructed flux ropes, as recorded in Fig. 14. The 31 km gridcell and 3G cases experience a prolonged period of transverse oscillations with larger amplitude as their magnetic field is relatively weak in comparison with the 10G case that displays much smaller amplitude oscillations that are quickly damped.

The first case explored in Sect. 3.1 was the 31 km grid step study under symmetric conditions and assuming the piecewise cooling curve of Hildner (1974). The resulting condensations appear symmetrically on either side of x = 0 and with a velocity vector oriented towards the concave-up dips of the magnetic topology. In line with the conclusions of Kippenhahn & Schlüter (1957), the cooling plasma appears to seek equilibrium where it can be supported against gravity by the tension of the magnetic field. The appearance of fine structure at this time, which was not reported in the study of Kaneko & Yokoyama (2015), is likely due to our more realistic handling of the low temperature limit. Indeed, they assumed a floor temperature of 0.3MK, far above the temperatures measured within prominences and more indicative of the cool yet diffuse plasma of the lower corona. The gradual evolution of this condensing, finely structured plasma towards the Kippenhahn & Schlüter (1957) equilibrium is actually influenced by a plethora of plasma instabilities. The interesting question is then of course which is (or which are) the most relevant.

Figures 3–7 focus on the physical processes driving the evolution. According to Fig. 3, the formation of the prominence is the end result of various interacting physical processes. Once the flux rope has formed, the density of the plasma within the flux rope can be considered a flux function i.e., constant along a given flux surface but varied with distance from the flux rope axis. Small density variations along a given flux surface caused by e.g., perturbations are then said to be CCI(-like) unstable, causing a redistribution of plasma to restore equilbrium. This is well-captured in the maps of $N_{\mathrm{m},\mathrm{p}}^2<0$ prior to the appearance of discrete condensations. The maps of isochoric C presented in Fig. 4 indicate the plasma in the core of the flux rope then undergoes the TI, driven by the temperature in this region decreasing due to ‘runaway’ radiative losses (see Fig. 2 and cf. Parker 1953; Field 1965; Klimchuk 2019; Antolin 2020). The subsequent formation of dense and cool regions additionally enforces the TI due to the ρ² dependence of ρℒ, but the now-large variations of plasma density along a given flux surface halts any valid comparison to the CCI. Nevertheless, the $N_{\mathrm{m},\mathrm{p}}^2<0$ metric remains valid for evaluating regions within the flux rope that are evolving as a consequence of being out of pressure balance. As expected, the field-projected pressure scale height responsible for enforcing hydrostatic equilibrium along a flux surface appears to shrink as the material dramatically cools (leading to the appearance of $N_{\mathrm{m},\mathrm{p}}^2<0$ regions within the flux rope as seen in Fig. 9). This out-of-equilibrium plasma then falls along field lines in the direction of the shrinking scale height i.e., towards the concave-up portions of the magnetic topology. At time t = 7264.6 s, the material within the outermost regions of the condensing region (compare with the distribution of optically thin radiative losses in Fig. 2) appears to undergo rotational motions as expected to occur when local baroclinicity is analysed. As the simulation evolves, this baroclinicity is then increasingly concentrated at the underside of the monolithic structure. The same may be said for the regions of $N_{\mathrm{m},\mathrm{p}}^2<0$, however additional regions remain present elsewhere and at higher heights within the condensing region. Therefore, it appears that within this first setup: The CCI is largely responsible for the evolution of plasma prior to the appearance of condensations, the TI drives the stable formation of condensations, and these condensations fall to the concave-up portions of the magnetic topology as a consequence of being out of pressure balance on a given flux surface. The enhanced signatures of baroclinicity located within the condensing material then indicates that the material evolving due to being out of pressure balance is also rotating. Although baroclinic signatures are expected in the wake of material evolving under the RTI, a direct comparison is not possible here due to the orientation of the 2.5D domain with respect to the magnetic field vector. However, the explicit coincidence between baroclinic signatures and the condensed material indicates such rotations may be identifiable within particularly high-resolution observations of prominences. The final result, as shown in the t = 7607.8 s snapshot of Fig. 3, is the formation of a monolothic prominence embedded within a density-depleted and cold coronal cavity.

Following the formation of the large-scale, coherent prominence, discrete blobs of plasma i.e., plasmoids are observed to fall from its underside, eventually developing pinchpoints in their wake and splitting off entirely after several episodes of reconnection (see the top panels of Fig. 5 and cf. Liu et al. 2012a). In order for such dynamics to occur, the gravitational forces have to exceed the magnetic tension forces so as to sufficiently drag the host field lines towards the bottom of the simulation domain. Indeed, we see from the second row of Fig. 5 that the values of β, taken as a first order approximation to this effect, within the falling plasmoids were of the order, or greater, than one. The bottom row of the same figure then details the distribution of baroclinicity wherein the negative-left and positive-right signatures at the edges of the monolithic structure and falling plasmoids can be attributed to the shape of the general pressure gradient within the flux rope. Then, for the falling plasmoids, theory suggests that clockwise (positive) and anticlockwise (negative) rotations are expected to be present on the left and right sides of structures falling and evolving under the influence of gravity. Indeed, smaller structures can be identified at the very periphery of the falling blobs. This suggests that although this 2.5D simulation does not contain distinct signs of the large-scale ‘roll-up’ of plasma involved in the RTI, there exists rotational motions at the edges of these falling plasmoids (cf. Hillier et al. 2012b).

Figures 6 and 7 focus on the distribution and evolution of plasma and magnetic parameters within these plasmoids in relation to the theory proposed by Low et al. (2012b). At all times explored in Figs. 6 and 7, the leading edge of the monolith and Blob (1), respectively, display a current enhancement consistent with the compression of field lines ahead of them. As this current enhancement grows, we see the localisation of a negative velocity enhancement that appears to bridge those cold flux surfaces containing the high-density prominence plasma with those below that are both hot and of low density. By contrast, the 1D study of this process by Low & Egan (2014) predicted the current enhancement would be co-spatial with the velocity enhancement, whereas here there is a clear displacement between the two signatures, with the velocity signature always leading the current signature. Blobs (2) and (3) then contain the same current enhancement but the velocity signature is significantly positive in comparison. Blobs (2) and (3) each exist above a reconnection point and in the wake of the blob below it. As such, it is conceivable that the reconnection process would provide a non-negligible positive velocity component between each blob as the field associated with the reconnection relaxes upwards. Hence, the cross-field diffusion described by Low et al. (2012a,b), if occuring for Blobs (2) and (3), may be supressed in these regions or simply be of a lengthscale below the resolution of the numerical grid (< 62 km) in this case. Crucially, however, the magnitude of the negative velocity measured ahead of the falling monolithic structure and Blob (1) is of the order and less than the velocity of each structure itself.

To further explore the initial results and conclusions of Sect. 3.1, in Sect. 3.2 we moved to consider two higher-resolution (5.7 km resolution), continuous cooling curve, symmetry-broken simulations of prominence formation and evolution proccesses above either a 3 or 10 G photospheric magnetic field. The general evolution of these 3G and 10G cases are presented in Figs. 9 and 11, respectively, with (H-α) visualisations presented in Fig. 14.

The initial formation process for the 3G case is driven in the same way as the lower-resolution case from Sect. 3.1, and indeed the resulting evolution is nearly identical with the exception of the ‘flux-feeding’ structures of Fig. 8. The global flux rope structure also behaves in much the same way, experiencing long period, weakly damped vertical and additional horizontal oscillations in response to both the formation dynamics and the symmetry-breaking perturbations (see Fig. 14 and cf. Molowny-Horas et al. 1999; Terradas et al. 2002). Furthermore, the catastrophic cooling within the 3G case is also concentrated first at the core of the flux rope before spreading out to the outer edges with cool condensations forming in its wake. As such, these condensations initially form at almost all locations within the rope before falling to the concave-up portions of the magnetic topology. It is important to note here that the combination of the realistic cooling curve (Schure et al. 2009) and the symmetry-breaking perturbation lead to the formation of these condensations at a much earlier time than the lower resolution, Hildner (1974) cooling curve, and symmetric setup.

The evolution of the prominence in the 10G case differs significantly from the two previous cases involving a 3 G photospheric field. First and foremost, the increase in magnetic field leads to an order of magnitude increase in the magnetic energy within the simulation. As such, the formation of the host flux rope occurs faster, and the associated oscillations induced by this process have both a much smaller amplitude and are damped much faster, see the O-point evolution in the bottom-right panel of Fig. 14. Then, by direct comparison with the 3G case, the condensations appear delayed in the 10G case. Indeed, thermal conduction is more efficient along those field lines of higher magnitude and as such any temperature perturbations that may therein lead to condensations are more easily stabilised. In the absence of any additional physical drivers of these condensations, a significant portion of an entire flux surface would be required to cool enough for catastrophic runaway cooling to initiate locally. Despite this, the size and shape of the initial condensations in the 10G case are far smaller and more localised than for the 3G case as shown in the TI maps of Fig. 13.

Fig. 13. Locations of TI within the 10G prominence just prior to the two main episodes of condensation formation (left: 5134.34 s and right: 6629.68 s). Top: ratio of thermal to density influence on the varying energy loss 𝒬 term in Eq. (3). Bottom: spatial map of the isochoric linear stability metric C (Eq. (23)) for the thermal mode instability.

Upon the eventual formation of cool condensations within the 10G prominence, the characteristic size and shape of the initial condensations are far smaller and more localised than those found in the 3G case. Figure 12 presents a zoomed-in view of the first condensation to form in the 10G case, in which it is immediately clear that multiple flux surfaces are involved across the length of this highly localised condensation. Indeed, this figure facilitates a direct comparison with those filamentary condensations previously found in the studies by Claes & Keppens (2019), Claes et al. (2020). Although these authors adopted a highly idealised setup, their work demonstrated that the selective excitation of the thermal MHD mode results in the formation of these structures. Furthermore and more generally, the authors suggested that these features may be formed as a natural consequence of radiative cooling + any perturbation e.g., the perturbation experienced by the entire flux rope shortly after formation, or the one applied to deliberately break the symmetry at t = 1945 s. Of course, the real Sun is highly dynamic and perturbations from MHD waves occur ubiquitously, therein suggesting the ubiquitous formation of thermally driven condensations.

We find that the CCI, TI, $N_{\mathrm{m},\mathrm{p}}^2<0$, and baroclinicity metrics, when applied to the 3G and 10G cases, indicate the identical relationship between the individual metrics and the initial stages of prominences’ evolution as was found in the 31 km gridstep setup. The CCI maps regions within the flux ropes wherein plasma is evolving slowly in response to changing pressure gradients induced by the global oscillations traced in Fig. 14. The isochoric C value of the TI is negative and large in locations where condensations are in the process of forming, particularly striking examples are presented in Fig. 13 for the highly localised condensations of the 10G case. After the appearance of distinct and discrete condensations, the flux ropes become filled with localised regions of $N_{\mathrm{m},\mathrm{p}}^2<0$ indicating plasma evolving throughout the flux rope due to varying temperature and associated pressure gradients. The condensations themselves are then cospatial with strong baroclinic signatures revealing that even those condensations on the smallest scales are rotating in tandem with evolving towards the concave-up dips of the magnetic topology.

Fig. 14. Filament (left) and prominence (H-α) visualisations (right) of the 3G (top) and 10G (bottom) cases using the method outlined in Heinzel et al. (2015). Right: evolution in the flux rope O-point and prominence mass are overplotted on the prominence evolution as the solid black and blue lines respectively. The large jump in altitude experienced by the 10G case is caused by the injection of the larger of two subordinate ropes into the main forming rope (cf. Fig. 8).

As time progresses for the higher-resolution setups, the locations of the $N_{\mathrm{m},\mathrm{p}}^2<0$ plasma become more evenly distributed, consistent with the observation of plasma cooling throughout the flux rope and falling towards the concave-up dips of the magnetic topology. Once material fully condenses to the concave-up dips we see it no longer highlighted as a region of $N_{\mathrm{m},\mathrm{p}}^2<0$. The material within the prominence has once more regained pressure equilibrium (Blokland & Keppens 2011a,b). This material does, however, appear to continue to undergo rotations once within the prominence as indicated by the persistent concentrations of strong baroclinicity at the edges of the formed prominence. On the other hand, the localisation of the signatures to the periphery of the prominence could indicate a purely local rotation induced by the ongoing condensation of material onto the already-condensed prominence.

Finally, Fig. 10 presented the profile of plasma density, current density, vertical velocity, and temperature along a vertical cut along x = 0 at a time within the 3G case when the oscillation of the formed monolithic structure was minimal. Previously, it was noted that the compression of the field lines is much smaller than in the previous case of Sect. 3.1, nevertheless the signature consistent with cross-field mass diffusion is still found (Low et al. 2012a,b). Specifically, a local current and negative velocity enhancement is present at the leading edge of the falling monolith as predicted in the 1D steady-state study of Low & Egan (2014). Once again, in comparison with the 1D steady-state model of Low & Egan (2014), we see that the velocity and current enhancements are not cospatial; the velocity enhancement is located below the current enhancement as was the case in the lower-resolution simulation. By definition, however, the authors state that the amount of material diffused from one flux surface to another in this way is quantitatively finite but as-yet unbounded, and so an explicit comparison is both non-trivial and outside the scope of this study. These results are, nevertheless, consistent with the hypothesis of Low et al. (2012a,b) and Low & Egan (2014) and indicate that this physical process should be studied in more detail.

4.2. Additional features and dynamics

The primary focus of this paper was in the physical processes involved in the formation and evolution of prominence condensations formed by catastrophic radiative losses. These simulations contain, however, a number of additional and interesting features that we will briefly summarise here.

The lower-resolution case explored in Sect. 3.1 may be directly compared to the previous study of Kaneko & Yokoyama (2015), wherein the main differences are (1) a different treatment of the upper boundary (K&Y: ‘free’ boundary; here: fixed density, pressure, and velocity + ‘extrapolation’ of magnetic quantities); (2) a choice of background heating function that resembles the assumed photospheric source for coronal heating (cf. Fang et al. 2013), and (3) the relaxation of their limit on a minimum temperature to below 0.3 MK. The difference in the choice of the upper boundary is unlikely to have a large impact on the resulting dynamics within the forming condensations. The different heating functions, on the other hand, are able to explain the delay in the formation encountered in the case of Sect. 3.1 in comparison with this previous work. Then, the relaxation of the lower limit demonstrates why the condensations of Kaneko & Yokoyama (2015) were unable to collapse into a monolith but, instead, remained fairly diffuse. As explored in this work, the hydrostatic equilibrium within prominence condensations may have a significant role to play in the internal structuring and evolution of these condensations. Indeed, the pressure scale height (Eq. (18)) for an atmosphere of hydrogen at 0.3 MK is ≈9 Mm versus 200 km for T ≈ 6000 K – the temperature of the monolithic structure in Figs. 3 and 5.

Figure 8 details the formation of subordinate flux ropes during the formation of the main 3G and 10G flux ropes. The exaggerated size of these structures has already been attributed to the large value of η adopted within the higher-resolution cases. Nevertheless, the previous studies of Zhao et al. (2017, 2019) simulated how chromospheric plasma can be levitated in the bottom part of a forming flux rope, which was formed in a similar process as studied here. Indeed, this is also further corroborated by the observational studies of Liu et al. (2012b) and Zhang et al. (2014, 2020). In this work we find that the subordinate flux rope axes are relatively void of material and note the decrease in density in each case. Two relatively minor enhancements in density are, however, located at both the bow and in the wake of the rising ropes (an observational feature commonly associated with flux ropes erupting from the solar atmosphere and forming coronal mass ejections (CMEs), cf. ‘three-part structure’, Dere et al. 1997). Such evolving CMEs are commonly likened to snow plows in how they appear to scoop up material in the path of their eruption. The material that trails behind the core of the subordinate ropes is then akin to the prominence material in the core of an erupting or stable rope. Since it is much hotter material, it would lack sufficient density or pressure to produce a characteristic absorption (a feature rectified with the inclusion of a cool chromosphere as in Zhao et al. 2017, 2019).

Upon the condensation of prominence plasma at locations away from the topological dips, the localised density enhancements start sliding downwards, triggered by $N_{\mathrm{m},\mathrm{p}}^2<0$ regions in each case. Figure 14 details the motion of those condensations that contain sufficient density and pressure (therein opacity) to produce an absorption or emission signature according to the (H-α) visualisation method of Heinzel et al. (2015). In each case, the filament visualisations indicate significant oscillatory motions that are not a consequence of some global MHD wave mode e.g., kink or sausage but instead simple oscillatory motions along a relatively rigid magnetic field topology. The filament visualisations also indicate these condensations originate on either side of the flux rope axis and flow along flux surfaces irrespective of helicity. According to the prominence visualisations with O-point overlaid, in many cases these condensations can have initial heights above the central flux rope axis before ultimately collecting in the bottom of the flux rope. Furthermore, the motions of the smaller, more concentrated condensations in the latter stages of the 3G case are clearly visible within the (H-α) visualisations. This suggests these TI-formed structures of dimensions ≈0.07″ will potentially be observable using the unprecedented spatial and temporal observations available with the upcoming Daniel K. Inouye Solar Telescope (DKIST) instrument suite (diffraction-limited to 0.011″ per pixel).

Each of the prominence formation cases explored in this paper ends up in a different state from the others. The lower-resolution case explored in Sect. 3.1 describes a monolithic structure that was slowly losing mass from the main body via isolated but repeated instances of reconnection. The 3G case does not display these small-scale islands falling from the underside but instead suggests the entire prominence itself may be falling through the host flux rope. In the wake of the falling prominence we see signatures of reconnection i.e., topology changes, in addition to localised heating that propagates throughout the core of the flux rope, heating it to temperatures in excess of 4 MK. The 10G case was not advanced in time to the same evolutionary state as the previous two but the final state of the 10G simulation suggested that the plasma was having a minimal effect on the topology of its host magnetic field. Consequently we do not see the core of the flux rope undergo sudden heating, instead the core is already at a temperature above 4 MK and the slight increase we do observe over time is most likely due to the dissipation of concentrated electric current as Joule heating. This hot, rather than cool, prominence cavity is consistent with the observations of hot X-ray coronal cavities in the literature, including a recent statistical study by Bąk-Stęślicka et al. (2019; see also e.g., Hudson et al. 1999; Fan & Gibson 2006; Reeves et al. 2012 although we note for the comparison with Fan & Gibson 2006 that the simulations presented here do not contain the bald patch separatrix surface (BPSS; Bungey et al. 1996) configuration, but the lowest concave-up dips of the flux rope are indeed close to the bottom of the simulation domain).

The aforementioned difference between the 3G and 10G cases is also well captured in Fig. 14. We see that for the 3G case, the O-point of the flux rope initially oscillates as a result of the formation process before the condensation of the prominence interrupts this bulk oscillation, as indicated with the clear displacement of the O-point to lower heights (cf. Shafranov Shift, Blokland & Keppens 2011a; Jenkins et al. 2019). Strikingly, the position of O-point and prominence material in time appear completely correlated between 2800 ⪅ t ⪅ 4500 s before the initiation of the internal reconnection that allows both the prominence material to continue falling and the O-point to return to a height of ≈15 Mm (an evolution that may prove related and/or comparable to eruptions that involve the ejection of a flux rope but leave the prominence behind e.g., Gilbert et al. 2001; Jenkins et al. 2018). The 10G case, on the other hand, also displays oscillations related to the formation process but upon the appearance of condensations the O-point appears comparatively unshifted. Considering the masses of all cases, the lower-resolution example attains a mass of ≈2650 g cm⁻¹, in comparison with ≈5000 and 3500 g cm⁻¹ for the 3G and 10G cases, respectively. Assuming a prominence with a characteristic length of 100 Mm such values would equate to total prominence masses of the order 10¹³ g – an order of magnitude less than typically observed, two orders of magnitude less than some of the larger observational examples (Labrosse et al. 2010; Vial & Engvold 2015). Hence, although the perturbation to the O-point is small in the 10G case, the prominence supplying this perturbation is correspondingly light.

The small masses of the prominences formed within this study are to be expected as it has long been known that the mass of a prominence cannot solely be supplied by the condensation of material already residing within a flux rope formed within the corona (e.g., Saito & Tandberg-Hanssen 1973). Instead, some additional ongoing mass-supply process is required for prominences of the order 10^{14 − 15} g to form. Zhao et al. (2017) presented a numerical study of the formation of a flux rope + prominence system wherein the cool chromospheric plasma was scooped up during the formation (reconnection) process (cf. the flux-feeding of Fig. 8). These authors report that their prominence attained a mass almost one order of magnitude higher than we find here. Furthermore, the study by Xia & Keppens (2016) described the highly dynamic and ongoing formation and draining of plasma across a fully 3D flux-rope + prominence system. These authors, too, were only able to create a prominence of order 10¹³ g. The studies by Fan (2017, 2018, 2020) have been able to produce prominences of order 10¹⁵ g but the authors concede that this is due to the use of an enhanced heating term implemented in the photosphere that drives a significant amount of plasma into a very large (> 500 Mm) flux rope.

And herein lies the particular problem. In order to obtain large, but observed, prominence masses, simulations are so far required to use empirical and unconstrained modifications to the physical evaporation processes to drive their formation (Fan 2017, 2018, 2020). Adopting physical processes that are more readily constrained, the community find large masses but in the absence of particularly exotic dynamics - in particular with regards to the vertical motions observed within limb-projected prominences (Zhao et al. 2017). And finally, in the cases explored here and by Xia & Keppens (2016), the reproduction of the dynamic motions within prominences then limits the mass budget but also results in visualisations that are saturated in observations according to the methods of Heinzel et al. (2015). As such, it is clear that to address these inconsistencies, we must continue to pursue the understanding of the physical plasma and magnetic processes occuring within solar prominences, and do so with the combination of high-resolution and fully 3D simulations (Gibson 2018).

5. Summary of conclusions

Solar prominences often appear within the solar corona with little indication of how the material got there. Multiple models have been proposed for this formation process, in this paper we have focussed on the spontaneous condensation model and characterised those processes responsible for the initial formation and subsequent evolution. Herein we have used a combination of advances provided by modern grid-adaptive simulations, to also describe evolutions at unprecedented resolutions. The CCI has been shown to map well onto those flux surfaces that are slowly redistributing plasma in response to an evolving pressure gradient. Shortly thereafter, the TI – driven by temperature decreases associated with optically thin radiative losses – proves responsible for the initial formation of the cool condensations. Maps of the field-projected Brunt-Väisälä (BV) frequency $N_{\mathrm{m},\mathrm{p}}^2<0$ indicate these condensations form out of pressure balance with their surroundings and subsequently evolve towards the concave-up topological dips in the magnetic field. In the process of falling towards equilibrium, local misalignments between the pressure and density gradients within the condensations induce rotations as indicated with the maps of the baroclinicity metric. Therefore, we conclude that the initial formation and evolution processes for prominence condensations follows the ordering of CCI – TI – BV and baroclinicity. Upon the complete formation of prominence condensations, their evolution appears to be a consequence of the BV and baroclinicity, with an apparent consistency with the mass-slippage theory of Low et al. (2012a,b), Low & Egan (2014). In future work, we plan to revisit this scenario in full 3D, as already done at modest resolutions by Kaneko & Yokoyama (2018), where the extra dimension will suddenly allow for additional interchange and RTI routes to dynamical rearrangements. It will be of interest to then see how the condensations (mis)align with the overall helical field line structure within the flux rope, as this possibility to misalign density structures with the magnetic field was found to be unavoidable in idealised 2D and 3D local box simulations (Claes & Keppens 2019; Claes et al. 2020).

Movies

Movie of Fig. 3 Access here

Movie of Fig. 9 Access here

Movie of Fig. 11 Access here

Acknowledgments

We wish to thank the anonymous referee for their constructive comments that improved the paper. We acknowledge the open source software that made possible the data visualisations presented within this work (Paraview;Python;yt-project;matplotlib). RK and JJ are supported by the ERC Advanced Grant PROMINENT and a joint FWO-NSFC grant G0E9619N. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 833251 PROMINENT ERC-ADG 2018). This research is further supported by 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.

References

All Tables

All Figures

	Fig. 1. Distribution of plasma density (top row) and temperature (bottom row) within the forming flux rope. Magnetic field lines have been overplotted in white to outline the isolated flux region commonly referred to as a magnetic flux rope.
In the text

	Fig. 2. Comparison between both the distribution and the magnitude of radiative losses within the flux rope immediately after the formation of the flux rope (left) and just prior to the appearance of condensing material (right).
In the text

Fig. 3. Evolution of the condensing material within the formed flux rope of the 31 km resolution case. The first row details the evolution of the density, the second details the temperature, the third details the Brunt-Väisälä instability ($N_{\mathrm{m},\mathrm{p}}^2$), and the fourth details the baroclinicity (positive → clockwise rotation). Velocity quivers of relative magnitude are overplotted as white arrows to indicate the direction of, and increase in, condensation rate with time. A movie of this figure is available online.

In the text

	Fig. 4. Locations of TI within the 31 km gridstep prominence. Top: ratio of thermal to density influence on the varying energy loss 𝒬 term in Eq. (3). Bottom: spatial map of the isochoric linear stability metric C (Eq. (23)) for the thermal mode instability.
In the text

	Fig. 5. Evolution of the monolithic prominence once formed by catastrophic, optically thin radiative cooling. Top row: evolution of the density with magnetic field lines indicated using white streamlines. Middle row: value of plasma β throughout the flux rope and prominence. Bottom row: evolution of the baroclinicity.
In the text

Fig. 6. Magnetic field and density map, and the evolution of plasma density, current density, vertical velocity, and temperature along a cut taken vertically through x = 0. Left: zoom-in of the region containing the bottom of the prominence condensation. Magnetic field lines are plotted in red with density contours overlaid according to the colourbar. Right: four panels detailing the evolution of the aforementioned physical parameters across the leading edge of the falling condensation, times are indicated as titles to each plot.

In the text

Fig. 7. Magnetic field and density map, and stratification of plasma density, current density, vertical velocity, and temperature along a cut taken vertically through x = 0. Left: a zoom-in of the region containing the multiple prominence condensations. Magnetic field lines are plotted in red with density contours overlaid according to the colourbar. Top: stratification of the aforementioned plasma parameters between 2–8 Mm, isolating the three prominence consensations. Bottom: three panels focussing on the distribution of parameters across the cusp of each blob.

In the text

Fig. 8. Production of subordinate flux ropes below the main rope and their involvement with the supply of material to the prominence. Top panels: distribution of Jz current density. Bottom panels: current and plasma density distribution in a vertical cut through the current sheet at x = 0. Left and right: construction of two additional flux ropes (t = 915.4, 1215.7 and t = 735.2, 1001.2 s) within the 3G and 10G cases, respectively.

In the text

	Fig. 9. Evolution of the condensing material within the 3G formed flux rope. First row: evolution of the density, second: temperature. Velocity quivers of relative magnitude are overplotted as before. A movie of this figure, including $N_{\mathrm{m},\mathrm{p}}^2$ and baroclinicity maps, is available online.
In the text

Fig. 10. Magnetic field and density map, and stratification of plasma density, current density, vertical velocity, and temperature within the 3G case along a cut taken vertically through x = 0 at t = 5136.76 s. Left: a zoom in of the bottom of the prominence condensation. Magnetic field lines are plotted in red with density contours overlaid at thresholds according to the colourbar. Middle: stratification of the aforementioned plasma parameters between 0–12.5 Mm, isolating the entire prominence condensation. Right: a zoom-in of the cusp of the condensation, showing the behaviour of the plasma parameters across the boundary.

In the text

	Fig. 11. Evolution of the condensing material within the 10G formed flux rope. First row: evolution of the density, second: temperature. Velocity quivers of relative magnitude are overplotted as before. A movie of this figure, including $N_{\mathrm{m},\mathrm{p}}^2$ and baroclinicity maps, is available online.
In the text

Fig. 12. Density, temperature, and pressure cross-section of the initial filamentary condensation within the 10G prominence formation case. Top panel: distribution of density within the forming filamentary condensation. Velocity quivers of relative magnitude are overlaid in white. Magnetic field lines are overlaid as solid-red lines. Bottom panel: density, temperature, and pressure profiles across the condensation cross-section, taken along the dashed-black line in the top panel. Profiles computed using yt (yt-project: Turk et al. 2011).

In the text

	Fig. 13. Locations of TI within the 10G prominence just prior to the two main episodes of condensation formation (left: 5134.34 s and right: 6629.68 s). Top: ratio of thermal to density influence on the varying energy loss 𝒬 term in Eq. (3). Bottom: spatial map of the isochoric linear stability metric C (Eq. (23)) for the thermal mode instability.
In the text

Fig. 14. Filament (left) and prominence (H-α) visualisations (right) of the 3G (top) and 10G (bottom) cases using the method outlined in Heinzel et al. (2015). Right: evolution in the flux rope O-point and prominence mass are overplotted on the prominence evolution as the solid black and blue lines respectively. The large jump in altitude experienced by the 10G case is caused by the injection of the larger of two subordinate ropes into the main forming rope (cf. Fig. 8).

In the text