© ESO, 2017

1. Introduction

Coronal loop plasma can be thermally unstable and subject to formation of cool and dense condensations (Field 1965). Thermal instability is likely to occur in loops with heating concentrated at the footpoints. If the thermal conduction along the loop is not efficient enough, the radiation losses near the loop top overcome the heating input, resulting in the onset of a thermally unstable regime (Antiochos et al. 1999; Karpen et al. 2001). Coronal rain is a phenomenon occurring in such footpoint-heated coronal loops as a result of catastrophic cooling. It consists of cool plasma condensations that are formed near the loop top and fall towards the solar surface, guided by the magnetic field lines. High-resolution solar observations show that the motion of coronal rain blobs is significantly sub-ballistic (e.g. Antolin & Verwichte 2011; Antolin & Rouppe van der Voort 2012; Kohutova & Verwichte 2016), suggesting that forces other than gravity have an important effect on its dynamics and evolution.

The formation and evolution of plasma condensations have been addressed by a number of numerical studies using 1D hydrodynamic simulations. The thermal instability onset and coronal rain formation in a coronal loop with footpoint-concentrated heating typically depends on the spatial distribution of the heating input and often occurs in periodically repeating limit cycles (Müller et al. 2003, 2004). Using a 1D approach, the pressure effects are found to have a strong influence on the motion of the individual coronal rain blobs, often counteracting the effect of gravity (Antolin et al. 2010; Oliver et al. 2014) in the case of a compressed plasma below the condensation, providing net upward pressure gradient force. Conversely, if a plasma condensation is moving in a low-pressure region, such as in the wake of a previously formed condensation, this can result in motion that is faster than free-fall (Müller et al. 2005).

The 1D hydrodynamic simulations modelling the evolution along a single field line neglect the effect of the finite magnetic field, however, as all of the plasma is confined below the condensation, and therefore 1D simulations likely overestimate the decelerating effects of coronal loop plasma pressure gradients. The effect of the finite magnetic field on the coronal rain evolution is only properly accounted for when using multidimensional magnetohydrodynamical (MHD) models.

Siphon flows due to local pressure variations can also have strong effect on the motion and morphology of coronal rain condensations, as shown by 2D MHD studies of coronal rain formation and evolution (Fang et al. 2013, 2015).

Using 2D MHD simulations, Mackay & Galsgaard (2001) investigated the evolution of a density enhancement in the context of cool prominence material. The density enhancement was found to rebound multiple times in this setup, which was explained as a two-step process: a deceleration phase caused by the pressure build-up below the enhancement, and a rebound phase caused by the restoring action of the Lorentz force, stressing the importance of the effect a finite magnetic field can have on the evolution of plasma condensations. Similar longitudinal oscillations of cool condensations can be seen in coronal rain observations (Kohutova & Verwichte 2016; Verwichte et al. 2017) as well as in prominences (e.g. Jing et al. 2003; Zhang et al. 2012; Luna et al. 2014).

We model the evolution of a cool plasma condensation in a realistically stratified atmosphere that includes a cool high-density chromosphere, a transition region layer, and a hot corona. We choose our problem setup to be representative for small coronal rain condensations that are typically formed in thermally unstable loops as a result of catastrophic cooling. The geometry of the problem is therefore set up to reflect the coronal loop geometry, accounting for the reduced effective gravity due to the semicircular shape of the loop. We furthermore analyse condensation trajectories, velocities, and accelerations in order to be able to compare them to recent high-resolution coronal rain observations. Finally, we propose an analytical model for the condensation dynamics in order to explain the oscillatory behaviour of the plasma condensations and compare it with the numerical findings.

2. Numerical model

We solve the nonlinear MHD equations using Lare2d (Arber et al. 2001) assuming neutral fully conductive plasma and including gravity and shock viscosity. We use the ideal equation of state. Thermal conduction and radiative transport are not included in the energy equation. We introduce a rectangular simulation domain with the extent −30 Mm Mm in horizontal direction and −120 Mm Mm in the vertical direction with 512 × 2048 resolution. The coronal loop is modelled as a straight slab along the y-direction. We adopt the variable s to describe the position along the loop from one footpoint to the other, that is, s = y + 120 Mm. Thus, the loop has a length of 240 Mm. The density variation between the loop and the background medium in the x-direction is given by the symmetric Epstein profile (Nakariakov & Roberts 1995): $\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{e}}\mathrm{+}\mathrm{(}{\mathit{\rho }}_{\mathrm{i}}\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}{\mathrm{sech}}^{\mathrm{2}}\left(\frac{\mathit{x}}{\mathit{a}}\right)\mathit{,}$(1)where ρ_e and ρ_i are external and internal densities, respectively, and a = 3 Mm is the loop scale width. We assume a constant density contrast ρ_i/ρ_e = 10 throughout the whole domain. In order for the setup to be representative of a semicircular coronal loop with both footpoints anchored to the photosphere, the effective gravity g_eff is determined assuming a semicircular coronal loop of length L and varies with the coordinate along the loop s as ${\mathit{g}}_{\mathrm{eff}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{\mathit{g}}_{\mathrm{\odot }}\cos \left(\frac{\mathit{\pi s}}{\mathit{L}}\right)\mathit{,}$(2)such that it equals zero at the loop top in the centre of the domain and g_⊙ = 274 m s^-2 at the loop footpoints at the top and bottom domain boundaries (Fig. 1). The temperature is constant in the x-direction. In the vertical direction we create a realistic temperature variation representative of an atmosphere consisting of a cool chromosphere, transition region layer and hot corona by adopting a smoothed step function temperature profile (Cargill et al. 1997): $\mathit{T}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{T}}_{\mathrm{cor}}\mathrm{+}{\mathit{T}}_{\mathrm{ph}}\mathrm{)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{T}}_{\mathrm{cor}}\mathrm{-}{\mathit{T}}_{\mathrm{ph}}\mathrm{)}\tanh \left(\frac{\mathit{h}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{-}{\mathit{s}}_{\mathrm{t}}}{\mathrm{\Delta }\mathit{s}}\right)\mathit{,}$(3)with photospheric temperature T_ph = 6 × 10³ K, coronal temperature T_cor = 10⁶ K, s_t = 4 Mm, Δs = 1 Mm and $\mathit{h}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathit{L}}{\mathit{\pi }}\sin \frac{\mathit{\pi s}}{\mathit{L}}$. The temperature variation controls the pressure scale height H(s): $\mathit{H}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{{\mathit{k}}_{\mathit{b}}\mathit{T}\mathrm{\left(}\mathit{s}\mathrm{\right)}}{\mathit{m}{\mathit{g}}_{\mathrm{eff}}}\mathrm{·}$(4)The density profile for the non-isothermal stratified atmosphere is then determined by numerically solving for a hydrostatic pressure balance:

$\begin{array}{ccc}& & \mathit{p}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\mathrm{0}}\exp \left(\mathrm{-}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{h}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\frac{\cos \mathrm{(}{\frac{\mathit{\pi }{\mathit{s}}^{\mathrm{\prime }}}{\mathit{L}}}^{\mathrm{)}}\mathrm{d}{\mathit{s}}^{\mathrm{\prime }}}{\mathit{H}\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}}\right)\mathit{,}\\ & & \mathit{\rho }\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}\frac{\mathit{mp}\mathrm{\left(}\mathit{s}\mathrm{\right)}}{{\mathit{k}}_{\mathit{b}}\mathit{T}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{·}\end{array}$The density stratification in the initial configuration is calculated using a footpoint density of ρ₀ = 6.5 × 10^-7 kg m^-3. This results in the densities in the upper region of the loop of the order of 10^-11 kg m^-3, which are representative of real coronal densities. The top and bottom boundaries are fixed to create a line-tied loop, and the boundary conditions along the vertical direction are symmetric (i.e. gradients are set to 0). The plasma condensation is superimposed on the background density and temperature profiles as follows. A 2D Gaussian enhancement is added to the equilibrium density profile and is positioned at x₀ = 0 Mm and s₀ = 100 Mm, that is, below the loop apex, of width σ = 0.5 Mm and height ρ_blob = r_bcρ_bg(x₀,s₀), with r_bc being the density contrast between the peak blob density and the density of the background loop plasma ρ_bg at the same position. We surround the condensation with a low-temperature region to maintain the plasma pressure balance and to prevent rapid initial expansion of the condensation in the vertical direction (the expansion in the transverse direction is counteracted by the magnetic Lorentz force). A grid convergence study using a grid with 1024 × 4096 resolution has been carried out in order to check the convergence of the numerical results.

Fig. 1 Initial effective gravity, density, and temperature profiles as a function of s at x = 0.

3. Blob evolution and kinematics

The evolution and kinematics of the plasma condensation, or blob, is analysed in detail for a magnetic field strength ranging from 20 G to 100 G and for three values of the initial density contrast between the condensation and the coronal loop plasma r_bc = 10², 10³, and 10⁴ corresponding to peak blob densities of 2.5 × 10^-9, 2.5 × 10^-8, and 2.5 × 10^-7 kg m^-3, respectively. The corresponding plasma-β ranges from 0.01 to 0.3.

The evolution of the density profile of the blob during the first 5600 s is shown in Fig. 2 for the B = 100 G, r_bc = 10³ case. This leads to an initial decrease of the blob density and to the emission of sound waves, most clearly seen in the case of r_bc = 10⁴. As the condensation falls, it develops an elongated tail as a result of the differential component of the effective gravity acting along the finite longitudinal extent of the blob. This elongation of the plasma condensations is in line with recent high-resolution coronal rain observations (e.g. Antolin & Verwichte 2011; Antolin & Rouppe van der Voort 2012; Kohutova & Verwichte 2016). The build-up of the density near the leading edge of the blob is further enhanced during the deceleration phase that occurs as the blob approaches the transition region (Fig. 2). Here the blob can be seen to rebound multiple times. When the blob hits the transition region, a rebound shock occurs that results in further sound wave emission. For low magnetic field strengths, the impact of the blob is accompanied by the ejection of chromospheric material since the finite plasma-β in the transition region and below does not restrict the transverse motion of the plasma. Except for heating the plasma along the blob edges, the overall temperature of the blob stays approximately constant during its downward motion. The temperature of the plasma below the blob increases as it is being compressed, whereas the plasma in the wake of the blob cools down. After the first rebound, the blob temperature slightly increases as a result of rebound shock dissipation.

Fig. 2 Evolution of the density profile along the bottom half of the loop plotted every 280 s during the first 500 time steps.

Fig. 3 Time-distance plot of the density along s at x = 0 for the rbc = 103, B = 100 G case.

The trajectory of the plasma condensation is determined by finding the position of the maximum in blob density along the vertical direction at each time step. This is subsequently used to deduce the evolution of the vertical velocity and acceleration. Two types of motion depending on the magnetic field strength and blob density are observed: a purely downward motion with the blob hitting the transition region, or damped oscillatory motion with the blob rebounding multiple times and eventually settling in an equilibrium position in the corona (Fig. 4). Higher magnetic field strengths lead to greater heights of the rebound points and greater heights of the equilibrium positions around which the blob oscillates. In addition, increasing the blob density leads to a decrease of the rebound point height and to a greater number of condensations reaching the surface. Similarly, Fig. 5 shows that the maximum downward velocity increases with increasing blob density and decreasing magnetic field strength. For r_bc = 10², the rebound motion occurs for all values of magnetic field strengths. For r_bc = 10³, purely downward motion occurs at low magnetic field strengths, while for r_bc = 10⁴, no rebound motion is observed.

Fig. 4 Height (top), velocity (middle), and acceleration (bottom) profiles of the condensation for different values of blob density contrast and magnitude of the magnetic field strength. The dotted line shows the free-fall height and velocity profiles.

Fig. 5 Dependence of the equilibrium height (left), maximum downward velocity (centre), and average period of blob oscillations (right) on the magnitude of the magnetic field strength for various blob densities.

For the lowest blob density, the blob motion shows distinct acceleration and deceleration phases: during the first ~300 s the blob accelerates downwards, afterwards it decelerates to t ~ 1000 s, followed by another acceleration phase lasting up to t ~ 1500 s and so forth, with the maximum values of the downward velocity ranging from 23 km s^-1 to 45 km s^-1 depending on the magnetic field strength.

The distinction between acceleration and deceleration phases is similarly clear for higher blob densities. There the maximum values of downward velocities are much higher, ranging from 107 km s^-1 to 130 km s^-1. The motion is sub-ballistic only in the case of lowest blob density, in the other two cases before the rebound, the blob falls approximately with free-fall speed. For the highest blob density, the effect of the varying magnetic field strength on the motion of the blob is negligible.

The motion of the coronal rain blobs deduced from high-resolution solar observations is mostly sub-ballistic with only few extreme cases (e.g. Antolin & Verwichte 2011; Antolin & Rouppe van der Voort 2012; Kohutova & Verwichte 2016). When we consider the significant effect that the peak blob density was found to have on its motion, the broad distribution of the blob velocities typically seen in the observations is therefore likely due to variations in masses of individual condensations. The extreme cases of observed velocities are likely caused by the variations in the plasma pressure across the coronal loop, for instance, when one blob travels in a wake of another, it can be siphoned into the region of the low pressure left behind by the first blob, which results in a motion that is faster than free fall (Müller et al. 2005).

4. Force balance analysis

In order to determine the relative influence of the individual forces on the motion and evolution of the plasma condensation, the vertical components of the gravitational force ρg_eff, pressure gradient force −∇p, magnetic pressure force −∇B²/ 2μ₀, and magnetic tension force (B·∇)B/μ₀ were calculated inside the loop and averaged in the transverse direction to obtain the longitudinal dependence. The evolution of the force balance during the first 5600 s is shown for the r_bc = 10³, B = 100 G case in Fig. 6.

For a rebounding blob, the force balance evolves as expected according to the conceptual model proposed by Mackay & Galsgaard (2001). As the condensation falls, it compresses the coronal loop plasma below it, leading to a build-up of the pressure gradient. For a high magnetic field strength (low plasma-β), the plasma is confined by the magnetic field, and as it expands below the blob, it pulls the magnetic field lines with it. The magnitude of the magnetic field strength therefore decreases and plasma moves away from the centre of the loop, as demonstrated by the positive divergence of the plasma velocity (Fig. 7). This results in a net upward magnetic tension force as shown in Fig. 6 before the first rebound at t~ 1500 s. The blob rebounds upwards, and as a result of the lower pressure below the blob, the plasma now moves towards the centre of the loop, which leads to a negative plasma velocity divergence. The field lines return to their original position, and the magnitude of the magnetic field strength increases (Fig. 7). This occurs multiple times until the blob eventually settles in an equilibrium position. The bending of the magnetic field lines is shown in Fig. 8.

In the case of weak magnetic field, the plasma-β in the transition region is high enough to allow the plasma below the blob to be displaced sideways, which prevents the pressure build-up. Hence no rebound motion occurs, and the blob falls directly towards the solar surface.

It should be noted that the initial uniform magnetic field configuration used here leads to a zero magnetic pressure gradient. However, in an expanding flux tube configuration, this is no longer the case. In this case, the magnetic pressure gradient force would have an additional decelerating effect and would therefore lead to lower downward velocities and greater equilibrium height of the oscillating blobs. Similarly, neglecting thermal conduction very likely affects the morphology of the cool condensations without significantly affecting the condensation dynamics. Given that the thermal conduction acts predominantly along the magnetic field lines, the only relevant exchange of energy will occur in the vertical direction. It is therefore not sufficient to remove significant amounts of thermal energy from the compressed underlying plasma, which would lead to large changes in plasma pressure.

Fig. 6 Evolution of the force balance along the bottom half of the loop plotted every 280 s during the first 500 time steps. We show the vertical components of the plasma pressure force (blue), gravity (green), magnetic pressure force (red), magnetic tension force (turquoise), and total net force in the vertical direction (black dashed line).

Fig. 7 Evolution of the magnitude of the magnetic field strength (top) and velocity divergence (bottom) averaged over the region below the plasma blob for the rbc = 103, B = 100 G case.

We further focus on the damped oscillatory motion of the plasma blob. The period of the individual rebound phases varies strongly with the blob density, while the dependence on the magnetic field magnitude is weak (Fig. 4). We therefore propose an analytical model for the period of the blob oscillations assuming a high β limit when the transverse motion of the plasma is prevented by the strong magnetic field in the vertical direction. We model a falling rain blob as a piston problem, where the rain blob is a piston compressing gas below it. We use a 1D model with s as the spatial coordinate along the loop. The rain blob has a fixed length 2Δs and its centre of mass is located at the position s_b(t). Its equation of motion is $\mathit{m}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{s}}_{\mathrm{b}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{mg}\mathrm{\left(}{\mathit{s}}_{\mathrm{b}}\mathrm{\right)}\mathrm{-}\left[{\mathit{p}}_{\mathrm{u}}\mathrm{(}{\mathit{s}}_{\mathrm{b}}\mathrm{+}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathrm{-}{\mathit{p}}_{\mathrm{d}}\mathrm{(}{\mathit{s}}_{\mathrm{b}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\right]\hspace{0.17em}\mathit{A,}$(7)where m is the blob mass, A is the effective cross-section of the blob, and g(s) is the solar gravitational acceleration along a semi-circular loop, that is, g(s) = g_⊙ cos(πs/L). It is measured at the blob’s centre of mass. In order to be able to solve the problem analytically, in the subsequent steps we neglect the stratification of the plasma. p_d(s_b−Δs) is the plasma pressure of the plasma below the blob, measured at the lower interface with the blob at s = s_b−Δs; p_u(s_b + Δs) is the plasma pressure of the plasma above the blob, measured at the upper interface with the blob at s = s_b + Δs.

Fig. 8 Bending of the magnetic field lines below the condensation at t = 1523 s for the rbc = 103, B = 100 G case.

Fig. 9 Solutions for the equilibrium position normalised by the loop length (left), and for the real and imaginary part of the normalised frequency of the blob oscillations (centre and right) as a function of loop-to-rain plasma mass ratios.

The blob is in equilibrium at position s₀. We assume that there is no exchange of mass between the background plasma. Subsequently, the plasma masses above and below are conserved, and we may write the equilibrium densities as ${\mathit{\rho }}_{\mathrm{u}\mathrm{0}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{M}}_{\mathrm{u}}}{\mathrm{(}\mathit{L}\mathrm{-}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{A}}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{\rho }}_{\mathrm{d}\mathrm{0}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{M}}_{\mathrm{d}}}{\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{A}}\mathrm{·}$(8)This allows us to rewrite the equilibrium pressures as ${\mathit{p}}_{\mathrm{d}\mathrm{0}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{d}\mathrm{0}}{\mathit{M}}_{\mathrm{d}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{A}}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{p}}_{\mathrm{u}\mathrm{0}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{u}\mathrm{0}}{\mathit{M}}_{\mathrm{u}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}\mathrm{(}\mathit{L}\mathrm{-}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{A}}\mathrm{·}$(9)Then, the equilibrium position of the blob is the solution of the transcendental equation ${\mathit{g}}_{\mathrm{\odot }}\cos \left(\frac{\mathit{\pi }{\mathit{s}}_{\mathrm{0}}}{\mathit{L}}\right)\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{u}\mathrm{0}}{\mathit{M}}_{\mathrm{d}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}\mathrm{(}\mathit{L}\mathrm{-}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{m}}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{d}\mathrm{0}}{\mathit{M}}_{\mathrm{u}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{m}}\mathrm{·}$(10)We introduce the lower density scale height H and sound speed C_Sd defined as $\mathit{H}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{d}\mathrm{0}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}{\mathit{g}}_{\mathrm{\odot }}}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{C}}_{\mathrm{Sd}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\sqrt{\frac{\mathit{\gamma }{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{d}\mathrm{0}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}}}\mathit{,}$(11)and the following dimensionless variables

$\begin{array}{ccc}& & {\mathit{\theta }}_{\mathrm{0}}\mathrm{=}\frac{\mathit{\pi }{\mathit{s}}_{\mathrm{0}}}{\mathit{L}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }\mathit{\theta }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathit{\pi }\mathrm{\Delta }\mathit{s}}{\mathit{L}}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{\theta }}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\pi }\mathrm{-}\mathrm{\Delta }\mathit{\theta ,}\\ & & \mathit{h}\mathrm{=}\frac{\mathit{\pi H}}{\mathit{L}}\mathit{,}\hspace{0.17em}\hspace{0.17em}\mathit{\nu }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{T}}_{\mathrm{u}\mathrm{0}}}{{\mathit{T}}_{\mathrm{d}\mathrm{0}}}\mathrm{·}\end{array}$(12)Subsequently, Eq. (10) may be rewritten in dimensionless form as $\cos {\mathit{\theta }}_{\mathrm{0}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{{\mathit{M}}_{\mathit{d}}}{\mathit{m}}\hspace{0.17em}\frac{\mathit{h}}{\mathrm{(}{\mathit{\theta }}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{\theta }\mathrm{)}}\mathrm{+}\frac{{\mathit{M}}_{\mathit{u}}}{\mathit{m}}\hspace{0.17em}\frac{\mathit{\nu h}}{\mathrm{(}{\mathit{\theta }}_{\mathrm{0}}\mathrm{-}{\mathit{\theta }}_{\mathrm{1}}\mathrm{)}}\mathrm{·}$(13)We linearise the equation of motion by considering small amplitude oscillations around the equilibrium position s₀, found by solving Eq. (10), such that s = s₀ + s₁ with | s₁ | ≪ L. Hence, $\begin{array}{ccc}\mathit{m}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{s}}_{\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}& \mathrm{=}& \mathrm{-}\mathit{m}\hspace{0.17em}\frac{\mathrm{d}\mathit{g}}{\mathrm{d}\mathit{s}}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right)}\hspace{0.17em}{\mathit{s}}_{\mathrm{1}}\\ & & \mathrm{-}\left[\frac{\mathit{\partial }{\mathit{p}}_{\mathrm{u}\mathrm{0}}}{\mathit{\partial s}}\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{+}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathrm{-}\frac{\mathit{\partial }{\mathit{p}}_{\mathrm{d}\mathrm{0}}}{\mathit{\partial s}}\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\right]\hspace{0.17em}{\mathit{s}}_{\mathrm{1}}\hspace{0.17em}\mathit{A}\\ & & \end{array}$(14)where p₁ is the linear perturbation of plasma pressure. The linear plasma displacement parallel to the equilibrium magnetic field ξ is governed by${\mathit{\rho }}_{\mathrm{0}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\xi }}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\frac{\mathit{\partial }{\mathit{p}}_{\mathrm{1}}}{\mathit{\partial s}}\mathit{,}$(15)where ρ₀ is the equilibrium plasma density, which we assumed to be uniform. Furthermore, this displacement satisfies the boundary conditions $\mathit{\xi }\mathrm{\left(}\mathrm{0}\mathit{,t}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\hspace{0.17em}\mathit{\xi }\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s,t}\mathrm{)}\mathrm{=}{\mathit{s}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{,}\hspace{0.17em}\mathit{\xi }\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right)}\mathrm{+}\mathrm{\Delta }\mathit{s,t}\mathrm{)}\mathrm{=}{\mathit{s}}_{\mathrm{1}}\mathrm{(}\mathit{t}\mathrm{)}\mathit{,}$(16)and is allowed to propagate in the upper region. We find $\mathit{\xi }\mathrm{\left(}\mathit{s,t}\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\{\begin{array}{c}\\ {\mathit{s}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}\frac{\sin \mathrm{\left(}{\mathit{k}}_{\mathrm{d}}\mathit{s}\mathrm{\right)}}{\sin \mathrm{\left(}{\mathit{k}}_{\mathrm{d}}\mathrm{\right(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}}& \mathrm{0}\mathrm{\le }\mathit{s}\mathrm{\le }{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\\ {\mathit{s}}_{\mathrm{1}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}\frac{\exp \mathrm{\left(}\mathrm{i}{\mathit{k}}_{\mathrm{u}}\mathit{s}\mathrm{\right)}}{\exp \mathrm{\left(}\mathrm{i}{\mathit{k}}_{\mathrm{u}}\mathrm{\right(}{\mathit{s}}_{\mathrm{0}}\mathrm{+}\mathrm{\Delta }\mathit{s}\mathrm{\left)}\mathrm{\right)}}& \mathit{s}\mathrm{\ge }{\mathit{s}}_{\mathrm{0}}\mathrm{+}\mathrm{\Delta }\mathit{s,}\end{array}$(17)where k_d and k_u are the wave numbers in the lower and upper regions, respectively. The corresponding pressure perturbation is found from Eq. (15): ${\mathit{p}}_{\mathrm{1}}\mathrm{\left(}\mathit{s,t}\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\{\begin{array}{c}\\ \frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{s}}_{\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathrm{1}}{{\mathit{k}}_{\mathrm{d}}}\hspace{0.17em}\frac{{\mathit{M}}_{\mathrm{d}}}{\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{A}}\frac{\cos \mathrm{\left(}{\mathit{k}}_{\mathrm{d}}\mathit{s}\mathrm{\right)}}{\sin \mathrm{\left(}{\mathit{k}}_{\mathrm{d}}\mathrm{\right(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{\left)}\mathrm{\right)}}\\ \frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{s}}_{\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathrm{i}}{{\mathit{k}}_{\mathrm{u}}}\hspace{0.17em}\frac{{\mathit{M}}_{\mathrm{u}}}{\mathrm{(}\mathit{L}\mathrm{-}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}\mathit{A}}\frac{\exp \mathrm{\left(}\mathrm{i}{\mathit{k}}_{\mathrm{u}}\mathit{s}\mathrm{\right)}}{\exp \mathrm{\left(}\mathrm{i}{\mathit{k}}_{\mathrm{u}}\mathrm{\right(}{\mathit{s}}_{\mathrm{0}}\mathrm{+}\mathrm{\Delta }\mathit{s}\mathrm{\left)}\mathrm{\right)}}\mathrm{·}\end{array}$(18)Equations (9) and (18) are substituted into Eq. (14): $\begin{array}{ccc}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{s}}_{\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}& \mathrm{=}& \frac{\mathit{\pi }{\mathit{g}}_{\mathrm{\odot }}}{\mathit{L}}\hspace{0.17em}\sin \left(\frac{\mathit{\pi }{\mathit{s}}_{\mathrm{0}}}{\mathit{L}}\right)\hspace{0.17em}{\mathit{s}}_{\mathrm{1}}\\ & & \mathrm{-}\frac{{\mathit{M}}_{\mathrm{d}}}{\mathit{m}}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{d}\mathrm{0}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}}\hspace{0.17em}\frac{{\mathit{s}}_{\mathrm{1}}}{\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{+}\frac{{\mathit{M}}_{\mathrm{d}}}{\mathit{m}}\hspace{0.17em}\frac{\mathrm{1}}{{\mathit{k}}_{\mathrm{d}}}\hspace{0.17em}\frac{\mathrm{1}}{\mathrm{(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{)}}\cot \mathrm{\left(}{\mathit{k}}_{\mathrm{d}}\mathrm{\right(}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{s}}_{\mathrm{1}}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\\ & & \mathrm{-}\frac{{\mathit{M}}_{\mathrm{u}}}{\mathit{m}}\hspace{0.17em}\frac{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{u}\mathrm{0}}}{\begin{array}{c}\mathrm{˜}\\ \mathit{m}\end{array}}\hspace{0.17em}\frac{{\mathit{s}}_{\mathrm{1}}}{\mathrm{(}\mathit{L}\mathrm{-}{\mathit{s}}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{s}{\mathrm{)}}^{\mathrm{2}}}\\ & & \end{array}$(19)for which normal mode solutions of the form s₁(t) ~ exp(−iωt) are sought. Equation (19) then turns into a dispersion relation for ω. We introduce the additional dimensionless variables $\mathrm{\Omega }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathit{\omega H}}{{\mathit{C}}_{\mathrm{Sd}}}\mathit{,}\hspace{0.17em}\hspace{0.17em}{\mathit{K}}_{\mathrm{d}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathit{k}}_{\mathrm{d}}\mathit{H,}\hspace{0.17em}\hspace{0.17em}{\mathit{K}}_{\mathrm{u}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathit{k}}_{\mathrm{u}}\mathit{H}\hspace{0.17em}\sqrt{\mathit{\nu }}\mathit{.}$(20)Equation (19) becomes in dimensionless form $\begin{array}{ccc}{\mathrm{\Omega }}^{\mathrm{2}}& \mathrm{=}& \frac{\mathit{h}}{\mathit{\gamma }}\left[\mathrm{-}\sin {\mathit{\theta }}_{\mathrm{0}}\mathrm{+}\frac{{\mathit{M}}_{\mathit{d}}}{\mathit{m}}\frac{\mathit{h}}{\mathrm{(}{\mathit{\theta }}_{\mathrm{0}}\mathrm{-}\mathrm{\Delta }\mathit{\theta }{\mathrm{)}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{M}}_{\mathit{u}}}{\mathit{m}}\frac{\mathit{\nu h}}{\mathrm{(}{\mathit{\theta }}_{\mathrm{0}}\mathrm{-}{\mathit{\theta }}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{2}}}\right]\\ & & \end{array}$(21)

Here ℜe(ω) determines the angular frequency of the blob oscillations and −1/ℑm(ω) sets the e-folding time for the damping due to wave radiation. Lastly, a dispersion relation is required in the two plasma regions to be able to connect K_d and K_u with ω. For a slow magnetoacoustic sausage mode, we find k in each region, defined by the solution to the dispersion relation (Edwin & Roberts 1983): $\frac{{\mathit{\kappa }}_{\mathrm{i}\mathit{s}}}{\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{v}}_{\mathrm{Ai}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{)}}\frac{{\mathrm{ℐ}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\kappa }}_{\mathrm{i}\mathit{s}}\mathit{a}\mathrm{\right)}}{{\mathrm{ℐ}}_{\mathrm{0}}\mathrm{\left(}{\mathit{\kappa }}_{\mathrm{i}\mathit{s}}\mathit{a}\mathrm{\right)}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\left(\frac{{\mathit{\rho }}_{\mathrm{0}\mathrm{e}\mathit{s}}}{{\mathit{\rho }}_{\mathrm{0}\mathrm{i}\mathit{s}}}\right)\frac{{\mathit{\kappa }}_{\mathrm{e}\mathit{s}}}{\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{v}}_{\mathrm{Ae}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{)}}\frac{{\mathrm{𝒦}}_{\mathrm{0}}^{\mathrm{\prime }}\mathrm{\left(}{\mathit{\kappa }}_{\mathrm{e}\mathit{s}}\mathit{a}\mathrm{\right)}}{{\mathrm{𝒦}}_{\mathrm{0}}\mathrm{\left(}{\mathit{\kappa }}_{\mathrm{e}\mathit{s}}\mathit{a}\mathrm{\right)}}\mathit{,}$(22)with the squared radial wave number ${\mathit{\kappa }}_{\mathit{ps}}^{\mathrm{2}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{C}}_{\mathrm{S}\mathit{ps}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{v}}_{\mathrm{A}\mathit{p}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{)}}{\mathrm{(}{\mathit{C}}_{\mathrm{S}\mathit{ps}}^{\mathrm{2}}\mathrm{+}{\mathit{v}}_{\mathrm{A}\mathit{p}}^{\mathrm{2}}\mathrm{)}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{c}}_{\mathrm{T}\mathit{ps}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{)}}\mathit{,}$(23)where $\mathit{a}\mathrm{=}\sqrt{\mathit{A}}\mathit{/}\mathit{\pi }$ is the loop cross-section radius, v_A the Alfvén speed, c_T the tube speed, and s ∈ { d,u }, p ∈ { i,e }. index i (e) refers to internal (external) conditions to the loop. The density contrast and Alfvén speed are assumed to be identical in the lower and upper regions. ℐ₀(x) and are the modified Bessel functions of the first and second kind, respectively. We further note that the oscillation seen in the simulations is essentially a slow magnetoacoustic sausage mode. We solve for the fundamental radial mode with the phase speed in the interval [c_T,c_S]. For ka ≪ 1 the solution is approximately k ≈ ω/C_Ti. This is also true for a slab geometry as used in the numerical simulations. Furthermore, for the range of values of density, temperature, and magnetic field strength, the tube speed varies from the sound speed by less than 10%. Therefore, it is reasonable to describe the mode with the dispersion relation of a one-dimensional acoustic mode with k = ω/C_Si. Then, K_d = K_u = Ω. Equation (21) is solved numerically together with the equilibrium Eq. (13) and the corresponding dispersion relation. We solve Eqs. (13) and (21) for a range of loop-to-rain mass ratios. There is a discontinuity in the solution for the equilibrium position of blob from one loop leg to the other in the case of a high mass of plasma that is confined below the blob and a low mass of the plasma in the rest of the loop (Fig. 9). This discontinuity then further propagates into the solutions for ℜe(ω) and ℑm(ω), resulting in discontinuity in the gradient. This is not likely to occur in a real gravitationally stratified loop that is initially symmetric, however, unless there is a direct mass injection occurring into one loop leg alone. We further focus on the behaviour of the oscillation parameters in limit cases. In high M/m limit (no coronal rain) both ℜe(ω) and ℑm(ω) increase linearly with M/m. In low M/m limit (no coronal plasma) ℜe(ω) decreases with $\sqrt{\mathit{M}\mathit{/}\mathit{m}}$ while ℑm(ω) remains constant (Fig. 10). Assuming realistic values of the loop-to-rain mass ratio are of the order of 1–10, the corresponding solutions for ℜe(ω) and ℑm(ω) are of the order of 0.001, or equivalently 1000 s for the period and damping scaling time.

Fig. 10 Horizontal (left) and vertical (right) cuts trough plots of the dependence of the angular frequency on the loop-to-rain mass ratio. Dashed lines mark the limit cases.

The comparison of blob oscillation periods determined from the simulation and periods predicted by the analytical model is shown in Table 1. We determine the loop-to-rain mass ratios (serving as an input for the analytical model) from the final simulation snapshot at t = 11 200 s, which we assume to be the best representation of the equilibrium state. The blob mass is determined by integrating the plasma density between blob boundaries, and the mass of the coronal loop plasma above and below the blob is determined by integrating the density between the blob boundary and loop footpoints while excluding the high-density chromosphere layer. Estimates of uncertainties in the parameters predicted by the analytical model are determined assuming 20% uncertainty on the position of blob boundaries. The agreement between the two is best for higher blob density and for high values of the magnetic field strength. This is as expected given the limitations of the analytical model. It should be noted that the analytical model considerably overestimates the height of the equilibrium position in the case of the lowest blob density. The equilibrium position predicted by the analytical model is heavily dependent on the input loop-to-rain mass ratios. In the low-density case, these are inherently more difficult to determine accurately from the simulation because we lack a well-defined upper boundary of the plasma blob. Here the elongated tail of the blob accounts for a higher fraction of the total blob mass than in the higher density cases. The low-density blob is also more sensitive to sound waves that are reflected from the boundaries (the analytical model neglects the presence of the upper domain boundary and assumes a radiating solution above the blob). The validity of the analytical periods and damping times in the r_bc = 10² case is therefore also limited.

Similarly, the agreement between the analytical model and the simulation is worse for cases with lower magnetic field strength when the plasma below the blob is less constrained in the transverse direction and hence allowed to expand, whereas the model explicitly assumes a constant loop cross-section. This also means that while the blob oscillates, transfer of plasma can occur from the lower loop leg to the region above the blob, thus invalidating the assumption of a piston-like behaviour. However, when we use the values corresponding to the case with B = 100 G and r_bc = 10³, which best adheres to assumptions made by the analytical model, this results in the predicted period of 1393 s and a damping scaling time of 1695 s, which corresponds to aboutthree clearly observable oscillation periods. This is in good agreement with the simulation results.

5. Discussion and conclusions

We carried out 2.5 MHD simulations of the dynamics of cool plasma condensations in a gravitationally stratified coronal loop. The motion and evolution of plasma condensations were found to be strongly affected by the pressure of the coronal loop plasma, and the pressure gradients can be high enough to account for the lower-than-free-fall speed of the coronal rain even in finite magnetic field cases. The fastest downward velocities are in agreement with recent coronal rain observations. High coronal magnetic field strength or a low mass of the condensations can lead to oscillatory motion consisting of multiple rebounds damped through sound wave emission, with the condensation eventually settling in an equilibrium position supported by the pressure of the underlying plasma. Rebounding of the condensation is due to a combined effect of the pressure gradient force and the magnetic tension force that results from bending of the magnetic field lines in the lower part of the coronal loop. The period and damping scaling time of the oscillatory motion are consistent with values determined using an analytical model for the balance of forces that act on the condensation.

Although the majority of coronal rain condensations are observed to fall directly towards the solar surface, the individual blobs are sometimes observed to longitudinally oscillate up and down before falling (Kohutova & Verwichte 2016). This loss of equilibrium has not been accounted for in our simulations and could be due to the change in mass of the coronal loop plasma that supports the blob or due to presence of other condensations. It should further be noted that in the non-equilibrium scenario, siphon flows caused by pressure differences in the loop can significantly affect the motion of the condensations, sometimes completely overriding the effects of the plasma pressure gradient and magnetic tension force addressed here.

The analytical model also highlights the fact that the dynamics of the plasma condensations (i.e. presence or lack of oscillatory motion and oscillation parameters) is determined by the loop-to-rain mass ratio. There is still considerable uncertainty about what fraction of the total mass of the coronal loop plasma condenses into coronal rain after catastrophic cooling takes place; current estimates of the loop-to-rain mass ratio from observations are in the order of 1–10 (Antolin et al. 2015). These estimates are subject to the spatial resolution limits of the instruments, however, it is therefore likely that a significant fraction of the condensation mass remains undetected.

The longitudinal oscillations of the coronal rain blobs were typically observed in transversely oscillating coronal loops (Kohutova & Verwichte 2016; Verwichte et al. 2017). This means that the action of the ponderomotive force that is due to transverse oscillations should be taken into account (Terradas & Ofman 2004). It has been proposed that the ponderomotive force can in fact affect the motion of the coronal rain condensations; however, typical amplitudes of the transverse loop oscillations are not sufficient to fully explain the observed oscillatory motion and sub-ballistic fall rates of coronal rain on their own (Verwichte et al. 2017). The ponderomotive force may still play a non-negligible role in the condensation dynamics, however, in addition to the effects of the coronal plasma pressure gradients and magnetic field effects addressed in this work. This is further supported by the fact that the oscillatory behaviour of coronal rain is usually observed near the loop top, and it suggests that the force that counteracts the motion under the gravity has a maximum near the loop apex, whereas the pressure gradient force was found to have greatest effect on the condensations in the lower part of the loop legs. The dynamics of plasma condensations in a transversally oscillating loop will be addressed in detail in future work.

Acknowledgments

P.K. acknowledges the support of the UK STFC Ph.D. studentship. E.V. acknowledges financial support from the UK STFC on the Warwick STFC Consolidated Grant ST/L000733/I.

References

All Tables

All Figures

	Fig. 1 Initial effective gravity, density, and temperature profiles as a function of s at x = 0.
In the text

	Fig. 2 Evolution of the density profile along the bottom half of the loop plotted every 280 s during the first 500 time steps.
In the text

	Fig. 3 Time-distance plot of the density along s at x = 0 for the rbc = 103, B = 100 G case.
In the text

	Fig. 4 Height (top), velocity (middle), and acceleration (bottom) profiles of the condensation for different values of blob density contrast and magnitude of the magnetic field strength. The dotted line shows the free-fall height and velocity profiles.
In the text

	Fig. 5 Dependence of the equilibrium height (left), maximum downward velocity (centre), and average period of blob oscillations (right) on the magnitude of the magnetic field strength for various blob densities.
In the text

	Fig. 6 Evolution of the force balance along the bottom half of the loop plotted every 280 s during the first 500 time steps. We show the vertical components of the plasma pressure force (blue), gravity (green), magnetic pressure force (red), magnetic tension force (turquoise), and total net force in the vertical direction (black dashed line).
In the text

	Fig. 7 Evolution of the magnitude of the magnetic field strength (top) and velocity divergence (bottom) averaged over the region below the plasma blob for the rbc = 103, B = 100 G case.
In the text

	Fig. 8 Bending of the magnetic field lines below the condensation at t = 1523 s for the rbc = 103, B = 100 G case.
In the text

	Fig. 9 Solutions for the equilibrium position normalised by the loop length (left), and for the real and imaginary part of the normalised frequency of the blob oscillations (centre and right) as a function of loop-to-rain plasma mass ratios.
In the text

	Fig. 10 Horizontal (left) and vertical (right) cuts trough plots of the dependence of the angular frequency on the loop-to-rain mass ratio. Dashed lines mark the limit cases.
In the text