Dynamics of coronal rain and descending plasma blobs in solar prominences: I. Fully ionised case.
Observations of active regions and limb prominences often show cold, dense blobs descending with an acceleration smaller than that of free fall. The dynamics of these condensations falling in the solar corona is investigated in this paper using a simple fully ionised plasma model. We find that the presence of a heavy condensation gives rise to a dynamical rearrangement of the coronal pressure that results in the formation of a large pressure gradient that opposes gravity. Eventually this pressure gradient becomes so large that the blob acceleration vanishes or even points upwards. Then, the blob descent is characterised by an initial acceleration phase followed by an essentially constant velocity phase. These two stages can be identified in published time-distance diagrams of coronal rain events. Both the duration of the first stage and the velocity attained by the blob increase for larger values of the ratio of blob to coronal density, for larger blob mass, and for smaller coronal temperature. Dense blobs are characterised by a detectable density growth (up to 60% in our calculations) and by a steepening of the density in their lower part, that could lead to the formation of a shock. They also emit sound waves that could be detected as small intensity changes with periods of the order of 100 s and lasting between a few and about ten periods. Finally, the curvature of the falling path is only relevant when a very dense blob falls along inclined magnetic field lines.
Coronal rain is a well-known example of cool material (with chromospheric to transition region temperature) falling down in the much hotter solar corona at heights of some tens of Mm. Coronal rain consists of cold, elongated blobs condensing at an active region loop and then falling down along its leg, with the loop magnetic field guiding the blobs descent (Beckers, 1962). The magnetic field structure has a much longer lifetime than the time taken for a condensation to fall (Beckers, 1962) and a given loop often supports a few falling blobs. For example, Antolin et al. (2010) tracked motions along 30 loops and detected more than one hundred descending condensations in a time span of 71 minutes. Moreover, Antolin & Rouppe van der Voort (2012) used an 84-minute long observational sequence to study the trajectories of 2552 condensations along 242 different paths. These two examples also reveal that coronal rain events take place continuously in an active region.
Regarding the temperature of these condensations, coronal rain is usually observed in cold “chromospheric” lines, both of neutral atoms (Ly and H) and of ionised elements (Ca ii and He ii). In addition, De Groof et al. (2004, 2005) analyzed several coronal rain blobs by combining simultaneous images in H, He ii 304 Å, and 171 Å and concluded that there is no hot “coronal” counterpart to this falling material because the blobs only showed up in H and He ii, but not in 171 Å. The work by De Groof et al. (2005) also provides evidence of coupling between the neutral and ionised blobs (detected through the H and He ii emission, respectively) because their kinematics was similar.
The first observational investigations of coronal rain gave falling speeds ranging from 40 to 60 km s (Beckers, 1962; Bruzek, 1969; Foukal, 1978; Engvold et al., 1979), but no details about the blobs kinematics were given in these works. For this purpose a monitoring of their position or speed as a function of height or time is necessary. As soon as this kind of information was gathered, it became clear that the cool knots acceleration was smaller than that of free fall (Loughhead & Bray, 1984; Xu, 1987; Heinzel et al., 1992). A detailed tracking of many coronal rain condensations provides the following picture of their kinematics (Wiik et al., 1996; De Groof et al., 2004, 2005; Antolin et al., 2010; Antolin & Rouppe van der Voort, 2012; Antolin et al., 2012): blobs display a free fall behavior only at the loop top, with velocities around 30–50 km s, whereas further down the loop leg the downward velocity becomes considerably smaller than that expected in a free fall regime and reaches 100–150 km s. Another description of the blobs kinematics is also given by Zhang & Li (2009), who also obtained moderate accelerations and found that, in general, the blobs speed increased during the first few Mm of their trajectory, then remained roughly constant for Mm, and finally decreased abruptly when the blobs were between 5 and 10 Mm above the solar surface. Two types of descending blobs were found in this study: a “fast” and a “slow” group, whose average falling speeds were 72 km s and 37 km s. Some authors have even computed the acceleration of these falling condensations and values smaller than the acceleration of gravity have always been derived: using data from a large ensemble of coronal rain blobs Antolin et al. (2010), Antolin & Rouppe van der Voort (2012), and Antolin et al. (2012) obtained average accelerations of 100–150 m s, 83.5 m s, and 137 m s, respectively. Furthermore, Schrijver (2001) analyzed a set of 45 coronal rain blobs and concluded that their acceleration was also smaller than that of free fall and that it had a very weak dependence on the loop inclination: it had an average value of 85 m s for loops inclined less than 30 from the vertical direction, 70 m s for loops with an average inclination between 30 and 47, and 71 m s for a highly inclined loop with an angle of 72 with the radial direction.
This observational background has often been acquired with the help of either spectra or narrow band images. These techniques have the drawback that they only provide the line-of-sight projection of the blob velocity or its position on the plane of the sky, so that its kinematics is not well described. In a few works, however, spectra and images have been combined together to derive the full velocity vector. In the analysis of Wiik et al. (1996), this was done from a combination of the loop shape reconstruction with the blobs Doppler shift and displacement in the plane of the sky. In the case of Antolin & Rouppe van der Voort (2012); Antolin et al. (2012), the velocity vector of descending blobs was obtained by combining their apparent motion and their Doppler velocity. The important result here is that the kinematics of coronal rain blobs, namely their typical speeds and acceleration and the variation of velocity with height, do not change when the full velocity vector is considered.
Several attempts have been made to understand the formation and dynamics of coronal rain material. Numerical simulations by Müller et al. (2004, 2005, see also references therein) indicate that a constant loop heating that decreases exponentially with height results in a thermal instability that leads to catastrophic cooling at the loop apex (see also Schrijver, 2001). A cold condensation then forms and falls down along the loop legs. Additional refinements were introduced by Antolin et al. (2010) by taking into account three different effects, namely the variation of the loop cross sectional area, a loop heating concentrated at the base of the structure and that changes episodically in time (so as to simulate nanoflare heating events), and Alfvén wave dissipation at the loop feet. An important conclusion of this work is that the structure and dynamics of the coronal rain blobs are more sensitive to the pressure variations arising from catastrophic cooling than to gravity itself. This is in agreement with Schrijver (2001), who suggested that the internal pressure evolution of the loops, rather than gravity, determines the speed of condensations. Fang et al. (2013) conducted two-dimensional numerical simulations of a sheared coronal arcade with a constant heating localised near the chromosphere. About 4000 coronal rain condensations formed during 80 minutes, although the authors found that only one fourth of them would be detectable at the best present observational resolution. This suggests that a large number of coronal rain events can go unnoticed in current observations. Blob descending speeds obtained in this work range from a few to 60 km s and in all cases are consistent with accelerations smaller than 180 m s. Finally, Murawski et al. (2011) put forward a completely different idea for the formation of coronal rain blobs. They showed that a pressure pulse at a magnetic null point can excite an entropy mode that creates a dense blob. Exhaustive information about coronal rain observations and numerical simulations leading to the formation of cold condensations in the corona can be found in Antolin et al. (2010); Antolin & Rouppe van der Voort (2012).
Another type of cool descending blobs is also present in limb prominences, a radically different environment from active region coronal loops in which coronal rain is observed. One of the first detailed reports of this phenomenon is that of Engvold (1976), who used time series of H filtergrams of 24 prominences located away from sunspots. He reported on bright knots that frequently appeared stationary for 2–10 min at the location where they were first detected and later descended with a speed in the plane of the sky between 15 and 35 km s. These events seem to be different from the thin vertical prominence threads in which material flows down at small speeds, of the order of 10 km s (Engvold, 1981; Berger et al., 2008, 2010).
Instruments onboard Hinode and SDO have been recently used to provide detailed information about bright descending knots in prominences. With the help of H and Ca ii H images obtained with SOT on Hinode, Chae (2010) tracked moving emission features in a large hedgerow prominence. An important result of this work is that blobs detectable in H images could always be detected in their Ca ii H counterparts, their descending motions showing no perceptible difference. As Chae (2010) pointed out, this implies that Ca ii and neutral H are dynamically coupled. Regarding the knots kinematics, this author found that their paths are almost vertical and measured speeds in the plane-of-the-sky between 10 and 30 km s, with an acceleration that is not constant in time, although it is always below 100 m s. Plasma knots ejected in a limb prominence and then undergoing ballistic motion were studied by Hillier et al. (2011). The kinematics of the descending motions were not investigated in the same detail as in Chae (2010), but analogous results were obtained: cool blobs had downward velocities smaller than 60 km s and their acceleration did not remain constant during their motion.
Liu et al. (2012) described the formation of a prominence by the condensation of coronal mass and concluded that the prominence under study was not a static object and that at a given time during the formation process, its total mass came from the balance between mass condensation and drainage through vertical downflows. Most of the condensed mass drained down along vertical paths, the mass drainage rate being about 96% the condensation rate. These downflows occured in the form of cool mass blobs descending with a roughly constant acceleration that is much smaller than that of free fall. A total of 874 downflowing trajectories were studied by Liu et al. (2012), who concluded that a typical downflowing event lasted between a few minutes and half an hour. The descending mass blob started at a height between 20 Mm and 40 Mm and traveled a distance in the range 10–30 Mm. The derived accelerations have a mean value of 46 m s and cover a wide range from 10 to 200 m s. Furthermore, Liu et al. (2012) measured the speed of blobs 25” above the surface and values around 30 km s were found.
Explaining vertical motions in prominences is a challenging task. Most observational determinations of the magnetic field vector in quiescent prominences support the idea that it is roughly horizontal inside these objects; see, for example, Leroy et al. (1983, 1984) and also the reviews by López Ariste & Aulanier (2007) and Mackay et al. (2010). And because very few observations (such as that by Merenda et al., 2006) contradict this conventional belief, most theoretical efforts to model the prominence internal magnetic configuration are based on horizontal or dipped magnetic fields. Under the conditions of the frozen-free theorem, ions moving vertically must deform the roughly horizontal magnetic field lines, which means the generation of electric currents and the creation of a Lorentz force to compensate the action of ions; the situation is different for neutrals, that can slip across the magnetic field. The different dynamics of ions and neutrals in a steady-state prominence was investigated by Gilbert et al. (2002), who derived falling velocities of the order of 4 m s and 80 m s for neutral H and He, respectively. The vertical drainage of neutrals (especially He) was proposed as the cause of the He deficit in the upper parts of several filaments (Gilbert et al., 2002, 2007). The descending speeds obtained in these steady-state calculations are thus much smaller than those of falling cool blobs in prominences.
Haerendel & Berger (2011) made order-of-magnitude calculations to understand the motions of dense blobs and vertical threads in prominences. They proposed that dense plasma elements condense at the top of prominences and then sink into the mainly horizontal magnetic field. These authors postulated that the plasma next becomes squeezed vertically by the action of gravity and so a current sheet forms just above them, where magnetic field reconnection subsequently takes place. The consequence of this process is that the plasma element becomes detached from the prominence magnetic field and quickly accelerates under the action of gravity. As the plasma knot makes its way through the ambient magnetic field, it looses momentum and energy by the emission of Alfvén waves, and so a more or less steady vertical velocity is attained. Magnetic reconnection has also been invoked as an essential process in falling prominence condensations by Chae (2010). A deeper insight into this phenomenon was obtained by Low et al. (2012a), who represented the prominence as a Kippenhahn-Schlütter plasma slab in thermal balance between optically thin radiation, heating, and conduction. Low et al. (2012a) suggested that the frozen-in condition may recurrently break down and so a downward resistive flow across the local horizontal magnetic field is produced, which would be detected as cold falling material in a prominence. Low et al. (2012a) predicted that these prominence structures may have cores much denser and less ionised than usually considered. Additional support for this hypothesis was provided by Low et al. (2012b), who argued that high electrical conductivity leads to the spontaneous formation and resistive dissipation of discrete currents, and that this may lead to the dynamical situation observed in the interior of some quiescent prominences.
In this series of papers we use a very simple model corona to address several issues associated to cool descending knots both in active region loops and in quiescent prominences. First of all, we want to understand the interplay between the various forces (especially gravity and the pressure gradient, but also the inertial term) in the dynamics of these features. We deliberately ignore the formation process of these condensations and so discard thermal effects in our numerical simulations. Second, we try to understand what conditions lead to falling blobs with either constant velocity or constant acceleration. Two scenarios are possible: one in which a blob with zero initial speed is left free under the action of the various forces, and a second one in which a continuous plasma injection takes place at a certain height in the corona so that a mass condensation forms and starts to fall. While the first situation is similar to most coronal rain and descending prominence knot events, the second one resembles the observations of Liu et al. (2012) with a continuous condensation and downward plasma drainage. Finally, we want to assess the strength of the coupling through friction between neutrals and ions in a descending blob. Although these two components can in principle move independently (as in the steady-state model of Gilbert et al., 2002), observations made simultaneously in suitable spectral lines show that their kinematics are the same (De Groof et al., 2005; Chae, 2010).
This paper is organised as follows: in § 2 the equations that describe the motions of dense, falling plasma blobs are presented. Here we concentrate in the fully ionised plasma case, while the dynamics of a partially ionised plasma blob are left for future work. In addition, only the first scenario (dynamics of a fully formed blob) is studied. § 4 and 5 respectively contain the results for a vertical and a circular descending path, and conclusions are drawn in § 6.
2.1 Governing equations
To study the temporal evolution of a fully ionised blob falling under the influence of gravity and pressure forces, only the dynamics in the vertical direction is taken into account. Hence, the horizontal structuring is ignored to facilitate our study. Moreover, the effect of magnetic fields is ignored. This corresponds either to a non-existent magnetic field or to motions along vertical magnetic field lines. The case of low-, high electric conductivity plasma blobs falling along curved magnetic field lines is discussed in § 5; the plasma is the ratio of the plasma pressure to the magnetic pressure.
A Cartesian coordinate system with the -axis pointing in the vertical direction is considered. In addition, the plasma is assumed to be fully ionised hydrogen with a density () that is approximately that of protons and a pressure () equal to that of protons plus electrons. The velocity of protons is denoted by and it is assumed that there are only vertical mass motions, that is, . The expressions describing the temporal evolution of these variables (, , ) are the mass balance, momentum, and energy equations,
with the acceleration of gravity. Joule heating and the divergence of the heat flux in the energy equation have also been neglected. All other terms missing in the energy equation vanish in the present configuration.
We supplement the above system of equations with the expression for the Lagrangian displacement of a plasma element, denoted by ,
Thus, if a particle is initially at position , at time its position is .
Now, we assume that variables depend on and only. Next we insert , , , and into Equations (1)–(4) and substitute . The horizontal components of the momentum equation and of the expression for the Lagrangian displacement are identically zero and so we end up with four non-linear partial differential equations for the four unknowns,
These expressions are supplemented with the ideal gas law for a fully ionised plasma,
with the ideal gas constant; and are Boltzmann’s constant and the proton mass. The reason for the presence of the factor 2 in Equation (9) is that the partial pressures of ions and electrons are summed up under the conditions of same temperature and same number density.
2.2 Static equilibrium
Here we adopt the following exponentially stratified solution to Equation (10)
where the density and pressure at the coronal base () satisfy the ideal gas law,
with the equilibrium temperature, assumed uniform. Furthermore, the vertical scale height () depends on the equilibrium temperature as follows,
2.3 Mass condensation
A fully ionised dense blob is superimposed to the previous static atmosphere and its temporal evolution is then investigated. Since the condensation process is ignored, a fully formed blob is initially at rest at a certain height and is left to evolve under the action of gravity, the pressure gradient, and the inertial term according to Equations (5)–(8). At the blob has a density described by
where is the central mass density of the condensation, is its initial position, and its length in the vertical direction is of the order of .
3 Numerical procedure
To summarise the expressions of § 2, the temporal evolution of the system is described by Equations (5) to (8), where the unknowns are the density, vertical velocity component, pressure, and Lagrangian displacement (, , , and ). The initial state is given by Equations (11) and (14), with the velocity and the plasma displacement both initially set to zero (). We focus our attention on the vertical dynamics of descending matter in the range of heights Mm, with the blob released at an initial height Mm. As we show later, inserting the dense condensation in an otherwise equilibrium corona results in the generation of two sound waves that propagate up and down at the coronal sound speed from the blob initial position. To avoid the influence of wave reflections coming from the boundaries, these are placed well away from our region of interest. Thus, the numerical integration is carried out in the range Mm. The numerical mesh consists of three regions with uniform spacing: in the ranges Mm and Mm the grid spacing is 100 km, whereas in the range Mm it is 10 km. The PDE2D code (Sewell, 2005) with Galerkin’s method and linear finite elements is used for the numerical integration. Finally, during the whole simulation , , and at the bottom and top boundaries are fixed at their initial values, while the -derivative of is set equal to zero there. This choice of boundary conditions provides satisfactory numerical stability.
In § 4.1 we consider a set of “reference” parameter values to study the blob dynamics with a numerical simulation lasting 1000 s. In the initial state the coronal gas has temperature K (so that the vertical scale height is Mm) and density at the coronal base kg m (corresponding to a particle density cm). The base pressure is Pa. The initial mass condensation density is kg m, its initial position is Mm, and the parameter , that determines the blob length, is 0.5 Mm. In § 4.2 we present the results of a similar numerical simulation with a blob ten times denser ( kg m) and with all other parameter values unchanged.
4.1 Blob dynamics
Figure 1a presents several snapshots of the temporal evolution of the density. The mass blob starts to fall under the action of gravity and apparently maintains its shape during its descending motion, although its maximum density increases slightly. From Figure 1a it is straightforward to appreciate that the blob does not fall with the acceleration of gravity: if this was the case, then the blob would reach at s and it would be outside the horizontal range of this figure. Next we draw our attention to the variation of the blob peak density in Figure 1a, which could be caused by the fact that, as the condensation falls, it moves in a denser environment. To check this possibility, in Figure 1b we represent the maximum blob density versus time (solid line) and the same quantity with the coronal density subtracted (dashed line). Both quantities display the same temporal increase and so the enhancement of the blob density is a real effect not caused by the change in the background coronal density with height. This issue is later discussed in more detail. The change of blob height with time is presented in Figure 1c. The blob trajectory is convex in the first s and then becomes roughly linear for s. Hence, the mass condensation accelerates downwards in the initial 200 s but then its acceleration seems to become practically zero and so the blob appears to fall with a constant velocity. This behavior can be found in time-distance diagrams of coronal rain events, as for example Figure 6 of De Groof et al. (2004), Figure 5 of De Groof et al. (2005), Figure 3 of Antolin et al. (2010), and Figure 7 of Antolin & Rouppe van der Voort (2012), and of descending prominence knots during the constant velocity phase (see Figures 3–5 of Chae, 2010).
To understand these two phases in the blob kinematics, we consider the path of the blob maximum density position, i.e., the pairs of values of the dashed line in Figure 1c. The blob falling speed is just the ions velocity obtained in the numerical simulation for these specific and ; this quantity is shown in Figure 1d. One can appreciate a rapid velocity increase (in absolute value) in the first 250 s followed by a slow deceleration that lasts for the rest of the simulation: these are the two phases mentioned before. The blob descending speed in the first stage reaches the maximum (unsigned) value km s, while in the second stage the blob speed decreases from 15 km s to 12 km s and so it is essentially constant although with a small upward acceleration.
Next, Equation (6) is divided by and so the terms on the right-hand side correspond to the acceleration caused by the inertial term, the pressure gradient, and gravity, respectively. These three accelerations, together with their sum, are plotted at the position of maximum blob density in Figure 2a. We see that the inertial term causes a negligible acceleration and that the pressure gradient results in no acceleration at . Nevertheless, the pressure gradient force increases in time and counteracts the blob weight at s. After this time the pressure gradient force overcomes the blob weight and the blob acceleration becomes slightly positive, as anticipated from the blob velocity in Figure 1d.
We now turn our attention to the plasma pressure and plot it as a function of height for various times at the beginning of the numerical simulation (Figure 2b). The presence of a mass condensation out of mechanical equilibrium in the corona generates a sound wave that propagates up and down at the (constant) sound speed, that in a fully ionised plasma at K is around 235 km s. This sound wave produces a significant modification of the pressure structure such that a steep pressure gradient builds up at the blob position. As the blob travels down, so does this pressure gradient, such as can be seen in the online movie accompanying Figure 2b. Then, we conclude that the vertical acceleration provided by the pressure gradient comes from a dynamical rearrangement of the plasma pressure in the corona.
4.2 Blob dynamics for higher blob density
In this section the calculations of § 4.1 are repeated for a blob ten times denser ( kg m), while keeping all other parameters unchanged. This heavier blob suffers both a stronger acceleration and a higher density increase. In the previous case (with kg m) the blob descends from Mm at to Mm at s, while in the present case the blob position after the same time is Mm (Figure 3a). In addition, now the blob density increases by more than 60% during the whole simulation (cf. Figure 3b), whereas in the previous (lower density) case the density increase was only about 6%. The two dynamical stages present in Figure 1c (initial acceleration followed by almost constant speed) are not so apparent in the time variation of the density (Figure 3c), but they become clear when we examine its velocity (Figure 3d). In this last plot we see that the acceleration stage lasts much longer in the present case (until s compared to s for the less dense blob) and for this reason the maximum descending speed is also much larger ( km s compared to km s). After this phase, the second stage with a moderate velocity increase is also found.
Now, the accelerations caused by the forces in the momentum equation (inertial term, pressure gradient, and gravity are represented (Figure 4a). In the present case the blob falls at the acceleration of gravity during the initial s. There are three main differences between Figure 4a and Figure 2a. Perhaps the most outstanding one is the presence, during the whole simulation, of wiggles in the inertial term acceleration. The origin of these oscillations is the emission of sound waves by the blob, that is explored in § 4.5. The second difference also concerns the inertial term acceleration, that increases by a factor when is varied from kg m to kg m. The reason for this behaviour can be found not only in the increase of the blob velocity (as seen from Figures 1d and 3d), but also in the increase of its -derivative. And the third difference between Figures 2a and 4a is that the rate of increase of the pressure gradient is much smaller now and so the total acceleration is negative during a longer time (cf. the two green lines). There are two reasons for this slower build up of the acceleration produced by the pressure gradient. First, to obtain the acceleration caused by this gradient is divided by , which is at least ten times larger for kg m than for kg m. Second, it is more difficult to generate a steep enough pressure variation at the blob position if it is denser. To explain this point, let us take, for example, the curve in Figure 4b that corresponds to s. The corona under the blob has been disturbed down to a height Mm by the sound wave excited at the initial blob position ( Mm), while the pressure distribution below this point is the initial one. Then, the pressure between the blob position and the lowest perturbed height cannot exceed that of this height, which is about 0.124 Pa for s. A similar argument gives the minimum pressure above the blob, Pa. The pressure jump between the top and bottom of the blob is limited by the range of heights that has been perturbed at any time, which in turn gives an upper limit to the pressure gradient. The conclusion is that denser blobs require the generation of larger pressure gradients to counteract their weight and so this process (to which we refer as the acceleration phase) lasts longer.
4.3 Blob temperature
Having omitted non-ideal effects in our equations, the blob temperature variation in time is determined purely by the ideal gas law once the plasma pressure and density have been computed. Then, the spatial variation of the temperature in the blob obtained here would probably be quite different had heat conduction (that tends to smear out temperature variations) and radiation (that tends to cool the plasma) been taken into account.
The temperature of the two simulations of § 4.1 and 4.2 are presented in Figure 5 for several times. In the case of the less dense blob (Figure 5a) the coronal temperature does not depart much from its initial value, K. The largest differences from the initial temperature are found below the falling condensation and arise because of the pressure restructuring caused by the initial sound wave emission (Figure 2a). Moreover, the blob temperature displays a Gaussian shape and shows a minimum value around K that remains essentially constant in time. This temperature is higher than expected for coronal rain and descending prominence knots, with observed emission in cool spectral lines. Regarding the simulation with the higher initial blob density (Figure 5b), now the coronal temperature presents a slight increase below (decrease above) the blob position, caused by a similar variation in the ambient pressure. The blob temperature has a shape similar to that of Figure 5a and, despite the strong increase of the blob density (see Figure 3a), its minimum temperature only increases slightly during the simulation, namely from K to K.
4.4 Blob shape
Here we examine the change in blob shape during the numerical simulations using the results of § 4.1 and 4.2. In Figure 6a we plot together the blob density at s with a common height origin. To do this plot, for each time () a new height variable () is defined by displacing so that the height of maximum density is . The increase of the maximum density noticed in Figures 1a and b is rather evident here and it is accompanied by a slight distortion of the Gaussian shape: after 1000 s the blob density has become a bit steeper in its bottom side () and a bit smoother in its upper side (). Since numerical dissipation is negligible in our computations, this seems a real effect. Now, we calculate a proxy for the total blob mass by integrating the density of Figure 6a with respect to in the range Mm Mm. Given that the environment density increases as the blob descends, the background is subtracted before obtaining this integral. The results display a very weak monotonic decrease of the blob mass in time. After 1000 s this mass reduction is about 0.8% and so this effect does not seem too relevant; we conclude that the blob mass is practically constant during its motion.
A similar plot is presented in Figure 6b for kg m. A comparison of Figures 6a and b reveals that a denser blob undergoes a larger maximum density increase (as discussed in § 4.2) and also a stronger shape distortion, the density at the bottom part of the blob displaying a steep profile at s. One may wonder whether shocks can develop at the bottom of falling blobs with high enough density. We have found no indication of this phenomenon in observations of rapidly falling cold, dense blobs. In addition, the total blob mass also remains essentially constant in time for the present initial blob density, with a reduction of only 0.5% with respect to the initial value after 1000 s.
4.5 Sound wave emission by the blob
We next pay attention to the outstanding oscillations of the inertial term acceleration (that is equal to ) seen in Figure 4a. Note that this acceleration (like all others in this figure) is taken at the position of the blob maximum density and so these oscillations must come from periodic variations in and/or at the maximum blob position. To investigate the presence of periodicities in these signals we have made use of the wavelet analysis, that has also been applied to the blob maximum density (Figure 3b). Before applying the wavelet analysis, the trends of and have been removed by fitting and subtracting a third degree polynomial. We start with , that has the strongest oscillatory power. Its wavelet diagram (Figure 7a) displays three periodicities with periods s, s, and s together with power at the full range of periods between 10 and 30 s. We concentrate in the three features with well-defined periods. They correspond to the oscillations visible in Figure 4a and have three properties: their period remains essentially constant in time, they have small amplitude, and the oscillatory power decays in time, with shorter periods displaying a quicker attenuation. The vertical velocity (Figure 7b) only shows the first of these periodicities, while the blob density contains the first and second ones (Figure 7c). Hence, the blob is subject to some phenomenon that produces at least two or three periodic, damped variations of , , and .
To understand these damped oscillations we notice that the blob has smaller sound speed than its environment. This sound speed depression is similar to the Alfvén speed depression of a coronal magnetic slab, a structure that supports both trapped and leaky fast waves (e.g. Terradas et al., 2005). We can then establish an analogy between the two cases and so we hypothesise that the falling blob can support small amplitude, leaky sound waves. This analogy is developed in Appendix A. The first three leaky sound waves whose transverse velocity is even about the blob center have periods s, 42.6 s, and 29.2 s, that agree quite well with those in the numerical simulations (see horizontal lines in Figure 7). Moreover, the corresponding exponential damping times of these modes are s, 182 s, and 135 s, so that our finding that shorter periodicities damp faster is also in agreement with these theoretical leaky mode results.
The calculations carried out in Appendix A give us more information than just the wave frequency. They indicate that solutions with even velocity about the blob center have non-zero velocity together with vanishing density perturbation and derivative of the velocity at the blob center (i.e., they are vertical blob oscillations that do not change its shape). This is in contradiction with our finding that the three quantities oscillate during the course of the blob descent. There are two reasons for this discrepancy: first, there is a continuous flow of material from the lower part of the blob into it and out from the blob above it. This issue is presented in § 4.6 and could be incorporated in Appendix A by assuming an equilibrium flow, i.e., a non-zero velocity in Equation (A1). In these circumstances, we do no longer have even or odd solutions about the blob center and all three variables (, , ) are non-zero there. This symmetry is also broken by the change of the blob shape, that is more important for larger blob densities. And finally, the simple equilibrium used in Appendix A is just a coarse approximation to the dynamical situation of a falling blob.
The results discussed in this section correspond to the numerical simulation with the densest blob. In the case of the numerical simulation with a smaller initial blob density the inertial term does not seem to have the same type of oscillations (cf. Figure 2a); we nevertheless have repeated the wavelet analysis with the same three variables (, , and at the blob maximum). The results are analogous for the two simulations, although those of the simulation with the largest initial density ( kg m) display clearer oscillatory features. So we conclude that the larger the blob density, the stronger the sound wave emission during its fall. In addition, the initial blob density is also crucial both in the period and damping time of the oscillations, that are both smaller for the less dense blob: s, 13.9 s, 9.54 s and s, 18.7 s, 13.6 s. The ratio of damping time to period determines the number of periods required to achieve a given attenuation. This quantity is much smaller for the simulation with initial blob density kg m than for that with kg m. In other words, the less dense blob has oscillations with smaller amplitude that attenuate much faster, so the detection of this sound wave emission during the blob fall is a signature of a high density. The reason for the smaller periods and shorter attenuation times in the less dense blob must be found in the larger value of the sound speed in the blob: this reduces the transit time inside the blob (and hence the period) and makes the sound speed depression shallower (and so worsens the blob capability for trapping sound waves).
4.6 Motion of plasma elements
The motion of each fluid element can be followed in time thanks to the Lagrangian displacement, . An element of plasma initially at a height will be at a height at time . Here the results for kg m are used to discuss the behavior of the plasma element initially at the blob center, i.e., with initial height Mm. Figure 8a reveals that this element moves upwards relative to the blob maximum and that after 250 s it is even outside the blob. Thus, as the blob falls down, plasma below it is engulfed by the condensation (and therefore highly compressed) and, at the same time, plasma in the blob is left above it and so must expand to meet the local, coronal density. The continuous separation between the plasma element initially at Mm and the blob maximum position is better seen by representing their respective positions (Figure 8b), that slowly diverge in time. The total acceleration of these two plasma elements is shown in Figure 8c. During most of the time, the blob maximum position has a downward acceleration that is larger by a few tens of m s, and this causes the departure displayed in Figures 8a and b. All the effects discussed in this paragraph are also found for the smaller blob initial density ( kg m), although with a lower magnitude.
4.7 Parametric study
In the numerical simulation presented above with kg m the blob speed is considerably smaller than that of many coronal rain blobs and descending prominence knots. On the other hand, the falling speed attained by the blob with initial mximum density kg m is considerably larger and surpasses that of descending prominence knots and even some coronal rain events. Our purpose now is to determine the influence of the parameter values on the mass condensation dynamics and on its maximum descending speed. We consider as reference values those used in § 4.1 and 4.2 and change only one parameter at a time.
First the blob maximum density, , is varied in the range kg m. The blob height versus time is presented in Figure 9a. In all cases the blob displays the initial acceleration stage followed by a period with a more or less straight trajectory, that is, a more or less constant speed. This figure illustrates well our previous finding that denser blobs fall faster, such as can be confirmed from the blob velocity, see Figure 9b. We note that consecutive values of differ more or less by a factor of two and that in these two panels the curves for kg m and kg m are closer than other adjacent pairs of curves. Hence, it seems that increasing the initial blob density even more will lead to a more moderate modification of the blob trajectory and maximum descending speed. Both Figures 9b and c confirm that the initial acceleration phase is present for all blob densities, although it lasts longer for larger densities of the mass condensation. In addition, in all cases except that of the densest condensation, the blob always has an upward acceleration at the end of the simulation. Finally, Figure 9c also confirms that denser blobs are characterised by larger oscillations in their acceleration, i.e., stronger sound wave emission.
Since the height dependence on time gives enough information on the blob kinematics, now we do not show the blob velocity and acceleration. First the influence of the background temperature is studied: along with the reference value K, we have also considered K and K (Figure 10a). Varying the coronal temperature does not modify the shape of the trajectory, but a reduction in leads to a descent with higher downward acceleration. This effect is explained as follows: smaller values of lead to a decrease of the vertical stratification scale height and, given that the base density has not been altered, this implies that the blob moves in an ambient plasma with reduced density. Thus, the blob to environment density ratio is larger. This is analogous to maintaining the coronal structure and increasing the blob density, such as in Figure 9, and so both a blob with higher density and a corona with smaller temperature result in a stronger acceleration and a larger descending speed. As a conclusion, the ratio of blob to background density controls the blob descent.
Now, a similar effect is obtained if the background temperature is maintained and the base coronal density (and hence plasma pressure) is decreased: the blob also falls in a rarer environment and achieves greater speeds (see triangles in Figure 10b). Moreover, if the blob length is doubled, while the rest of parameters are held fixed, then its mass doubles and once more it falls at a larger speed (Figure 10b). In fact, this figure shows that halving the coronal density or doubling the blob size produce the same height versus time variation. These two blob trajectories are also identical to that of the case with double initial blob density (green line in Figure 9a).
Thus, the results presented in this section show that the blob dynamics is governed both by the density contrast of the blob with respect to the ambient plasma and by its mass. Since the observations of cold descending condensations provide both the blob height and speed versus time, here we also plot the maximum (i.e., most negative) descending velocity, . Figure 11a shows that the strongest variation of the maximum speed occurs for the smaller values of in this plot, whereas the rate of change of is smaller on the right side of the plot. This confirms the prediction we made before. Figure 11a emphasises the relevance of the blob density (or the density ratio) in the blob dynamics. On the other hand (Figure 11b), changing the temperature in the range K has not such a big influence, although for small values of the maximum blob density it can cause to vary by a factor 1.5–2.
5 Path curvature
In this section we are concerned with a fully ionised plasma blob falling along circular magnetic field lines of radius . We assume a very large electric conductivity, thus ensuring that the plasma cannot move across magnetic field lines. In addition, we assume a low- situation so that the plasma cannot drag the magnetic field with it. Instead, the magnetic field maintains its shape and guides plasma motions.
We impose that magnetic field lines make an angle with the vertical direction at and so the height of a point in the falling trajectory is given by
with the distance from the level to this point along the field line. The paths used in this section are shown in Figure 12.
were the unknowns , , and here depend on and and is the component of the acceleration of gravity parallel to the circular path. Note that in Equations (5) to (7) represents the vertical velocity, but in Equations (16) to (18) it represents the velocity component parallel to the curved path.
Using the ideal gas law for charged particles, their pressure and density distribution along the circular trajectory at is
where the vertical scale height, , has the same value as for the vertical trajectory and is given by Equation (15). This formula is formally identical to its counterpart for a vertical path (Equation (11)).
Regarding the mass condensation at , we use the following expression (which is analogous to Equation (14))
Hence, the initial distance of the blob center from along the curved magnetic field lines is , although its height is obviously smaller than this value; cf. Equation (15).
This initial plasma configuration (vertically stratified coronal atmosphere and embedded mass condensation) is identical to the one we have considered so far if one substitutes for and for . Equations (16) to (18) are also formally identical to their counterparts for vertical motions (i.e., Equations (5), (6), and (7)), except for the projected acceleration of gravity along the circular magnetic field lines, . We can thus expect the effect of curvature to be noticeable only if the path is inclined from the vertical direction, so that departs significantly from . Moreover, if the condensation travels a long distance along the magnetic field then the effective acceleration of gravity will increase during this motion and this will also influence its dynamics.
We have solved Equations (16) to (18) using the parameter values of § 4.1 and 4.2 and for two different radii of magnetic field lines ( Mm and 250 Mm) and two different inclinations of magnetic field lines at (namely and ). To explore the blob dynamics, in Figure 13a we represent the blob height (for a vertical path) or its distance from the level along magnetic field lines (for a circular trajectory). The time-distance variation of the four cases is rather similar and, except for the one with inclined field lines at the coronal base, the trajectories of the other three are almost undistinguishable. A small departure can only be appreciated at the end of the numerical simulation. The variation of the blob acceleration with time, Figure 13b, displays an initial value of equal to the projection of on field lines. For the case Mm, this initial acceleration is of the order of 200 m s, whereas for the curved field lines that are vertical at the coronal base it is much closer to the vertical acceleration of gravity. In spite of these differences in the starting value of , in a very short time ( s and s for the cases and , respectively) the blob moving along a curved trajectory displays an acceleration that closely matches that of the vertically falling condensation.
These results confirm our previous speculations about the effect of curvature. We hypothesised that the path inclination would have some effect only if is different enough from . The two cases in Figures 13a and b that coincide with the vertical motion, i.e., the curved lines with Mm and 250 Mm, , have that differs 4% from at most along the blob track. On the other hand, the case with displays a maximum difference between and of 25%. It must be emphasised, nevertheless, that the interplay between the pressure gradient and gravity leads in the four cases to a similar blob acceleration after 100 s at most.
To complete our study of the effect of the path curvature we make use of the higher blob density considered above, kg m. The variation of blob position and the acceleration along the descending path versus time are shown in Figures 13c and d. Just as happens with the case kg m only the blob falling along the curved magnetic field with an angle at descends with a different dynamics, although its departure from the other three curves is much larger in Figure 13c than in Figure 13a. Moreover, the acceleration for is far from that of the other three cases because having a larger density ratio implies that the pressure structure takes much longer to readjust in order for its gradient to compensate the gravity force. An additional effect that appears for is that the inclination of field lines changes during the blob motion so as to increase and this requires a continuous readjustment of the pressure configuration in order to reach the balance with gravity.
In this work we have investigated some physical processes present in dense downfalling plasma blobs such as coronal rain and dense prominence knots. To simplify our task several important ingredients have been ignored: the formation process of the blobs, the presence of neutrals (that will be included in the future), non-ideal effects, and the effect of magnetic fields (apart from behaving as rigid guides for plasma motions). Moreover, the influence of the chromosphere and photosphere has also been omitted by putting the numerical boundaries far from the range of heights of interest. Our model thus contemplates the dynamical evolution of a fully ionised plasma condensation in an isothermal, vertically stratified corona, with motions that are either vertical or guided along rigid magnetic field lines (with assumed circular shape).
We have first focused our attention on a blob falling vertically. The mass condensation is released at rest and starts to accelerate, showing two distinct phases. In the first one, gravity is the dominant force and causes the downward acceleration of the blob. At a sound wave is emitted at the initial blob position. This wave results in a considerable modification of the coronal pressure, by which a strong pressure gradient is generated at the blob position. At the end of the first phase this gradient, that moves with the mass condensation as it falls down, becomes large enough to cancel the acceleration of gravity and the blob enters the second phase, characterised by a roughly constant velocity. This two-stage dynamics is very often found in time-distance diagrams of coronal rain events, as for example Figure 6 of De Groof et al. (2004), Figure 5 of De Groof et al. (2005), Figure 3 of Antolin et al. (2010), and Figure 7 of Antolin & Rouppe van der Voort (2012). Moreover, Schrijver (2001) has represented the temporal variation of the descending speed for 45 coronal rain features and many of them clearly display the initial acceleration phase. This remarkable similarity between our numerical simulations and observations points out that the simulations take into account some of the relevant physics for the dynamics of descending blobs.
A larger blob mass requires a larger pressure gradient to balance gravity and this leads to a lengthening of the downward acceleration phase, that can last between 100 s and 1000 s. A longer acceleration phase in turn allows the blob to achieve higher velocity and so we have found an increase of the maximum falling speed with the initial blob density (). For the parameter values used in this work, this maximum speed is more sensitive to changes of around kg m than to changes around kg m. Another important conclusion refers to the density contrast, that for a fixed coronal temperature is determined by the initial blob density and the density at the base of the corona. If these two parameters are modified simultaneously in such a way that the density contrast is held fixed, then the blob dynamics does not change. And finally, we have also established that the coronal temperature also has its influence because it determines the coronal vertical scale height. We have seen that a colder corona has a smaller vertical scale height and so, for a given density at , the mass condensation moves in a rarer environment compared to that of a hotter corona. Hence, this situation is similar to that of a blob with a higher density contrast, so that the condensation descends faster in a cold corona than in a hotter one.
Blob speeds in the range 10–150 km s have been obtained in this work. This velocity range covers all observations of coronal rain events and falling prominence knots, but a detailed comparison with observations is not pertinent because, as discussed in the previous paragraph, a given blob velocity can be obtained with different combinations of our parameters. Nevertheless, other results in our simulations could be compared to observations. For example, although the blob mass is conserved in time, its maximum density has been found to build up. In particular, for a plasma condensation with density of the order of kg m this increase is about 6% after 1000 s, while for larger initial blob densities (of the order of kg m) it is above 60%. In this last case the condensation density also presents a strong distortion with respect to its initial profile: the density variation at the blob bottom becomes steeper with time and so a large enough gradient could develop such that a shock might appear. Another feature of descending condensations is that they emit small amplitude sound waves as they fall. This wave emission is completely different from the one observed at , that distorts the background corona, and is present during the whole simulation. We have performed a simple perturbation analysis of this phenomenon and have derived leaky wave periods that are in good agreement with those found in the simulations. The sound wave emission during the blob lifetime is more evident for larger blob densities and might be observed as small periodic changes of the blob intensity. A trend subtraction could be necessary to achieve this detection.
Coronal rain blobs are observed to move along curved trajectories with a variety of shapes and inclinations (see e.g., Schrijver, 2001; Antolin & Rouppe van der Voort, 2012). It has been found that loop curvature is irrelevant for the blob dynamics if magnetic field lines are vertical at the coronal base. Noticeable differences are only present if magnetic field lines make an angle with the vertical direction at (here the value has been used) and if the condensation to corona density ratio is large enough. In view of these results, one could argue that the coronal rain condensations detected by Schrijver (2001), whose acceleration depends very little on the path inclination, had a moderate density and/or small path inclination. It is of great interest to solve the two-dimensional problem of partially ionised mass condensations falling along curved magnetic fields to see if the coupling between ions and neutrals is strong enough to drag neutrals along the magnetic field and to prevent them from falling vertically. In fact, observations seem to indicate that this must be the case since coronal rain blobs detected in the H line usually display curved trajectories. For a numerical study of fully ionised plasma blobs moving in a magnetic arcade the reader is referred to Fang et al. (2013).
The motion of individual plasma elements has also been investigated and it has been found that the position of blob maximum density moves faster than nearby particles. Thus, the background material penetrates inside the blob from below while the blob material can leave it at its top and become part of the background plasma. This is more than just a curiosity. In a partially ionised plasma it may lead to the recombination (ionisation) of material entering (leaving) the blob and so this may add more complexity to the interaction between neutrals and charged particles. This will be investigated in a forthcoming article.
The blob temperature remains constant during the whole simulation regardless of its initial density. It is worth mentioning that our model does not include thermal conduction, that could smear out the large temperature difference between the plasma condensation and the environment. However, radiative losses are very efficient at low temperatures and not only lead to the blob formation by catastrophic cooling, but they also help the condensation to keep its temperature after being formed (Müller et al., 2004, 2005; Antolin et al., 2010; Fang et al., 2013).
We finish with a comment on our one-dimensional assumption. If motions are guided by vertical or curved magnetic fields, then this assumption is reasonable because the material moving inside a magnetic tube is decoupled from that outside it. Nevertheless, in the absence of a magnetic field the blob behaves like a piston moving in a liquid. The piston presses the liquid, that can move around the piston, so that the strong pressure gradient we found in our simulations could actually be substantially smaller and blob speeds substantially larger. The second phase of blob dynamics, namely that with roughly constant speed, would also be modified or could even disappear. This is a motivation for undertaking two- or three-dimensional simulations.
Appendix A Small amplitude leaky sound waves in a dense blob
The usual procedure to study small amplitude trapped and leaky waves is to start from an equilibrium state and next to consider linear perturbations about this equilibrium. The problem here is that we want to study waves emitted by the blob, that is not in an equilibrium. For this reason, our results can only provide an approximation to the actual wave leakage found in the numerical simulations.
We neglect gravity and consider a static equilibrium with uniform pressure () and a density consisting of a uniform background () plus a condensation of the form of Equation (14). This equilibrium state is defined by
where now gives the central, unperturbed blob position. The next step is to take perturbations about this equilibrium and, given that we are interested in periodic variations, their time dependence is of the form . The pressure, density, and velocity perturbations are respectively given by the functions , , and multiplied by . With these assumptions, the linearised form of Equations (5)–(7) is
where the sound speed squared is defined as . The sound speed has a minimum value at the blob center and achieves its maximum value as . These two values of are denoted by and , respectively.
The two possible types of solutions to Equation (A5) are trapped waves, for which as and is real, and leaky waves, for which oscillates and grows as and is complex with positive imaginary part. Moreover, both trapped and leaky solutions can either have even and , odd about the blob center or odd and , even about . We refer to these solutions simply as even and odd solutions.
Solutions to Equation (A5) are obtained numerically using the well-known shooting and matching technique, with the boundary conditions obtained as follows. At a large distance from the plasma condensation the medium can be considered uniform and with a sound speed equal to . Then, the plasma velocity is . This solution represents the superposition of two waves propagating in opposite directions. If is real this solution is oscillatory and so the present equilibrium configuration does not support trapped sound waves. If is complex Equation (A5) has leaky solutions, characterised by a temporal attenuation and a spatial growth at large distances from the condensation. To find leaky solutions we write , where both and are real. In order to have perturbations that damp in time, we require to be positive. We also take with no loss of generality. Then, considering in the former expression for , the term with amplitude () corresponds to a wave travelling towards the blob (away from the blob), whereas the opposite applies for . To study wave emission by the blob we must only retain waves emanating from it and so we have for and for . We now proceed as follows: at a large distance from the blob ( Mm, say) we impose that the plasma velocity and its -derivative can be computed from , where a guess value of must be provided. Note that since we deal with a linear problem, we are allowed to take . Starting from these initial conditions, Equation (A5) is integrated to the position Mm, which is symmetrically placed about with respect to the starting point of the numerical integration. If the guess value of the frequency were correct, then we would obtain the right values of and . For even solutions they can be computed from (i.e., ), whereas for odd solutions they come from (i.e., ). Nevertheless, our initial value of will not be correct in general, so there will be a discrepancy between the numerical and analytical solutions at Mm. An iterative method is then applied by which and are varied so that this discrepancy is reduced until convergence of the numerical method is achieved. At this point we have acceptable approximations to the frequency and .
|Odd solutions||Even solutions|
|( s)||( s)||(s)||(s)||( s)||( s)||(s)||(s)|
To make a comparison with the simulation of § 4.2 we take the equilibrium density and pressure equal to those of the numerical simulation at the initial blob position, i.e., Pa, kg m, and kg m. Then, km s and km s. The results are summarised in Table 1, that contains , , the period (), and the exponential damping time () of the first three even and odd leaky waves. The period of the even solutions are shown in Figure 7 as horizontal lines and are used in § 4.5 to determine that the oscillations detected in Figure 4 are caused by the emission of leaky sound waves during the blob descent.
- affiliation: Also at Abastumani Astrophysical Observatory at Ilia State University, Cholokashvili Ave. 3/5, Tbilisi, Georgia
- affiliation: Also at Institute of Nuclear Physics, Moscow State University, Leninskie Gory, 119992, Moscow, Russia
- Antolin, P., & Rouppe van der Voort, L. 2012, ApJ, 745, 152
- Antolin, P., Shibata, K., & Vissers, G. 2010, ApJ, 716, 154
- Antolin, P., Vissers, G., & Rouppe van der Voort, L. 2012, Sol. Phys., 280, 457
- Beckers, J. M. 1962, Australian Journal of Physics, 15, 327
- Berger, T. E., Shine, R. A., Slater, G. L., et al. 2008, ApJ, 676, L89
- Berger, T. E., Slater, G., Hurlburt, N., et al. 2010, ApJ, 716, 1288
- Bruzek, A. 1969, Sol. Phys., 8, 29
- Chae, J. 2010, ApJ, 714, 618
- De Groof, A., Bastiaensen, C., Müller, D. A. N., Berghmans, D., & Poedts, S. 2005, A&A, 443, 319
- De Groof, A., Berghmans, D., van Driel-Gesztelyi, L., & Poedts, S. 2004, A&A, 415, 1141
- Engvold, O. 1976, Sol. Phys., 49, 283
- —. 1981, Sol. Phys., 70, 315
- Engvold, O., Jensen, E., & Andersen, B. N. 1979, Sol. Phys., 62, 331
- Fang, X., Xia, C., & Keppens, R. 2013, ApJ, 771, L29
- Foukal, P. 1978, ApJ, 223, 1046
- Gilbert, H., Kilper, G., & Alexander, D. 2007, ApJ, 671, 978
- Gilbert, H. R., Hansteen, V. H., & Holzer, T. E. 2002, ApJ, 577, 464
- Haerendel, G., & Berger, T. 2011, ApJ, 731, 82
- Heinzel, P., Schmieder, B., & Mein, P. 1992, Sol. Phys., 139, 81
- Hillier, A., Isobe, H., & Watanabe, H. 2011, PASJ, 63, L19
- Leroy, J. L., Bommier, V., & Sahal-Brechot, S. 1983, Sol. Phys., 83, 135
- —. 1984, A&A, 131, 33
- Liu, W., Berger, T. E., & Low, B. C. 2012, ApJ, 745, L21
- López Ariste, A., & Aulanier, G. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 368, The Physics of Chromospheric Plasmas, ed. P. Heinzel, I. Dorotovič, & R. J. Rutten, 291
- Loughhead, R. E., & Bray, R. J. 1984, ApJ, 283, 392
- Low, B. C., Berger, T., Casini, R., & Liu, W. 2012a, ApJ, 755, 34
- Low, B. C., Liu, W., Berger, T., & Casini, R. 2012b, ApJ, 757, 21
- Mackay, D. H., Karpen, J. T., Ballester, J. L., Schmieder, B., & Aulanier, G. 2010, Space Sci. Rev., 151, 333
- Merenda, L., Trujillo Bueno, J., Landi Degl’Innocenti, E., & Collados, M. 2006, ApJ, 642, 554
- Müller, D. A. N., De Groof, A., Hansteen, V. H., & Peter, H. 2005, A&A, 436, 1067
- Müller, D. A. N., Peter, H., & Hansteen, V. H. 2004, A&A, 424, 289
- Murawski, K., Zaqarashvili, T. V., & Nakariakov, V. M. 2011, A&A, 533, A18
- Schrijver, C. J. 2001, Sol. Phys., 198, 325
- Sewell, G. 2005, The Numerical Solution of Ordinary and Partial Differential Equations, Pure and Applied Mathematics Series (John Wiley & Sons)
- Terradas, J., Oliver, R., & Ballester, J. L. 2005, A&A, 441, 371
- Torrence, C., & Compo, G. P. 1998, Bulletin of the American Meteorological Society, 79, 61
- Wiik, J. E., Schmieder, B., Heinzel, P., & Roudier, T. 1996, Sol. Phys., 166, 89
- Xu, A.-A. 1987, in Solar Maximum Analysis, 147–155
- Zhang, J., & Li, L.-P. 2009, Research in Astronomy and Astrophysics, 9, 1368