© The Authors 2025

1. Introduction

Coronal rain is primarily observed in post-flare loops and active regions of the solar corona and plays a crucial role in understanding coronal heating, plasma dynamics, and mass transport within the solar atmosphere (Jejčič et al. 2018; Heinzel et al. 2020). Chromospheric evaporation initiates hot plasma inflow during a solar flare, while coronal rain marks the cooling and condensation outflow, completing the mass-energy cycle (Şahin & Antolin 2024). Single fluid magnetohydrodynamic (MHD) simulations of post-flare rain only recently progressed from 2D (Ruan et al. 2021) to 3D settings (Ruan et al. 2024), while more quiescent rain in loops and arcades started in multidimensional MHD in chromosphere to coronal settings (Fang et al. 2013, 2015; Moschou et al. 2015; Xia et al. 2017; Li et al. 2022). Active region rain was recently demonstrated to be linked to thermal instability in self-consistent radiative MHD, with sub-photosphere to coronal evolutions (Lu et al. 2024). Observations showed that coronal rain is a multithermal phenomenon, appearing as clumped structures organized in strands (Antolin et al. 2015). The density contrast across the blob interface is around 100 (Scullion et al. 2016). The coronal rain blobs have widths of 300–700 km and a broad velocity distribution, peaking below 50 km/h (Antolin et al. 2023). Higher velocities of the falling blobs, of up to 50–100 km/s, have also been observed (Mason et al. 2019). A two-fluid simulation of the evolution of coronal rain blobs showed a large decoupling in velocity due to the elastic collisions between ions and neutrals and associated frictional heating (Martínez-Gómez et al. 2022).

The coronal rain blobs appear as the nonlinear outcome of the thermal instability caused by radiative cooling of the plasma. The thermal instability was introduced conceptually by Parker (1953) and studied theoretically in the linear MHD assumption in the seminal paper by Field (1965). Local coronal volume simulations confirmed the thermal instability as the trigger for the formation of condensations in the MHD models with optically thin radiative cooling curves corresponding to the solar corona (Claes & Keppens 2019; Claes et al. 2020; Hermans & Keppens 2021; Jerčić & Keppens 2023). As the rain blobs represent a typically hundredfold increase in density and decrease in temperature, it is clear that plasma-neutral effects can become important as the thermal instability proceeds. None of the in situ formation studies on coronal rain have addressed this aspect in detail, although a recent 1D setup initiated this in fixed loop settings (Jerčić et al. 2025).

Analytical studies in the two-fluid plasma-neutral approach show that the presence of neutrals decreases the linear growth rate of the thermal instability when there is a larger fraction of neutrals or when they are coupled more strongly to the charges (Popescu Braileanu & Keppens 2024). The decrease in temperature inside the condensed regions triggers recombination of the charges, and the fraction of neutrals increases so much so that the contrast in the neutral density is orders of magnitude larger than that in the charged density (Popescu Braileanu & Keppens 2024). Similar important ionization-recombination processes have been identified in reconnection scenarios, or in the presence of shocks. Indeed, they stabilize the current sheets and decrease the growth of the coalescence instability in the nonlinear phase (Murtas et al. 2022). The inclusion of inelastic processes and the ionization potential in a two-fluid approximation can lead to a significant departure from MHD results in slow shock formation (Snow & Hillier 2021). In this paper, we perform 3D two-fluid simulations in which reconnection and thermal instability interplay, as is relevant for coronal rain.

Coronal rain was observed to be formed as a consequence of the impulsive heating following reconnection events in the solar corona (Kohutova et al. 2019) and it is frequently observed in post-flare loops (Mason et al. 2019). Unlike 2D reconnection, which happens only at “X”-type null points, 3D reconnection happens in regions of electric field concentrations and at null points, but also at separators and quasi-separators (Li et al. 2021; Pontin & Priest 2022). The reconnection rates in a 3D setup are usually smaller than in 2.5D and they converge in the limit of weak reconnection with a strong guide field (Daldorff et al. 2022).

Null points are ubiquitous in the solar corona, where small patches of magnetic field emerge through the photosphere and have an opposite polarity to the main surrounding field. In 3D the magnetic field around null points include isolated field lines, called “spines” and a surface of “fan” field lines, both structures starting or ending near the null point (Li et al. 2021; Pontin et al. 2016). The 3D null points have been studied in a linear field topology approximation around the null point, where the magnetic field vector is proportional to the position vector measured from the null point (Pontin et al. 2007; Priest & Pontin 2009). The simplest null point in this approximation, usually considered in the analytical studies, is the “proper radial null point,” where the spine is perpendicular to the fan and the field lines in the fan are radial, a particular case of the general linear null point when the spine is not perpendicular to the fan and the field lines in the fan can be curved (Li et al. 2021; Priest & Pontin 2009; Pontin et al. 2007, 2016).

Small incompressible shearing flows perturbing the spine of the radial null point create pure fan reconnection (Pontin et al. 2016), but generally the shearing driving produces spine-fan reconnection (Priest & Pontin 2009), which is the most common 3D null point reconnection (Li et al. 2021). The Lorentz force tends to increase the displacement created by the driving, producing a collapse of the spine toward the fan. As a result, the current sheet formed near the null point is inclined with respect to the fan plane with an angle that increases with increasing driving velocity; therefore, the current sheet contains parts of both spine and fan (Pontin et al. 2007). The current sheet, formed in the vicinity of the null point, becomes more localized and more intense while the null point collapses (Priest & Pontin 2009).

In many cases the spine is an open field line, and this situation usually occurs in coronal hole regions (Mason et al. 2019). The jets produced by reconnection in this spine-fan configuration have sizes varying from less than 10 Mm to hundreds of megameters and condensations have been observed for almost all the intermediate size jets of ≈100 Mm (Mason et al. 2019). The smaller-scale condensations (of sizes less than 10 Mm) are difficult to observe due to the instrument resolution and the inherent line-of-sight ambiguity for determining both velocity and size (Mason et al. 2019). Simulations using a 3D spine-fan magnetic reconnection setup provide a realistic representation of solar magnetic activity. Velocity driving of the bottom boundary in these simulations produces an accumulation of current density around the null point (Pontin & Galsgaard 2007; Galsgaard & Pontin 2011; Pontin et al. 2013; Yadav et al. 2022). The magnetic field reconnects in the region around the null point, leading to an intense localized energy release. Plasma flows develop along the field lines, guided by the spine and fan structures, transporting energy and momentum away from the reconnection site (Pontin et al. 2013). Observations showed slippage along flare ribbons following flare events moving in opposite directions, with velocities usually around 20–40 km/s, but also reaching 400–450 km/s (Li et al. 2021). The reconnection rate, quantified by the parallel electric field integrated along the field lines in the diffusion region around the null point, has a peak value for the field lines situated in the fan plane (Galsgaard & Pontin 2011). The current density accumulated near the null point is sensitive to the driving (Galsgaard & Pontin 2011) The reconnection rate time evolution follows the pattern of the current density evolution: both decrease when the driving stops, and the system relaxes (Al-Hachami & Pontin 2010; Pontin et al. 2013).

In this paper, we study the formation of coronal rain during the late evolution of spine-fan 3D magnetic reconnection. We use a two-fluid approach with ionization recombination processes taken into account. Section 2 describes the numerical setup. The results are presented in Section 3 and the conclusions in Section 4.

2. Numerical setup

The 3D computational domain corresponds to a coronal-only region of size −6.8 Mm ≤x≤3 Mm, −3 Mm ≤y≤3 Mm, and 0≤z≤3 Mm and is covered by a grid with a base resolution of 144×80×64 and three levels of refinement. Therefore, the effective resolution is 576×320×256 points and the cell sizes are dx=17.014 km, dy=18.75 km, and dz=11.719 km. Hence, we effectively resolve the structure down to about a few tens of kilometers, compared to the observationally inferred blob sizes on the order of a few hundred kilometers (Vashalomidze et al. 2015; Antolin et al. 2023; Scullion et al. 2016).

The spine-fan magnetic field configuration is very similar to the setup used by Pontin et al. (2013). The magnetic field is a superposition of five magnetic point charges situated below the bottom boundary, resulting in a divergence-free potential magnetic field:

$$\begin{array}{cc}\hfill N_p& =5,\hfill \\ \hfill x_p& =\{-4.1,-3.9,-1.9,0.1,0.3\},\hfill \\ \hfill {\mathit{y}}_p& =\{0.0,0.0,0.0,0.0,0.0\},\hfill \\ \hfill z_p& =\{-0.4,-1.0,-1.0,-1.0,-0.4\},\hfill \\ \hfill b_p& =\{0.3,-1.0,1.0,-1.0,0.3\},\hfill \\ \hfill {B_0}_x(x,\mathit{y},z)& =B_{00}{\displaystyle \sum _{j=1}^{N_p}}{b_p}_j\frac{x-{x_p}_j}{r_j^3},\hfill \\ \hfill {B_0}_{\mathit{y}}(x,\mathit{y},z)& =B_{00}{\displaystyle \sum _{j=1}^{N_p}}{b_p}_j\frac{\mathit{y}-{{\mathit{y}}_p}_j}{r_j^3},\hfill \\ \hfill {B_0}_z(x,\mathit{y},z)& =B_{00}{\displaystyle \sum _{j=1}^{N_p}}{b_p}_j\frac{z-{z_p}_j}{r_j^3},\hfill \\ \hfill \text{where}\hspace{0.5em}{r}_j& =\sqrt{(x-{x_p}_j{)}^2+(\mathit{y}-{{\mathit{y}}_p}_j{)}^2+(z-{z_p}_j{)}^2}.\hfill \end{array}$$(1)

We used the value B₀₀=100 G so we have a typical active region field at the bottom of our coronal box. The magnetic field configuration is symmetric with respect to x=−1.9, with two magnetic null points located at (x_N1,y_N,z_N) and (x_N2,y_N,z_N), where x_N1=0.123 ,y_N=0 ,z_N=0.42,x_N2=−3.923. The configuration of the magnetic field is shown in Figure 1, with the two null points indicated by two small spheres. These nulls are preferential sites for reconnection events, and we enhanced this possibility by adopting an anomalous, spatially localized resistivity.

Fig. 1. Initial configuration of the magnetic field. The vertical component of the magnetic field is plotted at the bottom of the domain. The two null points are indicated by small spheres. The field lines are plotted using a gray color map, which shows the magnitude of the magnetic field. The orange circle below one of the null points indicates the location of the velocity perturbation defined by Eq. (3), which has a center at the projection of the perturbed null point on the bottom boundary and a radius equal to $2\sqrt{{\sigma }_v/2}=0.447$ Mm. The isocontours of βplasma=1 are shown by a semiopaque violet mesh of two regions, one encapsulating the two null points and an outer region.

The resistivity localized around the two null points has a profile very similar to the one used in Popescu Braileanu & Keppens (2023):

$$\begin{array}{cc}\hfill \eta (x,\mathit{y},z)& =({\eta }_0-{\eta }_1)[\text{exp}\left(-\frac{(x-x_{\mathrm{N}1}{)}^2+(\mathit{y}-{\mathit{y}}_N{)}^2+(z-z_N{)}^2}{2{\sigma }_r}\right)+\hfill \\ \hfill & \text{exp}\left(-\frac{(x-x_{\mathrm{N}2}{)}^2+(\mathit{y}-{\mathit{y}}_N{)}^2+(z-z_N{)}^2}{2{\sigma }_r}\right)]+{\eta }_1,\hfill \end{array}$$(2)

where η₀=8 Ωm, η₁=0.8 Ωm, and σ_r=0.2 Mm².

We used open boundary conditions in the x and y directions and for the top boundary in the z direction. For the bottom z boundary, we used a driver during the whole period of the simulation for the x component of the velocity, centered around the projection of the first null point:

$$V_x=V_0{\text{sin}}^2\left(\frac{2\pi t}{P}\right)\text{exp}\left(-\frac{(x-x_{\mathrm{N}1}{)}^2+(\mathit{y}-{\mathit{y}}_N{)}^2}{{\sigma }_{\mathit{v}}}\right),$$(3)

where V₀=10 km/s, P=100 s, and σ_v=0.1 Mm². The location of the driver is indicated by an orange circle in Figure 1. This driving is adopted for both neutrals and charges. The y and z components of the velocities on the bottom boundary are antisymmetric, which implies that they vanish on the boundary and all the other variables are symmetric. This bottom boundary condition is similar to the one used by Pontin et al. (2013). The simulation was run for 598 s.

We solved numerically the nonlinear two-fluid equations, also shown in Popescu Braileanu & Keppens (2022), which in this case (without gravity and Hall term, but including the radiative losses) are:

$$\frac{\partial {\rho }_{\mathrm{n}}}{\partial t}+\nabla ·\left({\rho }_{\mathrm{n}}{\text{v}}_{\mathrm{n}}\right)=S_{\mathrm{n}},$$(4)

$$\frac{\partial {\rho }_{\mathrm{c}}}{\partial t}+\nabla ·\left({\rho }_{\mathrm{c}}{\text{v}}_{\mathrm{c}}\right)=-S_{\mathrm{n}},$$(5)

$$\frac{\partial ({\rho }_{\mathrm{n}}{\text{v}}_{\text{n}})}{\partial t}+\nabla ·\left({\rho }_{\mathrm{n}}{\text{v}}_{\text{n}}{\text{v}}_{\text{n}}+p_{\mathrm{n}}I\right)={\text{R}}_{\mathrm{n}},$$(6)

$$\frac{\partial ({\rho }_{\mathrm{c}}{\text{v}}_{\text{c}})}{\partial t}+\nabla ·\left[{\rho }_{\mathrm{c}}{\text{v}}_{\mathrm{c}}{\text{v}}_{\mathrm{c}}+\left(p_{\mathrm{c}}+\frac{1}{2}{B}^2\right)I-\text{B}\text{B}\right]=-{\text{R}}_{\mathrm{n}},$$(7)

$$\frac{\partial e_{\mathrm{n}}^{\mathrm{tot}}}{\partial t}+\nabla ·\left[{\text{v}}_{\mathrm{n}}\left(e_{\mathrm{n}}^{\mathrm{tot}}+p_{\mathrm{n}}\right)\right]=M_{\mathrm{n}},$$(8)

$$\frac{\partial e_{\mathrm{c}}^{\mathrm{tot}}}{\partial t}+\nabla ·[{\text{v}}_{\mathrm{c}}\left(e_{\mathrm{c}}^{\mathrm{tot}}+p_{\mathrm{c}}+\frac{1}{2}{B}^2\right)-\text{B}({\text{v}}_{\mathrm{c}}·\text{B})]$$

$$=\eta J^2-M_{\mathrm{n}}-{\rho }_c^2\mathrm{\Lambda }+{\rho }_{\mathrm{c}0}^2{\mathrm{\Lambda }}_0,$$(9)

$$\frac{\partial \text{B}}{\partial t}+\nabla ·\left({\text{v}}_{\mathrm{c}}\text{B}-\text{B}{\text{v}}_{\mathrm{c}}\right)=\eta \text{J}.$$(10)

The above set of equations (4)–(10) was written for conserved variables: mass densities, momentum, and total energies:

$$e_{\mathrm{c}}^{\mathrm{tot}}=e_{\mathrm{c}}+\frac{1}{2}{\rho }_{\mathrm{c}}{v}_{\mathrm{c}}^2+\frac{1}{2}{B}^2;e_{\mathrm{n}}^{\mathrm{tot}}=e_{\mathrm{n}}+\frac{1}{2}{\rho }_{\mathrm{n}}{v}_{\mathrm{n}}^2.$$(11)

The pressure was prescribed using the ideal equation of state,

$$p_i=(\gamma -1)e_i,\text{for}\hspace{0.5em}i=\mathrm{n},\mathrm{c}.$$(12)

The collisional terms are:

$$\begin{array}{cc}\hfill S_{\mathrm{n}}& ={\rho }_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}},\hfill \\ \hfill {\text{R}}_{\mathrm{n}}& ={\rho }_{\mathrm{c}}{\text{v}}_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{\text{v}}_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}}+{\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}\alpha ({\text{v}}_{\mathrm{c}}-{\text{v}}_{\mathrm{n}}),\hfill \\ \hfill M_{\mathrm{n}}& =\frac{1}{2}{\rho }_{\mathrm{c}}{\mathit{v}}_{\mathrm{c}}^2{\mathrm{\Gamma }}^{\mathrm{rec}}-\frac{1}{2}{\rho }_{\mathrm{n}}{\mathit{v}}_{\mathrm{n}}^2{\mathrm{\Gamma }}^{\mathrm{ion}}+\frac{1}{2}({{\mathit{v}}_{\mathrm{c}}}^2-{{\mathit{v}}_{\mathrm{n}}}^2){\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}\alpha \hfill \\ \hfill & +\frac{1}{\gamma -1}\left({\rho }_{\mathrm{c}}{T}_{\mathrm{c}}{\mathrm{\Gamma }}^{\mathrm{rec}}-{\rho }_{\mathrm{n}}{T}_{\mathrm{n}}{\mathrm{\Gamma }}^{\mathrm{ion}}\right)\hfill \\ \hfill & +\frac{1}{\gamma -1}(T_{\mathrm{c}}-T_{\mathrm{n}}){\rho }_{\mathrm{n}}{\rho }_{\mathrm{c}}\alpha .\hfill \end{array}$$(13)

The elastic collisional parameter is

$$\alpha =\frac{2}{{m_H}^{3/2}\sqrt{\pi }}\sqrt{k_B{T}_{\mathit{cn}}}{\mathrm{\Sigma }}_{\mathit{in}},$$(14)

with T_cn=(T_c+T_n)/2, and the ionization-recombination rates are

$${\mathrm{\Gamma }}^{\mathrm{rec}}\approx \frac{n_e}{\sqrt{T_e^{*}}}2.6·{10}^{-19}{\mathrm{s}}^{-1},$$(15)

$${\mathrm{\Gamma }}^{\mathrm{ion}}\approx n_{e}A\frac{1}{X+{\varphi }_{\mathrm{ion}}/T_e^{*}}{\left(\frac{{\varphi }_{\mathrm{ion}}}{T_e^{*}}\right)}^K{e}^{-{\varphi }_{\mathrm{ion}}/T_e^{*}}{\mathrm{s}}^{-1},$$(16)

where ϕ_ion=13.6eV, $T_e^{*}$ is the electron temperature in electronvolts, and constants have values of A=2.91·10⁻¹⁴, K = 0.39, and X = 0.232 in the SI unit system. The radiative cooling term, $-{\rho }_c^2\mathrm{\Lambda }$, is compensated for by a background heating: ${\rho }_{\mathrm{c}0}^2{\mathrm{\Lambda }}_0$, which cancels the cooling exactly at time t=0, the simulation starting from thermal equilibrium. For numerical reasons, the variables were split into time-dependent and time-independent quantities, the equivalent equations solved by the code being Eqs. (2.7)–(2.14) from Popescu Braileanu & Keppens (2024).

We used initially uniform thermodynamic background plasma parameters, characteristic to the solar corona, as was done in Popescu Braileanu & Keppens (2024): the coronal plasma pressure was set to the constant p₀=0.032 Pa, while the plasma density was given the constant value ρ₀=2.3×10⁻¹² kg/m³. This makes the coronal temperature have a value of around 10⁶ K. Since we started with a solar coronal volume, the ionization fraction corresponding to ionization-recombination equilibrium is very high, and has an almost negligible initial fraction of neutrals: 1−ξ_i=1.019×10⁻⁶. Therefore, we used the same initial conditions for thermodynamic variables and the same equations as in the simulation presented in Section 4a in Popescu Braileanu & Keppens (2024), but in a 3D setup and a different configuration of the force-free background magnetic field. As the magnetic field varies in 3D, the magnitude of the field defined in Eq. (1) has the maximum value of 99.55 G on the bottom boundary at the points (–4.1112,0,0) and (0.3112,0,0) and the minimum value of 0 at the two null points, resulting in an initial range of plasma beta from 8.1×10⁻⁴ to ∞. The isocontours of β_plasma=1 are shown in Figure 1 by a semiopaque isosurface revealing an inner region around the null points of about 200 km in size and an outer region, further away from the magnetic point sources.

The cooling curve used for the optically thin radiative losses is “Colgan_DM.” The equations were solved numerically with the code MPI-AMRVAC (Keppens et al. 2023; Popescu Braileanu & Keppens 2022). We used a finite volume approach for spatial discretization, employing an HLL Riemann solver (Toro 1997) with a “koren” limiter for variable reconstruction. We used the semi-implicit IMEX_ARS3 scheme for the time integration. The refinement criterium was based on the charged fluid density.

3. Results

3.1. Overall time evolution

Figure 2 shows the time evolution of the current density at the perturbed null point, the maximum values in the whole domain of the density of charges and neutrals and the decoupling in velocity (|v_n−v_c|) between the two species. Although the time cadence in Figure 2 is relatively coarse, the trends are obvious: the current density rises sharply at the null, and after several hundred seconds we witness the formation of coronal condensations in both charges and neutrals. Due to the bottom driving, the current density increases at the null point (panel a) shortly after the beginning of the simulation. The onset of thermal instability can be seen at around 300 s with a sudden increase in the maximum density of charges (panel b) and neutrals (panel d), when the decoupling in velocity (panel c) also reaches the maximum value. The maximum decoupling in velocity shows a pattern similar to the current density at the null point, suggesting its relationship to the magnetic activity. The maximum density of neutrals shows a much steeper increase by seven orders of magnitude, compared to almost two orders of magnitude that can be seen for the increase in the charge density. The timescale related to the periodicity of the velocity driver (50 s, being P/2 used in Eq. (3)) can be seen in the graphs of the evolution of the current density at the null point and the evolution of the velocity decoupling.

Fig. 2. Time evolution of: (a) current density at the null point; (b) logarithmic value of the maximum density of charges; (c) maximum value of the magnitude of the decoupling in velocity; (d) logarithmic value of the maximum density of neutrals.

3.2. Coronal rain blob

Since the magnetic null point is contained in the y=0 plane and the velocity driven at the bottom boundary only for the x component has the largest magnitude at y=0 (and x=x_N1), we expect most dynamics due to the magnetic reconnection to happen in the y=0 plane. Figure 3 shows a slice, y=0, of the charges density (left panel) and neutral density (right panel) at three times after the onset of the thermal instability: t=442 s, t=489.3 s, and t=532.8 s, capturing the formation and evolution of a condensation blob. The location where we encounter the maximum density for both neutrals and charges, located at x=−1.25 Mm and z=0.6 Mm, is shown by a blue and green dot in the charges and neutral panels, respectively, for the time t=489.3 s, when the fully formed blob is traveling through the domain.

Fig. 3. Slice at y=0 for the times t=442 s (top row), t=489.3 s (middle row), and t=532.8 s (bottom row), showing the formation and evolution of the rain blob. Left: Density of charges. The red arrows show the in-plane velocity of charges. Right: Density of neutrals. The red arrows show the in-plane velocity of neutrals. The in-plane magnetic field lines are overplotted by blue lines in all the panels. The perturbed null point is shown by a red dot, the other null point being outside the domain plotted in these panels. At the time t=489.3 s, when the blob is formed fully (middle row), the point with a high density, situated at x=−1.25 Mm and z=0.6 Mm, is marked in the charges and neutral density panels by a blue and green dot, respectively. This point is the same point marked in Figure 4 and Figure 6 below for the panels on the second row, for a better visualization. The line along the cut used in Figure 8 below is indicated by a dashed line: blue in the charge density panel and green in the neutral density panel.

In the left panels, we can observe field-aligned structures due to the reconnection at the null point. The velocities shown by red arrows in the left panels indicate the movement of the blob along the field lines from the reconnection point. At a subsequent time, t=532.8 s, the blob seems to have moved along the field lines and reaches the bottom boundary at x≈−1.75 Mm, where its further movement is slowed down by the zero-velocity bottom boundary conditions, similarly to a previously formed blob that can be seen at x≈1.25 Mm in all the panels. The two blobs have similar sizes and shapes.

Before the blob formation, at time t=442 s, we can observe a higher density above the location of the blob, indicating that most of the condensed material will come from above. The reconnection jets travel along the current sheets formed in the vicinity of the null point, and there is a spatial correlation between the charged density and the current density, shown in Figure 5 below for the time t=489.3 s. The density of charges seems to be higher along the spine; however, the recombination of charges show a depletion on the path traveled by the previously formed blob along the dome surface. We can also observe the spine becoming more inclined with time, consistent to the collapse of the null point during the shearing driving (Priest & Pontin 2009; Pontin et al. 2007).

The spatial contrast in the density of neutrals, much larger than that in the charged density, suggested by the limits of the color bars, indicates the important role of the recombination. Because of the recombination, the neutrals slow down and the charges speed up. We shall see how this effect, and the fact that the magnetic forces act on the charges only, affect the decoupling at different scales.

To infer the 3D shape of the blob, we looked at a x slice when it was formed. Figure 4 shows a slice of the charged density at x=−1.25 at time t=489.3 s, indicating that the condensation has sizes of the same order of magnitude of several hundred kilometers for all dimensions, consistent with the observations (Vashalomidze et al. 2015; Antolin et al. 2023). The condensations align perpendicularly with the field lines, especially when the magnetic field is strong (Hermans & Keppens 2021; Popescu Braileanu & Keppens 2024), shrinking in time in the direction parallel to the field lines. The blob is further deformed and it appears partially along the field lines, as is seen in the y=0 slice because of the high-velocity reconnection outflows in this plane, indicated by large red arrows, corresponding to velocities of around 50–60 km/s.

Fig. 4. Slice at x=−1.25 Mm for the time t=489.3 s, showing the density of charges. The in-plane magnetic field lines are overplotted by blue lines. The point with a high density, situated at y=0 and z=0.6 Mm, is marked by a blue dot, being the same point marked in Figure 3 above, and in Figure 6 below for the panels on the second row, for a better visualization. The line along the cut used in Figure 8 below is indicated by a dashed orange line.

The linear theory shows that the thermal instability growth rate is weakly sensitive to the wavelength of the perturbation (Field 1965); therefore, we expect the blobs to appear in the inverse order of their initial perturbation amplitude. The condensation blobs appear in the nonlinear stage of the thermal instability, as a result of the growth of the density perturbations created by reconnection. Therefore, the sizes of the initial perturbation are related to the sizes of the reconnection jets, which in turn are related to the width of the current sheets. Figure 5 shows the magnitude of the current density, in a base 2 logarithm, a decrease by 1 meaning that the current density magnitude dropped to half, making it easier to estimate the full width at half maximum, an indicator of the current sheet width. The current sheet contains part of the fan and spine, being stronger at the dome boundary along the separatrix surface than along the spine, similarly to the current density structures observed in Figure 9 in Pontin et al. (2013). The current sheet width seems to be similar to the size of the condensed structures in the direction perpendicular to the field lines. While the thermal instability evolves, the blobs shrink in the direction parallel to the field lines.

Fig. 5. Slice at y=0 for the time t=489.3 s of the magnitude of the current density. The quantity plotted is the base 2 logarithm of |J|, and the limits of the color bar of –13.5 and –9.5 correspond to |J|=8.6×10−5 A/m2 and |J|=1.38×10−3 A/m2, respectively. The in-plane magnetic field lines are overplotted by blue lines and the orange arrows indicate the in-plane decoupling velocity, the projection in the xz plane of the quantity vn−vc. The perturbed null point is shown by a red dot. The isocontours of the decoupling in temperature, Tn−Tc, are shown by green lines, indicating three values of 10, 30, and 50 kK (darker green for a larger value).

We overplot the in-plane decoupling in velocity (v_n−v_c projected in the xz plane) with orange lines and the isocontours of the decoupling in temperature (T_n−T_c) with green lines, showing three values of 10, 30, and 50 kK with progressively darker colors. The decoupling in velocity shows greater values around the blob in the bottom part of the domain, where the magnetic field is stronger. The neutral particles go inside the condensed structure across the magnetic field lines, unaffected by the Lorentz force. At the same time, we observe the highest decoupling in velocity above the previously formed blob (x≈1.25 Mm), on its path along the field lines toward the bottom. The large velocity decoupling, correlated spatially with a high current density (associated with a large Lorentz force) and with the depletion of the charged density, seen in the corresponding left panel of Figure 3, suggests that the decoupling is related to the recombination of the charges and the slippage of neutrals across the field lines. The decoupling in velocity produces frictional heating, which increases the temperature of the neutrals more than that of the charges, because of their lower density. This results in a difference of more than 50 kK at the location where the decoupling in velocity is maximum.

A better view of the 3D density structures can be seen in Figure 6, which shows the density of charges (left panels) and neutrals (right panels), integrated along the y direction for the three times shown before in Figure 3, around the time when the condensation formed in both charges and neutrals. The two null points are indicated by two red dots, and the location where we find the maximum density in the slice y=0, for both neutrals and charges at time t=489.3 s, situated at x=−1.25 Mm and z=0.6 Mm, is shown by a blue and green dot in the charges and neutral panels, respectively, similarly to Figures 3 and 4 above. The high-density rain blob surrounds this maximal density location. The charged density structures can be seen at the second null point as well, but are much less dense and have no condensations. The difference in the contrast of density between charges and neutrals can also be seen spatially. The integrated density of neutrals seems to have a much higher contrast, with sharp structures, compared to the charged density. No structure in the neutral density can be visually seen around the second null point. The long and thin 3D structures in the charged density, from which the blob condenses, can be better seen in a zoom-in around the high density point, shown in Figure 7. They are consistent with their 2D projection, the flows along the spine, seen in Figure 3.

Fig. 6. Evolution of the density integrated along the y direction around the time of the simulation when a condensation appear in the simulation. Left: Density of charges. Right: Density of neutrals. The magnetic null points are highlighted with small red circles. The point located at (–1.25,0,0.6) Mm, along the high density structure in the plane y=0, is marked for reference with green and blue, the same points highlighted in Figures 3 and 4 above.

Fig. 7. 3D integrated density of charges zoomed around the condensation point at time t=489.3 s, shown from a different angle than in Figure 6, as is indicated by the colored axes. The center of the condensation, located at (–1.25,0.0.6) Mm, is marked with a blue dot, and the perturbed null point is indicated by a red dot, consistent with Figures 3, 4, and 6.

3.3. Scales of ion-neutral decoupling

In order to better quantify the density contrast and the two-fluid effects, we looked at a cut along the line x=−1.25 Mm and y=0 for 0≤z≤1.2 Mm, containing the point with the highest density located at z=0.6 Mm, already mentioned and indicated in the previous figures. Panel (a) of Figure 8 shows the density of charges and neutrals, and the total density along this line. We can observe that the contrast in neutral density is seven orders of magnitude, indicating a major increase compared to its original coronal and negligible values. It is also more pronounced when compared to the contrast in the charges density, of two orders of magnitude. The latter is the typical density contrast between the corona and the interior regions of prominences and coronal rain. Still, throughout the blob we always have plasma at least one order of magnitude higher than neutral density. Therefore, the neutral density remains small compared to the charges’ density so that the total density contrast is that of the charges.

Fig. 8. Quantities along a vertical line cut at x=−1.25 Mm, y=0 for z between 0 and 1.2 Mm at time t=489.3 s. (a) Densities of charges (solid blue line), neutrals (dashed red line), and total (dotted black line). (b) Temperature of charges (solid blue line) and neutrals (dashed red line). (c)Vertical velocity of charges (solid blue line) and neutrals (dashed red line). (d) Decoupling in velocity, vnz−vcz (red line on the left axis), and decoupling in temperature, Tn−Tc (black line on the right axis). The three thin vertical reference lines correspond to: maximum density for all the three densities plotted: charges, neutrals, and total (green), maximum decoupling in both velocity and temperature (gray), and minimum temperature for charges and neutrals (violet).

The temperatures of charges and neutrals have a similar profile (panel b), with no visible difference, so plasma-neutral thermal coupling is near-perfect and the blob represents a hundredfold colder material than the corona, the typical observed contrast. The minimum peak in the temperatures (the vertical violet line) is slightly shifted up in the vertical direction from the maximum peak in the densities (the vertical green line). The distance between these two peaks is around 12 km, which is similar to the value of the collisional mean free path between ions and neutrals (as from Eq. (37) in Popescu Braileanu & Keppens 2022, where we used the charges fast magneto-acoustic speed), calculated at the center of the blob, suggesting that this misalignment could be explained by the two-fluid effects. Outside the blob the density of neutrals is very small, and therefore the collisional mean free path between ions and neutrals is very large (on the order 10⁹ km), orders of magnitude larger than the simulation domain. The higher density of charges makes the collisional mean free path between neutrals and ions (Eq. (37) in Popescu Braileanu & Keppens 2022, where we used the neutrals sound speed) smaller than that between ions and neutrals, being about 80 km outside the blob and reaching a value smaller than 1 km in the center of the condensation.

When the condensation forms, the charged and neutral particles accumulate inside the blob with velocities of about 10 km/s, which can be seen in the velocity profile, shown in panel (c). The whole structure moves downward, with the charges moving slightly faster everywhere, except for a small layer in the center of the structure. The downward velocity of about 10 km/s is similar to observations (Vashalomidze et al. 2015; Antolin et al. 2023).

The decoupling in velocity and temperature, shown in panel (d), has the same peak in magnitude (the vertical gray line), located toward the bottom, compared to the peak in the densities. The negative peak in the decoupling in velocity appears due to the movement of the neutrals across the field lines, when both neutrals and charges go toward the center of the moving blob, mostly from above. However, we can observe two smaller positive peaks in the decoupling at the sides, occurring at larger scales, when the neutrals are slowed by the recombination in the downward-moving blob. The frictional heating, proportional to the square of the decoupling in velocity, is much more efficient at heating the neutrals because of their lower density compared to that of the charges. This difference becomes reflected in the decoupling in temperature with the maximum located at the same point as the maximum decoupling in velocity, and the neutrals having higher temperatures than the charges. The maximum decoupling in both velocity and temperature is larger, over 10% of the values of velocity and temperature, respectively: about 1 km/second in velocity decoupling and a few tens of thousands of Kelvin in temperature. Because of the increase in temperature due to the frictional heating, the minimum temperature peak is shifted slightly toward the upper part. The radiative cooling also becomes more efficient for slightly higher temperatures, because of the dependence of the cooling curve on the temperature.

4. Conclusions

In this paper, we have demonstrated, for the first time in a full two-fluid setting, the reconnection-induced formation of coronal rain, in a 3D realistic setup of spine-fan reconnection. The shearing motions at the boundary will produce spine-fan reconnection. The spine collapses toward the fan and a strong current sheet forms in the vicinity of the null point, containing parts of both the fan and the spine. Even though the magnetic field configuration does not represent a “proper radial null point,” which is usually considered in analytical studies (Priest & Pontin 2009; Pontin et al. 2007), and the spine is inclined with respect to the plane tangential to the fan, the evolution is similar. The spine becomes progressively more inclined, similar to a collapse of the radial null point.

The thermal instability linear growth rate is not very sensitive to the scale size (Field 1965; Popescu Braileanu & Keppens 2024), the shape, size, and amplitude of the perturbations leading to the rain blob being determined by the reconnection process. The perturbations created by this reconnection scenario make the condensations appear earlier than in other studies of coronal rain (Li et al. 2022), and we witness blobs formed away from the perturbed null region after several minutes in the process. The blob formation time is inversely proportional to the amplitude of the perturbation, which in turn depends on many quantities, such as the magnetic field, amplitude of the driving, resistivity, and plasma beta. The growth of the thermal instability depends mainly on the cooling curve, which is a function of the temperature. The blob has a shape, size, and speed consistent with the observations. Another blob, formed previously, has a similar size and shape. We did not have gravity in our simulation; therefore, the speed of the blobs in our simulation depends mostly on the Lorentz force. The previous blob that traveled along the fan surface might have had a higher velocity because of the higher current density. The resemblance to observations is due to the fact that coronal rain is usually observed as a result of 3D null point reconnection (Mason et al. 2019), and we used values for density and temperature, magnetic field configuration, and driver velocity similar to the ones from observations of the solar corona.

The current density evolution at the null point is sensitive to the boundary driving and higher currents generally correlate to a higher decoupling in velocity. A steep increase in time in the maximum density of neutrals and charges occurs when the structure condenses to smaller scales, and the decoupling in velocity reaches its maximum. The perturbations around the secondary null point, which is not directly driven, are small and do not form condensations during the simulation time. The temperature inside the condensing blobs drops and triggers the recombination. The density of the neutrals throughout the blob remains smaller than that of the charges, and the overall density contrast is that of the charges, which is about two orders of magnitude, consistent with the observations. However, the contrast in the neutral density is almost five orders of magnitude larger than that in charges. Even if the overall density contrast is not modified by the presence of neutrals, the decoupling in velocity produces frictional heating, which can be significant, being proportional to the square of the difference in velocity. The decoupling between neutrals and charges can be observed at two different scales. Recombination slows the neutrals down and accelerates the charges, seen in a higher downward velocity of charges than of neutrals for the overall moving blob, producing decoupling at larger scales. However, in a very small layer inside the blob, the neutrals entering the structure mainly from above (because the blob is moving downward) move faster across the field lines compared to the charges, and the decoupling at these smaller scales has the opposite sign to that at larger scales.

The two fluid effects create small-scale structures due to very localized frictional heating, further modifying the temperature profile and consequent radiative cooling, taking into account that the cooling curves are very sensitive to the temperature. The two-fluid effects become more important inside the condensed blobs where the neutral density increases fast due to the recombination and there is a large decoupling in velocity and temperature.

The width of the small-scale structures created by the two-fluid effects due to elastic collisions is comparable to the mean free path between ions and neutrals, which is resolved in our simulation. However, the mean free path between neutrals and ions is not resolved at the center of the condensation, underestimating the collisional effect of the charges on the neutrals’ evolution. A better resolution would improve our understanding of this effect; however, a resolution of 1 km in a 3D setup is computationally challenging. We focused on the formation and evolution of only one blob, mentioning another previously formed blob. We expect two-fluid effects to be general for all blobs in this simulation, and simulations of a much longer duration can address statistical properties of blob sizes, typical characteristics, and overall dynamics.

The simulation does not include gravity, mainly because it is difficult to construct a stratified atmosphere in both a mechanical and ionization-recombination equilibrium for both neutrals and charges, which should be overcome in future work. In such more realistic, stratified settings, it will be relevant to look at rain cycles, while this work serves as a proof of concept that thermal instability can occur in a spine-fan setup, accounting for relevant plasma-neutral effects.

Acknowledgments

This work was supported by the FWO grant 1232122N and European Union (ERC, DynaMIT, 101086985). RK acknowledges the KU Leuven C1 project C16/24/010 UnderRadioSun and the FWO project G0B9923N Helioskill. The resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation – Flanders (FWO) and the Flemish Government.

References

All Figures

Fig. 1. Initial configuration of the magnetic field. The vertical component of the magnetic field is plotted at the bottom of the domain. The two null points are indicated by small spheres. The field lines are plotted using a gray color map, which shows the magnitude of the magnetic field. The orange circle below one of the null points indicates the location of the velocity perturbation defined by Eq. (3), which has a center at the projection of the perturbed null point on the bottom boundary and a radius equal to $2\sqrt{{\sigma }_v/2}=0.447$ Mm. The isocontours of βplasma=1 are shown by a semiopaque violet mesh of two regions, one encapsulating the two null points and an outer region.

In the text

	Fig. 2. Time evolution of: (a) current density at the null point; (b) logarithmic value of the maximum density of charges; (c) maximum value of the magnitude of the decoupling in velocity; (d) logarithmic value of the maximum density of neutrals.
In the text

Fig. 3. Slice at y=0 for the times t=442 s (top row), t=489.3 s (middle row), and t=532.8 s (bottom row), showing the formation and evolution of the rain blob. Left: Density of charges. The red arrows show the in-plane velocity of charges. Right: Density of neutrals. The red arrows show the in-plane velocity of neutrals. The in-plane magnetic field lines are overplotted by blue lines in all the panels. The perturbed null point is shown by a red dot, the other null point being outside the domain plotted in these panels. At the time t=489.3 s, when the blob is formed fully (middle row), the point with a high density, situated at x=−1.25 Mm and z=0.6 Mm, is marked in the charges and neutral density panels by a blue and green dot, respectively. This point is the same point marked in Figure 4 and Figure 6 below for the panels on the second row, for a better visualization. The line along the cut used in Figure 8 below is indicated by a dashed line: blue in the charge density panel and green in the neutral density panel.

In the text

Fig. 4. Slice at x=−1.25 Mm for the time t=489.3 s, showing the density of charges. The in-plane magnetic field lines are overplotted by blue lines. The point with a high density, situated at y=0 and z=0.6 Mm, is marked by a blue dot, being the same point marked in Figure 3 above, and in Figure 6 below for the panels on the second row, for a better visualization. The line along the cut used in Figure 8 below is indicated by a dashed orange line.

In the text

Fig. 5. Slice at y=0 for the time t=489.3 s of the magnitude of the current density. The quantity plotted is the base 2 logarithm of |J|, and the limits of the color bar of –13.5 and –9.5 correspond to |J|=8.6×10−5 A/m2 and |J|=1.38×10−3 A/m2, respectively. The in-plane magnetic field lines are overplotted by blue lines and the orange arrows indicate the in-plane decoupling velocity, the projection in the xz plane of the quantity vn−vc. The perturbed null point is shown by a red dot. The isocontours of the decoupling in temperature, Tn−Tc, are shown by green lines, indicating three values of 10, 30, and 50 kK (darker green for a larger value).

In the text

Fig. 6. Evolution of the density integrated along the y direction around the time of the simulation when a condensation appear in the simulation. Left: Density of charges. Right: Density of neutrals. The magnetic null points are highlighted with small red circles. The point located at (–1.25,0,0.6) Mm, along the high density structure in the plane y=0, is marked for reference with green and blue, the same points highlighted in Figures 3 and 4 above.

In the text

	Fig. 7. 3D integrated density of charges zoomed around the condensation point at time t=489.3 s, shown from a different angle than in Figure 6, as is indicated by the colored axes. The center of the condensation, located at (–1.25,0.0.6) Mm, is marked with a blue dot, and the perturbed null point is indicated by a red dot, consistent with Figures 3, 4, and 6.
In the text

Fig. 8. Quantities along a vertical line cut at x=−1.25 Mm, y=0 for z between 0 and 1.2 Mm at time t=489.3 s. (a) Densities of charges (solid blue line), neutrals (dashed red line), and total (dotted black line). (b) Temperature of charges (solid blue line) and neutrals (dashed red line). (c)Vertical velocity of charges (solid blue line) and neutrals (dashed red line). (d) Decoupling in velocity, vnz−vcz (red line on the left axis), and decoupling in temperature, Tn−Tc (black line on the right axis). The three thin vertical reference lines correspond to: maximum density for all the three densities plotted: charges, neutrals, and total (green), maximum decoupling in both velocity and temperature (gray), and minimum temperature for charges and neutrals (violet).

In the text