Issue 
A&A
Volume 634, February 2020



Article Number  A36  
Number of page(s)  13  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201937078  
Published online  03 February 2020 
Twodimensional simulations of coronal rain dynamics
I. Model consisting of a vertical magnetic field and an unbounded atmosphere^{⋆}
^{1}
Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain
email: dmartinez@iac.es
^{2}
Departament de Física, Universitat de les Illes Balears, 07122 Palma de Mallorca, Spain
^{3}
Institute of Applied Computing & Community Code (IAC3), UIB, Spain
^{4}
Departamento de Astrofísica, Universidad de La Laguna, 38205 La Laguna, Tenerife, Spain
Received:
7
November
2019
Accepted:
14
December
2019
Context. Coronal rain often comes about as the final product of evaporation and condensation cycles that occur in active regions. Observations show that the condensed plasma falls with an acceleration that is less than that of free fall.
Aims. We aim to improve the understanding of the physical mechanisms behind the slower than freefall motion and the twostage evolution (an initial phase of acceleration followed by an almost constant velocity phase) detected in coronal rain events.
Methods. Using the MANCHA3D code, we solve the 2D ideal magnetohydrodynamic equations. We represent the solar corona as an isothermal vertically stratified atmosphere with a uniform vertical magnetic field. We represent the plasma condensation as a density enhancement described by a 2D Gaussian profile. We analyse the temporal evolution of the descending plasma and study its dependence on such parameters as density and magnetic field strength.
Results. We confirm previous findings that indicate that the pressure gradient is the main force that opposes the action of gravity and slows down the blob descent, and that larger densities require larger pressure gradients to reach the constant speed phase. We find that the shape of a condensation with a horizontal variation of density is distorted during its fall because the denser parts of the blob fall faster than the lighter ones. This is explained by the fact that the duration of the initial acceleration phase and, therefore, the maximum falling speed attained by the plasma, increases with the ratio of blob to coronal density. We also find that the magnetic field plays a fundamental role in the evolution of the descending condensations. A strong enough magnetic field (greater than 10 G in our simulations) forces each plasma element to follow the path given by a particular field line, which allows for the description of the evolution of each vertical slice of the blob in terms of 1D dynamics, without the influence of the adjacent slices. In addition, under the typical conditions of the coronal rain events, the magnetic field prevents the development of KelvinHelmholtz instability.
Key words: Sun: corona / Sun: filaments / prominences
Movies associated to Figs. 1, 8 and 10 are available at https://www.aanda.org
© ESO 2020
1. Introduction
Coronal rain is a phenomenon that occurs in the coronal loops of active regions. The hot and dim coronal plasma condenses into cold elongated blobs that then fall towards the solar surface following the leg of the loop. It was already over fifty years ago when the first observational evidence of this phenomenon was presented by Beckers (1962). This, along with subsequent works such as Bruzek (1969), Foukal (1978) or Engvold et al. (1979), reported falling speeds ranging from 40 to 60 km s^{−1} but did not study the kinematics of the condensations. Later Loughhead & Bray (1984), Xu (1987), and Heinzel et al. (1992) performed a tracking of the positions and speeds of the blobs as functions of height and determined that they followed a slower than freefall motion.
More recent investigations measured velocities of around 30–40 km s^{−1} near the top of the loops, but values of up to 150 km s^{−1} further down the loop leg (Wiik et al. 1996; De Groof et al. 2004; Müller et al. 2005; Antolin et al. 2010; Antolin & Rouppe van der Voort 2012). These observations, in addition to the works of Schrijver (2001), Zhang & Li (2009), and Antolin et al. (2012), revealed that two stages can be distinguished in the motion of the condensations. The plasma starts falling with an acceleration that is less than that of gravity and decreases until an almost constant speed phase is reached. Therefore, there must be a physical mechanism that opposes the action of gravity and slows down the descent of the cold blobs considerably. This twostage feature is not unique for coronal rain events but a similar behaviour has also been detected in the descent of plasma after a prominence eruption (Xue et al. 2014).
The evidence that provides information about the low temperature of these blobs comes from the fact that they are typically observed in lines associated with cold plasma, such as Lyα, Hα, or He II 304 Å but do not appear in hot coronal lines like 171 Å (De Groof et al. 2004; Müller et al. 2005). At these temperatures, which are usually found in the chromosphere and the transition region, the plasma is not expected to be fully ionised but partially ionised, that is, it is expected to contain neutral species. The presence of this neutral component has been demonstrated by the observations of Chae (2010) and Ahn et al. (2014). When it comes to how the coronal plasma cools to such low temperatures and the condensation forms, several explanations have been developed. For instance, Antolin & Verwichte (2011) proposed that a large density enhancement can appear in a transversally oscillating coronal loop due to the ponderomotive force associated with large amplitude oscillations (Terradas & Ofman 2004). Some alternative mechanisms involve the presence of magnetic null points and explain the creation of the dense blob either by a pressure pulse that excites an entropy mode (Murawski et al. 2011) or by interchange reconnection (Mason et al. 2019).
However, most models that try to replicate the formation of coronal rain take the catastrophic cooling due to a thermal instability (Parker 1953; Field 1965) as the main ingredient for their recipe. In these scenarios, heating concentrated at the loops base produces the evaporation of the plasma from the solar surface that then starts filling the magnetic structures. The heating can be impulsive or steady. In the former case, an intense burst of heating, like that caused by a flare event, stops abruptly and due to the radiative losses, the plasma cools down to chromospheric temperatures (Antiochos 1980). This mechanism explains coronal rain events in postflare loops (Liu et al. 2015; Scullion et al. 2016), in which blobs form and fall only once. Yet it does not explain how an active region can undergo several cycles of coronal rain. Observations by Antolin & Rouppe van der Voort (2012), Antolin et al. (2015), Froment et al. (2015, 2017), and Auchère et al. (2018) show that the lifetime of the magnetic field structure is larger than the time taken for a condensation to form and fall, and that numerous blobs can consecutively descend along the same loop path. Such behaviour is better understood if the heating is assumed to be steady (or with a high frequency) and to decrease with height. In this case, the hot evaporated plasma starts to accumulate at the apex of the loop, where the heat input is not enough to maintain the energy balance with the radiative losses. This situation is referred to as a state of thermal nonequilibrium (TNE; Antiochos et al. 1999). Radiative losses increase with decreasing temperatures and increasing densities, leading to a catastrophic cooling and the fast condensation of the plasma. Then, since the condensation is not in mechanical equilibrium, it starts falling down under the action of gravity. The validity of this TNE scenario, which can also be applied to explain the origin of solar prominences (Antiochos & Klimchuk 1991), has been demonstrated by both hydrodynamic (Kuin & Martens 1982; Martens & Kuin 1983; Karpen et al. 2001; Müller et al. 2004) and magnetohydrodynamic (Fang et al. 2013, 2015; Moschou et al. 2015) simulations.
Regarding the reasons behind the slower than freefall motion, Schrijver (2001) suggested that the speed of the condensations is not determined by gravity but by the evolution of the internal pressure of the loops. This suggestion was later supported by the works of Antolin et al. (2010), who reached the conclusion that the structure and dynamics of the blobs are more sensitive to variations in pressure than to gravity, and Oliver et al. (2014), who reproduced the twostage evolution with a 1D numerical model and demonstrated the role of the pressure gradient in slowing down the plasma fall.
The previous paragraphs outlined several reasons for why the study of coronal rain is such a relevant topic in the field of solar physics. It can improve the understanding of the processes involved in coronal heating and thermal instabilities. Then coronal rain events can be used as magnetic field tracers and help to better describe the multistranded structure of the coronal magnetic field (Harra et al. 2014; Scullion et al. 2014). Furthermore, the low temperatures of the descending blobs allow to explore the interaction between the ionised and the neutral components of partially ionised plasmas (Oliver et al. 2016).
The present work is the first of a series in which we aim at obtaining a better comprehension of the dynamics of coronal rain. This is a quite complex task, as evidenced by the amount of past research devoted to the study of this phenomenon. So, instead of trying to reproduce every phase of the coronal rain cycle, here we follow the same approach as in Oliver et al. (2014) of focusing only on its final stages, that is, on the descent of the condensation towards the solar surface. This will allow us to explore in great detail the forces that determine the evolution of the falling condensations. In this series, we improve the numerical 1D model used in Oliver et al. (2014, 2016) by taking into account the 2D evolution of the falling blobs and including the effects of the presence of the magnetic field. In the first installment of the series, we start by assuming that the plasma is fully ionised and we use a singlefluid numerical code to compute its temporal evolution. We apply the simple model of an unbounded isothermal atmosphere with a vertical magnetic field and investigate the dependence of the blob motion on parameters such as its density and the magnetic field strength. We also analyse the conditions that could lead to development of the KelvinHelmholtz instability (Chandrasekhar 1961) during the plasma descent. In forthcoming works, we will add new features to this model, such as more complex configurations of the background magnetic field or the inclusion of a large density enhancement at the lower boundary to represent the chromosphere. In addition, we will use a twofluid code to take into account the fact that the plasma in the condensations is partially ionised.
This paper is organised as follows. In Sect. 2, we present the equations that govern the motion of the plasma and the numerical setup we use in our simulations. Section 3 contains the analysis of the simulations we performed, which includes the study of the kinematics and dynamics of the falling blobs, in addition to the research of the possible development of instabilities. Finally, the conclusions of this investigation are summarised in Sect. 4.
2. Model
2.1. Governing equations
For the present research, we assumed that the plasma is fully ionised. To analyse its dynamics, we considered the singlefluid approximation and neglected nonideal effects. Therefore, the temporal evolution of the density (ρ), velocity (V), energy (E) and magnetic field (B) is described by the set of ideal magnetohydrodynamics (MHD) equations:
where P is the pressure, μ_{0} is the magnetic permeability, 𝕀 is the identity tensor and g is the gravitational acceleration. The total energy is defined as
where γ is the adiabatic index.
In addition, to close the system of equations we used the perfect gas equation, which is given by
where m_{p} is the proton mass, k_{B} is the Boltzmann constant, T is the temperature, and is the mean atomic weight. In this work, we assume that the plasma is only composed of hydrogen, which corresponds to a value of .
2.2. Numerical setup
All the simulations analysed in the present work were carried out with the MANCHA3D code (Khomenko & Collados 2006; Felipe et al. 2010; GonzálezMorales et al. 2018). This numerical code solves the nonideal nonlinear equations of singlefluid MHD in three dimensions. The nonideal terms include the ambipolar diffusion due to the presence of neutrals in the plasma, the Hall effect, the battery effect and radiative losses. However, for the present investigation we neglected the dependence on the ydirection, so our simulations were twodimensional, and we did not take into account those nonideal effects but we retained only the terms that correspond to the ideal MHD description.
The MANCHA3D code treats each variable in the system of Eqs. (1)–(4) as a combination of a static equilibrium value plus a perturbation (that can be nonlinear), and computes the temporal evolution of the perturbations. Therefore, for instance, the variables density and pressure can be written as
where ρ_{eq} = ρ_{eq}(x, z) and P_{eq} = P_{eq}(x, z) are the equilibrium values, and ρ_{1} = ρ_{1}(x, z, t) and P_{1} = P_{1}(x, z, t) are the perturbations, respectively.
The code uses a centered sixth order scheme for spatial integration and the temporal evolution is computed by means of an explicit fourstep Runge–Kutta method. To prevent the exponential growth of highfrequency noise related to the small scales of the numerical grid, artificial diffusion terms are added to the righthand side of Eqs. (1)–(4). This artificial diffusion is formed by a shock resolving term and a hyperdiffusivity part, as described in Vögler et al. (2005). The coefficients of these diffusion terms are kept as low as possible, so that the actual dynamics of the studied system is not noticeably affected while, at the same time, the stability of the simulations is guaranteed.
2.2.1. Equilibrium state
For the equilibrium state in our simulations, we considered the solar corona as an isothermal vertically stratified atmosphere, with a temperature given by T_{0} and a uniform vertical magnetic field given by B_{0} = (0, 0, B_{0}), where the zaxis points upward. Mathematically, this equilibrium is described as follows. From Eq. (2), if we ignore the temporal variations, assume V = 0 and consider the variations along the zdirection, we get
which, using Eq. (6), transforms into
Solving the previous equation, we get that the equilibrium state is described by
where ρ_{0} and P_{0} are the values of density and pressure at z = 0, respectively, and satisfy the relation , and H is the vertical scale height, given by
In all the simulations of the present study, the values of density and pressure at the base of the corona were ρ_{0} = 5 × 10^{−12} kg m^{−3} and P_{0} = 0.165 Pa. The chosen temperature was T_{0} = 2 × 10^{6} K, which, assuming the typical value of g ≃ 274 m s^{−2} for the solar surface gravity, gives a vertical scale height H ≃ 120 Mm.
2.2.2. Initial and boundary conditions
At the initial step of our simulations, the blob was represented by a density enhancement, described by the following twodimensional Gaussian profile:
where ρ_{b0} is the density at the center of the blob, x_{0} and z_{0} represent its initial position and Δ is proportional to its width. This enhancement was inserted into the code through the variable corresponding to the perturbation of density, ρ_{1}. Therefore, the initial condition for the total density was given by ρ(x, z, t = 0) = ρ_{eq} + ρ_{1}(x, z, t = 0) = ρ_{eq} + ρ_{b}(x, z).
In the xdirection, we considered a physical domain of length L_{x} = 10 Mm, that goes from x_{left} = −5 Mm to x_{right} = 5 Mm. The height of the domain, L_{z}, and its bottom and top limits varied from one simulation to another, as we explain later. The simulation with the largest height used a value of L_{z} = 100 Mm. With a resolution of 10 km, this domain was represented by a mesh of 10^{3} × 10^{4} points. We used periodic boundary conditions in the xdirection while in the vertical direction we used Perfectly Matched Layers (PML; Berenger 1994, 1996; Parchevsky & Kosovichev 2009). These layers prevented the reflection of perturbations at the boundaries, which may introduce unphysical behaviors in the system, and allowed us to use a much smaller verticaldomain in comparison with that employed in Oliver et al. (2014).
3. Results
According to the observations, once a dense blob is formed in the corona, it starts to fall towards the solar surface under the action of gravity. However, it follows a slower than freefall motion. With our simulations, we wanted to replicate this behaviour to determine which mechanisms act against the gravity to slow down the fall and to study how the motion of the blob depends on parameters such as its density or the strength of the background magnetic field. For this reason, we performed a series of simulations combining different values of the peak density of the blob, ρ_{b0}, and the magnetic field strength, B_{0}. Other parameters such as the width of the blob or its initial position, were kept fixed throughout this series of simulations, with Δ = 500 km and x_{0} = 0, z_{0} = 50 Mm, respectively.
In Sects. 3.1 and 3.2, we analyse the results of a series of simulations in which we combined three values of the blob density (ρ_{b0} = 1, 5, 10 × 10^{−10} kg m^{−3}) and three magnetic field strengths (B_{0} = 5, 10, 20 G). In Sect. 3.3, we also include results from simulations with B_{0} = 1 G and without magnetic field. The set of parameters used in the different simulations are summarised in Table 1. For the sake of simplicity, when discussing the results of a given simulation, instead of repeating the values of the chosen parameters, we will refer to it by the abbreviated name presented in Table 1.
Parameters of simulations.
3.1. 2D blob kinematics
In the first place, we focus on the output of the simulation denoted as I5, which has low values of the blob density and the magnetic field. The height of the physical domain in this simulation was L_{z} = 80 Mm, with the bottom boundary at z_{bot} = 0.
Figure 1 displays several snapshots of the temporal evolution of the density. In addition, a horizontal red dashed line has been drawn to mark the position that the blob would have if it followed a freefall motion. The blob is initially situated at the position z_{0} = 50 Mm. At this height, the density ratio between the condensation and the surrounding hotter corona is ∼30. Then, since the blob is not in mechanical equilibrium, it starts to fall under the action of gravity. However, it can be clearly seen that its descent is much slower than what is expected for a freefalling object: as shown by the rightmost panel, after a time of 400 s the central part of the blob has only fallen a distance of ∼7 Mm, while a freefalling plasma would be out of the limits of the plot.
Another noteworthy feature can be appreciated in Fig. 1. It can be seen even at the initial instants of the simulation that the blob is not able to keep its original shape as it falls. The central denser part moves faster than the lighter sides, which leads to the development of a characteristic Vshape. We later explain the reason why this shape appears.
To better understand the behaviour described in the previous paragraphs, we selected different points of the blob and investigated their kinematics during the whole simulation. Thus, we performed several vertical cuts along the numerical domain and computed the position and velocity of the peak of density in each cut for each snapshot of the simulation. The positions of these peaks are marked as black asterisks in Fig. 1. The results of this study are shown in Fig. 2, which includes plots of the position (left panel) and falling speed (right panel) of the selected points as functions of time. All the points start at the same height, but as time advances, it can be seen that those associated with a smaller density start to lag behind the denser ones.
Fig. 1. Density at different times of the simulation I5. The horizontal red dashed lines represent the position that would correspond to a freefalling object. The black asterisks show the position of the points whose kinematics is represented in Fig. 2. An animation of this figure is available online. 

Open with DEXTER 
Fig. 2. Left panel: position as a function of time of selected points from simulation I5. Right panel: velocity as a function of time. Vertical dashed lines represent the time of the snapshots displayed in Fig. 1. The colours of the lines are related to the density scale shown in the previous figure. 

Open with DEXTER 
The different lines in the plots of position and falling speed show that, qualitatively, all the points have the same behaviour: at first, they are found in an acceleration phase, that is then followed by a stage in which the falling speeds become constant and even suffer a small decrease. The disparities appear in the duration of the acceleration phase, the values of the maximum falling speeds, and the heights reached by each point. It can be corroborated that larger densities lead to longer acceleration phases and larger falling speeds. As a consequence, the points with a larger density reach lower heights faster. The twostage behaviour represented by the plots of position as a function of time is in good agreement with the timedistance diagrams obtained from observations of coronal rain events (De Groof et al. 2004; Antolin et al. 2010).
To better study the relation between the density of the descending plasma and its maximum falling speed, we also performed simulations with larger values of the peak density: ρ_{b0} = 5 × 10^{−10} kg m^{−3} and ρ_{b0} = 10^{−9} kg m^{−3}, which corresponded to density ratios at the blob peak of ∼150 and ∼300, respectively. For these simulations, the height of the domain was L_{z} = 100 Mm, starting at z_{bot} = −40 Mm. The reason for using this different height was that the blob fell much faster in these simulations than in the previous one and arrived at positions below z = 0 (this means that, in a more realistic model, the blob would impact against the chromosphere). We also repeated all these simulations but increasing the magnetic field strength to 10 and 20 G, and we again selected different points of the condensation to investigate their kinematics. The results of this research are illustrated by Fig. 3, where the maximum falling speed is represented as a function of the density ratio between the plasma of the blob and the surrounding corona at t = 0. The black circles show the results obtained from the 1D simulations analysed in Oliver et al. (2014, 2016). The remaining symbols correspond to the 2D simulations described in the present work: triangles, squares and diamonds represent the cases with ρ_{b0} = 1, 5, 10 × 10^{−10} kg m^{−3}, respectively, and red, green, and blue correspond to simulations with B_{0} = 5, 10, 20 G, respectively (the specific colour and symbol for each simulation is shown in Table 1).
Fig. 3. Maximum falling speed as a function of density ratio. Black circles correspond to the results from 1D simulations. Triangles, squares and diamonds show the results of 2D simulations with ρ_{b0} = 1, 5 and 10 × 10^{−10} kg m^{−3}, respectively. Red, green and blue colour correspond to the cases with B_{0} = 5, 10 and 20 G, respectively. 

Open with DEXTER 
A clear trend is found in Fig. 3: an increase in the density ratio leads to a larger falling speed. The black dashed line represents the fitting of the numerical data to the analytical expression
where Θ is the density ratio. We found that the observed trend is well described by the fitting parameters a = 2.56 and b = 0.64.
It is remarkable to observe the good agreement between the 1D and the 2D results displayed in Fig. 3. Yet it is also unexpected. Here we used a twodimensional instead of a onedimensional approach to investigate coronal rain and we included some effects, like the presence of a magnetic field, that were not taken into account in Oliver et al. (2014, 2016). So, why is there such a good agreement in the results in this large range of parameters when the models are, in principle, clearly different?
Figure 3 also shows that there are some points corresponding to 2D simulations that do not conform to the overall trend. These “misbehaving” points are obtained from simulations with large densities and small values of magnetic field, while those with higher values of the initial magnetic field strength follow well the 1D general trend.
The results presented in Figs. 1–3 allow us to describe the motion of the falling blobs and draw a relation between some parameters but they do not explain the reasons of the observed behaviours. For instance, up to this point of the research we can say that for small enough magnetic fields the 2D results depart from the 1D predictions. On the other hand, the agreement between the 1D and 2D results for larger values of B_{0} indicates that the dynamics of each vertical slice x = const. is independent of the neighbouring slices. To explain all this, we need to look into the dynamics of the system. This is done in the following section.
3.2. Blob dynamics
From the equation of momentum conservation we can determine the different accelerations the plasma is subject to during its descending motion. Equation (2) can be recast as
where is the current density. In this form, we can identify the first term on the righthand side as the inertial term and the second one as the acceleration due to the pressure gradient. The third one is the Lorentz term related to the presence of the magnetic field and the last one is the acceleration due to gravity.
In Fig. 4 we represent the different components of the acceleration at the density peak of the condensation as functions of time, that is, following the motion of the black cross at x = 0 in Fig. 1. The results are taken from simulation I5. The top panel displays the accelerations along the horizontal direction, while the bottom panel shows the accelerations in the vertical direction.
Fig. 4. Dynamics of the density peak in the simulation I5. Top: accelerations in the xdirection. Bottom: acceleration in the zdirection. 

Open with DEXTER 
First, we focus on the horizontal accelerations. The top panel of Fig. 4 shows that the inertial component of the acceleration (represented by the solid red line) is negligible and we do not have an acceleration due to gravity because it is directed along the zdirection. Thus, the dynamics in the horizontal axis are mainly determined by the pressure gradient (green line) and the Lorentz term (black line). In absolute values, both terms reach very large magnitudes. Nevertheless, they are oriented in opposite directions and they balance each other almost perfectly, which leads to a total horizontal acceleration (orange line) that is almost zero: as shown in the inset, it varies in the range [ − 1, 1] m s^{−2}. The same fact is observed at different points away from the x = 0 axis. Horizontal accelerations at these outofaxis points are slightly larger than those for x = 0, with absolute values of up to 5 m s^{−2}, which are still negligible in comparison to the vertical accelerations. This means that the motion of the blob is mainly directed along the zaxis.
On the other hand, the bottom panel of Fig. 4 shows that the inertial term and the Lorentz term are both negligible in the vertical direction. Consequently, the descending motion depends on the relation between the pressure gradient and gravity. At t = 0, there is no pressure gradient and the only acting force is gravity. Thus, plasma starts falling with an acceleration of ∼274 m s^{−2}. But immediately a pressure gradient starts to build up, which opposes the action of gravity. This explains the acceleration phase seen in the observations and explains why the condensation does not follow a freefall motion. The same conclusion was reached in Oliver et al. (2014, 2016). As time advances, the pressure gradient keeps increasing until it completely balances the action of gravity, which would correspond to the constant velocity phase. According to the plot, this occurs at a time around t = 200 s. From this time on, it can be seen that the total vertical acceleration is not exactly zero but becomes positive and, in the last steps of the simulation, reaches values of up to 10 m s^{−2}. This means that the total acceleration points upwards and the falling speed is slightly reduced, in agreement with the velocity profiles shown in Fig. 2.
The time it takes the pressure gradient to balance gravity, that is, the duration of the acceleration phase, depends on the density, a dependence that is then reflected in the maximum falling speed reached by the plasma, as shown in Fig. 3. According to Eq. (14), the acceleration related to the pressure gradient is inversely proportional to the density. So, points of larger density require a larger pressure gradient to compensate the acceleration of gravity (see Appendix A for more details on the pressure gradient). The fact that horizontal forces cancel out implies that each vertical slice x = const. does not interact with its neighbours, so that its dynamics is only governed by the density ratio at this value of x. This explains why the central part of the blob falls faster than its sides and the Vshape develops: it happens due to its horizontal variation of density.
The analysis detailed in the paragraphs above focused on specific points in the numerical domain of the simulations. Now, we want to give a broader perspective of the dynamics. For this reason, in Fig. 5 we display 2D maps of the xcomponent (top), zcomponent (middle) and absolute value (bottom) of the magnetic field at three different times of the same simulation discussed before. The purple and green colours represent the strength of the magnetic field and the black lines are contours of density that show the position and shape of the condensation at each time.
Fig. 5. Colour maps of B_{x} (top), B_{z} (middle) and B (bottom) from simulation I5. Black lines are density contours. 

Open with DEXTER 
The left column of Fig. 5 shows that initially the magnetic field does not have a component in the xdirection, but it has only a zcomponent, with a uniform strength of B_{0} = 5 G. As time advances and the plasma starts to fall, perturbations appear in both components. The amplitude of these perturbations is less than 5% of the initial magnetic field strength for the case of B_{z}, and even less for B_{x}. The middle and bottom panels show that the perturbations mainly propagate in the vertical direction and form a column in front of and behind the descending blob, while the rest of the domain remains almost unaffected. However, there is also a small spread along the horizontal direction, which is more evident in the top panels.
A similar behaviour is found in other variables, as shown by Fig. 6, where we represent the plasmaβ (top) and the pressure (bottom). The black lines correspond to density contours. Again the perturbations are seen to propagate mainly along a vertical column whose horizontal extension is given by the size of the blob. Although the plasmaβ is larger at the front of the blob than at its back, its magnitude stays around unity, which means that neither the magnetic field nor the gas dynamics completely dominates the motion of the plasma.
Fig. 6. Colour maps of plasmaβ (top) and pressure (bottom) from simulation I5. Black lines are density contours. 

Open with DEXTER 
These results correspond to the case of small peak density and small magnetic field. However, as mentioned before, we also performed simulations with larger values of these parameters. Since β is proportional to the gas pressure (which is proportional to the density) but is inversely proportional to the square of the magnetic field strength, we expect very different values of β in these simulations. Increasing the density yields an increase of β, meaning that the gas dynamics becomes more important than the magnetic field. However, an increase of the magnetic field produces a much smaller β and, therefore, the motion of the plasma becomes more dominated by the magnetic field.
In Fig. 7, we show snapshots from simulations with a larger density than before. The panel on the left corresponds to the simulation V5 while the panel on the right corresponds to the simulation V20. The colour map represents the plasmaβ and the black lines are the magnetic field lines. The position of the density peak is marked by the horizontal red dashed lines. The left panel, associated with β > 1, shows that the field lines are not completely vertical but they are distorted and have some curvature around the position of the blob. On the other hand, for the case of a stronger magnetic field, the lines are almost completely straight. They are hardly distorted by the presence of the plasma. Moreover, it can be seen that the blob is not at the same height in both simulations although the snapshots correspond to the same time. The plasma has fallen faster in the case with a smaller magnetic field.
Fig. 7. Snapshots from simulations V5 (left) and V20 (right) at time t = 250 s. The colour maps represent the plasmaβ, magnetic field lines are displayed in black and the horizontal red dashed lines mark the height of the density peak of the blob. 

Open with DEXTER 
Now we can fully understand the results summarised by Fig. 3. The effect of the magnetic field is to try to confine the motion of the plasma and the propagation of perturbations into the vertical direction and reduce or prevent the motion and propagation along the horizontal direction. With a strong enough magnetic field, each plasma element falls following the path determined by a field line, with no influence from the adjacent elements. In addition, perturbations such as that of pressure or velocity can only propagate along the same line. Thus, the motion of each vertical cut of the descending blob can be described in terms of the 1D dynamics analysed in Oliver et al. (2014, 2016). On the other hand, with a weaker B_{0}, the dynamics is not completely dominated by the magnetic field and the perturbations can spread along the horizontal direction. For instance, the variation of pressure below the falling blob extends along a larger horizontal length, which leads to a smaller value of the pressure gradient in each vertical slice. A smaller magnitude of the force opposing the action of gravity results in a longer acceleration phase and a larger maximum falling speed. Thus, we have found the explanation for the outlier points in Fig. 3.
3.3. KelvinHelmholtz instability
All the simulations analysed in the previous section corresponded to cases in which the magnetic field has a considerable contribution to the dynamics of the descending condensations. So, up to this point, we left unexplored the regime in which the influence of the magnetic field is much smaller or can be completely neglected. Now we want to investigate this regime. In this section, we start by analysing the results of a simulation in which no magnetic field has been included in the equilibrium state. Therefore, this simulation is purely hydrodynamic. Then we turn back to the magneothydrodynamic case but with a weak magnetic field.
Figure 8 displays several snapshots of the simulation V0. Except for the magnetic field, it has the same parameters as those used to generate Fig. 7. The top and bottom panels show the temporal evolution of density and pressure, respectively. Looking at the top panels, it can be seen that some ripples appear at the sides of the blob as it falls. As time advances, these ripples turn into vortexes, features with smaller scales appear and the shape of the blob becomes much more distorted than in the cases that took into account the presence of the magnetic field. This behaviour reminds of the one expected when the conditions for the development of the KelvinHelmholtz instability (KHI) are fulfilled. Another remarkable fact can be seen in the bottom panels: in contrast with the magnetohydrodynamic simulations analysed before, here the perturbations of pressure clearly propagate in all directions and are not confined to a vertical column.
Fig. 8. Colour maps of density (top) and pressure (bottom) from simulation V0. An animation of this figure is available online. 

Open with DEXTER 
The KHI is a phenomenon that occurs when there is a shear flow at the interface between two different fluids. However, the existence of the shear flow does not always guarantee the emergence of this instability. There are some physical effects, such as gravity or the surface tension of the fluids, that prevent its development or greatly reduce its growth rate (Chandrasekhar 1961).
The most simple way to study the KHI is to consider two fluids of different but uniform densities, ρ_{1} and ρ_{2}, and a shear flow represented by a discontinuous jump at the interface between the two fluids. In this configuration, according to Chandrasekhar (1961), the presence of a magnetic field oriented in the same direction as the shear flow has an stabilising effect on fully ionised plasmas. Then, the plasma only becomes unstable if the velocity of the flow is larger than a certain threshold. The value of this threshold is given by
which assumes that the magnetic field felt by each fluid may be different.
In the simulations analysed in Sects. 3.1 and 3.2, we found no signs of the development of the KHI. The smallest value of the magnetic field strength used in that series of simulations was B_{0} = 5 G, which may produce a large enough velocity threshold.
To compute an approximate expression for the velocity threshold associated with that B_{0}, we set B_{1} = B_{2} = B_{0}, because we have a uniform magnetic field over all the domain, and ρ_{1} = ρ_{eq} and ρ_{2} = ρ_{b0}. We also take into account that the equilibrium density varies with height but we assume that each height is an independent layer not affected by the layers above or below. This means that we neglect the effects related to the vertical gradient of density. Thus, the velocity threshold can be expressed as
Moreover, due to the density of the blob being much larger than the density of the surrounding corona, we can obtain the following approximate expression:
where is the Alfvén speed of the corona. This expression tells us that superAlfvénic speeds are needed for the development of the KHI. Note that in this approximation the velocity threshold does not depend on the density of the blob.
The work of Chandrasekhar (1961) also provides the growth rate of the instability, which in our case can be written as
where k_{z} is the wavenumber associated with the perturbation that may become unstable and ΔU is the shear flow velocity.
Now to check if the KHI may develop during the fall of the plasma, we can compare the falling speeds of the blob with the velocity thresholds provided by the previous equations. For this comparison we defined V_{fall} as the zcomponent of the velocity at the density maximum of the blob, whose time dependence is shown for a particular simulation at the right panel of Fig. 2. Since the plasma surrounding the blob remains approximately static, the shear flow velocity is equal to V_{fall} and the KHI only emerges if V_{fall} > V_{KHI}.
In Fig. 9, we represent V_{KHI} as a function of height for two different magnetic field strengths: the red line corresponds to B_{0} = 5 G and the black line corresponds to B_{0} = 1 G. We also display the values of V_{fall} obtained from two different simulations with a large density blob: red crosses and black diamonds are the results from simulations X5 and X1, respectively.
Fig. 9. V_{fall} and V_{KHI} as functions of height. Lines represent V_{KHI} and symbols represent the falling speed of blobs in simulations with ρ_{b0} = 10^{−9} kg m^{−3}. Red colour corresponds to the case with B_{0} = 5 G while black corresponds to B_{0} = 1 G. 

Open with DEXTER 
In the first place, we pay attention to the case with B_{0} = 5 G. We can see that, although the blob reaches falling speeds of ∼120 km s^{−1}, V_{fall} always stays well below V_{KHI}, which has a minimum value of ∼280 km s^{−1}. This means that the KHI does not appear under these conditions. On the contrary, for B_{0} = 1 G, the velocity threshold is small enough so that at a certain height V_{fall} becomes larger than it. In this simulation, the condition for the development of the KHI is fulfilled once the blob falls to a height of ∼42 Mm. Then, the falling speeds keep increasing while the threshold reduces, which leads to larger growth rates of the instability.
The direct comparison of the results from simulations X5 and X1 without taking into account their relation with the KHI thresholds may give the impression that both cases have a similar behaviour and that the KHI does not actually develop in the latter. The dependence of each respective V_{fall} on height follows an almost identical trend. However, examining them in more detail, it can be checked that the condensation in X1 falls slightly faster than in X5. Moreover, we remind that V_{fall} only provides information about the motion along the x = 0 axis. Therefore, we cannot infer the evolution of the whole condensation by only looking into this parameter.
To check if the prediction from Fig. 9 that the KHI should appear in simulation X1 is correct, we need to examine its 2D behaviour. In Fig. 10 we display colour maps of density (top panels) and plasmaβ (bottom panels) at different snapshots of this simulation. The bottom panels also include a representation of the magnetic field lines. The panels for times t = 100 s and t = 200 s show the formation of the Vwings behind the blob. However, they do not emerge from the sides, as in Fig. 1, but from a more centered position. This is a consequence of the magnetic field being very weak, so it cannot prevent the horizontal motions as in the cases discussed in Sect. 3.1. The bottom panels show how the gas dynamics dominates and the magnetic field lines become very distorted.
Fig. 10. Colour maps of density (top) and plasmaβ (bottom) from simulation X1. Magnetic field lines are plotted in black colour in the bottom panel. An animation of this figure is available online. 

Open with DEXTER 
Regarding the KHI, no sign of the ripples typically associated with this instability is seen in the leftmost panels of Fig. 10. Nevertheless, the panels in the center, in which the condensation has already fallen below the height z = 40 Mm, display a larger deformation of the blob. And at later times, the ripples and vortexes of the instability can be clearly seen. According to Fig. 9, the KHI may start to develop at heights below z ∼ 42 Mm, but since the growth rate of the instability depends on the difference between the falling speed and the velocity threshold, as shown by Eq. (18), it would take some time until the perturbations on the plasma are large enough that they greatly alter the shape of the blob. Therefore, the results presented in Fig. 10 seem to agree with the prediction from Fig. 9.
It must be noted that Fig. 9 displays the results of two of the most unstable cases we studied with our simulations. They have the largest blob densities of all the series, which lead to the largest falling speeds (see Fig. 3), and the smallest magnetic field strengths, that correspond to the smallest velocity thresholds as shown by Eq. (15).
In addition, this analysis assumed a simplified model of the KHI, which considers that there is a discontinuous change of density and velocity between the two media and corresponds to the most unstable configuration. However, in our simulations we do not have that discontinuity but a continuous variation of density, which increases the stability of the system.
Therefore, from our analysis we conclude that, due to the stabilising effect of the magnetic field, the coronal rain blobs are not expected to develop KH instabilities during their falling motion.
4. Conclusions
In the present work, we used a 2D numerical model to investigate the dynamics of descending condensations in the solar corona. We assumed that the plasma is fully ionised and did not take into account nonideal effects (e.g. the ambipolar diffusion due to the presence of neutral particles or thermal conduction). We represented the solar corona as an isothermal vertically stratified atmosphere with no boundaries, which means that we neglected the effects related to the presence of the chromosphere and photosphere at the bottom boundary. Since we wanted to focus only on the descending phase of the coronal rain cycle, we did not include the formation process of the blobs in our study. Therefore, we assumed that the condensations are already formed at the initial time of our simulations. We represented them by a density enhancement described by a twodimensional Gaussian profile. In addition, we took into account the presence of a uniform vertical magnetic field. In such a configuration, the condensation is not in equilibrium and starts to fall down under the action of gravity. In the present work, we analysed the features of this falling motion.
First, we described the general evolution of the descending condensation, which does not behave as a freefalling object. The denser central parts fall faster than the lighter sides, which leads to the deformation of the initial Gaussian profile and the appearance of the Vshape displayed in Fig. 1. Then, we studied the kinematics of several selected points along the horizontal axis of the condensation. Our numerical results show a good agreement with the timedistance diagrams elaborated from different observations of coronal rain events (Schrijver 2001; De Groof et al. 2004; Antolin et al. 2010). Those diagrams reveal that the motion of the descending plasma can be separated in two phases: a first phase of acceleration and a second phase of almost constant falling speed.
To understand the reasons behind the two phases of the descending motion and the development of the Vshape, we analysed the different forces that determine the evolution of the condensation. We confirmed the findings of Oliver et al. (2014) that the reason for the slower than freefall descent is the appearance of a pressure gradient that opposes the action of gravity. As time advances, the pressure gradient builds up and, in some cases, is able to fully balance the acceleration due to gravity, which explains the existence of the constant speed phase. The magnitude of the pressure gradient required to balance gravity depends on the density of the plasma. Larger densities need larger pressure gradients, which leads to longer accelerations phases and larger maximum falling speeds. Therefore, the shape of a condensation with a horizontal variation of density will become distorted as the denser parts fall faster than the lighter ones. This explains the Vshape found in our simulations.
The maximum falling speeds found in our study are below 150 km s^{−1}, which is in good agreement with the observations (Wiik et al. 1996; Antolin & Rouppe van der Voort 2012). Moreover, our 2D results show a remarkable agreement with the 1D results from Oliver et al. (2014, 2016), except for the cases with a weak magnetic field. For weak fields, the maximum falling speeds are larger than in the 1D model. The analysis of the accelerations due to the presence of the magnetic field shows that the effect of the Lorentz term in the vertical direction is negligible in comparison with those of the pressure gradient or gravity. On the other hand, it becomes fundamental in the motions along the horizontal direction. It almost perfectly balances the effect of the pressure gradient, leading to a negligible acceleration along the xdirection. As a consequence, the evolution of each vertical slice of the condensation becomes almost independent of its neighbours and its evolution can be described by the 1D model. Another effect of the presence of a strong enough magnetic field is that the propagation of the perturbations becomes confined to the vertical direction. For instance, the perturbation of pressure cannot spread along the xdirection, which results in a larger vertical pressure gradient and, consequently, a smaller falling speed.
We also simulated cases in which the magnetic field strength is very small or there is no magnetic field at all. In the purely hydrodynamic configuration, that is, when the effect of the magnetic field is neglected, the evolution of the descending blob is very different from that described in the previous paragraphs. As the plasma falls, ripples and vortexes related to the KHI appear and the shape of the blob is greatly distorted. In these circumstances, the 2D evolution of the condensation cannot longer be treated as a combination of independent vertical slices that follow the 1D dynamics.
It is known that the presence of a magnetic field oriented along the direction of the shear flow has a stabilising effect on the KHI (Chandrasekhar 1961). So we performed a study to determine under which conditions the KHI may develop during coronal rain events. Using a simple analytical model, we compared the theoretical velocity thresholds for the emergence of the instability with the falling speeds computed from several simulations. Such a comparison shows that, for the majority of cases, the maximum falling speeds are well below the superAlfvénic velocity thresholds. Only for very small values of the magnetic field strength are the falling speeds comparable to or larger than the theoretical thresholds. This means that under the typical conditions in which coronal rain occurs, the KHI is not expected to develop.
However, some circumstances that might modify the previous conclusion should be mentioned here. We assumed a vertical uniform magnetic field. In a more realistic model, B_{0} should have a decreasing dependence on height, which may lead to smaller values in the upper layers of the solar corona than the ones used in our simulations, thus reducing the instability velocity thresholds. So, it may be possible for the KHI to appear in blobs falling from higher layers of the corona. Furthermore, it is known that the presence of neutrals in a plasma can have a dramatic effect on the development of the instability (Watson et al. 2004). The neutral component of the plasma is always unstable, even for subAlfvénic speeds, although the collisional coupling with ions reduces the growth rate of the instability (Soler et al. 2012). The influence of the dependence of the magnetic field on height and the presence of a neutral component in the condensation should be studied in the future.
We should also remark that the Vshape of the blob found in our simulations does not agree with what has been actually detected. The picture we get from observations is that of a condensation that is elongated along the direction of its falling motion, with no wings at its sides (Wiik et al. 1996; Antolin & Rouppe van der Voort 2012). This discrepancy may be explained by assuming that the catastrophic cooling leading to the formation of the blob is a continuous process that is still active while the plasma starts to fall. In this way, the coronal plasma can continue to condense and fall along the same path. On the contrary, our simulations start when the blob has already been completely formed and the cooling has already stopped. It is also possible that the wings are actually present in the coronal rain events but the instruments cannot detect them because their emission is much smaller than that of the central part of the blob, due to the great difference in density. Another alternative explanation comes from extrapolating our 2D results to a 3D scenario: if we assume that the blobs have axial symmetry, the 2D Vshape transforms into a funnel whose density decreases with height. This funnelshape would have a greater resemblance to what the observations show. In any case, more complex simulations including the formation process of the condensations and forward modeling tools should be used to solve the discrepancy.
Finally, as mentioned previously, this is the first paper of a series. Hence, we left out from the present study some effects that we would like to include in forthcoming works. For instance, we plan to perform simulations that take into account the influence of the chromosphere. The presence of this much denser layer can contribute to slowing down the blob descent and even produce a rebounding of the plasma, as shown by Mackay & Galsgaard (2001). In addition, here we only considered a completely vertical magnetic field but observations show that the plasma in coronal rain moves along curved paths, following the magnetic structure of the coronal loops (Tripathi et al. 2006, 2007; Antolin & Rouppe van der Voort 2012; Antolin et al. 2012). We want to study the evolution of the condensation when the magnetic field lines are inclined or follow a curved path. In such configurations, we expect to find smaller falling speeds due to the effect of the projection of gravity along the path (Petralia et al. 2016; Xia et al. 2017). The influence of the inclination or the curvature of the field lines may be enhanced when the condensation is assumed to be partially ionised. The neutral component of the plasma is not directly impacted by the guiding effect of the magnetic field and can drag the ionised component through the magnetic field lines. However, it is also possible that the collisional coupling between ions and neutrals is strong enough to prevent the latter from falling vertically (Terradas et al. 2015). It would prove interesting to investigate how the dynamics of partially ionised blobs depend on the ionisation degree in addition to the inclination or curvature of the field lines.
Acknowledgments
D.M., E.K., and M.C. acknowledge support from the European Research Council through the Consolidator Grant ERC2017CoG771310PI2FA. R. O. acknowledges support from the Spanish Ministry of Economy, Industry and Competitiveness through the grant AYA201785465P. R.O. also acknowledges the travel support received from the International Space Science Institute (Bern, Switzerland) as well as discussions with members of the ISSI team on “Observed multiscale variability of coronal loops as a probe of coronal heating”, led by C. Froment and P. Antolin. The authors thankfully acknowledge the technical expertise and assistance provided by the Spanish Supercomputing Network (Red Espanola de Supercomputacion), as well as the computer resources used: the LaPalma Supercomputer, located at the Instituto de Astrofisica de Canarias.
References
 Ahn, K., Chae, J., Cho, K.S., et al. 2014, Sol. Phys., 289, 4117 [NASA ADS] [CrossRef] [Google Scholar]
 Antiochos, S. K. 1980, ApJ, 241, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Antiochos, S. K., & Klimchuk, J. A. 1991, ApJ, 378, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Antiochos, S. K., MacNeice, P. J., Spicer, D. S., & Klimchuk, J. A. 1999, ApJ, 512, 985 [NASA ADS] [CrossRef] [Google Scholar]
 Antolin, P., & Rouppe van der Voort, L. 2012, ApJ, 745, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Antolin, P., & Verwichte, E. 2011, ApJ, 736, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Antolin, P., Shibata, K., & Vissers, G. 2010, ApJ, 716, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Antolin, P., Vissers, G., & Rouppe van der Voort, L. 2012, Sol. Phys., 280, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Antolin, P., Vissers, G., Pereira, T. M. D., Rouppe van der Voort, L., & Scullion, E. 2015, ApJ, 806, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Auchère, F., Froment, C., Soubrié, E., et al. 2018, ApJ, 853, 176 [NASA ADS] [CrossRef] [Google Scholar]
 Beckers, J. M. 1962, Aust. J. Phys., 15, 327 [NASA ADS] [CrossRef] [Google Scholar]
 Berenger, J.P. 1994, J. Comput. Phys., 114, 185 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Berenger, J.P. 1996, J. Comput. Phys., 127, 363 [NASA ADS] [CrossRef] [Google Scholar]
 Bruzek, A. 1969, Sol. Phys., 8, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Chae, J. 2010, ApJ, 714, 618 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Oxford University Press) [Google Scholar]
 De Groof, A., Berghmans, D., van DrielGesztelyi, L., & Poedts, S. 2004, A&A, 415, 1141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Engvold, O., Jensen, E., & Andersen, B. N. 1979, Sol. Phys., 62, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, X., Xia, C., & Keppens, R. 2013, ApJ, 771, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, X., Xia, C., Keppens, R., & Van Doorsselaere, T. 2015, ApJ, 807, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Field, G. B. 1965, ApJ, 142, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Foukal, P. 1978, ApJ, 223, 1046 [NASA ADS] [CrossRef] [Google Scholar]
 Froment, C., Auchère, F., Bocchialini, K., et al. 2015, ApJ, 807, 158 [NASA ADS] [CrossRef] [Google Scholar]
 Froment, C., Auchère, F., Aulanier, G., et al. 2017, ApJ, 835, 272 [NASA ADS] [CrossRef] [Google Scholar]
 GonzálezMorales, P. A., Khomenko, E., Downes, T. P., & de Vicente, A. 2018, A&A, 615, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harra, L. K., Matthews, S. A., Long, D. M., Doschek, G. A., & De Pontieu, B. 2014, ApJ, 792, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Heinzel, P., Schmieder, B., & Mein, P. 1992, Sol. Phys., 139, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Karpen, J. T., Antiochos, S. K., Hohensee, M., Klimchuk, J. A., & MacNeice, P. J. 2001, ApJ, 553, L85 [NASA ADS] [CrossRef] [Google Scholar]
 Khomenko, E., & Collados, M. 2006, ApJ, 653, 739 [NASA ADS] [CrossRef] [Google Scholar]
 Kuin, N. P. M., & Martens, P. C. H. 1982, A&A, 108, L1 [NASA ADS] [Google Scholar]
 Liu, W., De Pontieu, B., Vial, J.C., et al. 2015, ApJ, 803, 85 [NASA ADS] [CrossRef] [Google Scholar]
 Loughhead, R. E., & Bray, R. J. 1984, ApJ, 283, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Mackay, D. H., & Galsgaard, K. 2001, Sol. Phys., 198, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Martens, P. C. H., & Kuin, N. P. M. 1983, A&A, 123, 216 [NASA ADS] [Google Scholar]
 Mason, E. I., Antiochos, S. K., & Viall, N. M. 2019, ApJ, 874, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Moschou, S. P., Keppens, R., Xia, C., & Fang, X. 2015, AdSpR, 56, 2738 [NASA ADS] [CrossRef] [Google Scholar]
 Müller, D. A. N., Peter, H., & Hansteen, V. H. 2004, A&A, 424, 289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Müller, D. A. N., De Groof, A., Hansteen, V. H., & Peter, H. 2005, A&A, 436, 1067 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murawski, K., Zaqarashvili, T. V., & Nakariakov, V. M. 2011, A&A, 533, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oliver, R., Soler, R., Terradas, J., Zaqarashvili, T. V., & Khodachenko, M. L. 2014, ApJ, 784, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Oliver, R., Soler, R., Terradas, J., & Zaqarashvili, T. V. 2016, ApJ, 818, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Parchevsky, K. V., & Kosovichev, A. G. 2009, ApJ, 694, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Parker, E. N. 1953, ApJ, 117, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Petralia, A., Reale, F., Orlando, S., & Testa, P. 2016, ApJ, 832, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Schrijver, C. J. 2001, Sol. Phys., 198, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Scullion, E., Rouppe van der Voort, L., Wedemeyer, S., & Antolin, P. 2014, ApJ, 797, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Scullion, E., Rouppe van der Voort, L., Antolin, P., et al. 2016, ApJ, 833, 184 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Díaz, A. J., Ballester, J. L., & Goossens, M. 2012, ApJ, 749, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., & Ofman, L. 2004, ApJ, 610, 523 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., Soler, R., Oliver, R., & Ballester, J. L. 2015, ApJ, 802, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Tripathi, D., Solanki, S. K., Schwenn, R., et al. 2006, A&A, 449, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tripathi, D., Solanki, S. K., Mason, H. E., & Webb, D. F. 2007, A&A, 472, 633 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Watson, C., Zweibel, E. G., Heitsch, F., & Churchwell, E. 2004, ApJ, 608, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Wiik, J. E., Schmieder, B., Heinzel, P., & Roudier, T. 1996, Sol. Phys., 166, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Xia, C., Keppens, R., & Fang, X. 2017, A&A, 603, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Xu, A. A. 1987, in Solar Maximum Analysis, eds. V. E. Stepanov, & V. N. Obridko, 147 [Google Scholar]
 Xue, Z., Yan, X., Qu, Z., & Zhao, L. 2014, New Astron., 26, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, J., & Li, L.P. 2009, Res. Astron. Astrophys., 9, 1368 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Profile of the pressure gradient
In Sect. 3.2, we show that the inertial and the Lorentz terms of the momentum equation do not have an important effect on the vertical motion of the descending blob. Therefore, neglecting those terms, Eq. (14) can be rewritten as
After the initial acceleration phase, the plasma keeps descending with a constant velocity, so ∂V_{z}/∂t = 0 and ∂P/∂z = −ρg.
The initial twodimensional profile of the density enhancement is given by Eq. (12), but if we only take now into account its dependence with height, we have that
Thus, the profile of the perturbation of pressure that balances the action of gravity is given by
which results in
where erf(z) is the Gauss error function. Note that here we used z_{peak} instead of z_{0} to take into account the actual position of the blob during its fall.
Then, to obtain the total pressure required for a constantspeed phase, we have to combine the expression from Eq. (A.4) with the one corresponding to the pressure in the equilibrium state, given by Eq. (10). Hence,
In Fig. A.1, we compare the analytical profile derived above with vertical cuts of the pressure at x = 0 from two different simulations in which the magnetic field strength is B_{0} = 5 G. Both panels show the results at a time t = 500 s, but the top panel corresponds to a simulation with ρ_{b0} = 10^{−10} kg m^{−3}, while in the bottom panel the density is 10 times larger, ρ_{b0} = 10^{−9} kg m^{−3}. The black symbols represent the data from the simulations and the red lines represent the analytical results. The dashed vertical lines mark the position of the condensation.
For the case of smaller density, the numerical and the analytical results overlap, meaning that the condition for the constantspeed phase is fulfilled. The magnitude of the pressure gradient is enough to balance the action of gravity. On the contrary, the bottom panel shows that a larger density requires a larger pressure gradient to reach the constantspeed phase. The pressure gradient computed from the simulation is larger than in the previous case but it is still not enough to completely compensate the acceleration of gravity. Therefore, the blob is still gaining speed at this time of the simulation. Consequently, it has reached a lower height than the lighter condensation of the top panel.
Fig. A.1. Pressure as a function of height. Black symbols represent the numerical results from simulations with B_{0} = 5 G. Red lines show the analytical results given by Eq. (A.5). Top: case with ρ_{b0} = 10^{−10} kg m^{−3}. Bottom: case with ρ_{b0} = 10^{−9} kg m^{−3}. The vertical black dashed lines mark the positions of the condensations. 

Open with DEXTER 
Movies
Movie of Fig. 1 (Access here)
Movie of Fig. 8 (Access here)
Movie of Fig. 10 (Access here)
All Tables
All Figures
Fig. 1. Density at different times of the simulation I5. The horizontal red dashed lines represent the position that would correspond to a freefalling object. The black asterisks show the position of the points whose kinematics is represented in Fig. 2. An animation of this figure is available online. 

Open with DEXTER  
In the text 
Fig. 2. Left panel: position as a function of time of selected points from simulation I5. Right panel: velocity as a function of time. Vertical dashed lines represent the time of the snapshots displayed in Fig. 1. The colours of the lines are related to the density scale shown in the previous figure. 

Open with DEXTER  
In the text 
Fig. 3. Maximum falling speed as a function of density ratio. Black circles correspond to the results from 1D simulations. Triangles, squares and diamonds show the results of 2D simulations with ρ_{b0} = 1, 5 and 10 × 10^{−10} kg m^{−3}, respectively. Red, green and blue colour correspond to the cases with B_{0} = 5, 10 and 20 G, respectively. 

Open with DEXTER  
In the text 
Fig. 4. Dynamics of the density peak in the simulation I5. Top: accelerations in the xdirection. Bottom: acceleration in the zdirection. 

Open with DEXTER  
In the text 
Fig. 5. Colour maps of B_{x} (top), B_{z} (middle) and B (bottom) from simulation I5. Black lines are density contours. 

Open with DEXTER  
In the text 
Fig. 6. Colour maps of plasmaβ (top) and pressure (bottom) from simulation I5. Black lines are density contours. 

Open with DEXTER  
In the text 
Fig. 7. Snapshots from simulations V5 (left) and V20 (right) at time t = 250 s. The colour maps represent the plasmaβ, magnetic field lines are displayed in black and the horizontal red dashed lines mark the height of the density peak of the blob. 

Open with DEXTER  
In the text 
Fig. 8. Colour maps of density (top) and pressure (bottom) from simulation V0. An animation of this figure is available online. 

Open with DEXTER  
In the text 
Fig. 9. V_{fall} and V_{KHI} as functions of height. Lines represent V_{KHI} and symbols represent the falling speed of blobs in simulations with ρ_{b0} = 10^{−9} kg m^{−3}. Red colour corresponds to the case with B_{0} = 5 G while black corresponds to B_{0} = 1 G. 

Open with DEXTER  
In the text 
Fig. 10. Colour maps of density (top) and plasmaβ (bottom) from simulation X1. Magnetic field lines are plotted in black colour in the bottom panel. An animation of this figure is available online. 

Open with DEXTER  
In the text 
Fig. A.1. Pressure as a function of height. Black symbols represent the numerical results from simulations with B_{0} = 5 G. Red lines show the analytical results given by Eq. (A.5). Top: case with ρ_{b0} = 10^{−10} kg m^{−3}. Bottom: case with ρ_{b0} = 10^{−9} kg m^{−3}. The vertical black dashed lines mark the positions of the condensations. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.