© The Authors 2025

1. Introduction

Quiescent prominences are large, cool, and dense structures of plasma that exist stably in the hot solar corona. They have their origins in filament channels and manifest themselves above magnetic polarity inversion lines (PILs), which are the demarcation between mostly positive and negative magnetic fields at photospheric heights. These prominences are typically observed in quiet areas within the Sun’s atmosphere, spanning various latitudes and occupying the most extensive scales which range from 60–600 Mm in length, 15–100 Mm in height, and 5–15 Mm in thickness (Tandberg-Hanssen 1995). Early studies of their fundamental properties deduced characteristic temperatures and densities of ∼10⁴ K, densities ∼10⁻¹³ g cm⁻³, respectively, with initial magnetic field strength estimates around 3–10 G (Tandberg-Hanssen 1995). Such ranges insist that prominences contain dynamics that are low plasma β in nature, embedded within the equivalently low plasma beta solar corona. Furthermore, Rust (1967) calculated prominence β values of order ≈0.02, as derived from average measurements of the magnetic field ≈10 G, and the pressure ≈0.1 dyne cm⁻². In more recent observations, Casini et al. (2003) have reported the existence of localised regions within quiescent prominences that exhibit significantly higher magnetic field strengths, ranging from 60–80 G.

Early studies of prominences focused on the velocity field, in particular noting the presence of upflows and downflows of the order of a few km s ⁻¹ (Engvold 1981; Kubota & Uesugi 1986). Stellmacher & Wiehr (1973) reported the development of a cavity within a quiescent prominence, as seen in Hα slit-yaw images, and linked it with a fast magnetohydrodynamic (MHD) wave signal. Engvold (1981) initially interpreted observed downflows as occurring along magnetic field lines that are driven by the force of gravity. Motions observed in prominences include the formation of shear flows, strongly suggesting the presence of instabilities arising from this shear. Furthermore, Liggett & Zirin (1984) observed the presence of rotating regions within non-eruptive prominences (five out of 51 observed cases), with sizes between 3000 × 7000 km² and up to 30 × 104 Mm² and exhibiting rotation rates of around 30 km s⁻¹. Shortly thereafter, however, early observations of magnetic fields in prominences revealed a predominantly horizontal magnetic field within the developing vertical plasma structure (Leroy 1989). On first inspection, it was not at all clear how shear flows can develop perpendicular to the magnetic field under low-beta conditions – in 2D for a magnetic field oriented across the plane, this is prohibited without significant fieldline deformation, and thus is in clear contradiction with the observations.

The high-resolution HinodeSolar Optical Telescope (Hinode; Kosugi et al. 2007, SOT; Tsuneta et al. 2008) has improved our understanding of the nonlinear behaviour of prominence internal dynamics (Chae 2010; Hillier et al. 2012a). Chae (2010) found from the observations that the vertical fine substructures, termed knots, were impulsively accelerated in the downward direction. The authors related this to magnetic reconnection and the interchange of the magnetic configuration. The SOT observations of Berger et al. (2008, 2010, 2011) simultaneously linked these dynamics to the Rayleigh-Taylor (RT) instability (Hillier et al. 2012a,b; Keppens et al. 2015; Xia & Keppens 2016); in a 3D volume, or a carefully constructed 2D plane, such cross-field-like behaviour is permitted even in low-beta plasmas. In turn, the Kelvin-Helmholtz (KH) instability is a well-known shear flow instability that leads to the fragmentation of coherent vorticity structures into individual vortices during its nonlinear development. In relation to solar prominences, and in tandem with the RT instability, the KH instability has subsequently been investigated in numerical works where good qualitative and quantitative agreement has been found with the observations (Berger et al. 2017; Hillier & Polito 2018). Ryutova et al. (2010) combined SOT observations with theoretical criteria to categorise observed plasma instabilities in quiescent prominences, identifying RT and thermal instability (TI) processes, as well as how prominence cavities form.

To summarise, quiescent prominences are seen to exhibit spatio-temporal evolution characterised by high variability. It is particularly interesting that the dynamic nature and the large prevailing Reynolds numbers then indicate the fluctuations present within prominences to be turbulent in nature. For example, Berger et al. (2017) used the SOT to study and analyse the characteristics of a quiescent prominence bubble and the associated instability dynamics. The characteristic speed of the largest rising bubble was found to be 1.3 km s ⁻¹. The prominence downflows deposited plasma onto the bubble, forming an increasingly thick boundary layer with a speed of 20–35 km s ⁻¹. This led to a strong shear flow of the order of 100 km s ⁻¹ across the bubble boundary and in turn the coupled KH-RT instability. The nonlinear evolution of the KH instability is known to play a significant role in the generation of turbulent flows by means of reconnection and secondary instabilities (Matsumoto & Hoshino 2004).

Solar prominences are dictated by complex interactions between gravitational, thermal, and magnetic forces. In Changmai et al. (2023) (hereafter MC23), we modelled this complex evolution as RT instability that later developed fully turbulent characteristics under the ideal-MHD approximation. We performed high-resolution simulations where a stratified atmosphere containing a prominence structure evolves due to Rayleigh-Taylor dynamics. This RT initiated at the density inversion at the bottom of the prominence, eventually deforming the prominence body into falling fingers and rising bubbles. Our study addressed the far-nonlinear turbulent state of the entire prominence body, as fingers got reflected upwards at the transition region and interacted with still-falling prominence matter (cf. Rees-Crockford et al. 2024). The truly mixed and turbulent state obtained was analysed statistically, performing structure-function analysis to quantify the degree of intermittency. A fair agreement with a purely observational study on observed RT turbulence in prominences by Leonardis et al. (2012) was obtained. Recently, Hillier & Polito (2021) used Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014) slit-jaw images at Si IV and Mg II lines, finding evidence for bi-directional jets that would signal internal magnetic reconnections within quiescent prominence bodies. They interpret these reconnections as resulting from shearing, buoyant flows. In addition to the general theory, we use this observation to motivate this follow-up numerical study, where we will investigate secondary instability processes in particular, within RT deforming prominences. We run highly resolved, resistive MHD models to focus on secondary shear flow effects, KH vortices, and induced reconnections.

This paper is organised as follows: we discuss about the numerical setup in Section 2, and in Section 3 we analyse and characterise the energetic events due to the secondary instabilities. In Sections 4 and 5, we present our summary and conclusions, respectively.

2. Numerical setup

To investigate small-scale dynamics in a quiescent prominence due to initial RT instability, we performed 2.5D high-resolution resistive MHD simulations of a prominence inserted in the solar atmosphere using the parallelised, open-source Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al. 2012, 2021, 2023; Porth et al. 2014; Xia et al. 2018). This work employs the same methodology as in MC23. Using Cartesian geometry, we establish a large simulation domain encompassing 30 Mm horizontally along the solar surface and 30 Mm vertically from the photosphere – low chromosphere into the solar corona. We use a maximum refinement level of seven and a base grid resolution of 40 × 40 to obtain a final grid resolution of 2560 × 2560. We achieve a spatial resolution of ∼11.7 km and run the simulation for about 10 minutes of solar time. Compared to MC23, we added another AMR level to achieve this higher resolution (MC23 had cells of 23 km size), and we add explicitly resolved resistivity.

We solved the resistive MHD equations for the conservation of mass, momentum, total energy density, and induction. The equations are given by

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

$$\begin{array}{cc}\hfill {\partial }_t(\rho \mathbf{v})+\mathrm{\nabla }·(\mathbf{v}\rho \mathbf{v}-\mathbf{BB})+\mathrm{\nabla }{p}_{\mathrm{tot}}& =\rho \mathbf{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {\partial }_{t}e+\mathrm{\nabla }·(\mathbf{v}e-\mathbf{BB}·\mathbf{v}+\mathbf{v}{p}_{\mathrm{tot}})& =\rho \mathbf{g}·\mathbf{v}+\mathrm{\nabla }·(\mathbf{B}\times \eta \mathbf{J}),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\partial }_t\mathbf{B}+\mathrm{\nabla }·(\mathbf{vB}-\mathbf{Bv})& =-\mathrm{\nabla }\times (\eta \mathbf{J}),\hfill \end{array}$$(4)

where p = (γ − 1)(e − ρv²/2 − B²/2), p_tot = p + B²/2, and J = ∇ × B. The quantities ρ, v, p, p_tot, e, and B denote density, velocity vector, plasma pressure, total pressure, total energy density, and magnetic field vector, respectively, and γ (adiabatic index) is the ratio of specific heats, taken as 5/3 under the assumption of a mono-atomic ideal gas. A uniform resistivity, η = 10⁻⁴ (or 1.2 × 10¹³ cm² s⁻¹ in physical units) is taken throughout the entire simulation domain. The magnetic field vector is measured in units for which magnetic permeability is equal to μ₀ = 1. The plasma temperature, T, follows from the ideal gas law,

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

where R and μ are the gas constant and mean molecular mass, respectively. The normalisation of the quantities is done identically to MC23.

The Eqs. (1)–(4) are solved numerically on a dynamically changing grid structure, known as a hierarchical block-adaptive grid. The temporal integration strategy is a four−stage, third−order Runge-Kutta method, as described by Ruuth & Spiteri (2002). In the numerical flux computations, we employ the Harten-Lax-van-Leer approximate Riemann solver for multiple discontinuities (HLLD) following the approach proposed by Miyoshi & Kusano (2005). The HLLD method allows to accurately simulate high-resolution magnetohydrodynamic (MHD) evolutions while retaining positivity. Additionally, a third-order asymmetric slope limiter, originally introduced by Koren (1993), is used in cell-centre to cell-edge reconstruction. In our grid-adaptive simulations, refinement criteria are employed. These criteria are based on the instantaneous properties of the plasma and their gradients, following the methodology outlined by Lohner (1987), and use the refinement criteria based on density. In the AMR code, the structure block-tree indicates that the option to refine or coarsen must be taken for each individual block (Keppens et al. 2012). In our simulation, we used a block size of 10 × 10. The generalised Lagrange multiplier (GLM) method was employed to control erroneous numerical magnetic field divergence, as described by Dedner et al. (2002). A Courant value of 0.8 was selected in order to maintain temporal stability during the explicit time integration process. The initial setup and boundary conditions are identical to our previous 2D ideal MHD simulations from MC23. The introduction of uniform η extends that previous setup to a resistive MHD simulation. The mean magnetic field strength is taken to be around 10 G. The initial magnetic field makes an angle (θ) of ∼88° with the horizontal x-axis, so it is nearly perpendicular to the simulated plane. We use the same localised multi-mode perturbation as in MC23 to induce turbulence due to Rayleigh Taylor (RT) instability in the initial prominence setup and let the system evolve self-consistently. We note that without this perturbation, the prominence stays totally pressure-supported at its initial height. The boundary conditions are identical to our previous work, that is, periodic along the horizontal x-direction and a combination of hydrostatic and extrapolated in the vertical y-direction.

3. Results

3.1. Overall evolution of the prominence

The evolution of density, as well as the corresponding compressibility quantification through the sound speed (L/c_s, L = Δ_(x, y)) scaled divergence of the velocity field, for the resistive simulation is shown in Figure 1 for time = 85.87, 429.37, and 618.29 s. From this figure, we can see how the quiescent prominence evolves from the onset of instability which occurs as a result of gravity’s influence on an initial multi-mode perturbation and leads to the formation of vertical structures. The coherent structures, which manifest as pillars and bubbles, exhibit a progressive increase in fine-scale structures over time. Note that the presence and movement of (even very weak) shocks can be clearly identified in the velocity divergence.

Fig. 1. Combined plots of log ρ and $\frac{L}{c_{\mathrm{s}}}\mathrm{\nabla }·\mathbf{v}$ at time 85.87 s, 429.37 s, and 618.29 s (left, middle, right) showing the evolution of bow shocks in the evolution of the solar prominence. The horizontal red dashed line indicates the 2 Mm threshold relevant for Figure 2.

At time = 85.87 s, the deformation of the initial interface into vertical RT structures in the form of bubbles and pillars is in progress, with the nonlinear nature of the evolution already clear not only in the density distributions but also the presence of secondary compressive features within the prominence body, shown here in the bottom row of Figure 1. In particular, the ∇ ⋅ v panel describes the formation of bow shocks above each of the rising plumes and falling fingers, which are visible as arcs.

From time = 429.37 s, we see that the cool and dense material of the prominence falls and collides the chromospheric boundary, such that it is partly reflected back to higher heights. Portions of the fallen prominence material remain in contact with the chromosphere but do not successfully mix. This deflection leads to the formation of more horizontal coherent structures within the prominence which was otherwise dominated by vertical structures in the longitudinal direction. At about this time (cf. MC23), the higher regions become beta unity and above, which results in the rising of the buoyant matter through our top boundary and fine-scale mixing, as the thermal pressure rises accordingly. In Figure 1 (bottom-middle), shock waves which were seen as internal arch-shaped shocks in the initial stages, now appear to follow the large-scale deformations within the prominence material. The overall kink deformation of the prominence body, caused by the rising plumes, gradually displaces the upper prominence interface towards the upper boundary.

Finally, at time = 618.29 s, we see the complete development of the principle plumes and pillars of the nonlinear RT-unstable prominence. The falling fingers and rising bubbles are accompanied by a wealth of horizontal structures closer to the chromospheric boundary, filled with smaller scale structures that appear as discrete ‘cells’ – not to be confused with the individual grid cells of the simulation. We will explore these features more in Section 3.2. An equal number of two large scale pillars and plumes can be discerned, with the latter now in the process of lifting a significant portion of the initial prominence mass through the upper boundary. In our previous study (MC23), we statistically explored the far nonlinear evolution of the RT-deformed prominence and its turbulent properties, but here we will focus only on this earlier phase of its evolution, where the prominence structure is still largely identifiable as seen in these figures.

Figure 2 shows the evolution of the percentage of the cool dense prominence material within the simulation domain over time where the material below 2 Mm, the chromosphere, has been neglected. Cool and dense prominence material is defined as cells having a temperature below 5 × 10⁴ K, and a density higher than 10⁻¹¹ kg m⁻³. The dashed line in Figure 2 represents the first strong instance of mixing of the prominence material around time = 350 s when two large pillars collide with the chromosphere boundary and result in the initial drop of the mass of the cool and dense volume. This decrease is entirely explained as the material breaching under the 2 Mm threshold. The solid line in Figure 2, represents the phase of the prominence at time = 429.37 s shown in Figure 1 (middle). At this time, there is an ongoing exchange of material between the prominence and the chromosphere, with some being reflected back towards higher heights but the majority remaining below 2 Mm, hence the overall decrease. At the same time, the buoyancy of the rising bubbles is in the process of lifting the upper portions of the prominence towards the top boundary. The dotted line represents the prominence around time ∼620 s, when the initial horizontal prominence material is completely deformed due to the evolving shocks and rises towards the top boundary leading to an ∼28% decrease in the cool and dense volume of the prominence at that time. We have considered the evolution of the prominence due to the initial RT instability until time = 700 s, where the loss of the cool and dense volume is around ∼33% of the initial prominence material. As already discussed, we will therefore restrict our following analysis to this initial roughly 11-minute time frame.

Fig. 2. Time evolution of the cool and dense solar prominence matter above the chromosphere boundary (above 2 Mm from the bottom boundary) taken below a threshold temperature of 5 × 104 K, and denser than 10−11 kg m−3. The dashed line represents the phase of the initial decrease of mass of the cool and dense volume in the solar prominence due to the merging of the cool dense material into the chromosphere boundary. The dotted line represents the phase where the prominence material rises towards the top boundary with a significant loss of the prominence material. The solid line represents the evolution at time = 429.37 s (see Figure 1 (middle)).

Throughout the initial evolution of our model, we find the formation of nonlinear structures from the initial stratification down to the small fine structures within the aforementioned cells of Figure 1. The cold and dense material is primarily located in the pillars that descend under the effect of gravity and hit the chromosphere-to-corona transition region. At this point, there is a clear change and accumulation of kinetic energy in the horizontal direction attributed to the nonlinear stage development of the RT instability. Thereafter, vortices and swirling structures transfer energy among the scales and in both directions. In what follows, we will further explore the presence and characteristics of these complex plasma motions.

3.2. Presence of secondary processes

3.2.1. The KH instability

From the initial time when the setup is in a slightly perturbed equilibrium state, the dynamics quickly become nonlinear due to the multi-mode RT instability induced by the perturbation within the prominence-corona interface. In this early stage, we see the formation of falling fingers or pillars and rising bubble-shaped structures, in both cases there is obvious mixing at the interface layers. Under the effect of gravity, denser and cooler plasma precipitates towards the sun’s surface, while buoyancy drives the less dense, hotter material upwards. This differential motion sets up convective patterns within the prominence and its coronal environment, thereby establishing shear flow layers.

Under hydrodynamic conditions, shear layers can be susceptible to KH instability if the velocity jump is sufficient. The differential motion of magnetised plasma induces flow-related instabilities, such as the KH type, that are sensitive to the magnetic field orientation; if k ⋅ B ≠ 0, where k is the wave vector, magnetic tension can suppress the development of rotational motion. In our simulation, the main magnetic field component is largely perpendicular to the simulated plane, and so we approximately recover the hydrodynamic behaviour in the sheared regions in the (x, y) plane since k ⋅ B ≈ 0 and the restoring magnetic tension force is very weak. This is particular to our 2.5D setup, where invariance in the z-direction is adopted.

Figure 3 shows the evolution of a localised shear flow interaction region that begins as an approximately laminar shear interface before fragmenting into discrete vortices. We show the velocity field magnitudes over the whole domain for time = 274.79, 291.97, and 343.49 s in its middle row. The top panels show a zoomed region manifesting KH instability, tracked using a combination of vorticity as ω = ∇ × v on the left of each pair, and the so-called ‘R-criterion’ on the right (Canivete Cuissa & Steiner 2020, sometimes referred to as a ‘Rortex’). Here, following Wang et al. (2019) and Xu et al. (2019), $R=|\omega |-\sqrt{{|\omega |}^2-{\lambda }^2}$, with λ = 2λ_ci, where ${\lambda }_{\mathrm{ci}}=\frac{1}{2}\sqrt{-[(u{\prime }_{\mathrm{x}}-v{\prime }_{\mathrm{y}}]-4u{\prime }_{\mathrm{y}}v{\prime }_{\mathrm{x}}}$ for a negative discriminant and 0 otherwise, and $u{\prime }_{\mathrm{x}}=\frac{\delta u}{\delta x}$, etc., are the ∇u 2D tensor components (as in Canivete Cuissa & Steiner 2020). Although numerous formalisms exist in the literature to track the location, orientation, and strength of rotating fluids, such as ω, the Q-criterion, and λ_ci, R was constructed to be completely independent of the shear components that frequently complicate associated analysis (Tian et al. 2018). Locations of rotating fluid with very large radii of curvature, resembling sheared flows at the infinite limit, can then be excluded if the ratio of $\frac{{\lambda }_{\mathrm{cr}}}{{\lambda }_{\mathrm{ci}}}$, where λ_cr is the real solution to the associated eigenvalue analysis, is above some threshold, oftentimes 1. We set a threshold on this ratio to 0.1 to be even more restrictive. The bottom insets show the same zoomed region for time = 291.97 s, presenting the in-plane v_x and v_y component of the velocity field at that time. We find the vorticity to well describe the shear interface regions, with the additional R-criterion refining the spatial extent to just those features where rotation dominates over shear or general deformation. At this time, the utility of the R-criterion is clear as it refines the mapping onto portions of the discrete KH ‘rollups’ identifiable by eye between 8 and 11 Mm. We note that local-box, two-fluid plasma-neutral simulations by Snow & Hillier (2024) recently highlighted the role of recombination and ionisation processes in the thermodynamics of a KH shear layer. Future work could exploit this aspect in even higher resolution, global setups as realised here in single fluid settings.

Fig. 3. Evolution of secondary KH instability in the prominence at time = 274.79, 291.97, and 343.49 s as shown in the red insets. The top row shows the zoomed image for the KH instability regions, represented as both vorticity ∇ × v and Rortex R, identified using the velocity magnitude in the middle row for each time. The bottom insets show the zoom region for time = 291.97 s, representing the vy and vx component of the velocity field at that time. Together, they signal a trail of (rising) vortex structures being formed.

In Figure 4, we present a snapshot of the condition of the simulated plasma at time = 472.3 s in those locations identified by the R-criterion such that R > 5 × 10⁻² s⁻¹ where λ_ci > 10λ_cr. In the left panel we see that many discrete locations satisfying the selection criteria are present across the simulation domain. In particular, at higher temperatures and this is reflected in the right panel presenting the 2D density – temperature histogram of those regions satisfying the selection criteria. The combination of higher temperature and lower density agrees with the 2D representation on the left, where it is clear the majority of the rotating plasma is located within the more coronal material, with some but far less present in the cooler prominence material itself. At this time, we see that 6.1% of the 30 × 30 Mm domain is undergoing rotational motion. In the animation that accompanies the Figure, we see that this fraction increases with time, peaking at 7.3%, in agreement with the steady increase of secondary processes as the initially RT-unstable plasma develops mixing interface regions (cf. the prominence-corona-transition-region (PCTR)) that in turn interact. In the left panel of Figure 4 and the associated movie, it is clear that smaller, discrete, and often short lived or ‘flickering’ regions are identified by the R-criterion - in fact these features are present in all methods referred to above as a consequence of gradients on a discrete grid. Following authors such as Yadav et al. (2020), one may homogenise the signatures of rotating cells by smoothing the velocity field, and obtain a more robust estimate for the number of cells present within the simulation at a given time. The FWHM of the smoothing Gaussian is set assuming identification over four pixels at the IRIS platescale of 0.167″. We find the number of cells vary by a factor of two throughout the temporal evolution, peaking at ∼10. Such a statistic is heavily influenced by the selection criteria, of course, where a broader smoothing to match the platescale of the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) instrument onboard the Solar Dynamics Observatory (SDO; Pesnell et al. 2012), for example, would lead to less identified cells.

Fig. 4. Condition of the plasma undergoing rotational motion within the domain, as highlighted according to the R criterion. Left; the 2.5D domain of temperature with positions of R > 5 × 10−2 s−1 and where λci > 10λcr overlaid in white, see text for description of these terms. A few clear KH instability cells have been indicated with teal boxes. Right; the density – temperature phase space of these extracted regions. A broad distribution is present, with a clear peak at the high temperatures and low densities of the solar corona. The inset axis plots the number of distinct KH cells in time (s), current time indicated by the dashed-black line, explanation in the text. An animation of this Figure will be available online.

The observational work of Hillier & Polito (2018) shows the presence of KH instability vortices associated with downflows of the plasma in quiescent prominences. Different flow patterns were associated with the observed instabilities. Firstly, a coiling pattern as the blob falls and begins to wrap around the structure as well as a sinusoidal pattern that is mildly symmetric about the axis of the thread which moves upwards. The velocity associated with the downflow pattern was observed to be 9–16 km s⁻¹ whereas the velocity of the upflow sinusoidal pattern was 34 km s⁻¹. The width of the first pattern was 0.9 Mm and the width of the second pattern was 0.48 Mm. The length covered by the thread throughout the process of rolling up in the downflow was 3.2 Mm and that of the sinusoidal structure was 2 Mm. At time 291.97 s, the zoomed insets in the bottom row of Figure 3 show the individual velocity (v_x and v_y) components of our simulation. The v_y component shows the vertically dominated shear flows while the v_x component demonstrates the kind of coiling pattern found in observations. We find similar characteristics to Hillier & Polito (2018) for those vortices evolving due to the KH instability, matching in particular the length scales with velocity of the same order of magnitude or higher. In the regions immediately surrounding the roll-ups, however, there are very high velocities of several hundreds of km s⁻¹. The key difference here being that, and as generally summarised in Figure 4, the vortices present within our simulation are found in the much hotter > 10⁵ K plasma rather than the cooler prominence material. This will be explored more in Section 3.3.

As the system evolves, this velocity shearing effect within the plasma may lead to the generation of high-energy events and structures such as plasmoids, current sheets, and jets. Additionally, KH instability may also foster the onset of tearing mode instabilities, a resistive MHD instability. This typically occurs on pronounced and thin two-dimensional (2D) current sheets, which in turn can be subject to secondary instabilities such as the aforementioned tearing mode. In the next section, we highlight how we can identify current sheets and get statistical info on their properties from our simulations. The interplay of these primary and secondary instabilities contributes to the complex dynamical evolution of the solar prominence environment.

3.2.2. Current sheets

The dynamic evolution of the solar prominence environment establishes the complex formation of current sheets, which are critical for magnetic reconnection and energy release. As we have discussed in Section 3.2, the preceding instabilities may play a significant role in the formation of these thin, high current-density regions through the deformation of magnetic field lines. The pathway towards prominence turbulence, driven by factors such as gravity, buoyancy, instabilities, and shear flows, creates a challenging backdrop for isolating the signatures of current sheets generated by these instabilities. The characterisation of these instability induced current sheets is vital for understanding prominence stability, the triggering mechanisms of prominence eruptions, and the associated release of energy into the solar corona.

3.2.2.1. Methodology

The identification of current sheets in nonlinear plasmas is often nontrivial in full 3D, time-dependent settings. Here, we are in 2.5D, but still a rather nonlinear setting, and we are interested in locating fine-scale current structures that form due to evolving spatial eddies. These are formed due to the onset of KH instabilities which leads to a topological change of the magnetic fields. Identification of these current sheets is important to understand the energetic events associated with magnetic reconnection. Hillier & Polito (2021) showed the presence of current sheets and the observation of bi-directional jets in quiescent prominences. We will now evaluate whether we can use a clearly automatable procedure to identify and count current sheets in our snapshots.

We follow the methodology given by Zhdankin et al. (2013) to identify local maxima associated with the current sheet formation in the MHD simulations, which suits our inhomogeneous, compressible, non-periodic boundary and gravity-driven evolution. The initial procedure remains the same in determining the regions of local maxima in a designated snapshot, here we used time 680.12 s as a reference to first detect these current sheets. Figure 5 shows the methodology used to detect current sheets in the prominence. The current density for time 680.12 s is shown in Figure 5 (left). Given the distinctive features of current sheets, which are defined by extreme values in the current density profile, the task at hand is to identify local maxima in the magnitude of the current density. In order to accomplish this, the algorithm examines all data points that exceed a certain threshold current density, denoted as J_threshold. This threshold is set to be significantly higher than the average magnitude of current density across the whole dataset, we will discuss this later. The method then identifies the local maxima inside a cubic subarray around each of these selected locations. Each subarray is positioned at the candidate point with a radius of 3, a demonstration of the well-resolved peaks is shown in the middle panel of Figure 5 denoted by the red blobs. Each maximum is thereafter associated with a current sheet denoted by a current sheet index i, and the related current density is referred to as the peak current density, J_max, i. The right panel of Figure 5 shows the correlation of the local maxima events at time 680.12 s, to the magnitude current density of the simulation at that time. We will now need to group these points into individual current sheets, requiring us to cluster nearby points and then derive the sheet typical dimensions (width and length).

Fig. 5. Detection of current sheets using the modified algorithm based on current density Jz at time 680.12 s. The left plot shows the current density magnitude, the middle plot shows the correspondence of local maxima events using a threshold higher than the mean |Jz| with the current density (|Jz|), and the right plot shows the identification of clusters derived from density clustering methodology using local maxima events at time 680.12 s. 123 separate clusters could be identified.

The next step is to find clusters around the local maxima points iteratively and then identify the points belonging to each cluster that lie above a certain threshold. We augmented the previous identification algorithm to use a clustering technique called density-based spatial clustering of applications with noise (DBSCAN; Ester et al. 1996) to find the clusters representing the current sheets. The DBSCAN method is a density-based approach that identifies and groups data points inside a specified space based on their proximity to one another. It forms aggregates by collecting points that are closely clustered together. The DBSCAN algorithm utilises a distance function to determine the proximity of data points. In this implementation, a neighbourhood growth radius of four grid points is used to expand regions around points of interest. The DBSCAN parameter ϵ, set to three grid points, defines the clustering neighbourhood size, while the minimum number of points is set to one to permit small spatial jumps and distinguish a cluster from isolated outliers. We then group those points in the neighbourhood that are half of the maximum current density threshold of that given cluster, J_min, i = J_max, i/2, and further exclude point identifications and very short current sheets with an area below 3.3 Mm². We find 123 individual clusters which are the current sheets for the reference time as shown in Figure 5. The choice of a threshold parameter and number of neighbourhood cells can greatly affect the detection of individual current sheets. Using the parameters above, a visual comparison demonstrates that this approach accurately labels most of the extended |J_z| ridges as current sheets with few exceptions (cf. Lapenta et al. 2022).

The current sheets in the simulations are evolving structures that show varying current density intensities in an overall dynamic environment. To determine the properties of these evolving current sheets, we consider J_threshold = J_mean + nσ rather than a static value. In doing so, we maintain robustness and consistency in detecting the current sheets in time. We experimented with different settings for this parameter, a summary of which is shown in Figure 6. A choice of n = 1 shows a tendency to overfit in the sense that many local, often singular maxima are identified within the highly dynamic, evacuated regions in between descending fingers, while a value of n = 3 neglects many of the more elongated current sheets that drape around the fingers and plumes. For our study, we decided to work with the intermediate parameter value of n = 2, a trade-off which captures both elongated and localised current maxima satisfactorily, and importantly does so throughout the simulation run (not only the time shown in Figure 6 and evidenced in Figure 5). Hereafter, we make a distinction between positive and negative current sheets, as in Figure 6 (red and blue), to find the individual clusters and avoid biasing caused by |J_z| grouping closely co-located but oppositely oriented sheets as single sheets. Of course, the bias of choosing n = 2 remains.

Fig. 6. Comparison of different threshold parameters of Jthreshold at time 680.12 s in correlating the local maxima events with the current density (|Jz|). The red current sheets represent the positive current sheets whereas the blue represent the negative ones. Depending on the threshold used, the algorithm may overlook significant current sheet regions to fulfil the parameter (right) or overfill the local maxima regions (left).

To characterise the evolution of the current sheets quantitatively, we follow Zhdankin et al. (2013). The largest distance between any two pairs of points in the current sheet’s x − y cross-section, given by the maximum of the euclidean pairwise distance map, is how we define the length. We then use a method similar to Uritsky et al. (2010) to estimate the second dimension (width). By assuming that the current sheet is mostly uniform along its extent, i.e. no taper, and given our definition of a current sheet considers those adjacent cells for which J_min, i = J_max, i/2 is satisfied, dividing the cross-sectional area of a given cluster’s concave hull by its length yields an average width according to the full-width at half-maximum (FWHM) of the cluster (Barber et al. 2013).

3.2.2.2. Statistics

We characterised all extracted current sheets present at times t = 670.17, 680.12, and 685.60 s, shown in Figure 7 (top) as a scatter plot of width (the sheet shortest dimension) versus length (the sheet longest dimension). The positive current sheets are marked by ‘+’ and the negative ones by ‘−’. We found 73 and 65 clusters for t = 671.53 s, 75 and 62 clusters for t = 680.12 s, and 73 and 70 for t = 704.16 s for positive and negative sheets, respectively. Note that the distinction of positive and negative sheets understandably leads to a higher number of extracted current sheets than |J_z| in Figure 7. The majority of the clusters found by the algorithm have lengths ranging between between 631–1388 km in length and 90–152 km in width. The smallest current sheets appear to cluster in tight distributions within the simulation domain, and likely contribute more heavily to the dissipation of stressed magnetic energy (left panel of Figure 6. The largest current sheet found goes up to 11355 km in length and around 116 km in width.

Fig. 7. Top: Scatter plot of all current sheets characterised in length being the longest dimension (horizontal) versus width (vertical) for three times: 73 clusters (positive) and 65 clusters (negative) for t = 671.53 s, 75 clusters (positive), and 62 clusters (negative) for t = 680.12 s, and 73 clusters (positive) and 70 clusters (negative) for t = 704.16 s (blue, green, and red, respectively). The ‘+’ refers to the positive current sheets while ‘−’ refers to the negative current sheets. Most of the current sheets are smaller ranging between 631–1388 km in length and of 61–78 km in width. The largest current sheet analysed by the algorithm in the simulation goes up to ∼11 355 km in length and around ∼116 km in width. Bottom: Statistical distribution of current sheet lengths and widths at the same times, shown as violin plots for positive (red) and negative (blue) current sheets. The width of each violin indicates the probability density estimated by kernel density estimation, while the solid orange horizontal lines mark the medians. The ‘L’ and ‘W’ labels denote distributions for length and width, respectively. Shading is used to visually separate the three time instances. These distributions reveal asymmetries between positive and negative sheets in both length and width, highlighting temporal evolution and polarity dependence in current sheet morphology of the simulation.

The violin plots in Figure 7 (bottom) display the distributions of current sheet lengths (L) and widths (W) separated by polarity and time in the simulation. Each violin visualises a kernel density estimate of the probability density function of current sheet sizes.

At t = 671.53 s, positive polarity length distributions peak near a median of ∼861.8 km with an interquartile range (IQR) from ∼692.3–1352.3 km, while negative sheets exhibit slightly larger typical lengths with median ∼991.4 km (IQR 663.3–1544.8 km). However, the corresponding sheet widths are tightly clustered: for positive sheets, the median is ∼75.6 km (IQR: 67.2–78.0 km), and for negative sheets, ∼70.2 km (IQR: 65.5–77.5 km). The width distributions are notably narrow compared to the broader length distributions, indicating far less spread and the absence of extreme outliers, which reflects the asymmetrical temporal evolution of the current sheets.

At t = 680.12 s, the distributions shift towards slightly smaller and more uniform scales. The length medians shift slightly downwards: positive sheets show a median length of ∼752.3 km (IQR 573.7–1134.5 km), while negative lengths hold a median of ∼930.8 km (IQR 708.6–1228.0 km). Width remains tightly distributed with medians of ∼74.4 km (IQR: 65.5–78.6 km) for positive, and 70.2 km (IQR: 60.2–76.3 km) for negative sheets reflecting increased variability in sheet width. This temporal evolution suggests a gradual shortening and narrowing of positive sheets, while negative sheets remain comparably larger in size.

By t = 704.16 s, the typical sizes converge yet further: positive and negative sheet lengths have nearly equal medians (∼775.7 and ∼774.3 km, respectively) and overlapping IQRs. Width distributions also narrow and align closely around ∼74.9 km (IQR: 67.9–77.4 km) for positive and ∼70.0 km (IQR: 59.9–78.5 km) for negative sheets. This convergence likely reflects the dynamic balance reached in sheet formation and dissipation as the turbulent phase develops.

The shapes and widths of the violins additionally reflect the length distributions being broader and more skewed compared to width distributions which are more sharply peaked, as seen in the top pane of the same Figure. Lengths exhibit a tail towards larger values, consistent with elongated resistive layers, whereas widths remain relatively constrained and reflect typical properties governed by the chosen η.

Overall, these violin plots quantitatively establish the temporal evolution and polarity-dependent morphology of current sheets in the simulation, highlighting weak initially marked asymmetries – likely inherent to the specific initial condition and nonlinear development i.e. clockwise versus counterclockwise KH instability development – that approach statistical similarity as time progresses and the imprint of the linear – nonlinear development is lost. Such a mild trend should be explored at higher resolutions and at corresponding η.

3.2.3. Jets

One of the energetic events associated with the occurrence of current sheets that demonstrate magnetic reconnection are jets. Observational work of Hillier & Polito (2021), has revealed the presence of bi-directional jets in solar prominences using IRIS slit-jaw imagers. In the current study, the shearing motion of plasma due to the KH instability is one of the main reasons for forming current sheets and associated jets. In Figures 8D–F, we show the time evolution of jets which we identify based on the Alfvénic Mach number

$$\begin{array}{c}\hfill M_A=\frac{|v|}{v_a},\end{array}$$(6)

where v is the flow speed of the plasma and $v_a=\frac{B}{\sqrt{\mu \rho }}$ is the Alfvén velocity. We look at in-plane velocity, while B is the total magnetic field strength and ρ is the density of the plasma. If the value of M_A > 1, the flow is super-Alfvénic. Super-Alfvénic outflows are usually a clear indication of local reconnection events. These jets are found to be short-lived events in our simulations (Figures 8A–C) having a lifetime of around 18 s. They show sideways motion as they rise up from the points where current sheets deform due to tearing-like events at magnetic-reconnection sites in the simulations. Figures 8G–L shows that the jets originate from pronounced current sheet regions as quantified by the current density field (|J_z|). The evolution of these current sheets gives rise to reconnection regions that lead to the jet exhausts shown in Panels (A–C). The reconnection is essentially component-reconnection, noting that the main magnetic field component B_z stays unidirectional out of the plane, but is clearly rearranging in-plane magnetic field (B_x, B_y). In the examples shown here, the jets mainly occur near a pronounced, curved current sheet located at about x, y = (6, 7) Mm. For the duration that we have studied, we identify two main jet exhausts that follow in quick succession over a period of around one minute, each comprised of many components in the fine structure as the initial laminar flow becomes unstable.

Fig. 8. Time series of the occurrence of jets in the prominence identified by regions where MA > 1 at time = 671.53, 680.12, and 704.16 s. Panels (A–C) show the zoomed images of Panels (D–F) and the occurrence of jets, here essentially flowing upwards. Panels (G–L) show the time evolution of current density (|Jz|) and magnetic field (|Bplane|) for zoomed regions where the jets emerge. Both the current sheets and subsequent magnetic islands responsible for the jets are clearly visible.

3.3. Synthetic observations

Prominences are frequently observed in extreme ultraviolet (EUV) using the EUV instruments from the space-based SDO/AIA, and also in the hydrogen Hα spectral line using optical instruments from ground-based observatories. To compare the spatio-temporal dynamics found in our simulations with the motivating observations of, for example, Hillier & Polito (2021), we have made synthetic observations of our simulations. In order to do so, the radiative transfer equation for the relevant elements and ionisations is required. From Rybicki & Lightman (1986), the general solution to the radiative transfer equation is written as

$$\begin{array}{c}\hfill I_{\lambda }=I_{\lambda }(0)e^{-{\tau }_{\lambda }}+{\displaystyle {\int }_0^{{\tau }_{\lambda }}}{e}^{-({\tau }_{\lambda }-{\tau }_{\lambda }\prime )}{S}_{\lambda }({\tau }_{\lambda }\prime )d{\tau }_{\lambda }\prime ,\end{array}$$(7)

where I_λ(τ_λ) is the intensity of measured light, of wavelength λ, τ_λ is the total optical thickness along the chosen line of sight (LOS), and I_λ(0) is the intensity of any background illumination. The source function is defined as $S_{\lambda }=\frac{j_{\lambda }}{{\alpha }_{\lambda }}$, where j_λ and α_λ are the emission and absorption coefficients, respectively, dependent nonlinearly on temperature and density. They are the measure of the addition and removal of the intensity to and from the light due to the local plasma conditions having local optical thickness τ_λ′. The implementation here represents a similar approach to that of Jenkins & Keppens (2022).

To qualitatively compare this work against observations, we synthesise images in three of the passbands observed by AIA, namely 171, 193, and 094 Å. These passbands were chosen as the simulation contains temperatures between 10³–10⁷ K and their temperature response functions cover a range peaking at ≈0.8, 1.5, and 7 MK, respectively, as shown in Figure 9. Plasma emitting around these temperatures will contribute intensity to the synthetic observation accordingly. The physical variables of plasma density and temperature were converted into these spectroscopic observables under a combination of the ‘coronal approximation’ and local thermodynamic equilibrium (LTE) (Verner et al. 1996; Keady & Kilcrease 2000). Generally, coronal emission will appear bright in the resulting images, whereas prominence material will appear dark. Additionally, photoionisation of hydrogen and helium by such emission (j_λ) is responsible for active absorption (α_λ) (Kucera 2015). As these simulations are in 2.5D and hence contain no background, we assumed a geometrical integration depth of 500 km, and a background illumination of 1 × 10¹, 0.5 × 10¹, and 5 × 10⁻² erg s⁻¹ cm⁻² sr⁻¹ for 171, 193, and 094 Å, respectively. Figure 10 then shows the synthesised prominence appearance for times = 85.87, 274.79, 429.37, and 680.12 s for the three aforementioned AIA passbands.

Fig. 9. Temperature response function of the different AIA passband filters as in Barnes & Chen (2022).

Fig. 10. Specific intensity counterparts to the simulation, where we account for emission and absorption. Synthetic images for the broadband 171 (left), 193 (middle), and 094 Å (right) SDO/AIA filters are shown for time = 85.87, 274.79, 429.37, and 680.12 s from top to bottom respectively.

From time = 85.87 s, the early development of the simulation describes two primary domains of darker (cooler) and brighter (mixed) material. In all passbands, the clear absorption of the prominence body captures the early nonlinear RT instability with dark falling fingers developing along the previously plane-parallel interface. The emitting features are that of prominence material that has already mixed with plasma initially at coronal temperatures, leading to dense plasma at hotter temperatures than in the prominence interior appearing bright in each of the passbands (cf. Hillier et al. 2023, and the PCTR phenomenon). By time = 274.79 s, rotating blobs that have separated from the cool and dense prominence layer create bright shear interfaces throughout the region between the prominence and the chromosphere. Then, at time = 429.37 s, it is well-captured how the falling prominence material interacts with the underlying chromosphere, with a significant portion being reflected back upwards and sidewards at the same time as inducing additional rotation. At this time, the plasma mixing caused by the combined nonlinear developments of the RT and KH instabilities leads to the development of additional bright fine structures for each of the passbands. By time = 680.12 s, the dual action of the KH instability and magnetic reconnection has isolated structures into identifiable ‘cells’, which can be contrasted to density views as seen in Figure 1.

To facilitate a qualitative comparison against observations, we have convolved the synthesis with the 0.6″ platescale PSF of each of the AIA passbands (Grigis et al. 2012; Barnes et al. 2020). The result, shown in Figure 11, significantly complicates any further interpretation, in particular for the jet-like feature of Figure 8. The bottom-right panel of Figure 11 compares the synthesis at both the AMRVAC and AIA resolutions for the 094 Å passband, with contours of M_A > 1 overlaid in cyan. At the AMRVAC resolution, the jet is visible in the 094 Å synthesis whereas at AIA resolutions all trace of it has been lost. The jet is entirely not visible in the 171 or 193 Å panels of the same figure, at either AMRVAC or AIA resolutions, as the thermodynamics of the jet are incompatible with the formation temperatures of these passbands. In all panels, the hotter material that bounds the individual rotating ‘cells’ remains visible, clearly demarcating distinctly different regions within the nonlinearly evolving RT and KH instabilities, in particular the largest of such examples.

Fig. 11. Synthetic intensity counterparts to the simulation for AIA 171, 193, and 094 Å at 680.12 s, convolved and de-resolved with the 0.6″ platescale PSF of the AIA instrument on board SDO. The lower row zooms in on the jet indicated by a high MA in Figures 8D–F (overlaid and contoured in cyan); the left and right panels of each passband are at the AIA 0.6″ and AMRVAC 0.02″ resolution, respectively.

Ground-based observatories commonly observe, amongst others, the strong Hydrogen n = 3 → n = 2 (Hα) line at 6563 Å. To synthesise the Hα emission and absorption we follow Heinzel et al. (2015), who presented a method to convert the local pressure and temperature values within the plasma to Hα opacity (α_λ), a relationship obtained from a comprehensive array of 1.5D radiative transfer models. Whilst an approximation, we have already demonstrated its accuracy in comparison to non-local thermodynamic equilibrium (NLTE) synthesis applied to more complex models (Jenkins et al. 2023, see also Osborne & Sannikov 2024). A core consideration within these models is the assumption of a constant source function of Eq. (7) along a given LOS, yielding a simplified equation,

$$\begin{array}{c}\hfill I_{\lambda }=I_{\lambda }(0)e^{-{\tau }_{\lambda }}+S_{\lambda }(1-e^{-{\tau }_{\lambda }}).\end{array}$$(8)

The resulting emergent (specific) intensity of the Hα line is therefore found through a LOS integration of the approximate α_λ according to the tables of Heinzel et al. (2015), wherein a coarse height-dependent estimation to the value of S_λ is also provided.

As already discussed, the dark features that are observed in the EUV passband images are due to the photoionisation of the cooler Hydrogen and Helium material, and these are the prominence falling pillars. This cool plasma with temperatures ∼10⁴ K instead appears bright in the equivalent Hα images. Figure 12, shows the evolution of our prominence according to Hα synthesised images for times 85.87, 274.79, 429.37, and 680.12 s. To first order, the locations of the dark features within the EUV synthesis are cospatial with the brighter structures in Hα. In particular for the snapshot at time = 680.12 s, however, the relationship between the EUV and Hα features is not so trivial, with overlap found for even some of the brighter EUV pixels. Indeed, depending on the moment within the simulation, some of the material rotating due to shear flows is simultaneously visible in both the EUV and Hα syntheses. As already discussed, however, the hot jet feature is not present within the Hα panel at time = 680.12 s. Interestingly, the shocks of Figure 1 are present throughout the prominence, evidenced by the formation of density compression layers that are subsequently encoded within the Hα emission due to its sensitivity to the local plasma pressure quantity (Heinzel et al. 2015). Depending on the applicability of the reduced Hα synthesis approach, and the 2D versus 3D simulation setup, similar features may be visible within observations (Jenkins et al. 2023).

Fig. 12. Synthetic images for the narrowband Hydrogen-Hα filter, for a LOS view into the plane of the simulation at time = 85.87, 274.79, 429.37, and 680.12 s.

4. Discussion

In this paper, we have performed a resistive extension to the 2.5D simulation of MC23 and detailed the subsequent evolution considering both theory and observations. MC23 advocated for an immediate 3D extension to the prominence model, but we here maintain the 2.5D geometry so as to study the resistive behaviour in isolation from additional out-of-plane dynamics. The prominence here evolves from the same initial multimodal perturbation triggering the RT instability and leading to the development of distinct ‘falling fingers’ and ‘rising plumes’. The subsequent interaction of the falling fingers with the underlying chromosphere strongly influences the overall evolution, generating a significant horizontal momentum and forcing the collision of adjacent structures (in contrast to Rees-Crockford et al. 2024). Secondary processes are then formed, with both the KH instability and magnetic component reconnection driving additional dynamics. For an identical initial condition and numerical scheme, the markedly different distributions of plasma when compared to MC23 (see their Figure 3) is attributed to the inclusion of finite resistivity throughout the domain. For an equivalent timestep, the model presented here yields more coherent structures, whether in ‘fingers’ or ‘plumes’. This distinction is a function of the magnitude of η chosen; a necessarily large value of 10⁻⁴ (≫10⁻⁹) was set so as to resolve its influence on the discrete grid, consequently introducing a minor diffusive smoothing of the thermodynamics that becomes more noticeable over longer timescales (cf. Ripperda et al. 2019). We note that compared to MC23, which only had numerical resistivity at play, we here targeted actually resolved resistive evolutions (at unrealistically low magnetic Reynolds number), using simulations with a higher effective resolution.

In quiescent prominences, the primary RT instability initiates upflows and downflows, successfully reproduced here in our simulation as shown in Figure 2 (Hillier et al. 2012a,b; Jenkins & Keppens 2022; Donné & Keppens 2024). The transition to turbulence in prominence simulations is marked by the emergence of secondary instabilities, such as the KH instability, which facilitates the turbulent energy cascade (previously explored in observations Leonardis et al. 2012; Hillier et al. 2017). In our simulation, the interaction of gravity and buoyancy induces shear flows at the edges of vertically dominated flows, leading to the presence of complex motions including KH instability-driven shear flows as shown in Figures 3 & 4. As turbulence further develops and fine-scale structures form through secondary processes, short-lived jet exhausts are found within our simulation, shown in Figure 8, indicating that these shear flows distort the magnetic field and promote the formation of current sheets and magnetic reconnection Hillier & Polito (2021). In our simulations, the formation of these infrequent jets facilitates short lived energy transfer and release in the upper solar atmosphere, where the plasma β is significantly higher than in lower regions (as shown in MC23, Figure 3).

The time and lengthscales of the secondary instabilities recovered from our simulation appear to agree with those derived from observations (Berger et al. 2010; Hillier & Polito 2018; Yang et al. 2018). As shown in Figure 4, however, the majority of the structures or dynamics of interest are present in the hotter coronal or intermediate PCTR environments, rather than the cooler prominence material. Given that these structures have been predominantly recorded in observations of prominences in the cooler, chromospheric and transition region lines, this places our results in direct contrast with previous results from observations. This is further emphasised in the following synthetic analysis, after which we will suggest potential extensions to the base model to address these discrepancies.

The results presented in our paper correspond to a very high-resolution resistive MHD study. In Section 3.3, we presented the synthetic observation equivalents to the simulation. Herein, all of the considered passbands described clear signatures of the ‘cells’ demarcating the discrete bubbles containing rotational motions (as found in the observations of e.g. Berger et al. 2017; Rees-Crockford et al. 2024; Zhang et al. 2024). In general, the appearance of the prominence in the 171 and 193 Å passbands are qualitatively similar, in contrast to the 094 Å passband where the prominence absorption is expectedly weak. The jet feature of Section 3.2.3 is, nevertheless, well captured in the latter. However, when convolving these synthetic observations at a resolution of 0.02″ with the PSF and 0.6″ platescale of the AIA instrument, all of the passbands considered lose a significant amount of the finescale detail. Indeed, Figure 11 demonstrates that the resolution of the SDO/AIA instruments is insufficient to recover the smallest scales, yet the larger cells with radii of several Mm remain visible. Based on these syntheses, it appears that particularly bright, curvilinear intra-prominence sheets between the dark falling fingers are strongly correlated to locations of shear flows and/or current sheets. Geometric considerations between this 2.5D representation at that of the 3D Sun will complicate any direct comparison. Nevertheless, given the much higher spatial ≈100 km perihelion resolution of Solar Orbiter’s Extreme Ultraviolet Imager (EUI; Rochus et al. 2020) High Resolution Telescope (HRT) instrument, instances of this should be identifiable and verifiable (cf. Jenkins & Keppens 2022).

The approximately linear dependency of the hydrogen Hα absorption coefficient α_λ on pressure and density (Heinzel et al. 2015) leads the associated synthesis to match the primitive distributions almost 1–1, see Figure 1. Consequently, the roll-ups associated with the collision between the falling fingers and the underlying chromosphere is well captured. Once again, the geometry of the simulation at hand, the definition of the chromosphere, and the adopted magnetic field will influence how such an evolution would manifest within the 3D atmosphere. Despite these limitations, these results suggest a serendipitous observation may indeed be able to capture such an evolution if the material is not dispersed before reaching chromospheric heights (as also suggested by Kaneko & Yokoyama 2018).

The motivating observations of Hillier & Polito (2021) were recorded using IRIS, whose slit-jaw imager (SJI) has a spatial resolution higher than SDO/AIA but still 10× lower than the present simulation (MC23). The spatial extent of the simulated jet is once more too small to be directly imaged. More fundamentally, however, the formation temperature of either the Mg II or Si IV lines observed with SJI are several orders of magnitude lower (10^4 − 5 K) than recorded in SDO/AIA 094 Å (10^6 − 7 K). Hence, although a jet was self-consistently formed within our simulated solar prominence undergoing a combination of nonlinear RT and KH instabilities, it is clearly not of the same origin as those that were observed by Hillier & Polito (2021). Indeed, it is also completely absent from the cooler hydrogen Hα dynamics. To first order, the location of the jet with respect to dense plasma is responsible for the difference in appearance, that is, the jet occurred within more ‘coronal’ plasma already at > 10⁶ K. Should reconnection occur internal to a cool, dense plasma concentration then, even though the reconnection site is almost certainly too small to be directly imaged, the subsequent heating may be observed. Similar to how it appears the passage of shocks within the simulation are visible within the hydrogen Hα syntheses. Recently, Jerčić et al. (2024) performed 1.5D non-local thermodynamic equilibrium (NLTE) spectral synthesis to demonstrate that a reconnection event internal to a cool ∼10⁴ K prominence condensation yields a significant increase in the emission of the Mg II line core.

There remain a number of improvements that can be made to the current simulation. The 2.5D domain with a constant magnetic field such that k ⋅ B ≈ 0 has here enabled us to achieve a high spatial resolution, resolve the numerical resistivity, and hence reproduce RH and KH instability dynamics at those scales (cf. Kaneko & Yokoyama 2015; Jenkins & Keppens 2021). It is imperative, however, that we explore the linear and nonlinear developments of these instabilities, including subsequent secondary processes, under more physically realistic conditions.

Considering the thermodynamics, this current model neglects the non-adiabatic terms such as radiative losses that may lead to the formation of additional, cooler material via the aforementioned KH mixing (Hillier et al. 2023). This is of particular interest as the secondary processes we find here are located within these mixing regions. For the magnetic field, we should extend the domain to consider magnetic flux ropes. Previously, Popescu Braileanu et al. (2021) explored the linear development of the RT instability within a 2.5D domain that included magnetic shear so as to approximate the internal topology of flux ropes, but as shown in Jenkins & Keppens (2022) and Donné & Keppens (2024) prominences are anything but invariant in the third dimension. Both authors demonstrated clearly how plasma processes throughout a flux rope can influence the development in any given k ⋅ B ≈ 0 plane, as a simple consequence of condensations being localised but free to move along field lines. Furthermore, and as discussed in Section 3.2, the 2.5D description artificially suppresses any influence that the magnetic tension may have. This is certain to play a larger role in the less dense environments, currently dominating the KH unstable populations, than the more dense prominence fingers.

Finally, the syntheses presented in Section 3.3 consider a LOS perpendicular to the plane, that is, along the invariant direction. An actual observation contains structure in this direction and would invariably alter the appearance, even if the underlying plasma distribution resembled closely our simulation (cf. Gunár et al. 2018). Whether for the underlying physics or the subsequent observational synthesis and as similarly concluded in MC23, the computationally advantageous yet nevertheless simplified 2.5D domain setup commonly employed to study the internal dynamics of prominences is fast becoming insufficient. We are already actively engaged in a 3D extension to this work.

5. Conclusions

In this paper, we have used 2D high-resolution resistive MHD simulations to analyse the spatio-temporal dynamics of a quiescent prominence. The system starts from an initial state where the prominence boundary is set up within the box and the RT instability is set to self-consistently evolve the system from equilibrium. The initial perturbation leads to the formation of rising plumes and falling bubbles forming coherent structures commonly present in observations of quiescent prominences. The addition of resistive dynamics leads to the formation of secondary instabilities and the growth of nonlinear dynamics which in turn excites energetic jet-like events in the prominence environment. A combination of these events forming outside of relevant temperature domains and insufficient observational resolution, comparatively, currently renders a direct 1–1 comparison impossible. Further study is necessary to map the available parameter space, with future endeavours focused on the step-up to 3D so as to excite yet-more realistic evolutions.

Data availability

Movie associated to Fig. 4 is available at https://www.aanda.org

Acknowledgments

We acknowledge the comments and suggestions by the anonymous referee that significantly improved the accuracy, completeness, and thus relevance of this work. This work is supported by FWO grant G0B9923N Helioskill and by Internal funds KU Leuven, through the project C16/24/010 UnderRadioSun. JMJ is supported by a European Space Agency (ESA) Internal Research Fellowship.

References

All Figures

	Fig. 1. Combined plots of log ρ and $\frac{L}{c_{\mathrm{s}}}\mathrm{\nabla }·\mathbf{v}$ at time 85.87 s, 429.37 s, and 618.29 s (left, middle, right) showing the evolution of bow shocks in the evolution of the solar prominence. The horizontal red dashed line indicates the 2 Mm threshold relevant for Figure 2.
In the text

Fig. 2. Time evolution of the cool and dense solar prominence matter above the chromosphere boundary (above 2 Mm from the bottom boundary) taken below a threshold temperature of 5 × 104 K, and denser than 10−11 kg m−3. The dashed line represents the phase of the initial decrease of mass of the cool and dense volume in the solar prominence due to the merging of the cool dense material into the chromosphere boundary. The dotted line represents the phase where the prominence material rises towards the top boundary with a significant loss of the prominence material. The solid line represents the evolution at time = 429.37 s (see Figure 1 (middle)).

In the text

Fig. 3. Evolution of secondary KH instability in the prominence at time = 274.79, 291.97, and 343.49 s as shown in the red insets. The top row shows the zoomed image for the KH instability regions, represented as both vorticity ∇ × v and Rortex R, identified using the velocity magnitude in the middle row for each time. The bottom insets show the zoom region for time = 291.97 s, representing the vy and vx component of the velocity field at that time. Together, they signal a trail of (rising) vortex structures being formed.

In the text

Fig. 4. Condition of the plasma undergoing rotational motion within the domain, as highlighted according to the R criterion. Left; the 2.5D domain of temperature with positions of R > 5 × 10−2 s−1 and where λci > 10λcr overlaid in white, see text for description of these terms. A few clear KH instability cells have been indicated with teal boxes. Right; the density – temperature phase space of these extracted regions. A broad distribution is present, with a clear peak at the high temperatures and low densities of the solar corona. The inset axis plots the number of distinct KH cells in time (s), current time indicated by the dashed-black line, explanation in the text. An animation of this Figure will be available online.

In the text

Fig. 5. Detection of current sheets using the modified algorithm based on current density Jz at time 680.12 s. The left plot shows the current density magnitude, the middle plot shows the correspondence of local maxima events using a threshold higher than the mean |Jz| with the current density (|Jz|), and the right plot shows the identification of clusters derived from density clustering methodology using local maxima events at time 680.12 s. 123 separate clusters could be identified.

In the text

Fig. 6. Comparison of different threshold parameters of Jthreshold at time 680.12 s in correlating the local maxima events with the current density (|Jz|). The red current sheets represent the positive current sheets whereas the blue represent the negative ones. Depending on the threshold used, the algorithm may overlook significant current sheet regions to fulfil the parameter (right) or overfill the local maxima regions (left).

In the text

Fig. 7. Top: Scatter plot of all current sheets characterised in length being the longest dimension (horizontal) versus width (vertical) for three times: 73 clusters (positive) and 65 clusters (negative) for t = 671.53 s, 75 clusters (positive), and 62 clusters (negative) for t = 680.12 s, and 73 clusters (positive) and 70 clusters (negative) for t = 704.16 s (blue, green, and red, respectively). The ‘+’ refers to the positive current sheets while ‘−’ refers to the negative current sheets. Most of the current sheets are smaller ranging between 631–1388 km in length and of 61–78 km in width. The largest current sheet analysed by the algorithm in the simulation goes up to ∼11 355 km in length and around ∼116 km in width. Bottom: Statistical distribution of current sheet lengths and widths at the same times, shown as violin plots for positive (red) and negative (blue) current sheets. The width of each violin indicates the probability density estimated by kernel density estimation, while the solid orange horizontal lines mark the medians. The ‘L’ and ‘W’ labels denote distributions for length and width, respectively. Shading is used to visually separate the three time instances. These distributions reveal asymmetries between positive and negative sheets in both length and width, highlighting temporal evolution and polarity dependence in current sheet morphology of the simulation.

In the text

Fig. 8. Time series of the occurrence of jets in the prominence identified by regions where MA > 1 at time = 671.53, 680.12, and 704.16 s. Panels (A–C) show the zoomed images of Panels (D–F) and the occurrence of jets, here essentially flowing upwards. Panels (G–L) show the time evolution of current density (|Jz|) and magnetic field (|Bplane|) for zoomed regions where the jets emerge. Both the current sheets and subsequent magnetic islands responsible for the jets are clearly visible.

In the text

	Fig. 9. Temperature response function of the different AIA passband filters as in Barnes & Chen (2022).
In the text

	Fig. 10. Specific intensity counterparts to the simulation, where we account for emission and absorption. Synthetic images for the broadband 171 (left), 193 (middle), and 094 Å (right) SDO/AIA filters are shown for time = 85.87, 274.79, 429.37, and 680.12 s from top to bottom respectively.
In the text

Fig. 11. Synthetic intensity counterparts to the simulation for AIA 171, 193, and 094 Å at 680.12 s, convolved and de-resolved with the 0.6″ platescale PSF of the AIA instrument on board SDO. The lower row zooms in on the jet indicated by a high MA in Figures 8D–F (overlaid and contoured in cyan); the left and right panels of each passband are at the AIA 0.6″ and AMRVAC 0.02″ resolution, respectively.

In the text

	Fig. 12. Synthetic images for the narrowband Hydrogen-Hα filter, for a LOS view into the plane of the simulation at time = 85.87, 274.79, 429.37, and 680.12 s.
In the text