© ESO 2021

1. Introduction

Optically thin radiative emission is an important physical process. The energy loss due to recombination of ions can cause a significant decrease in the temperature of a plasma. Hence it plays a key role in the formation of condensations in various fields of astrophysics.

Examples of condensed structures in solar physics are prominences and coronal rain. Prominences are cool and dense structures that are suspended in the hot and tenuous corona by the magnetic field (Parenti 2014). Coronal rain is a more transient and recurrent phenomenon, with a lifetime of only minutes to days. After being formed at the top of magnetic loops, the dense blobs rain down along field lines in cycles (Antolin 2020; Froment et al. 2020). In the interstellar medium, turbulent flows can lead to the formation of dense molecular clouds, from which stars can form in turn (Audit & Hennebelle 2005). Condensations are also observed as neutral-atomic and molecular outflows in galaxies (Veilleux et al. 2020) and as filaments in the intergalactic medium (Sharma 2018).

The modelling of the thermodynamics of such phenomena is not straightforward. Solving the full radiative transfer equations alongside the magnetohydrodynamics (MHD) equations is possible. An example is the Bifrost code (Gudiksen et al. 2011), which has recently been used to study coronal rain formation (Kohutova et al. 2020). However, this is computationally demanding for multidimensional simulations.

Under several assumptions, especially when dealing only with optically thin radiation, the density and temperature dependencies of the radiative energy loss can be separated (Dopita & Sutherland 2003). The temperature dependence is typically calculated in advance and tabulated as a cooling curve. A large variety of cooling curves exists in the literature (see e.g. Cox & Tucker 1969; Hildner 1974; Rosner et al. 1978; Colgan et al. 2008; Dere et al. 2009; Schure et al. 2009), and in our appendix, we provide an overview of the most frequently employed cooling curves, which we have all implemented in our open-source software MPI-AMRVAC (Keppens et al. 2012, 2020; Porth et al. 2014; Xia et al. 2018). Differences arise due to different assumptions in the calculation, varying atomic physics parameter, and the solar abundances that are used (Rosner et al. 1978). The effect of improved atomic physics has been proven, for example, in the field of asteroseismology. Advances in opacity tables for elements such as iron and nickel have led to an explanation of the excitation of pulsations in B-type stars (Dziembowski 1994).

In this work we study the effect of varying the optically thin cooling curve on the formation of condensations. Various physical processes can lead to condensations. We form in situ condensations by thermal instability (TI), first described by Parker (1953) and Field (1965). The theory was recently revisited by Claes & Keppens (2019) and was simulated in an idealised setup. This setup of a local typical solar coronal volume showed that condensations arise naturally due to optically thin cooling and a perturbation, independent of the origin of the perturbation. Claes et al. (2020a) included the effects of anisotropic thermal conduction and showed that the fragmentation of the condensations and the magnetic field lines is misaligned. This complicates the interpretation of observed fine structure in solar prominences, for instance. Soler et al. (2012) investigated the stability of thermal modes in cool prominence plasmas for three typical cooling curves.

Thermal instability is a very likely mechanism to form solar prominences, as was recently shown in numerical simulations by Jenkins & Keppens (2021). A review by Mackay et al. (2010) discussed a variety of mechanisms to form prominences, not necessarily invoking thermal instability. The formation of other condensed structures, such as filaments in the interstellar and intergalatic medium are also being explained by thermal instability (see e.g. Cox 1972; Sharma et al. 2010). Recently, Falle et al. (2020) revisited thermal instability in the context of shock-induced condensation formation in the interstellar medium as a process to generate density inhomogeneities without relying on self-gravity. Jennings & Li (2021) studied cloud formation in the interstellar medium by thermal instability using Gaussian random fields as initial conditions and including magnetic fields, thermal conduction, and viscosity. They also discussed fragmentation processes for hydrodynamic clouds. Furthermore, Waters & Proga (2019a) investigated the possible regimes of thermal instability in hydrodynamics. They discussed the stability of the entropy and acoustic modes between the classical isobaric and isochoric limits. They also performed 1D numerical simulations of the non-linear evolution of the thermal instability and discussed fragmentation mechanisms. With respect to this, it has been shown that clumps in active galactic nucleus (AGN) winds can be explained by thermal instability (Dannen et al. 2020; Waters et al. 2021). It should be noted that the effect of typical AGN net cooling curves, which depend on the temperature and photoionisation parameters, on the stability of plasma has been studied in the recent past (see e.g. Chakravorty et al. 2008, 2009, 2012).

The results of this case study using thermal instability might also be applicable to other mechanisms that rely heavily on optically thin radiative cooling. In thermal non-equilibrium models (TNE), evaporation and condensation cycles can create coronal rain (Antolin 2020). A difference between TNE and thermal instability is that most studies on thermal instability start from a thermal equilibrium, which is then perturbed (Klimchuk 2019). In TNE studies, the focus lies instead on time-varying conditions that may never be in perfect thermal equilibrium, and the non-linear evolution shows runaway cooling, denoting the fast decrease in temperature due to optically thin cooling. In both processes the cooling, and thus the formation of condensations, is dependent on the cooling curves.

In this work we use a simplified setup of thermal instability with solar coronal conditions, following the method of Claes & Keppens (2019) and Claes et al. (2020a). The initial thermal equilibrium will be perturbed by slow MHD waves. Slow MHD waves are commonly observed in the solar corona (see Banerjee et al. 2020, and references therein) and are therefore a good choice for the perturbation.

In Sect. 2 we discuss the optically thin cooling curves. In Sect. 3 we present the setup, initial conditions, and numerical methods. The benchmark simulation is discussed in Sect. 4, and a comparison is made to the setup and results of Claes et al. (2020a). The effect of several numerical settings, such as the resolution and the detailed low-temperature treatment, are investigated in Sect. 5. We compare simulations with different cooling curves in Sect. 6. Using the insights gained in preceding sections, we follow the evolution of the thermal instability far into the non-linear regime in Sect. 7. We present a summary of the results and discussion in Sect. 8. Appendix A gives references to all cooling curves currently available in the MPI-AMRVAC software we use. Appendix B contrasts the low plasma beta MHD case studied in this manuscript with a related hydrodynamics case.

2. Optically thin radiative cooling curves

The energy loss due to optically thin radiative cooling is typically written as

$$\begin{array}{c}\hfill \frac{\partial e}{\partial t}=-n_{\mathrm{i}}{n}_{\mathrm{e}}\mathrm{\Lambda }(T),\end{array}$$(1)

where n_e, n_i, and T are the local electron and ion number density and the temperature of a plasma, respectively. The variable e denotes the internal or total energy density, and the function Λ(T) is the predefined cooling curve that is also called the cooling rate. This equation is also commonly written in normalised variables as

$$\begin{array}{c}\hfill \frac{\partial \overline{e}}{\partial \overline{t}}=-{\overline{\rho }}^2\overline{\mathrm{\Lambda }}(\overline{T}),\end{array}$$(2)

where $\overline{\rho }$ is the normalised local density. All variables in this equation are normalised, in contrast to the variables in Eq. (1), as denoted by the overlying bar. To normalise the equations, we use the unit values of the quantities of the initial state given in Table 1. All other unit quantities can be derived from the chosen three quantities and are also given in Table 1. As an example, the normalisation of the density from ρ to $\overline{\rho }$ and the definition of ρ_u is given by

$$\begin{array}{c}\hfill \overline{\rho }=\frac{\rho }{{\rho }_{\mathrm{u}}}=\frac{\rho }{(1+4A_{\mathrm{He}})m_{\mathrm{p}}{n}_{\mathrm{u}}},\end{array}$$(3)

where ρ_u, n_u, m_p, and A_He are the unit density, the unit number density from Table 1, the proton mass, and the assumed hydrogen-helium abundance ratio, respectively. We used a hydrogen-helium abundance ratio of 0.1. The factor four arises from mass conservation, that is, the fact that a helium ion has four times more mass than a hydrogen ion. The cooling curve can be normalised as

$$\begin{array}{c}\hfill \overline{\mathrm{\Lambda }}=\frac{\mathrm{\Lambda }}{{\mathrm{\Lambda }}_{\mathrm{u}}}=\mathrm{\Lambda }\frac{t_{\mathrm{u}}{n}_{\mathrm{u}}^2(1+2A_{\mathrm{He}})}{p_{\mathrm{u}}},\end{array}$$(4)

where Λ_u, t_u, and p_u are the unit cooling rate, unit time, and unit pressure, respectively. The last two can be derived similarly. The factor two arises from charge conservation. Unless stated otherwise, all equations are written in normalised variables, hence we drop the overlying bar.

The energy loss is due to emission, such as line and continuum radiation and bremsstrahlung. It is typically calculated under the assumption of collisional ionisation equilibrium (CIE), meaning that collisional ionisation balances the recombination processes. A necessary condition is that the plasma is optically thin, so that all radiation produced within is able to escape and/or is fully ionised so that photoionisation can be ignored. Furthermore, it is commonly assumed that all excited ions return to their ground state faster than the time between collisions is long, so that only the ground state of each ionic species is populated (Dopita & Sutherland 2003). However, this is not a necessary assumption. Cooling curves can also be constructed by solving the full collisional-radiative rate equation, that is, calculating the populations of all the excited levels (see e.g. Colgan et al. 2008). Under these assumptions, the density and temperature dependences can be separated. A cooling curve thus represents the dependence of the energy loss on the local temperature of the medium.

In addition to the specific method adopted within the computation, the accuracy of a cooling curve depends on two other things: the accuracy of the atomic physics parameters, and the accuracy of the assumed solar abundances (Rosner et al. 1978). The effect of several parameters, such as solar abundances, ion fractions, electron densities, and transition probabilities, has been investigated (Landi & Landini 1999). The abundances seem to be particularly important, although other parameters involved can also have a strong effect on a cooling curve.

We compared five different radiative cooling curves for optically thin cooling that were used in past or present solar physics research. We briefly describe these cooling curves in this section. In some works that studied solar phenomena, such as flares and prominences (see e.g. Forbes & Malherbe 1991; Kaneko & Yokoyama 2015) or thermal instability in molecular cloud formation (e.g. Sharma et al. 2010), cooling curves were modified for numerical reasons, such as to speed up the growth of the thermal instability or to increase the characteristic length scale of thermal conduction. We did not investigate these particular changes, but rather compare the effects of changing the entire cooling curves.

Before we discuss the five cooling curves, we address the numerical implementation of a cooling curve. The Λ(T) prescriptions typically come in two forms: interpolatable tables, or piece-wise power laws. As the name suggests, an interpolatable table is a set of data points, typically 50 to 100 points, that can be interpolated with a higher temperature-resolution to a more extended table at the start of a simulation. During runtime, the values of the cooling rate at a given temperature are selected from this table. The piece-wise power laws are analytical expressions. Their cooling rate is just calculated based on the given temperature. They are coarser than the interpolatable tables because they typically consist of fewer than ten segments. The piece-wise power laws have the advantage that they can easily be used in analytic calculations, such as in 0D models of coronal rain (Cargill et al. 1995). They have also been used in early studies of magneto-thermal instabilities in non-uniform coronal plasmas (van der Linden et al. 1992), which have recently been revisited with the linear MHD code Legolas (Claes et al. 2020b).

We now describe the five cooling curves. For details about their calculations, we refer to the accompanying papers. Appendix A gives an overview of the various available curves, with the original references. The first three cooling curves are interpolatable tables. The first cooling curve is the SPEX_DM table as described in Schure et al. (2009). The main temperature regime, 10^3.8 to 10^8.16 K, was calculated with the SPEX package (Kaastra & Mewe 2000). The solar abundances of Anders & Grevesse (1989) were used. For temperatures below 10 000 K, the radiative cooling is approximated by the cooling curve first described by Dalgarno & McCray (1972), but implemented as presented by Schure et al. (2009). The cooling of this _DM part is due to the collisional excitation of singly charged ions of O, C, N, Si, Fe, Ne, and S with thermal electrons and neutral hydrogen. There is also a contribution by collisions between neutral hydrogen with electrons. An ionisation fraction of 10⁻³ is assumed for the low-temperature cooling curve. The SPEX_DM table is also the reference cooling curve in this work.

The Colgan_DM cooling curve is given in Colgan et al. (2008) as their ‘DW ground’ model. It is constructed assuming that only the ground state of ions is populated. The atomic and collisional calculations were performed with the Los Alamos plasma kinetics code ATOMIC. The abundances of Grevesse & Sauval (1998), modified by Feldman & Widing (2003), were used. This table is extended for low temperatures by the _DM part described above.

The third cooling curve, Dere_corona_DM, is described and given in Dere et al. (2009). This table was constructed with version 6 of the CHIANTI atomic database (Dere et al. 1997). They used the abundances from Feldman (1992). This table is also extended by the _DM part for low temperatures.

The first piece-wise power law is the cooling curve by Rosner et al. (1978). It is actually an analytical fit to an unpublished interpolatable table by J. Raymond. It uses the solar coronal abundances that were accepted at that time, but the specific details are not known. We used a version that is extended for high and low temperatures by Priest (1982).

The last cooling curve is the piece-wise power law by Klimchuk et al. (2008). Although the cooling curve was described in 2008, it was based on much older physical data and calculations. The atomic physics calculations were provided by J. Raymond by means of priv. comm. in 1994. The coronal abundances of Meyer (1985) were modified by a factor of two.

These five cooling curves are depicted in Fig. 1 on a log–log scale. There are many differences between the curves. The two tables that differ the least are the Colgan_DM and the Dere_corona_DM. This is expected because they both use current values for the abundances. As discussed by Dere et al. (2009), this might be due to the different amount of atomic fine structure that is used. The piece-wise power laws are much coarser and hence do not represent all variations in this detail. A striking difference between the Colgan_DM and Dere_corona_DM and the other three curves is that the last three have a much lower cooling rate from 10⁶ K and up. This is mostly because these three cooling curves used older solar abundances. The difference is almost the same as between solar corona and photospheric abundances, as reported by Dere et al. (2009). In the 1980s, it became known that the quiet-Sun upper atmosphere can have a fourfold increase in abundance of low first-ionisation potential elements, such as Fe, Si, and Mg, with respect to the solar photosphere (Feldman 1992). It might therefore be argued that cooling curves based on older abundances are less representative of the optically thin radiative energy loss in the solar corona. Nevertheless, because all five curves have been (and still are) used in contemporary MHD simulations of the solar corona, we here set forth to study a simple idealised setup that truly isolates their effect on the formation of condensations.

Fig. 1. Five cooling curves that we compare.

Several comparisons can be made without multidimensional MHD simulations. One of them is the linear cooling timescale τ_r, also called radiative timescale. It describes the time it takes for a parcel of plasma of a given density and temperature to cool. For a given temperature T, it can be written as

$$\begin{array}{c}\hfill {\tau }_r(T)=\frac{T}{(\gamma -1)\rho \mathrm{\Lambda }(T)},\end{array}$$(5)

where ρ and Λ are the density and cooling rate, respectively. We take γ, the ratio of specific heats, to be 5/3 by assuming a monatomic ideal gas.

We can also define a non-linear radiative timescale as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{nl}}(T)=\frac{T_{\max }\mathrm{\Lambda }(T)}{T\mathrm{\Lambda }(t_{\max })}{\tau }_r(T),\end{array}$$(6)

where T_max is the temperature at which the cooling curve has a maximum (Forbes & Malherbe 1991). The physical interpretation of the timescales are depicted in Fig. 2, which shows the solution of Eq. (2) as a function of time for a fixed density and initial temperature. The radiative timescale is the time it takes to cool from a starting temperature, while the non-linear timescale is the time of the non-linear part of the cooling.

Fig. 2. Left: general thermal evolution of thermal instability and its timescales. Right: linear and non-linear timescales of five cooling curves. The linear timescale is depicted by the length of the full bar, and the non-linear part is represented by the red part.

For typical solar coronal conditions, that is, a temperature of 10⁶ K and a density of 2.34 × 10⁻¹⁵ g cm⁻³, these two timescales are shown in Fig. 2 for the five cooling curves. The full length of the bar denotes the radiative timescale, and the red part denotes the non-linear timescale. A first observation is that the length of the timescale corresponds to the cooling rate, as evidently follows from Eq. (5). There are large differences between the radiative timescales of the different cooling curves. The timescales for the Colgan_DM and Dere_corona_DM are quite similar, which is also seen in the similarity of their cooling curves. This suggests that they will give similar results in numerical simulations. However, in full MHD simulations, other non-linear and dynamical effects will be present. Therefore these values can only be used as a guideline. For example, from this we can assume that numerical simulations using the Rosner table will take much longer than for the other tables if we are interested in obtaining thermally driven condensations.

3. Numerical setup

We study the effect of using different cooling curves on condensation formation by thermal instability. To initiate the thermal instability, we follow the method from Claes & Keppens (2019) and Claes et al. (2020a).

3.1. Governing equations

We make use of the MPI-parallelised Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC)¹ (Keppens et al. 2012, 2020; Porth et al. 2014; Xia et al. 2018) to numerically solve the following set of normalised MHD equations:

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

$$\begin{array}{cc}& \frac{\partial \rho \mathit{v}}{\partial t}+\mathrm{\nabla }·(\rho \mathit{v}\mathit{v}+p_{\mathrm{tot}}\mathit{I}-\mathit{B}\mathit{B})=0,\hfill \end{array}$$(8)

$$\begin{array}{cc}& \frac{\partial e}{\partial t}+\mathrm{\nabla }·(e\mathit{v}+p_{\mathrm{tot}}\mathit{v}-\mathit{B}\mathit{B}·\mathit{v})=-{\rho }^2\mathrm{\Lambda }(T)+H,\hfill \end{array}$$(9)

$$\begin{array}{cc}& \frac{\partial \mathit{B}}{\partial t}+\mathrm{\nabla }·(\mathit{vB}-\mathit{Bv})=0,\mathrm{\nabla }·\mathit{B}=0,\hfill \end{array}$$(10)

where I is the unit tensor, with the magnetic field in units, where μ₀ = 1. The quantities ρ, T, v, and B are the density, temperature, velocity, and magnetic field, respectively. The total energy density e is given by

$$\begin{array}{c}\hfill e=\frac{p}{\gamma -1}+\frac{\rho v^2}{2}+\frac{B^2}{2},\end{array}$$(11)

and the total pressure is given by

$$\begin{array}{c}\hfill p_{\mathrm{tot}}=p+\frac{B^2}{2}·\end{array}$$(12)

The relation between the density and temperature is described by the ideal gas law, thus closing our set of equations. We considered a fully ionised plasma of hydrogen and helium with an abundance ratio of ten hydrogen atoms for each helium atom. This ratio was used in the normalising of the equations and variables.

To keep the setup as simple as possible, we neglected gravity, resistivity, conduction, and other non-adiabatic effects, except for optically thin radiative cooling and a constant background heating, denoted by H. The optically thin radiative cooling depends on the local density and cooling rate of the prescribed cooling curves. In order to have an initial state in thermal equilibrium, we took the background heating equal to the radiative loss of the initial state,

$$\begin{array}{c}\hfill H={\rho }_0^2\mathrm{\Lambda }(T_0),\end{array}$$(13)

where ρ₀ and T₀ are the equilibrium density and temperature, respectively. The background heating is therefore constant. However, it should be noted that this background heating H is not the heating per unit mass, which is part of the net heat-loss function, denoted by ℒ, as typically used in analytical research of thermal instability and as defined by Field (1965). The net heat-loss function is here given by

$$\begin{array}{c}\hfill \mathcal{L}=\rho \mathrm{\Lambda }(T)-\frac{{\rho }_0^2\mathrm{\Lambda }(T_0)}{\rho }·\end{array}$$(14)

Having the initial state in thermal equilibrium also allows direct comparisons between predicted growth rates of the thermal instability and the onset stages of the condensations, as demonstrated previously by Claes & Keppens (2019).

3.2. Initial conditions

To simulate the in situ condensation-forming process by thermal instability, we set up a local 2D coronal volume with periodic boundary conditions. We started from a thermal equilibrium state that was nearly uniform, except for a superposed wave solution. The equilibrium number density and temperature were used together with the characteristic length scale L of the simulation box to normalise the equations and variables. Their values are given in Table 1. The equilibrium values for the density and the temperature are thus equal to unity in normalised notation.

The thermal mode, which is the entropy mode, can be excited directly, but here we followed Claes & Keppens (2019) and used two interacting slow MHD modes to perturb the thermal equilibrium and facilitate the thermal instability. The slow waves become damped due to the radiative cooling, while the thermal mode is unstable for the physical conditions in Table 1, as was shown by Claes & Keppens (2019). The isobaric and isochoric instability criteria are both satisfied for the given initial conditions. Because the plasma has been changed by the interacting, damped slow wave propagation, it is not trivial to categorise the linear thermal mode as a pure isobaric or isochoric mode, although the observed growth rate matches analytic predictions reasonably well (Claes & Keppens 2019).

We excited one slow wave along each axis. The wave number, k, was equal to 2π divided by the domain length L to satisfy the periodic boundary conditions by having one full oscillation within the domain.

The direction of the background magnetic field was taken to be along the diagonal of the square domain. This was not the case in Claes et al. (2020a), who allowed it to vary for different simulations. The implications and differences are discussed in Sect. 4. The background magnetic field had a strength of 10 G, and the medium had a plasma beta of approximately 0.08, which are typical for prominences in the solar corona (Gibson 2018).

The eigenfrequency ω of the slow modes does not include the effect of the cooling curve on the initial condition. We used the adiabatic slow MHD eigenfrequency as can be found in any textbook on MHD (see e.g. Goedbloed et al. 2019). It can be calculated as follows:

$$\begin{array}{c}\hfill {\omega }_{\mathrm{s}}=k\sqrt{\frac{1}{2}(v_{\mathrm{A}}^2+c^2)-\frac{1}{2}\sqrt{{(v_{\mathrm{A}}^2+c^2)}^2-4c^2{v}_{\mathrm{A}}^2{cos}^2\theta }},\end{array}$$(15)

where the sound- and Alfvén speeds are defined by

$$\begin{array}{c}\hfill c^2=\frac{\gamma p_0}{{\rho }_0}\mathrm{and}{v}_{\mathrm{A}}^2=\frac{B_0^2}{{\rho }_0}·\end{array}$$(16)

The angle θ is defined between k and B and had a value of 45° for all simulations in this work because we set up the magnetic field along the diagonal of the square domain and then adopted one slow mode travelling along the x-axis and one slow mode travelling along the y-axis.

Lastly, the slow MHD waves were initiated by perturbation of the thermal equilibrium. The eigenfunctions were derived and described by Claes et al. (2020a). The initial velocity perturbation was taken as a plane wave and given by

$$\begin{array}{c}\hfill v_x=A[cos(k_{x}x)+isin(k_{x}x)],\end{array}$$(17)

with A an amplitude taken to be 10⁻³. The velocity perturbation along the y-axis is similar. The real part of complex eigenfunctions was used. A density view of the initial state, with the magnetic field lines superimposed, is given in Fig. 3.

Fig. 3. Density view of the initial condition of two interacting slow MHD waves perturbing the thermal equilibrium. Magnetic field lines are superimposed in blue.

3.3. Numerical methods

The set of equations, Eqs. (7)–(10), was solved on a 2D Cartesian grid using a strong stability preserving five-step fourth-order Runga-Kutta method (Spiteri & Ruuth 2002) and an HLL flux scheme, first described by Harten et al. (1983). To compute fluxes at the cell edges, cell centers values were reconstructed where we used a second-order symmetric total variation diminishing (TVD) limiter, denoted ‘woodward’ (van Leer 1977; Woodward & Colella 1984). We used a Courant number of 0.8 to limit the time step used by the explicit Runga-Kutta method according to the Courant–Friedrichs–Lewy condition (CFL; Courant et al. 1928). The divergence-free condition of the magnetic field, ∇ ⋅ B = 0, was numerically handled by the constrained transport method (Gardiner & Stone 2005). A lower limit to the size of the time step was set to ensure the feasibility of the simulations.

The square grid has a physical domain spanning from 0 to 10 Mm, as shown in Fig. 3. The base resolution is 100². These cells have a width of 100 km. With adaptive mesh refinement (AMR), the effective resolution of the smallest cells can be increased to 3200² for six levels of AMR, which corresponds to a physical width of 3.125 km. The number of levels of AMR we used, and hence the smallest resolvable length scale, depends on the simulations itself and is mentioned when appropriate. Our highest physical resolution was the same as in Claes et al. (2020a). A user-defined criterion based on the density was used to refine and coarsen the grid.

In Sect. 6 we compare the five cooling curves discussed in Sect. 2. A summary of all the cooling curves available in the MPI-AMRVAC framework and their properties is given in Appendix A. During a simulation, the lowest temperature of the adopted cooling curve is taken as the lowest possible temperature that the simulation can ever reach in order to avoid unphysical negative pressure values. An artificial Tfix parameter was set to ensure this. For the five cooling curves discussed above, this is not a problem because their lowest temperature is 1000 K at most and the lowest temperature in the simulations performed here does not reach to this value. However, in Sect. 5.2 we discuss the effect of the low-temperature treatment of the cooling curves. Radiative cooling was handled as an additional source term. It was computed last to ensure that no other physical or numerical process could modify the temperature below this limit during a time step. The radiative energy loss was solved to use the exact integration scheme developed by Townsend (2009). For piece-wise cooling curves, the analytic formulae were used. However, for interpolatable tables, the temporal evolution function, denoted by Y(T) in Townsend (2009), was calculated and tabulated in advance by numerically approximating the integral of the reciprocal of the cooling curve with respect to the temperature. This table was then used to calculate the energy loss by the optically thin radiative cooling cooling. More extensive details of the implementation for interpolatable tables are described in van Marle & Keppens (2011). The exact method is much faster than a traditional explicit scheme because it does not set a limit on the time step of the simulation. It is also more reliable than implicit schemes (van Marle & Keppens 2011). The constant heating term, which is needed to initialise the setup with a background in thermal equilibrium, is given by Eq. (13). It was handled as a user-set source term to the energy equation.

4. Benchmark simulation

The numerical setup is based on the work by Claes et al. (2020a). Before the effect of several parameters on the formation of condensations can be compared, a benchmark needs to be set. We used the same cooling table and effective resolution as were used for the highest-resolution run in Claes et al. (2020a), which are the SPEX_DM table and 3200². This resolution means that the smallest cells have a width of 3.125 km. We did not include thermal conduction, the effects of which are discussed later in this section. In addition to this difference, the largest difference is the orientation of the background magnetic field. In our case, the angle between the x-axis and the magnetic field is 45°, while for the highest-resolution run in Claes et al. (2020a), the angle was 54°. This small difference is significant, however, when in the benchmark simulation the magnetic field and combined perturbations of the two slow MHD modes travel in the same direction. This artificial symmetry ensures additional stability with respect to dynamical effects.

The evolution of the thermal instability is shown in the four panels in Fig. 4, where we plot the density view. The basic principles are the same as explained in detail in Claes et al. (2020a). The top left panel is 4.12 h after the initial state shown in Fig. 3. At this point, the interacting slow MHD waves are sufficiently damped due to the energy loss by radiative cooling. The thermal mode has started to grow because it is excited in an unstable regime. A region of increased density and decreased temperature starts to develop along the diagonal, that is, in the direction in which the interacting slow MHD waves were travelling. The lower temperature and hence thermal pressure causes an inflow of matter towards this spot due to a pressure gradient. The inflow occurs along the magnetic field lines because of flux freezing in the low plasma beta medium. After 4.12 h, this region has reached a density of twice the initial background value. The increased density leads to a larger contribution of radiative cooling, which is dependent on the density squared. The cooler plasma also facilitates a higher cooling rate because the cooling curve globally increases from 10⁶ to 10⁵ K. This causes the catastrophic cooling associated with the thermal instability. In a mere 0.04 h, or 2.4 min, the density increases to 20 times the initial background value, as we show in the top right panel, due to the counter-streaming inflows of matter along the magnetic field lines. This is the onset of the filament formation stage. Similar to Claes et al. (2020a), we wish to note that with the term ‘filament’, we mean a thin, high-density, low-temperature structure, not necessarily the solar feature related to prominences.

Fig. 4. Density view of the benchmark simulation using the SPEX_DM table and a resolution of 32002. The top left and right panels have a density of twice and twenty times the background value, respectively. The bottom left panel has a density of 95 times the background value, and the bottom right panel denotes the end of the simulation due to the small time-step restriction.

Another 1.2 min later, a long, thin filament has formed, and its density has increased to approximately 95 times the background value. At this point, the filament has reached a temperature below 10 000 K, at which the optically thin radiative cooling becomes less effective, see Fig. 1, ending the thermal instability. At the same time, a large part of the finite amount of mass in the strip perpendicular to the filament has been accumulated to the filament, limiting further growth in density. However, it should be noted that for the initial configuration typical for the solar corona, this density ratio matches the difference of two orders of magnitude that is observed between the solar corona and cool condensed prominences (Parenti 2014). The rapid increase in density is due to accretion of mass along the magnetic field lines from plasma that is sucked into the central region. In this process, rebound shocks are observed to form that move outwards against the infalling material, as discussed in Claes et al. (2020a). The depletion in the region perpendicular to the filament is apparent in the bottom left panel of Fig. 4. In the right panel of Fig. 5, the rebound shocks are shown in velocity view. As previously discussed in Claes et al. (2020a), these rebound shocks are formed by the collision of the counter-streaming inflow of matter along the field lines that created the high-density filament. Matter is fed to the filament from the region in between the filament and the shock. The uniform medium sound speed and the Alfvén speed are approximately 150 and 583 km s⁻¹, which causes the velocities to become subsonic in the simulated reference frame. Similar rebound shocks have been observed in numerical simulations of prominence formation (Xia et al. 2012) and coronal rain (Fang et al. 2015). Claes et al. (2020a) quantified the nature of these shocks as slow MHD shocks.

Fig. 5. Left: ram pressure view at the end of the benchmark simulation with insets zoomed in on the filament edges, where ram pressure differences fragment the filament. Right: total velocity magnitude after the filament has formed. It clearly shows the rebound shocks due to the formation of the filament.

From this point in our simulation, noticeable differences with the simulation of Claes et al. (2020a) start to arise. Owing to the lack of inclination between the magnetic field and the direction of the travelling waves, the filament has formed completely perpendicular to the magnetic field. The ram pressures at the front and backside of the filament are nearly equal, making it dynamically stable with respect to rotation. In the simulation by Claes et al. (2020a), the misalignment allowed for a ram pressure imbalance, which in turn caused a pinching effect. This caused the filament to stretch, effectively thinning it and making it more susceptible to small velocity and pressure imbalances. In our simulation there is no torque, and the filament remains straight.

The fourth panel of Fig. 4, after another half hour, shows the filament at its full extent with an inset showing fragmentation. A direct comparison of the timescale with the result of Claes et al. (2020a) is not possible because they used a different equilibrium temperature. At this point, the filament fragments. Unlike in Claes et al. (2020a), where the filament was first twisted and thinned, the filament does not fragment everywhere, but only at the edges where it is thin enough. Although the location of the fragmentation is not the same as in Claes et al. (2020a), the mechanism is. This is caused by small pressure or velocity imbalances, where small seeds may be numerical in origin (round-off errors). These imbalances grow due to the thin-shell instability. The best-known thin-shell instabilities are the linear (Vishniac 1983) and non-linear thin-shell instability (Vishniac 1994). The first occurs when a thin shell is caught between ram pressure and thermal pressure, and in the second case, the thin shell experiences ram pressure on both sides. In the left panel of Fig. 5 we show the ram pressure view with insets zoomed in on the filament edges of this final moment. The large ram pressure differences around the horizontal fragmenting strands of lower density is apparent. This type of fragmentation of thin shells and slabs has also been observed in simulations of star formation in molecular clouds (e.g. Hueckstaedt 2003) and in the interaction of stellar winds with circumstellar disks (e.g. van Marle et al. 2011; van Marle & Keppens 2012). The fragmentation in our simulation shows that the pinching effect described by Claes et al. (2020a) is not a necessary condition, but it helps by thinning the filament.

Within the fourth panel of Fig. 4, the inset shows a fragmenting edge with three perpendicular strands. This is the last frame of the evolution before our settings of the code, which warns the user when the CFL limit becomes too computationally demanding (time steps below 10⁻⁹ in normalised units) and causes the simulation to stop. This small time step is caused by the significant decrease in density in between the fragmenting strands, which in turn reduces the CFL condition. When no artificial measures are taken to handle near vacuum conditions, this is the farthest into the non-linear evolution we can reach. In Sect. 7 we explain how some bootstrap measures for low density can avoid this problem, and the simulation can continue.

As mentioned in Sect. 3.1, we only included optically thin radiative cooling and a constant heating term. This also means that we did not account for thermal conduction, in contrast to the simulations by Claes et al. (2020a). The effect of thermal conduction on a setup with optically thin radiative cooling, such as this thermal instability setup, is to facilitate a cut-off wavelength called the Field length. Perturbations with wavelengths shorter than this cut-off wavelength are stabilised, that is, they become damped because thermal conduction smooths these perturbations out. This was first described by Field (1965) and was recently studied in detail by Claes et al. (2020a). The Field length can be estimated as

$$\begin{array}{c}\hfill {\lambda }_F=\sqrt{\frac{\kappa (T)T}{n^2\mathrm{\Lambda }(T)}},\end{array}$$(18)

where n is the number density, T is the temperature, and Λ(T) and κ(T) are the optically thin radiative loss and conductivity at a temperature T, respectively (Field 1965; Begelman & McKee 1990). This equation is not normalised. For a plasma under solar conditions, we can use the Spitzer conductivity parallel to the magnetic field (Spitzer 2006), defined as

$$\begin{array}{c}\hfill {\kappa }_{\Vert }=8\times {10}^{-7}{T}^{5/2}\mathrm{erg}{\mathrm{cm}}^{-1}{\mathrm{s}}^{-1}{\mathrm{K}}^{-1}.\end{array}$$(19)

For the initial state and SPEX_DM cooling curve, as the benchmark simulation, the Field length is about 10¹⁰ cm and thus an order of magnitude larger than our perturbation wavelength. Hence the thermal mode will not be excited if we were to include thermal conduction in our setup. After the slow MHD modes are damped, the medium remains in balance. The idealised setup, without conduction, that we study here means that our Field length is zero and therefore leads to the scenario discussed in our benchmark case.

Thermal conduction also causes evaporation of cool, small parcels of plasma (Begelman & McKee 1990). Hence, when thermal conduction is not included or is not propely resolved, it can lead to numerical errors and non-convergence, as discussed in Koyama & Inutsuka (2004) and shown in Sharma et al. (2010). The critical length scale that needs to be resolved is again the Field length because it is the characteristic length scale of the thermal conduction and depends on the local physical parameters of the plasma, on the thermal conductivity, and the optically thin radiative cooling rate. In setups without conduction, such as we discuss here, the local Field length is always and everywhere effectively zero. Numerical instabilities can thus arise at the size of the grid, that is, the smallest size in the simulation, because they are not damped out by conduction. Evaporation of cold plasma is thus not possible in our setups. Evaporation might alter the non-linear evolution when the filament has fragmented into treads, but this is not included in our work. However, the same issues arise when thermal conduction is present but the local Field length is unresolved, as would be the case if we were to incorporate thermal conduction. Koyama & Inutsuka (2004) stated that at least three cells per local Field length are needed to obtain a physically correct solution.

A complementary view emphasises that this is because the transition region between the cold and dense condensation and the hot corona would otherwise not be correctly resolved, and hence thermal conduction would not be able to transport energy to the correct grid cells (Koyama & Inutsuka 2004). A higher resolution minimises this problem. Nevertheless, even for our highest-resolution simulation, the local Field length, corresponding to the condensations, would still be unresolved by several orders of magnitude. At the typical temperatures and densities of condensations in the solar corona, the local Field length is of the order of a few metres, while high-resolution multi-dimensional simulations can resolve up to a few kilometres. In the right panel of Fig. 6 we plot the Field length and energy loss by radiative cooling along a strip denoted by the red line in the left panel of the same figure at the end of the simulation. The selected strip is along the diagonal of the domain. It intersects the filament and thus contains information on the Field length in the filament and corona, as well as in the transition region between them. The Field length is several metres long in the filament and clearly unresolved. The Field length in the hot corona is well resolved. However, in the transition region, the Field length is even shorter than a few kilometres. This is due to the peak in radiative energy loss in the transition region. This peak is caused by the high cooling rates for temperatures between 10⁴ and 10⁶ K, in combination with a still high density with respect to the thin corona. The evaporation, driven by thermal conduction, in the transition region, that is, at intermediate temperatures, would not be resolved correctly. Even if we were able to resolve the local Field length everywhere in the domain along the magnetic field lines, the length scale of thermal conduction perpendicular to the field is still smaller because the perpendicular conductivity is 10¹² times lower (Braginskii 1965), and the associated resolution required would make a multi-dimensional MHD simulation computationally unfeasible. This limitation can be overcome if the conductivities are modified or isotropic conductivity is used, as discussed in Sharma et al. (2010). The TRAC method that was recently developed by Johnston & Bradshaw (2019) can patch this problem by artificially increasing the conductivity and decreasing the optically thin radiative cooling rate. This leads to a larger Field length and hence lifts the high-resolution restriction. The method was recently implemented in MPI-AMRVAC by Zhou et al. (2021), where it was extended to multi-dimensional MHD setups with anisotropic thermal conduction for the first time. However, this patch is typically used for the chromosphere-corona transition region and not for the transition region around in situ condensations.

Fig. 6. Field length and optically thin radiative energy loss along a selected line. Left: density view at the end of the benchmark simulation with the selected line, along the diagonal of the grid and through the filament and transition region, superimposed in red. Right: field length and optically thin radiative energy loss along the selected line.

5. Effect of numerical settings

In this section we investigate the effect of certain numerical, that is, unphysical, parameters or treatments of the evolution of the condensation formed by the thermal instability. First, we varied the resolution to understand differences between low- and high-resolution simulations. Second, we examined the low-temperature treatment of cooling curves. As mentioned in the preceding sections, there are several ways to handle the optically thin radiative cooling below 20 000 K. At this temperature, the plasma becomes opaque and the optically thin assumption breaks down. Some cooling curves stop at this temperature, leaving the simulation without information. We used a setup similar to the benchmark setup and highlight the differences.

5.1. Resolution

The benchmark simulation was run with six levels of AMR, leading to an effective resolution of 3200². The width of the smallest cells was 3.125 km. Density blobs had a width of several cells, corresponding to a few tens to 100 km. This is of the order of the observational limit of the finest instruments, such as CRisp Imaging Spectro-Polarimeter (CRISP; Scharmer et al. 2008). Scullion et al. (2014) showed statistically that the fine structure in post-flare loops might have a cross-sectional width smaller than 100 km. However, commonly used instruments, such as the Atmospheric Imaging Assembly on board the Solar Dynamics Observatory (SDO/AIA; Lemen et al. 2012) and the EUV Imaging Spectrometer (EIS) on board Hinode (Culhane et al. 2007; Kosugi et al. 2007), can only resolve down to around 1000 km (Lang et al. 2006; Grigis et al. 2012). It is worthwhile to investigate the effect of the resolutions on the simulations to determine how much of the small-scale dynamics of the condensation formation process will be revealed at a given resolution. This is particularly timely, with upcoming possibilities to reach several dozen kilometre resolution using the Daniel K. Inouye Solar Telescope (DKIST; Rimmele et al. 2020).

We used the SPEX_DM cooling curve and kept the base resolution constant to 100². The levels of AMR were varied from one to six. The MPI-AMRVAC framework has a block-based refinement scheme with a fixed refinement ratio between grid levels of two (Keppens et al. 2012). When refining in 2D, a single block becomes four blocks, which each have a physical dimension half of that of the original block. We denote the simulations by the effective resolution. We thus examine the cases 100², 200², 400², 800², 1600², and 3200². They correspond to a smallest width of 100, 50, 25, 12.5, 6.25, and 3.125 km, respectively.

In the panels of Fig. 7 we show the distribution of density after 4.22 h. This is during the fragmentation phase. We chose this moment for the comparison because at this time, the 800² run encountered the small time-step limit. For all resolutions, a filament is formed, and the filaments are of approximately the same length. However, the width is different. The filament is thinner for higher resolutions. This behaviour has also been observed in simulations of filamentary structure in galaxy clusters by Sharma et al. (2010). They discuss the need of resolving the Field length in every direction to obtain converged results. As mentioned before, the Field length here is effectively zero because we did not incorporate thermal conduction, and it is thus unresolved. The large gradients in quantities such as temperature or density necessarily induce small numerical errors. Strictly speaking, the 3200² simulation does not converge, but the results do not differ much from a resolution of 1600² to 3200². At the highest resolution, the filament has 20 cells across the diagonal to ensure a smooth transition region. This resolution hence truly resolves the internal filament structure.

Fig. 7. Density views of the simulations at different resolutions. The snapshots are all taken after 4.22 h.

From the first panel to the last, it is clear that the fragmentation occurs in each case, but demonstrates a different dominant wavelength along the filament. For a lower resolution, the wavelength with which the filament fragments is longer. In the case of the 100² run, there are only three wavelengths. The 800² run, in contrast, has approximately 12 wavelengths. This is because the wavelength of the thin-shell instabilities depend on the thickness of the shell. The wavelength is shorter for thinner shells. This has also been observed in simulations of stellar winds (Schure et al. 2009). In the 3200² run, no fragmentation has occurred. However, from the benchmark, we know that it will start to break up at a later time and only at the edges, as shown in Fig. 4.

For each of the runs, we tracked the evolution of the maximum density and minimum temperature within the entire grid. These values are important properties of the filament. The maximum density reflects the amount of mass that was able to accumulate before the fragmentation begins. The total filament can still accumulate mass during the fragmentation phase. The minimum temperature denotes the temporal evolution of the thermal instability, which causes the formation of the condensation in our simulations.

In the left panel of Fig. 8 we show the evolution of the maximum density. It only shows the later phases because the density changes slowly while the slow waves are being damped. Several things can be observed. First, the simulations with lowest resolution, 100² and 200², do not reach the same maximum density as the other simulations at the end of the non-linear cooling phase. This can be explained by the fact that the filaments of these simulations fragment earlier and in larger strands. Second, the maximum density continues to increase for the highest resolution. This is because it has not yet started to fragment. Lastly, even though the filaments of the simulations with lower resolution all fragment, their evolution in terms of the instantaneous maximum density is still slightly different. This is due to the erratic behaviour of the fragmentation process, which is seeded from the numerical round-off errors, as mentioned earlier.

Fig. 8. Evolution of the maximum density and minimum temperature for the simulations at different resolutions. Left: maximum density. Right: minimum temperature.

The minimum temperature within the grid is shown in the right panel of Fig. 8. The typical evolution of the thermal instability is very clear, where one can compare this panel to the left panel of Fig. 2. For the four simulations with the lowest resolution, the minimum temperature reaches a plateau. During the fragmentation, the density has dropped or stopped increasing. This causes the cooling by optically thin radiation to be less efficient because it depends on the density squared. This does not occur in the simulations at the highest two resolutions. For these, the minimum temperature continues to decrease, although more slowly than during the non-linear cooling phase. This can be explained by the increase in density and the non-zero cooling rate, despite the decrease in cooling rate for temperatures below 10 000 K for the SPEX_DM table.

The density of the filament increases due to the accretion of plasma along the magnetic field lines, as discussed previously. Two other interesting characteristics of the filament are its surface area and accumulated mass. The simulations have only two dimensions. The mass is the surface mass, defined as the density of the filament times its surface area. The grid has a non-uniform resolution due to the adaptive mesh refinement; the region around the high-density filament has a higher resolution due to the refinement criterion. The surface area and surface mass were determined after extrapolating the grid to the highest resolution within the simulation. We defined the filament as the cells with a temperature lower than or equal to 15 000 K. There are two reasons for this temperature limit. First of all, this temperature is near the end temperature reached by condensations formed by the non-linear phase of the thermal instability. Second, observations have shown that prominences and coronal rain formed in the solar corona reach these kinds of temperatures (see e.g. Parenti 2014; Antolin 2020).

The left panel of Fig. 9 clearly show that the surface area of the filament is smaller for simulations at higher resolutions. The reason is that at higher resolutions, thinner structures are formed. This behaviour was also present in the density views of Fig. 7. In the right panel of Fig. 9 we show the evolution of the surface mass. A lower resolution leads to more mass in the filament. This is likely related to the larger surface area that is covered by their filament. For these panels the exact values do not really matter because the final mass and surface are not reached due to the small time-step limit. For completeness, the total surface area and mass in the grid are 100 Mm² and 2340 g cm⁻¹. In all cases only a small fraction of the total surface area and mass are accumulated in the filament, while the conservative discretisation employed and the use of double periodic boundary conditions ensures that the total mass in the domain is perfectly conserved.

Fig. 9. Evolution of the surface area and surface mass of the filament for the simulations at different resolutions. Left: surface area. Right: surface mass.

5.2. Low-temperature treatment

At low temperatures and high densities, such as for the condensations formed in these types of simulations, the plasma has a high opacity. It cannot cool efficiently and is optically thick. The plasma is not fully ionised for temperatures below 20 000 K. Hence, self-absorption and photoionisation can occur within an optically thick plasma. This violates the assumptions of the collisional ionisation equilibrium and hence the optically thin radiative cooling (Dopita & Sutherland 2003).

To correctly determine the energy loss by radiation, we would need to solve the full radiative transfer equations. However, this quickly becomes prohibitive in combination with multidimensional MHD, creating a need for more pragmatic measures that approximate the actual physics and keep simulations running.

As an example, we considered the table of Colgan et al. (2008). In the first low-temperature treatment, no changes to the table were made, but the Tfix parameter was used. This has the effect that when the plasma reaches the lowest temperature of the table, 11 482 K for the Colgan table, instead of losing more energy, the plasma is kept at that temperature. Energy is artificially increased in the simulation. This is in contrast to the simulations without such a fix, where energy conservation is always in accord with the governing equation.

The second method is commonly used in the astrophysics community. The cooling curve is modified to vanish for low, approximately 20 000 K, temperature. The modified version of the Colgan table is the JCcorona table (Xia et al. 2011), and it is shown in the top left panel of Fig. 10 together with the cooling curves of the other two treatments. This figure shows that by vanishing at 20 000 K, a clear bump in the cooling rate is neglected. This bump is due to the Lyman-α radiation (Cox & Tucker 1969). Because the plasma is not fully ionised, it has been suggested in the past that neglecting the Lyman-α radiation is better than the optically thin radiative approximation at these temperatures (McClymont & Canfield 1983). Nevertheless, corrections for the radiative loss rate for Lyman-α have been attempted (Athay 1986).

Fig. 10. Three cooling curves with different low-temperature treatment and their density view at the end of the evolution. Top left: cooling curves. Top right to bottom left and bottom right: density views at the end of the evolution. They correspond to the Colgan_DM, JCcorona, and Colgan cooling tables, respectively.

The third low-temperature treatment is also a modification, or rather an extension. It is the Colgan_DM table, as discussed in Sect. 2. It consists of a high-temperature part of the Colgan table Colgan et al. (2008), augmented with cooling rates for temperatures below 10 000 K by Dalgarno & McCray (1972). A low-temperature extension is justified because even though the optically thin approximation fails, there will still be some radiative cooling. The cooling is due to collisional excitation of singly charged ions of O, C, N, Si, Fe, Ne, and S with thermal electrons and neutral hydrogen. There is also a contribution by collisions between neutral hydrogen with electrons, as described in Schure et al. (2009). The cooling rate at these low temperatures is several orders of magnitudes lower than at high temperatures, as shown in the top left panel of Fig. 10, but is never actually zero (as adopted without this extension of the cooling curve).

To show the effect of the low-temperature treatment, we ran the thermal instability setup for the three cooling curves discussed in the previous paragraphs. We used four levels of AMR, leading to an effective resolution of 800² and a width of the smallest cell of 12.5 km. From the previous subsection we have learned that this will yield fragmentation at larger wavelengths than for higher resolutions. We would like to emphasise that we studied a single cooling curve here, the Colgan table, and applied different treatments for the low-temperature regime. In Sect. 2 we discuss the comparison between different cooling curves that might have the same low-temperature treatment, such as the Colgan_DM and SPEX_DM tables.

The density views at the end of the simulations, the moment at which the low density limits our time step, are shown in the right and bottom panels of Fig. 10. The top right, bottom left, and bottom right panels correspond to the Colgan_DM, JCcorona, and Colgan cooling tables, respectively. In all cases a filament is formed, but there are differences in morphology. For the Colgan and Colgan_DM the filaments fragment, but for the JCcorona table, it becomes thicker and only the edges stretch. The difference with and without the _DM part is not so pronounced. In the case of the Colgan table, the fragmentation occurs in a similar way as for the Colgan_DM, but the filament remains unperturbed in the centre. We are led to conclude that the vanishing and the ceasing of the optically thin cooling undermines or at least delays the thin-shell instability. Another way of viewing this is that the radiative cooling becomes ineffective for the JCcorona table. The evolution of the filament is therefore governed by its dynamics, which appear to be near equilibrium.

As before, we tracked the evolution of the maximum density and minimum temperature in the grid for the three simulations. The results are presented in Fig. 11. We note a small difference in timescale of the thermal instability, although the cooling rate at 10⁶ K is exactly the same. Due to interpolation when constructing the interpolatable cooling curves numerically, very small differences in the cooling curves arise. This can lead to slightly more or less energy loss after many iterations. Hence the timescale of the thermal instability might be slightly longer or shorter, respectively. In all the cases the instability begins and a filament is created, as shown in the density views of Fig. 10.

Fig. 11. Evolution of the maximum density and minimum temperature for the simulations with the Colgan_DM, JCcorona, and Colgan cooling tables, which only differ in their low-temperature treatment. Left: maximum density. Right: minimum temperature.

In the right panel of Fig. 11 the different low temperature treatments are clear. The temperature remains constant for the Colgan table after cessation of the thermal instability. For the simulation with the JCcorona table, the thermal instability stops abruptly at 20 000 K, and its further temperature evolution is determined by the changes in pressure and density. For the Colgan_DM table, the evolution reaches lower temperatures and stops because of the fragmentation and small time-step limit.

The evolution of the maximum density is also different, as expected from the density views. This is shown in the left panel of Fig. 11. The results of the Colgan and Colgan_DM table are quite similar, but the simulation with the JCcorona table reaches a much lower maximum value. It also remains at this value afterwards because there is no apparent fragmentation to alter the density significantly. The density views also showed that the extent of the filament perpendicular to the magnetic field lines is different. This might be because the onset of fragmentation, which varies, appears to limit the perpendicular growth.

We can conclude that the low-temperature treatment has a strong effect on the density and temperature evolution of condensations. Indirectly, it has an effect on the fragmentation process, and hence on the morphology of condensations. In reality, this means that the true physical conditions of being non-optically thin when the condensation has formed should be accounted for in future simulations.

6. Effect of changing the cooling curve

In the previous section, the effect of numerical settings, such as the resolution and low-temperature treatment of the radiative cooling, were investigated. Here we discuss the effect on MHD simulations of completely changing the optically thin radiative cooling curve. To this end, we again study the effect on the condensation process of the thermal instability. The general setup is discussed in Sect. 4, with six levels of AMR corresponding to 3.125 km. We used the five cooling curves discussed in Sect. 2. The Colgan_DM, Dere_corona_DM, and SPEX_DM are interpolated tables, which are extended for low temperatures with the DM table. The Klimchuk and Rosner tables are piece-wise power laws.

The evolution of the thermal instability and formation of a condensation for the five cooling curves is shown in the panels of Fig. 12. Each row denotes a different phase. The first four rows are the same four times as discussed for the benchmark simulation of Sect. 4. The benchmark is the SPEX_DM run. The last row shows a zoom-in of the filament and has an inset that is zoomed into the low-density region that triggers the small time-step condition.

Fig. 12. Density views during the simulations with different cooling curves. Each cooling curve has its own column. The rows represent similar moments and are scaled to the same density colour bar. The first three rows are when the filament has a maximum density of approximately 2, 20, and 95 times the background density. The fourth row is the last moment, at which the time step becomes restrictive. The last row is a zoom-in onto the filament of the last moment, it also contains a zoomed-in inset of the low-density region.

For all five cooling curves, the first stages of the evolution are similar. The thermal instability is not obstructed by changing the cooling curve. A thin filament with a maximum density of 95 times the background value always forms. From the fourth and fifth row, it is clear that the final shape, that is, the morphology, can be different. The filaments of the simulations with the Colgan_DM, Dere_corona_DM, and SPEX_DM tables remain straight, with fragmentation only at the edges. They also have approximately the same length. The filament belonging to the Klimchuk table is much shorter and fragments at several places, but most prominently at the centre. The filament of the Rosner table is longer and becomes wavy. There is again a suggested link between the length of the filament and the amount or onset of fragmentation. The perpendicular growth appears to be hindered as soon as the filament fragments.

These global differences in morphology are due to the low-temperature treatment of the cooling curves, because for all cooling curves, the thermal instability creates a condensation. It is thus not surprising that the results for the first three cooling curves are similar. They have the same low-temperature treatment. There are still small differences between these three cooling curves because their earlier evolution was slightly different. This difference leads to small variations in density and pressure and hence in the fragmentation by the thin-shell instability.

Based on the results of Sect. 5.2, we can try to explain the behaviour of the other two curves. We have noted before that when the cooling curve vanishes, the filament did not fragment as easily. This is also the case for the Rosner table. Its cooling rate does not vanish completely at 20 000 K, but drops off very fast and does not reach a low-temperature plateau such as for the tables with _DM extension. The Klimchuk table is much coarser. Its cooling rate for low temperatures remains at least an order of magnitude higher than for the tables with the _DM extension, see Fig. 1. A possible explanation could be that a higher cooling rate at low temperatures leads to small temperature and hence pressure variations or blobs within the filament. This in turn leads to easier fragmentation by thin-shell instabilities due to the presence of small seeds.

We also wish to note that by changing the cooling curve, we inevitably also modify the heating of a given model. To recall the heating function, given by Eq. (13), it is constant, depends on the cooling curve, and was chosen this way to ensure thermal equilibrium of the initial state. Thus the differences observed here between models with different cooling curves are due to the combined effect of a different cooling rate and background heating.

The right panel of Fig. 13 shows the evolution of the minimum temperature. As previously discussed, the linear timescales are very different on the basis of the cooling curve itself (see Fig. 2). The numerical simulations indeed confirm this. The Rosner cooling curve has a much lower cooling rate at 10⁶ K. The thermal instability takes much longer to form a condensation. In addition to the difference in timescales, it is important to note that the final temperature is different. However, all temperatures reach below 20 000 K, at which optically thin cooling is no longer strictly valid. The maximum density reached is also different, as shown in the left panel of Fig. 13. These two effects are intertwined. The simulations stop because of the fragmentation, which is different even for cooling curves with the same low-temperature behaviour. Small differences in the cooling curves during the evolution of the filament can cause structural differences and hence small differences in the final result of certain parameters, although the global density view looks similar.

Fig. 13. Evolution of the maximum density and minimum temperature for the simulations with different cooling curves. Left: maximum density. The horizontal dotted grey lines denote the maximum densities of the rows of Fig. 12. Right: minimum temperature.

The Klimchuk and Rosner cooling curves give a totally different result. The minimum temperature of the Klimchuk simulation continues to decrease and causes its maximum density to increase. This is most likely because its cooling rate does not decrease much compared to the other cooling curves, and hence energy is still efficiently radiated away. This allows the temperature and pressure to drop and facilitates the accumulation of more matter. The simulation with the Rosner cooling curve reaches a more steady end phase, with both the maximum density and temperature quite stable for more than an hour. This behaviour can also be assigned to the low-temperature behaviour of its cooling curve, as discussed before.

It seems trivial that the timescales are different, but it has important consequences. The setup of the thermal instability we used is simplified. In the real world, many forces and effects occur at once, and many external perturbative events can take place. In the context of solar physics and prominence formation, a flare or wave might perturb the environment. If this occurs after a given time, then the condensation might have formed for one cooling curve, but not for another.

For the simulations with varying resolution, we track the evolution of the surface area and surface mass of the filament, defined by the cells with a temperature lower than 15 000 K. The results are shown in Fig. 14. The surface masses of the filaments of the three _DM tables are nearly equal, but their surface areas are not. Just like for their maximum density, small variations in the physical variables and small differences in the fragmentation lead to differences in surface area. The SPEX_DM table has the highest maximum density and the smallest surface area of the three curves. The opposite is true for the Colgan_DM table. This is because the slightly earlier fragmentation in the case of the SPEX_DM table leads to more contracted or compressed density clumps in the strands. The combination of the two effects leads to the equal surface mass.

Fig. 14. Evolution of the surface area and surface mass of the filament for the simulations with different cooling curves. Left: surface area. Right: surface mass.

For the Rosner table the filament did not fragment. It continued to grow in both directions, along and perpendicular to the magnetic field. This is clear from the evolution of the surface area and surface mass, and from the density views. The Klimchuk table produced distinctly different results. The surface area is at least four times smaller than for the other tables and twice the surface mass. Due to the strong fragmentation in this case, the growth is halted. Its perpendicular length is also much shorter than for the other curves, as shown in the bottom panels of Fig. 12. The smaller filament has thus a smaller surface area, and this leads to a lower surface mass. The difference in surface mass is smaller because the overall density of the filament is higher, counteracting the small surface area.

7. Extending far into the non-linear regime

We can now use the insights gained in the previous sections to study the evolution of a condensation formed by the thermal instability far into it non-linear regime. To do this, we study the fragmentation phase in more detail.

All the simulations we showed so far and the similar simulations performed by Claes et al. (2020a) were terminated due to a small time-step restriction. This small time step is caused by the Courant–Friedrichs–Lewy condition (CFL). This is a necessary condition for stability when any kind of explicit time-marching scheme is used for the essentially hyperbolic MHD equations that were first derived by Courant et al. (1928). For a given resolution normalised Δx, it is given by

$$\begin{array}{c}\hfill C=\frac{|v|\mathrm{\Delta }t}{\mathrm{\Delta }x},\end{array}$$(20)

where C is the Courant number, Δt is the normalised time step, and v is the fastest velocity at which information travels through the grid. In this simulation, it is the local velocity magnitude plus the local fast magnetosonic speed. From a physical point of view, this means that the explicit time step has to be shorter than the time required for the fastest wave to propagate from one cell to another. More information can be found in standard textbooks on magnetohydrodynamics and numerical modelling (see e.g. Goedbloed et al. 2019).

In the simulations presented here, the short time step was caused by the sudden drop in density in between the fragmented strands perpendicular to the filament and aligned with the magnetic field. For our simulations, the time step can be approximated in terms of the density and resolution as

$$\begin{array}{c}\hfill \mathrm{\Delta }t\approx 0.1\mathrm{\Delta }x\sqrt{\rho },\end{array}$$(21)

where we used a Courant number of 0.8 and used the fast MHD wave velocity. The actual restriction on the time step thus comes from the fact that the filament collects almost all matter on an initial field line, and our periodic boundaries then mean that no additional mass is available to counter the drop in density.

In order to overcome this problem, we made use of a numerical bootstrap technique. This amounts to introducing an artificial lower boundary to the allowed density, that is, if the density drops below a certain value, we adjust it. This is one of the different uses of the low value fixes, such as negative pressure handling, within MPI-AMRVAC. The framework has several methods to handle low values, such as replacing the triggered cell by a user-set lower bound or by averaging in a box around the cell. The average method is normally more physical, but in our case, the low density occurs in between the fragmented strands. This causes the averaged value to be very high compared to the original low density. We therefore made use of the replace method. This method violates the conservation of mass in the system by introducing additional plasma.

Because we aimed to extend a simulation for a long time, we needed a reasonable time step for a given iteration while still resolving the physical details. We chose to use four levels of AMR, that is, an effective resolution of 800². We used the SPEX_DM table as for the benchmark case. Our lower density bound, also denoted by ρ_min, was set to 0.01 in normalised units, which is about 10⁻¹⁷ g cm⁻³. This was sufficient to simulate the entire further non-linear evolution we discuss next.

In Fig. 15 we show the evolution of thermal instability and its fragmentation. In the top left panel we show the starting point. We started from the 800² simulation of Sect. 5.1 at the point where the condensation started to form by the thermal instability. The lower density bound is not yet reached here. In the top right panel we show the fragmenting filament at the moment in which the low-density treatment starts working. The white regions in between the strands have this low density value. The bottom left image depicts the density distribution roughly 22 min later. The highest density has decreased, but most of the mass in the filament has coagulated in the centre of the filament, while the outer parts become strands guided by the magnetic field. At this stage, 1% of the total number of cells have been affected by the treatment. These are the white areas between the strands. In the bottom right panel we show the density distribution after one hour and twenty minutes. We ended the simulation at this point because about 5% of additional mass was introduced to the simulation by the continuous bootstrap strategy. At this point, the filament is completely broken up and the strands are scattered throughout the grid. However, their orientation is still aligned with the magnetic field, which is to be expected because it is a low beta plasma with a strong magnetic field. There are extended regions with low density.

Fig. 15. Density view of the full evolution of the thermal instability. Top row: onset of the thermal instability and the moment in which the low-density bootstrap measure kicks in, in the left and right panel, respectively. Bottom left: moment in which 1% of the area is at the prescribed lowest allowed density. Bottom right: last moment of the simulation. The animation of this simulation is available online.

In Fig. 16 we show the evolution of the maximum density and minimum temperature in the grid. As mentioned previously, the starting point is the onset of the thermal instability, which can be seen as the rise and drop in density and temperature, respectively. In these figures we can study the full evolution in contrast to the results in the previous sections.

Fig. 16. Evolution of the maximum density and minimum temperature for the simulation of the full evolution of the thermal instability. Left: maximum density. Right: minimum temperature. The evolution is shown from the start of the simulation, that is, the onset of the thermal instability. The dotted grey and red lines correspond to the top right and bottom left density views of Fig. 15, respectively.

The temperature evolution is straightforward to interpret. Because of the thermal instability, the temperature drops non-linearly. The drop stops at around 10 000 K. This is to be expected because the optically thin radiative cooling becomes less efficient. During the fragmentation and redistribution of matter, the minimum temperature remains approximately constant. However, there appears to be a slight increase in temperature while the filament is redistributed. It is noted here, and further explained in our Appendix B, that the fragmentation and dynamics seen in this low plasma beta situation are very different from purely hydrodynamic settings.

The maximum density increases to its maximum value due to the pressure gradient. The gradient facilitates the flow of material towards the filament along the magnetic field lines. The accretion of matter is clearly shown in the top right panel of Fig. 15 by the lower-density regions perpendicular to the filament. Through collisions within the filament and fragmentation, the maximum density drops slightly. The grey line corresponds to the time in the top right panel of Fig. 15, while the red line corresponds to the bottom left panel. At the latter time, the filament breaks up completely, which causes larger drops in maximum density. During the redistribution, several strands of plasma might interact with one another, causing small bumps in the maximum density, as shown in Fig. 16 at about 18 000 to 18 500 s. However, the global trend is a decrease in maximum density, meaning a redistribution of the material over the box along the magnetic field lines. This behaviour is quite reminiscent of coronal rain or downflows in the legs of prominences. Moreover, despite the initial orientation of the filament, which was perfectly orthogonal to the magnetic field, the final state we obtain is one of many filaments that are clearly aligned with the background field.

To study this complete fragmentation phase, we used a numerical trick to prevent the time step from becoming impractically short. Nevertheless, we need to investigate the effect of this bootstrap approach on the physical behaviour of the numerical simulation. In Fig. 17 we show the evolution of the percentual surface mass in the filament, that is, temperature lower that 15 000 K, the total surface mass, the number of cells or amount of surface area with the artificial limit value, ρ_min, and the surface area in the filament. The dotted grey line indicates the introduction of ρ_min and corresponds to the top right panel of Fig. 15. The dotted red line indicates the moment in which 1% of the number of cells or the surface area has a density equal to ρ_min and corresponds to the bottom left panel of Fig. 15. In the top panel of Fig. 17 we show that after the condensation is formed, its growth and increase in mass is quite steady. At most 77% of the total mass is condensed into a cold dense filamentary structure. The fragmentation and introduction of a lower density bound does not substantially change the initial growth rate. After 1%of the cells are modified, the growth stagnates. However, this almost coincides with the peak of fragmentation. The second panel of Fig. 17 shows that the additional mass induced by this measure is negligible until the 1% mark. Afterwards, the additional mass increases to 5%, at which point the simulation was terminated. The total number of cells or area with ρ_min also increases rapidly after the 1% mark, reaching a maximum of 34%. The behaviour and effect of the measure after the 1% mark can be better understood by considering the bottom panels of Fig. 15. In the moment in which the filament starts to break up into large strands, the growth in mass is stopped. From the pockets of the strands emerge long, white, low-density, regions. From these regions matter has been transported to the filament until ρ_min was reached. After breaking up, these large regions inject much additional mass, causing a continued increase in the total mass in the grid. Through the redistribution and interaction of the strands, the number of cells with ρ_min stops increasing. The surface area of the filament and strands is also almost constant during the redistribution and has a maximum of only 24% of the grid. Hence, while the bootstrapping indeed allows us to continue far into the non-linear regime of the fragmentation of our filament, the effect of its continued mass addition is certainly noticeable.

Fig. 17. Top panel: evolution of the surface mass of the filament. Second panel: evolution of the total surface mass in the grid. Third panel: evolution of the percentage of the number of cells or area with the prescribed lowest allowed density. Bottom panel: evolution of the surface area of the filament with respect to the total area.

8. Conclusions and discussion

Optically thin radiative cooling is one of the important basic processes in physics. The precomputed, tabulated prescriptions or cooling curves vary widely due to differences in computations, atomic physics parameters, and assumed solar abundances (Rosner et al. 1978). In the literature they are mostly used as unquestioned truthful ingredients, but their effect on the simulations is not considered.

We investigated how the condensation formation process is affected by the choice of cooling curve, and hence by the physical details of the curve. We used thermal instability as the mechanism in a simplified setup. For all tested cooling curves, the non-linear phase of the thermal instability led to the formation of a condensation that became a filament due to ram pressure. However, the growth rate of thermal instability is different for the various cooling curves. The growth rate obviously depends on the cooling rate at the initial temperature of 10⁶ K, which is dependent on the abundances. Cooling curves based on older values for the abundances tend to underestimate the cooling rate because of the lower values for the abundances of low first-ionisation potential elements (Feldman 1992). The second difference is in the morphology of the filaments that form. In some cases, the filament is much smaller in the direction perpendicular to the magnetic field. This is because the onset of fragmentation, which varies with respect to its growth rate for the thermal instability, is a determining factor. The shapes during fragmentation are also different. This might be because the seeds for the thin-shell instabilities, that is, small variations in the physical parameters, are slightly different as a result of the earlier linear evolution. We also found that some cooling curves are more stable against these instabilities, and that this is affected by the low (lower than 20 000 K) temperature treatment. This is especially the case when the cooling curve is artificially caused to vanish at low temperatures. The fact that the choice of cooling curve and low-temperature treatment affects the dynamical stability might also be important in more complicated simulations of prominences and coronal rain, where it might be easier to form a stable prominence, but harder to form transient coronal rain using a certain cooling curve.

Furthermore, we have shown the need for high-resolution simulations by investigating the effect of varying the resolution. In our highest-resolution simulations, the dynamic behaviour of the condensation during formation and fragmentation occurred on length scales much finer than currently observable. Upcoming instruments such as DKIST (Rimmele et al. 2020), which are able to resolve down to several dozen kilometres, are a great step forward.

Unfortunately, we noted that even at our highest resolution, we still encountered extremely low-density values in between strands of the fragmenting filament, similar to what has been observed in simulations by Claes et al. (2020a). These cells become void because of the accretion onto the filament. The region heats up, but keeps the pressure steady. It can therefore not regain mass by dynamical means, that is, pressure gradients. Thermal conduction was not included in the simulations and might help solve this problem. In the cases of Claes et al. (2020a), thermal conduction is present. However, the typical length scale for thermal conduction, the Field length, was not resolved. The low-density values lead to unfeasible short time steps due to the CFL condition. We provided all details of a possible bootstrap procedure that allowed us to simulate far into the non-linear regime. It involves the sudden introduction of a lowest allowed density value. The filament breaks up and is redistributed as long strands throughout the domain. This phase is likely the most relevant for the ultimate shapes seen in (filamentary structured) prominences. It would be exciting to find observational evidence for the change in the filament orientation we predict from our simulations. The same effects have been discussed in 3D settings in Claes et al. (2020a), but continuing into the far non-linear regime will require a bootstrap approach as advocated here. The additional effects of anisotropic thermal conduction, and the use of a numerically cured TRAC (transition region thermal conduction), will also play a role in future 3D simulations.

Movies

Movie 1 associated with Fig. 15 (dens_extending) Access here

Movie 2 associated with Fig. B.2 (dens_HD) Access here

Acknowledgments

We wish to thank the anonymous referee, for the constructive comments that improved the paper, and the editorial office. JH would like to thank Niels Claes for providing the initial setup and Jack Jenkins for the useful discussions. The visualisations were achieved using the open source software ParaView (paraview.org), Python (python.org) and yt (yt-project.org). The 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. Both authors 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.

References

Appendix A: MPI-AMRVAC cooling curves

We present the optically thin radiative cooling curves implemented in MPI-AMRVAC. Table A.1 lists basic properties of the cooling curves and their references. There are four piece-wise power laws: ‘Hildner’, ‘FM’, ‘Klimchuk’, and ‘Rosner’. All other cooling curves are interpolatable tables. The tables are extended to include a cooling prescription at temperatures beyond the upper temperature limits. The piece-wise power laws are extended by using their last segment. For the interpolatable tables, bremsstrahlung is used. All tables with the ‘_DM’ part are modified with the same low-temperature treatment as prescribed in Schure et al. (2009), where an ionisation fraction of 10⁻³ is used. The parameter Λ_max is defined for temperatures lower than 10⁷ K because of the ever-increasing cooling rate due to bremsstrahlung at high temperatures. The radiative and non-linear radiative timescales τ_r and τ_nl are calculated using Eq. (5) and Eq. (6), respectively. A typical solar density of 2.34 10⁻¹⁵ g cm⁻³ is used in these calculations. Not all tables always assume collisional ionisation equilibrium. The ‘ML’ and ‘cloudy’ cooling curves do not. We refer to the accompanying papers for more details about their computations.

Appendix B: Hydrodynamics (HD) case

Some differences between thermal instability in hydrodynamics and magnetohydrodynamics are highlighted in this section. Fig. B.1 shows the isochoric growth rate, that is, the imaginary part of the complex eigenfrequency ω in our adopted exp−iωt temporal variation, of the thermal, slow, and fast MHD modes in their dependence on the background magnetic field strength. The plot was obtained by solving the isochoric dispersion relation for a uniform infinite medium, as can be found in Claes & Keppens (2019). The medium corresponds to the setup used in this manuscript with the background values given in Table 1. The wave number and the angle between the background magnetic field and the wave vector were kept fixed to 2π divided by the domain length L and π/4, respectively.

Fig. B.1. Isochoric linear growth rates as a function of the background magnetic field B0. The thermal, slow, and fast mode growth rates are denoted by the green, red, and blue lines, respectively.

We used the SPEX_DM cooling curve and the background heating H, given by Eq. (13). We recall that ℒ is density dependent and given by Eq. (14). The background magnetic field was varied from 0 to 30 G. At this angle and physical medium, the growth rate of the thermal mode remains constant for all background magnetic field strengths, as shown by the green line. This is not the case for the slow and fast MHD modes, which are shown by the red and blue lines, respectively. The slow MHD mode becomes less damped as the background magnetic field decreases. The mode completely vanishes in the hydrodynamical limit, B₀ = 0 G, as shown in the real part of the wave number in Eq. (15). The damping rate of the fast MHD mode increases as the background magnetic field decreases. This mode has a sound wave character in the hydrodynamic limit.

To qualitatively show the difference between the result of a hydrodynamics and a low plasma beta magnetohydrodynamics thermal instability case, we modified the setup described in Section 3. The values of the background physical parameters were kept the same, but we omitted the magnetic field completely. We used the SPEX_DM cooling curve. We also initiated the setup with two sound waves instead of slow MHD waves because there is no hydrodynamic analogue of the slow MHD wave. Fig. B.1 shows that the role of slow and fast MHD waves changes drastically when the plasma beta changes to values below 1 (typical for our solar corona application).

In Fig. B.2 the density views are shown. The top left panel depicts the initial density state, which is not very different from the MHD case. The top right panel shows the formation of the condensation after the non-linear growth of the density due to the thermal instability. The timescale at which this occurs is slightly different from the one of the MHD case because the sound waves are damped with a different damping rate than the slow MHD waves. The condensation that forms is nicely spherical symmetric because of the uniform inflow of matter, that is, there is no preferential direction because there is no magnetic field. There is only the fact that the combined sound wave pattern travels along the diagonal. Afterwards, the condensation fragments in all directions, as shown in the bottom panels. In the HD case the fragmentation is not the same as in the MHD case, and it is not caused by ram pressure differences. There is no consensus about the mechanism for fragmentation of high beta hydrodynamic gas clouds (see e.g. Waters & Proga 2019a,b; Gronke & Oh 2020; Jennings & Li 2021; Das et al. 2021). One of those processes is ‘shattering’ (McCourt et al. 2018), where a cloud fragments into small pieces that each cool isobarically, independently of each other. Because our result resembles the hydrodynamic simulation of Jennings & Li (2021), the fragmentation process is most likely similar, although the initial conditions are very different. They indicate two processes for which the condensation is shredded apart through high vorticity, and they are thus explained differently from the ‘shattering’ mechanism of McCourt et al. (2018). In the first process, vorticity is generated in a rapid expansion phase by a Richtmyer-Meshkov (Richtmyer 1960; Meshkov 1972) process before the condensation fragments into smaller pieces (Gronke & Oh 2020). In the second process, smaller condensations are fragmented by pressure gradients between merging large clouds (Waters & Proga 2019b). We also find an expansion phase before fragmentation in our simulation, as shown in Fig. B.2. The high vorticity is also apparent in Fig. B.3 at the same moment as in the bottom left panel of Fig. B.2. We did not include thermal conduction, and thus no evaporation, in our simulation.

Fig. B.2. Density view during the hydrodynamical case. Top left and right: initial state and onset of condensation formation, respectively. Bottom left and right: rapid expansion and fragmentation of the condensation, respectively. The animation of this simulation is available online.

Fig. B.3. Vorticity magnitude at the same moment as in the bottom left panel of Fig. B.2, zoomed onto the condensation.

Fig. B.4. Curve representing the solutions to the thermal equilibrium, analogously to the S-curves used by the AGN community. The black line represents all the solutions to the equilibrium ℒ = 0 for the SPEX_DM cooling curve and corresponding heating term. The red dot indicates the initial conditions.

Using equations 6 and 8 from Waters & Proga (2019a), we quantified the mode behaviour of our HD case, shown in Fig. B.2. Using their notation, we have R = 2.689, N_ρ = −0.575, and N_P = −2.575. The setup is in case 6, where the plasma is both isobaric and isochoric unstable. In this regime, we indeed find two damped and propagating sound waves, and the unstable entropy mode responsible for thermal instability. The transition from strong background magnetic field and low plasma beta regime to the hydrodynamic limit and high plasma beta regime is worth a more detailed investigation.

As an analogue to the S-curves used in the AGN and circumgalactic medium (CGM) community (Lepp et al. 1985), we can plot the information of the cooling and heating functions differently. We solved ℒ = 0 with ℒ given by Eq. (14) for different densities. Plotting the temperature against the reciprocal of the density yields Fig. B.4. This figure clearly shows that there is far more detailed variation than the S-curve used for example by Waters & Proga (2019a).

All Tables

All Figures

	Fig. 1. Five cooling curves that we compare.
In the text

	Fig. 2. Left: general thermal evolution of thermal instability and its timescales. Right: linear and non-linear timescales of five cooling curves. The linear timescale is depicted by the length of the full bar, and the non-linear part is represented by the red part.
In the text

	Fig. 3. Density view of the initial condition of two interacting slow MHD waves perturbing the thermal equilibrium. Magnetic field lines are superimposed in blue.
In the text

	Fig. 4. Density view of the benchmark simulation using the SPEX_DM table and a resolution of 32002. The top left and right panels have a density of twice and twenty times the background value, respectively. The bottom left panel has a density of 95 times the background value, and the bottom right panel denotes the end of the simulation due to the small time-step restriction.
In the text

	Fig. 5. Left: ram pressure view at the end of the benchmark simulation with insets zoomed in on the filament edges, where ram pressure differences fragment the filament. Right: total velocity magnitude after the filament has formed. It clearly shows the rebound shocks due to the formation of the filament.
In the text

	Fig. 6. Field length and optically thin radiative energy loss along a selected line. Left: density view at the end of the benchmark simulation with the selected line, along the diagonal of the grid and through the filament and transition region, superimposed in red. Right: field length and optically thin radiative energy loss along the selected line.
In the text

	Fig. 7. Density views of the simulations at different resolutions. The snapshots are all taken after 4.22 h.
In the text

	Fig. 8. Evolution of the maximum density and minimum temperature for the simulations at different resolutions. Left: maximum density. Right: minimum temperature.
In the text

	Fig. 9. Evolution of the surface area and surface mass of the filament for the simulations at different resolutions. Left: surface area. Right: surface mass.
In the text

	Fig. 10. Three cooling curves with different low-temperature treatment and their density view at the end of the evolution. Top left: cooling curves. Top right to bottom left and bottom right: density views at the end of the evolution. They correspond to the Colgan_DM, JCcorona, and Colgan cooling tables, respectively.
In the text

	Fig. 11. Evolution of the maximum density and minimum temperature for the simulations with the Colgan_DM, JCcorona, and Colgan cooling tables, which only differ in their low-temperature treatment. Left: maximum density. Right: minimum temperature.
In the text

Fig. 12. Density views during the simulations with different cooling curves. Each cooling curve has its own column. The rows represent similar moments and are scaled to the same density colour bar. The first three rows are when the filament has a maximum density of approximately 2, 20, and 95 times the background density. The fourth row is the last moment, at which the time step becomes restrictive. The last row is a zoom-in onto the filament of the last moment, it also contains a zoomed-in inset of the low-density region.

In the text

	Fig. 13. Evolution of the maximum density and minimum temperature for the simulations with different cooling curves. Left: maximum density. The horizontal dotted grey lines denote the maximum densities of the rows of Fig. 12. Right: minimum temperature.
In the text

	Fig. 14. Evolution of the surface area and surface mass of the filament for the simulations with different cooling curves. Left: surface area. Right: surface mass.
In the text

Fig. 15. Density view of the full evolution of the thermal instability. Top row: onset of the thermal instability and the moment in which the low-density bootstrap measure kicks in, in the left and right panel, respectively. Bottom left: moment in which 1% of the area is at the prescribed lowest allowed density. Bottom right: last moment of the simulation. The animation of this simulation is available online.

In the text

Fig. 16. Evolution of the maximum density and minimum temperature for the simulation of the full evolution of the thermal instability. Left: maximum density. Right: minimum temperature. The evolution is shown from the start of the simulation, that is, the onset of the thermal instability. The dotted grey and red lines correspond to the top right and bottom left density views of Fig. 15, respectively.

In the text

	Fig. 17. Top panel: evolution of the surface mass of the filament. Second panel: evolution of the total surface mass in the grid. Third panel: evolution of the percentage of the number of cells or area with the prescribed lowest allowed density. Bottom panel: evolution of the surface area of the filament with respect to the total area.
In the text

	Fig. B.1. Isochoric linear growth rates as a function of the background magnetic field B0. The thermal, slow, and fast mode growth rates are denoted by the green, red, and blue lines, respectively.
In the text

	Fig. B.2. Density view during the hydrodynamical case. Top left and right: initial state and onset of condensation formation, respectively. Bottom left and right: rapid expansion and fragmentation of the condensation, respectively. The animation of this simulation is available online.
In the text

	Fig. B.3. Vorticity magnitude at the same moment as in the bottom left panel of Fig. B.2, zoomed onto the condensation.
In the text

	Fig. B.4. Curve representing the solutions to the thermal equilibrium, analogously to the S-curves used by the AGN community. The black line represents all the solutions to the equilibrium ℒ = 0 for the SPEX_DM cooling curve and corresponding heating term. The red dot indicates the initial conditions.
In the text