When Hot Meets Cold: Post-flare Coronal Rain

Solar flares represent explosive phenomena in the solar atmosphere, where 10²⁸–10³² erg of energy originally stored in the solar magnetic field can suddenly be released via magnetic reconnection (Sweet 1958; Shibata & Magara 2011). The time development of a solar flare event can be divided into three phases: a pre-flare phase, a sudden impulsive phase, and a gradual or post-flare phase (Kane 1974). The magnetic energy is released rapidly in the impulsive phase within a typical timescale of (tens of) minutes. A large fraction of this released energy is transported from the tenuous and hot corona downward to the denser and colder solar chromosphere via thermal conduction and by energetic electrons (Antiochos & Sturrock 1978; Ruan et al. 2020). This deposition of energy in the chromosphere leads to a sudden heating of the local plasma, and causes upward evaporation of the plasma to form super hot (∼10 MK) and dense (∼10¹⁰ cm⁻³) arcade-like loop systems at coronal heights. The loops return to their usual coronal conditions (∼1 MK, ∼10⁸–10⁹ cm⁻³) in the following, gradual phase, where field-guided thermal conduction and radiative losses generally contribute to the cooling process (Cargill et al. 1995; Aschwanden & Alexander 2001). These loops, visible in extreme ultraviolet (EUV), but also in Hα images, are usually called post-flare loops (PFLs; Bruzek 1964).

Thanks to dramatically increased spatio-temporal resolutions in observations, this gradual phase of solar flares is now known to show PFLs that spontaneously develop fine-scale coronal rain (Martínez Oliveros et al. 2014; Scullion et al. 2014; Jing et al. 2016; Scullion et al. 2016). In the multi-thermal coronal rain events, cool and dense rain blobs form in situ in the hot corona, to fall to the chromosphere with speeds up to 100 km s⁻¹. Coronal rain in flaring coronal loops has been reproduced in 3D simulations by Cheung et al. (2019) and Chen et al. (2021), although the dynamics of this rain has not been analyzed in any detail. Coronal rain is also observed in non-flaring coronal loops, is frequently found in loops of active regions (Leroy 1972; Levine & Withbroe 1977; Schrijver 2001; O'Shea et al. 2007; Antolin & Rouppe van der Voort 2012; Ahn et al. 2014; Antolin et al. 2015), and this type of coronal rain has been studied previously using magnetohydrodynamic (MHD) simulations (Fang et al. 2013, 2015a; Moschou et al. 2015; Xia et al. 2017; Kohutova et al. 2020). Kohutova et al. (2020) successfully reproduce coronal rain in a self-consistent 3D radiative MHD simulation, in which the radiative loss in chromosphere and corona is offset by ohmic and viscous heating. It has been generally accepted that these rain blobs are generated in a catastrophic cooling process, essentially caused by thermal instability (Parker 1953; Field 1965; Claes & Keppens 2019; Claes et al. 2020). In a catastrophic cooling event, local temperatures drop from 1 MK to below 0.1 MK within one minute, while local densities can increase by orders of magnitude (Scullion et al. 2016). Observations of flare-driven coronal rain demonstrate that this catastrophic cooling can also happen in PFLs, but this has thus far never been modeled in detail. Another phenomenon thought to have a close relationship with these sudden condensations are the so-called dark post-flare loops (DPFLs), where some PFL loops suddenly vanish from specific EUV passbands, e.g., at 17.1, 30.4, and 21.1 nm (Song et al. 2016; Jejčič et al. 2018; Heinzel et al. 2020). In these DPFLs, an EUV loop that was bright for a while suddenly darkens for several minutes, so effectively disappears between the adjacent EUV loops seen at the same height. It has been suggested that the formation of cool and dense coronal rain may contribute to EUV emission and absorption, inducing DPFL formation (Jejčič et al. 2018; Heinzel et al. 2020). However, this suggestion must still be confirmed in an ab initio model.

Here we perform an MHD simulation of a flare event from its pre-flare phase all the way into the gradual phase. This simulation finally allows us to understand the complex thermodynamic evolutions of PFLs. We (1) reproduce post-flare coronal rain, (2) quantify the chromosphere–corona mass and energy cycles during PFLs, and (3) demonstrate the intricate relationship between condensations and the disappearing EUV loops, or DPFLs.

We perform the simulation with the open-source MPI-AMRVAC code (Xia et al. 2018; Keppens et al. 2021). The simulation is 2.5D, where the domain is 2D but all vector quantities have three components. The simulation domain is given by −75 Mm ≤ x ≤ 75 Mm and 0 ≤ y ≤ 100 Mm. This simulation box has an initial resolution of 96 × 64, but an equivalent high resolution of 6144 × 4096 is achieved with our block-adaptive mesh. The governing equations are the MHD equations with effects of gravity, thermal conduction, radiative loss, and magnetic-field dissipation due to resistivity included, also shown in Ruan et al. (2020; the source terms related to fast electrons have not been activated here). The new multi-dimensional field-line-based transition-region adaptive conduction (TRAC-L) method is adopted to properly handle the chromosphere–corona interaction at affordable resolution (Zhou et al. 2021).

A relaxation has been done to obtain a static pre-flare atmosphere before we perform the flare simulation. A background heating is required to offset the energy losses due to radiative cooling and thermal conduction. Inspired by Ye et al. (2020), this background heating is a function of initial and spatio-temporally evolving values, given by

$$\begin{eqnarray}&&\begin{array}{l}{H}_{b}(x,y,t)=0.5\{\tanh [(y-{h}_{\mathrm{tra}})/{h}_{a}]+1\}\\ \times \ {N}_{{\rm{e}},0}(y){N}_{{\rm{e}}}(x,y,t){\left(\displaystyle \frac{{T}_{0}(y)}{T(x,y,t)}\right)}^{2}{\rm{\Lambda }}({T}_{0}(y)),\end{array}\end{eqnarray} \tag{ 1 }$$

where h_tra = 3 Mm, h_a = 0.1 Mm, N_e indicates electron-number density, subscript 0 indicates the t = 0 initial value and Λ(T) is the radiative cooling curve adopted in our simulations. The cooling curve from Colgan et al. (2008) is used. Many other optically thin cooling curves, such as from Schure et al. (2009), are also available in MPI-AMRVAC. Hermans & Keppens (2021) showed that changing the cooling curve does not truly change whether condensations happen, but only impacts the growth rate of the thermal instability and the morphology of the formed condensations. The initial vertical temperature T₀(y) profile in Avrett & Loeser (2008, model C7) is employed. The number density at y = 40 Mm is set to 2 × 10⁹ cm⁻³ and the initial density profile is calculated based on hydrostatic equilibrium. In the relaxation stage, a uniform vertical magnetic field is adopted. A numerically static atmosphere is obtained after a relaxation time corresponding to 3.5 hr. The local number density at y = 40 Mm decreases to about 10⁹ cm⁻³ after this relaxation. The final instantaneous background heating rate of this relaxed stage is saved and then used in the subsequent flare simulation. Such a background heating (∼10⁻⁴ to ∼10⁻³ erg cm⁻³ s⁻¹) is crucial for keeping the temperature profile outside the flare region and the PFLs, but it is not impacting strongly the condensation process. We verified that condensation can still happen without this added background heating, with an additional simulation, which can be found on the website of the ERC PROMINENT project (https://erc-prominent.github.io//media/Ruan/).

After the relaxation, the magnetic field configuration is changed. The magnetic configuration from Ye et al. (2020) is employed then, given by

$$\begin{eqnarray}&&{B}_{x}=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}{B}_{y}=\left\{\begin{array}{ll}-{B}_{0}, & \quad x\lt -\lambda \\ {B}_{0}, & \quad x\gt \lambda \\ {B}_{0}\sin [\pi x/(2\lambda )], & \quad \mathrm{else}\end{array}\right.\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{B}_{z}=\sqrt{{B}_{0}^{2}-{B}_{y}^{2}},\end{eqnarray} \tag{ 4 }$$

where B₀ = 30 G is the new initial magnetic field strength and λ = 10 Mm. Such a configuration allows magnetic reconnection. A three-stage resistivity strategy is then activated. A spatially localized resistivity inside the initial current sheet triggers magnetic reconnection in the first stage. This localized resistivity is given by

$$\begin{eqnarray}\eta (x,y,t\lt {t}_{\eta 1})=\left\{\begin{array}{ll}{\eta }_{1}[2{\left(r/{r}_{\eta }\right)}^{3}-3{\left(r/{r}_{\eta }\right)}^{2}+1], & \quad r\leqslant {r}_{\eta }\\ 0, & \quad r\gt {r}_{\eta }\end{array}\right.\end{eqnarray} \tag{ 5 }$$

where η₁ = 0.1, $r=\sqrt{{x}^{2}+{\left(y-{h}_{\eta }\right)}^{2}}$, h_η = 40 Mm, r_η = 2.4 Mm, and t_η1 = 31 s. In the second stage, we use an anomalous resistivity given by

$$\begin{eqnarray}&&\begin{array}{l}\eta (x,y,{t}_{\eta 1}\lt t\lt {t}_{\eta 2})\\ \ \ =\left\{\begin{array}{ll}0, & \quad {v}_{d}\leqslant {v}_{c}\\ {\rm{\min }}\{{\alpha }_{\eta }({v}_{d}/{v}_{c}-1)\exp [{\left(y-{h}_{\eta }\right)}^{2}/{h}_{s}^{2}],1\}, & \quad {v}_{d}\gt {v}_{c}\end{array}\right.\end{array}\end{eqnarray} \tag{ 6 }$$

where α_η = 1 × 10⁻³, h_s = 10 Mm, t_η2 = 7.78 minutes, v_d (x, y, t) = J/(eN_e), and v_c = 128,000 km s⁻¹. The resistivity is set to zero in the third stage t > t_η2 to force the flare to enter the gradual phase, where only numerical dissipation happens. This means that the explicit influence of Joule heating on the cooling of PFLs and the triggering of the condensation is neglected. The resistivity strategy used in our first and second stage is similar to that in Yokoyama & Shibata (2001).

We employ symmetric boundary conditions for number density, pressure, and magnetic-field components, while anti-symmetric conditions are employed for velocity components at the left and right boundaries. At the upper and bottom boundaries, density and pressure are fixed to their initial values. The magnetic-field components at the bottom boundary are also fixed to the initial values. An anti-symmetric condition is applied for the x-component of the magnetic field at our upper boundary, while the other two components of the magnetic field employ symmetric conditions there. The velocity at the upper boundary is set to zero, as we are not interested in an erupting flux-rope structure here, but focus on the PFLs. Anti-symmetric conditions are employed for the velocity components at the bottom boundary.

The SXR emission is calculated with the method reported in Pinto et al. (2015). The EUV emissions are calculated with the contribution function provided by the CHIANTI database and the optically thin assumption (Del Zanna et al. 2015).

Figure 1 demonstrates the full evolution in magnetic topology, from the pre-flare current sheet to the formation of PFLs at the impulsive phase (Figure 1(b)), extended to the entire evolution of the PFLs through the gradual phase (Figures 1(c), (d)). There is a vertical current sheet separating regions of opposite field directions at the beginning of our simulation (Figure 1(a)). Magnetic reconnection inside this current sheet produces closed magnetic arcades below and a flux rope above a reconnection site in the impulsive phase (t ≲ 8 minutes). The magnetic energy is converted into heat at the reconnection site and then conducted into the chromosphere. Deposition of the thermal energy there leads to a sudden increase in local pressure and this produces upward evaporation flows. They fill the generated coronal loops with hot plasma (∼10 MK) and increase the electron-number density by one order of magnitude in the PFLs (Figures 1(b), (f)). Thereafter, the flare enters the gradual phase when we have a rapidly decreasing magnetic-reconnection speed (t ≳ 8 minutes). The PFL temperature decreases slowly due to thermal conduction and radiative losses in this gradual phase, but suddenly triggers thermal instability near t ≃ 35 minutes. The loop density also decreases in this period, but is still much higher than the external coronal density (Figure 1(c)). Energetic electron-beam deposition can also lead to chromospheric evaporations, but this mechanism is less efficient in weak flares where the beam energy flux is small (Fisher et al. 1985). Electron-beam heating has not been included in our simulation for this reason, though this heating has been successfully reproduced in a 2.5D flare simulation in our previous work (Ruan et al. 2020). The peak flux in our simulation is about 4 × 10⁻⁷ W m² when assuming that the loop depth in the third direction is 100 Mm (see Figure 1(i)), therefore the simulated flare is a B-level flare.

In the impulsive phase, the energy gain of the flare loop is determined by the reconnection process whereas the energy loss is mainly handled by thermal conduction, which overwhelms the radiative cooling by more than one order. Therefore, the flare temperature can be predicted from the coronal magnetic field strength, density, and flare-loop length or simply from coronal plasma β with the scaling laws provided in Yokoyama & Shibata (1998, 2001) and Takasao et al. (2015). This flaring temperature estimated from Equation (1) of Yokoyama & Shibata (1998) well matches our loop-top temperature in the impulsive phase (∼20 MK) where a combination of coronal magnetic field strength 30 G, coronal number density 10⁹ cm⁻³, and half-loop length 25 Mm gives an estimated temperature of 26 MK.

Catastrophic cooling driven by thermal instability condenses local plasma in situ and leads to high density (close to 10¹¹ cm⁻³) and cold (close to 0.01 MK) structures in the coronal PFLs (Figures 1(d), (h)). Figure 1(i) shows how the SXR flux reaches its peak value at the impulsive phase, to then gradually decrease as the loop temperature drops. The temporal evolution of the total coronal rain material (i.e., T_e < 0.1 MK, N_e > 10¹⁰ cm⁻³ and y > 5 Mm) is also illustrated in Figure 1(e). We note that sudden condensations happen in two successive events during the entire simulated period. These are located in different loop systems, with the second rain event appearing at a higher altitude.

Once condensations happen within PFLs, the formed cold and dense plasma structures will likely fall down from coronal heights due to gravity and hence appear as observed coronal rain blobs. This is fully reproduced in our simulation as demonstrated in Figure 2. Cold plasma is formed at a PFL loop top at the beginning of the runaway condensation (Figure 2(a)). Thereafter, the cold structure extends to lower and higher loops, meanwhile sliding down to one side along the magnetic post-flare arcade (Figures 2(b), (c)). This falling cold plasma gets accelerated to a speed of ∼100 km s⁻¹ by gravity before it enters the chromosphere (Figures 2(d)–(f)). Such a speed is close to that found in coronal rain observations (Martínez Oliveros et al. 2014). Considering an acceleration timescale of 10 minutes, the average acceleration rate is lower than the acceleration of gravity. A detailed analysis of the rain blob acceleration process for non-flaring (or quiescent) coronal rain was given in Fang et al. (2013). Different from Fang et al. (2013), the guide field B_z of the loops where condensation happens is nearly zero. But since coronal rain is formed through thermal instability, which depends only slightly on the magnetic field, whether or not there is a guide field should have little influence on our results.

Here we investigate the mass and energy cycles in our entire 100 minute simulation of the gradual phase and the role of condensation in it. To do so, we track mass and energy budgets into the coronal part (y > 5 Mm) of a loop section in which the first round of condensation happens. This loop section is always bound by (a) the evolving magnetic field line with a fixed footpoint at x = −25 Mm at our lower y = 0 boundary; and by (b) a similarly evolving field line with footpoint at x = −15 Mm (Figure 3(a)). As seen in Figures 2 and 3(a), this region gets emptied during the first round of condensation as matter collects along the field lines into localized rain blobs. This entire loop system continuously moves downward as a result of the above reconnection dynamics and the corresponding area of the selected region continuously decreases during the simulation as quantified by the solid line in Figure 3(b). The area-integrated total energy (kinetic, thermal, and magnetic combined; dashed line in Figure 3(b)), area-integrated mass (solid line in Figure 3(c)), and the area-averaged temperature (solid line in Figure 3(d)) of this region also continuously decrease before the first condensation (t ≃ 35 minutes). Plasma that evaporated upward into the corona in the previous impulsive phase now leaks back into the chromosphere in this period, as seen in the downward mass flux (dashed line) in Figure 3(c). The decrease in temperature leads to a decrease in the atmospheric scale height. However, this downward leakage of coronal plasma is severely reduced when the condensation happens after t ≳ 35 minutes since the gas pressure in the coronal part of the loop then decreases rapidly. Therefore, we see a drop of the downward mass flux in Figure 3(c) near t ≈ 38 minutes when condensations are fully formed. Later on, the downward mass flux experiences a sudden increase, exactly when the cool coronal rain material goes through the lower y = 5 Mm boundary of the studied region. The rain material is dense and heavy, so the material can still drop even if there is an upward pressure gradient below it. The area-integrated mass reaches its minimum value when all cool material leaves and enters the chromosphere. Thereafter, plasma from the chromosphere is injected to the loop again, due to the low pressure inside the loop, leading to an upward mass flux. The total mass then gradually returns to its value before condensation.

The changing coronal energy budget shows a similar tendency with changing coronal mass cycle. The energy also experiences a decrease in the pre-condensation and during the condensation phase to then experience an increase after the condensation merges into the chromosphere due to renewed plasma injection (Figure 3(b)). It has been suggested that thermal conduction determines the PFL's energy loss at the beginning of the gradual phase of a flare and that subsequently, radiative losses will become dominant (Cargill & Klimchuk 2004). Our simulation shows that this suggestion is correct before and also during condensation. The efficiencies of radiative cooling and of thermal conduction to the energy loss in the selected region are compared in Figure 3(d). The contribution of thermal conduction is greater than radiative losses for t ≲ 25 minutes, but conductive losses drop gradually owing to the decreasing temperature gradient. At t ≈ 25 minutes, the average temperature is about 3 MK, and the efficiency of radiative losses becomes most prominent. However, conductive losses become stronger than radiative losses again when the condensation vanishes from the loop system, as collisions between re-injected flows from both footpoints make the loop hot again and the radiative losses drop for a while due to the decrease in loop density.

The first round of condensation happens near t ≈ 35 minutes. The temporal evolutions of instantaneous maximum/minimum temperature/number density in the condensation region (the same as marked in Figure 3(a)) are illustrated in Figure 4. Triggering of thermal instability switches the radiative cooling process from linear to nonlinear, and then leads to a catastrophic cooling of local plasma. We get an average temperature decreasing rate of −9000 K s⁻¹ in the catastrophic cooling phase. In contrast, the cooling rate before catastrophic cooling is −3000 K s⁻¹. As a result of catastrophic cooling, the local temperature decreases from 0.2 to 0.02 MK within half a minute (Figure 4(a)), while the local number density increases by one order, from 10¹⁰ to 10¹¹ cm⁻³ (Figure 4(b)). Noting that the cooling rate is density-dependent, the cooling rate we get in the catastrophic cooling phase is on the same order as that observed in Scullion et al. (2016) (−22,700 K s⁻¹).

Zoom In Zoom Out Reset image size

A quasi-periodic pulsation (QPP) with a period of ∼3 minutes appears in the maximum density curve, just after the rain condensation disappears from the PFL system. This QPP is caused by the injected flows mentioned in the previous section, refilling and reheating the PFL. The density variation due to these flows along a field line is shown in Figure 4(c). Injected flows propagating from one footpoint to the other produce reflected slow-mode waves. Such a process has previously been studied in Fang et al. (2015b) for isolated loop systems. The sharp density changes in Figure 4(c) are shocks ahead of the injection flows and the wave fronts of the slow-mode waves. Such QPPs hence reflect density variations in the low corona due to combined flows and wave propagation. The period of our QPP is close to the time for the slow-mode wave to propagate from one footpoint to the other as the wave speed is about 300 km s⁻¹.

In synthesized EUV images of our simulation, coronal loop(s) appear in the gradual phase. An example is shown in Figure 5 at the 17.1 nm wave band. Interestingly, this loop disappears for about 10 minutes during the evolution, as demonstrated in panels (a)–(e) of Figure 5. The sudden darkening of the bright EUV loop in our simulation resembles the dark post-flare loop (DPFL) phenomenon, previously observed at the same 17.1 nm passband and with a similar timescale of 10 minutes (Song et al. 2016). Observed DPFLs and the disappearing EUV coronal loop in our simulation also share the same time evolution of integral EUV flux: the EUV flux reaches its minimum value when the darkening happens (compare Figure 3 in Song et al. 2016 and our Figure 5(f)). The formation of a darkened coronal EUV loop needs to satisfy one or both of the following conditions: (1) an emission drop in an existing bright EUV loop and (2) an absorption of the background EUV emission (Anzer & Heinzel 2005). Here we explain how these conditions can be satisfied based on our simulation results.

The role of cool and dense coronal plasma in EUV emission and absorption leading to DPFLs has been emphasized previously (Jejčič et al. 2018; Heinzel et al. 2020), with coronal rain observations (Martínez Oliveros et al. 2014; Scullion et al. 2014, 2016; Jing et al. 2016) and our simulation results showing how this cool and dense plasma can be generated in PFLs. The drop in the loop emission is understood from our simulation: loop temperature changes relate to nearby rain condensation. Indeed, the temperature of bright loops in this 17.1 passband (about 10^5.8 K) is not far from the critical temperature for the onset of catastrophic cooling (this is density- and temperature-dependent, but generally happens below 2 MK according to Cargill & Bradshaw 2013). Condensation can be triggered by thermal instability near bright loops and these suddenly formed structures grow fastest across magnetic field lines (which is counterintuitive due to the field-aligned thermal conduction, but see Fang et al. 2013; Claes & Keppens 2019; Claes et al. 2020; Zhou et al. 2020). Rain that forms near (Figures 5(a)–(b)), and ultimately inside (Figure 5(c)), the bright coronal loops is thus causing the darkening as illustrated in Figures 5(c)–(d). Once condensation happens, a lot of plasma will collect into a small region, so the plasma density elsewhere in the loop decreases. To maintain pressure balance, these evacuated loop regions will increase in temperature, so EUV brightness decreases due to these combined temperature and density changes. Ultimately, the loop refills and brightens once more (Figure 5(e)).

To fully understand coronal rain in PFLs and the mass and energy budget in the gradual phase of solar flares, we performed a flare simulation from onset all the way into the long-duration post-flare phase. Post-flare coronal rain successfully and, at least in our simulation, repeatedly forms. Flare-induced rain is a result of catastrophic cooling by thermal instabilities, and our simulation shows successive rain formation at increasing heights in the PFL configuration. Rain blobs falling into the chromosphere lead to sudden mass drops in PFLs, but their mass increases again due to spontaneously forming injection flows. Therefore, coronal rain events do not accelerate the PFL mass loss in the longer term. Such longer-term mass loss is more determined by the change of the gravity scale height due to the cooling of PFLs.

Both thermal conduction and radiative losses contribute to the energy budget in PFLs. Thermal conduction dominates the PFL energy loss at the beginning of the gradual phase. Thereafter, it becomes less efficient than radiative losses, owing to decreases in loop temperature and in temperature gradient. However, thermal conduction can efficiently recover again after a condensation falls to the chromosphere, as the loop again reaches a high temperature. In this phase, an emptied loop refills and can show a slow-wave-related QPP.

We showed that the formation of DPFLs can result from post-flare rain condensation. Condensation changes loop temperatures and can make existing bright EUV loops temporarily disappear for several minutes. This timescale of EUV loop darkening is identical to observed DPFLs.

A recent paper by Reep et al. (2020) suggests that an additional energy supplement (heating) is necessary in the formation of flare-driven coronal rain based on a 1D parameter study. In their simulations, a 1D flare loop is filled with hot/dense plasma by non-thermal electron-beam-driven chromospheric evaporations at the beginning and then gradually cools down due to radiative loss. The loop will directly cool down to chromospheric temperature without condensation happening if there is no extra heating provided during the cooling, in which the coronal density also drops by orders of magnitude during the cooling. It is interesting to note that condensation can still happen inside our 2D post-flare loop, even if the background heating is switched off (see the additional simulation mentioned in Section 2). Compression caused by thermal pressure of the surrounding environment or by the Lorentz force may serve as an additional energy supplement, as it has been shown in Figure 3(b) that the loop area decreases gradually during the cooling phase. When this interaction with the surrounding environment is taken into account, the coronal temperature and density are not likely to drop by orders of magnitude at the same time as shown in Reep et al. (2020). This difference between 1D simulations and 2D simulations of post-flare loop development will be studied in future work.

This work was supported by the ERC Advanced Grant PROMINENT, an FWO project G0B4521N 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 program (grant agreement No. 833251 PROMINENT ERC-ADG 2018). This research is further supported by Internal funds KU Leuven, through the project C14/19/089 TRACESpace. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government, department EWI.