© ESO 2021

1. Introduction

The low solar atmosphere is always in a dynamic state even in the so-called quiet Sun region. Different scales of magnetic flux emergence happen all the time. The observational results indicate that the emerged magnetic fields can be as strong as several thousand Gauss in the photosphere (e.g., Yan et al. 2017; Leenaarts et al. 2018; Getling & Buchnev 2019; Liu et al. 2020; Yan et al. 2020). Magnetic reconnection can happen between emerging magnetic fields with opposite directions, or when emerged magnetic fields interact with background magnetic fields. The frequently observed compact transient brightenings, the bidirectional flows in these brightenings, and the corresponding opposite magnetic fields approaching each other in the photosphere are all indicators of these reconnection events. Numerous high-resolution observations at different wavelengths have led to the discovery of different kinds of small-scale reconnection events (e.g., Xue et al. 2016; Zhao et al. 2017; Tian et al. 2018a; Huang et al. 2018, 2019; Yang et al. 2019; Romano et al. 2019) from the photosphere to the transition regions.

Ellerman bombs (EBs; e.g., Ellerman 1917; Ding et al. 1998; Georgoulis et al. 2002; Watanabe et al. 2011; Nelson et al. 2015; Vissers et al. 2013; Yang et al. 2013; Vissers et al. 2015) are considered as one type of these low-solar-atmosphere reconnection events, and usually appear around the solar temperature minimum region (TMR). They are seen as compact intense brightenings in images of the extended Hα wings, but leave no obvious signatures in Hα core images. The temperature increases in EBs are about 400 to 3000 K (e.g., Fang et al. 2006; Rutten 2016; Hong et al. 2017). The typical lifetime of an EB is about several minutes, and the typical size is about 1″.

In the past few years, the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014) has identified a type of “hot bomb (Peter et al. 2014; Grubecka et al. 2016; Tian et al. 2016, 2018a,b; Rouppe van der Voort et al. 2017; Chitta et al. 2017; Guglielmino et al. 2019; Chen et al. 2019a)” called a UV burst (e.g., Young et al. 2018). Observations indicate that UV bursts are also formed in the magnetic reconnection process in the cool, lower solar atmosphere. They share some characteristics with traditional EBs (e.g., similar lifetimes and sizes). The strong Si IV emission, enhanced emission in the Mg II wings but not in the core, Mn I absorption in the wings of Mg II, deep absorption in Ni II lines, and compact brightenings in the AIA 1700 are usually observed in UV bursts. Therefore, there are suggestions that some UV bursts might also be formed around the solar TMR region (e.g., Peter et al. 2014; Tian et al. 2016). However, the Si IV emission in the UV bursts requires a temperature increase of about 2 × 10⁴ K in the dense photosphere, or about 8 × 10⁴ K in the upper chromosphere with much smaller density (Rutten 2016). Though UV bursts have a signature in the UV continua at 1600 and 1700 observed by the Solar Dynamics Observatory (SDO) Atmospheric Imaging Assembly (AIA; Lemen et al. 2014), they remain invisible in its He II and higher temperature coronal channels. The spectral profiles of Si IV, C II and Mg II h and k emission lines in UV bursts are significantly enhanced and broadened (e.g., Peter et al. 2014; Tian et al. 2016), and the red and blue wings of the line profile of Si IV emission lines are usually about 100 km s⁻¹ away from line center. The commissioning observations with CHROMIS in Ca II K uncovered a whole new level of fine structures in UV bursts, which are highly dynamic blob-like substructure evolving on a timescale of seconds (Rouppe van der Voort et al. 2017). The non-local thermodynamic equilibrium (non-LTE) inversions of the low-atmosphere reconnection based on SST and IRIS observations also revealed these blob-like substructures (Vissers et al. 2019), which have been proposed as one of the mechanisms causing the non-Gaussian broadened spectral profiles. However, we should also point out that UV bursts with narrow line widths of 15−18 km s⁻¹ have also been observed (Hou et al. 2016).

Recently, the relationship between EBs and UV bursts was debated. Observations showed that the occurrence of co-spatial and co-temporal EBs and UV bursts is between 10% and 20% (e.g., Vissers et al. 2015; Tian et al. 2016; Chen et al. 2019b), which suggests that the UV bursts and their associated EBs are caused or modulated by a common physical process. However, Fang et al. (2017) pointed out that the observed Hα emission cannot be reproduced using non-LTE semi-empirical modeling if the temperature is above 10⁴ K. Therefore, the Hα and Si IVemissions in these events are likely to come from the plasmas with different temperatures. Chen et al. (2019b) identified 161 EBs from the 1.6 m Goode Solar Telescope (GST), and 20 of them reveal signatures of UV bursts in the IRIS images. These latter authors found that most of these UV bursts have a tendency to appear at the upper parts of their associated flame-like EBs. However, in about one-third of the events, the formation heights of the UV bursts are not found to be higher than their corresponding EBs.

The early numerical simulations of magnetic reconnection in the low solar atmosphere can qualitatively reproduce and explain some of the typical characteristics of the observed EBs (e.g., Chen et al. 2001; Isobe et al. 2007; Archontis & Hood 2009; Danilovic 2017). The large anomalous or numerical resistivities are usually applied to trigger the reconnection processes. Partial ionization effects are rarely considered in these simulations, and the flux density of the reconnecting magnetic fields is approximately dozens of Gauss, which causes the maximum temperature increase to be only about several thousand Kelvin. Ni et al. (2015, 2016) studied magnetic reconnection around the solar TMR with much stronger magnetic fields, and included the ambipolar diffusion effect resulting from the decoupling of ions and neutrals in their simulations. The authors include extremely high resolutions to ensure that the magnetic diffusion in the reconnection process are close to the realistic one (Ni et al. 2015, 2016). For the first time, the plasmoid cascading process has been shown and studied in the chromosphere magnetic reconnection by numerical simulations (Ni et al. 2015), and the plasmas can be heated from several thousand Kelvin to above tens of thousands of Kelvin by small-scale shocks inside the plasmoids even around the solar TMR (∼500 km above the solar surface; Ni et al. 2015, 2016). By including the iterations between ions and neutrals, the multi-fluid simulations further prove that the plasmas indeed can be heated above tens of thousands of Kelvin when the reconnection magnetic fields around the solar TMR are strong (Ni et al. 2018a,b; Ni & Lukin 2018). Nevertheless, the non-equilibrium ionization-recombination effect makes the temperature increase in the reconnection process more difficult (Ni et al. 2018a).

The plasmoids in the reconnection current sheets are suggested to correspond to the blob like structures in the observations. Recent simulation results show that the plasmoid instability in the reconnection region can significantly broaden the Si IV spectral line profiles (Innes et al. 2015; Rouppe van der Voort et al. 2017), and the strong enhancement of line cores and increased emission in the line wings can both be produced when plasmoid instability occurs.

The flux emergence process has been modeled many times and used to simulate magnetic reconnection events (e.g., flares, jets and surges) with different length scales from solar chromosphere to corona (e.g., Heyvaerts et al. 1977; Shibata et al. 1992; Yokoyama & Shibata 1995; Galsgaard et al. 2005; Archontis et al. 2005; Takasao et al. 2013; Nóbrega-Siverio et al. 2016, 2018). When the resolutions are high enough and the Lundquist number exceeds the critical value (e.g., Leake et al. 2013; Ni et al. 2015; Murphy & Lukin 2015), the plasmoids can always be identified in the interaction regions where the magnetic fields with opposite directions meet with each other even in the chromosphere. The recent flux emergence simulations (Nóbrega-Siverio et al. 2016; Rouppe van der Voort et al. 2017; Nóbrega-Siverio et al. 2018) with a length scale of 10 Mm have shown the formations of plasmoids above the upper chromosphere (1.5 Mm above the solar surface).

Hansteen et al. (2017) showed that magnetic reconnection between the emerged magnetic fields give rise to EBs in the photosphere and UV bursts in the middle chromosphere. However, the EBs and UV bursts appear at different reconnection process and no co-spatial or co-temporal EBs and UV bursts were produced. The same three dimensional (3D) radiative MHD code was then used by Hansteen et al. (2019) to study the UV bursts that connect with EBs. The numerical results of these latter authors show that a long-lasting current sheet that extends over different scale heights through the low solar atmosphere is formed. The part of such a long current sheet above 1 Mm is very hot and the temperature reaches about 10⁵ K, while the temperature of the lower part is below 10⁴ K. The plasmoid-like fine structures are not mentioned or shown in their work. Their simulation results indicate that EBs and UV bursts are occasionally formed at opposite ends of a long current sheet (Hansteen et al. 2019), which could explain the UV bursts that connect with the EBs.

In this work, we study the magnetic reconnection resulting from flux emergence in the stratificated low solar atmosphere. We propose a model to explain the UV bursts with multi-thermal characteristics and the UV bursts that appear to be associated with EBs. Simulations were performed using different resolutions, and the effects of radiative cooling, ambipolar diffusion, and heat conduction are discussed. The remaining part of this paper is structured as follows. Section 2 describes the numerical models, with their initial and boundary conditions. Section 3 presents our numerical results. A summary and discussions are presented in Sect. 4.

2. Numerical setup

2.1. Model equations

We performed 2.5D MHD simulations in Cartesian geometry using the single-fluid MHD code NIRVANA(version 3.6; Ziegler 2011). The solved MHD equations are as follows:

$$\begin{array}{cc}\hfill \frac{\partial \rho }{\partial t}& =-\mathrm{\nabla }·\left(\rho \mathit{v}\right)\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\partial \left(\rho \mathit{v}\right)}{\partial t}& =-\mathrm{\nabla }·[\rho \mathit{vv}+(p+\frac{1}{2{\mu }_0}{\left|\mathit{B}\right|}^2)I-\frac{1}{{\mu }_0}\mathit{BB}]\hfill \\ \hfill & +\mathrm{\nabla }·\tau +\rho \mathit{g}\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial e}{\partial t}& =-\mathrm{\nabla }·\left[(e+p+\frac{1}{2{\mu }_0}{\left|\mathit{B}\right|}^2)\mathit{v}\right]\hfill \\ \hfill & +\mathrm{\nabla }·\left[\frac{1}{{\mu }_0}(\mathit{v}·\mathit{B})\mathit{B}\right]\hfill \\ \hfill & +\mathrm{\nabla }·[\mathit{v}\tau +\frac{\eta }{{\mu }_0}\mathit{B}\times (\mathrm{\nabla }\times \mathit{B})]\hfill \\ \hfill & -\mathrm{\nabla }·[\frac{1}{{\mu }_0}\mathit{B}\times {\mathit{E}}_{\mathrm{AD}}+{\mathit{F}}_c]\hfill \\ \hfill & +\rho \mathit{g}·\mathit{v}+L_{\mathrm{rad}}+H\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\partial \mathit{B}}{\partial t}& =\mathrm{\nabla }\times (\mathit{v}\times \mathit{B}-\eta \mathrm{\nabla }\times \mathit{B}+{\mathit{E}}_{\mathrm{AD}})\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill e& =\frac{p}{\gamma -1}+\frac{1}{2}\rho {\left|\mathit{v}\right|}^2+\frac{1}{2{\mu }_0}{\left|\mathit{B}\right|}^2\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill p& =\frac{(1+Y_i)\rho }{m_i}{k}_{\mathrm{B}}T\hfill \end{array}$$(6)

where ρ, v, B, p, T, e, Y_i are mass density, fluid velocity, magnetic field, thermal pressure, temperature, total energy density, and the ionization fraction of the plasma, respectively. The gravitational acceleration of the sun is g = −273.9 m s⁻²e_y, and m_i is the mass of a proton. We set the ratio of specific heats as γ = 5/3. Dynamic viscosity was included, and $\tau =\nu [\mathrm{\nabla }\mathit{v}+{\left(\mathrm{\nabla }\mathit{v}\right)}^{\mathrm{T}}-\frac{2}{3}(\mathrm{\nabla }·\mathit{v})I]$ is the stress tensor, where ν is the dynamic viscosity coefficient and its unit is kg m⁻¹ s⁻¹.

We applied physical magnetic diffusion in our simulations. Similar to Khomenko & Collados (2012), collisions between electrons and ions and those between neutrals and electrons both contribute to the magnetic diffusivity η, which is given by:

$$\begin{array}{c}\hfill \eta =\frac{m_{\mathrm{e}}({\nu }_{\mathrm{ei}}+{\nu }_{\mathrm{en}})}{e^2{n}_{\mathrm{e}}{\mu }_0}\end{array}$$(7)

where ν_ei and ν_en are collisional frequencies between electrons and ions and between electrons and neutrals, respectively. According to Spitzer (1962) and Braginskii (1965), Eq. (7) can be simplified as:

$$\begin{array}{c}\hfill \eta ={10}^9{T}^{-\frac{3}{2}}+1.762\times {10}^{-3}(\frac{1}{Y_i}-1)\sqrt{T}\end{array}$$(8)

We also used the classic form of anisotropic heat conduction flux:

$$\begin{array}{c}\hfill {\mathit{F}}_c=-{\kappa }_{\Vert }(\mathrm{\nabla }T·\widehat{\mathit{B}})\widehat{\mathit{B}}-{\kappa }_{\perp }[\mathrm{\nabla }T-(\mathrm{\nabla }T·\widehat{\mathit{B}})\widehat{\mathit{B}}]\end{array}$$(9)

where $\widehat{\mathit{B}}=\mathit{B}/\left|\mathit{B}\right|$ is the unit vector in the direction of magnetic field. According to Orrall & Zirker (1961) and Ni et al. (2016), the parallel and perpendicular conductivity coefficients, κ_∥ and κ_⊥, are given by:

$$\begin{array}{cc}\hfill {\kappa }_{\Vert }& =1.96\times {10}^{-11}{T}^{\frac{32}{13}}+{10}^8{T}^{-2.5}\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\kappa }_{\perp }& ={10}^8{T}^{-2.5}\hfill \end{array}$$(11)

in which 10⁸T^−2.5 and $1.96\times {10}^{-11}{T}^{\frac{32}{13}}$ are the terms associated with neutral and ionized particles, respectively.

The radiative cooling function L_rad was taken from Gan & Fang (1990):

$$\begin{array}{c}\hfill L_{\mathrm{rad}}=-1.547\times {10}^{-42}{Y}_i{\left(\frac{\rho }{m_i}\right)}^2\alpha T^{1.5}\end{array}$$(12)

and the heating function H is the same as that in Ni et al. (2016):

$$\begin{array}{c}\hfill H=1.547\times {10}^{-42}{Y}_i\frac{\rho {\rho }_0}{m_i^2}\alpha T_0^{1.5}\end{array}$$(13)

where ρ₀ and T₀ are the initial mass density and temperature.

The ambipolar diffusion field E_AD in the energy Eq. (3) and induction Eq. (4) was given by:

$$\begin{array}{c}\hfill {\mathit{E}}_{\mathrm{AD}}=\frac{1}{{\mu }_0}{\eta }_{\mathrm{AD}}[(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}]\times \mathit{B}\end{array}$$(14)

where η_AD is the ambipolar diffusion coefficient. We took the same η_AD given in Ni et al. (2015):

$$\begin{array}{c}\hfill {\eta }_{\mathrm{AD}}=1.65\times {10}^{-11}(\frac{1}{Y_i}-1)\frac{1}{{\rho }^2\sqrt{T}}\end{array}$$(15)

2.2. Initial conditions

The simulation domain extends from −8 to 8 Mm in the horizontal (x) direction and from 0 to 2 Mm in the vertical (y) direction, and the photosphere is at y = 0. We applied adaptive mesh refinement (AMR) in our simulations, and the simulations started from a base-level grid of 768 × 384. We used different highest refinement levels of 1, 3, 5, and 7 in different runs referred to Cases 1, 2, 3, and 4 to investigate the effects of increasing spatial resolution. The minimum grid size is about 45 m when the highest refinement level reaches 7. The output raw data with nonuniform grids based on the AMR skill are converted to uniform IDL data, which are then used to plot the figures in this work. We can choose to produce the IDL data with different grid sizes.

Based on the C7 atmosphere model (Avrett & Loeser 2008), the function of the initial temperature profile is as follows:

$$\begin{array}{cc}\hfill T& =-1228{\left(\frac{y}{L_0}\right)}^5+9749{\left(\frac{y}{L_0}\right)}^4-2.843\times {10}^4{\left(\frac{y}{L_0}\right)}^3\hfill \\ \hfill & +3.621\times {10}^4{\left(\frac{y}{L_0}\right)}^2-1.717\times {10}^4\left(\frac{y}{L_0}\right)+6912\hfill \end{array}$$(16)

where L₀ = 10⁶ m, y is in meters (m), and T is in Kelvin (K). We set the ionization fraction (Y_i) as:

$$\begin{array}{c}\hfill Y_i=3\times {10}^{-5}\frac{exp(3.486\times y_1^{1.2}+2.4\times y_1^{0.9}+y_1^{0.3})}{17\times y_1^{2.2}+2.639\times y_1^{1.8}+0.006}\end{array}$$(17)

where y₁ = y/L₀ + 0.12. It is worth mentioning that the Y_i does not change with time in our simulations. We chose the mass density at the photosphere in the C7 model as the mass density at our bottom boundary, and then we numerically solved the hydrostatic equilibrium equation ∇p = −ρg to obtain the initial mass density in our simulations. The initial temperature, ionization fraction, and mass density of the atmosphere are presented in Fig. 1. We mimicked a temperature minimum region at the heights around 400 km, that is, the temperature drops from ∼6900 K to ∼4200 K and then rises to ∼6600 K. The ionization fraction varies from ∼10⁻⁴ to ∼0.6, and the resulting mass density drops from ∼10⁻⁴ kg m⁻³ at the bottom boundary to ∼10⁻⁹ kg m⁻³ at the top boundary.

Fig. 1. Initial state of the simulations. Panel a: temperature (black), density (blue), and ionization fraction (red) of the atmosphere. Panel b: initial Lundquist number varying with height. The corresponding dashed lines represent the exact values in C7 model.

We also assumed a uniform initial magnetic field:

$$\begin{array}{cc}\hfill B_{x0}& =-0.85b_0\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill B_{y0}& =-0.35b_0\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill B_{z0}& =c_0{b}_0\hfill \end{array}$$(20)

where c₀ and b₀ are constants. We chose c₀ = 0.5 and b₀ = 600 G for Cases 5, 5a, 5b, and c₀ = 1 and b₀ = 500 G for other cases.

Based on the initial temperature and ionization fraction of our simulations, we find that the initial η is in the order of 10⁴ m² s⁻¹ in the whole simulation domain. We also estimated the corresponding Lundquist number Lv_A/η, where L = 1 Mm is the typical length scale of the system and ${\mathit{v}}_A=B/\sqrt{{\mu }_0\rho }$ is the Alfvén speed. The initial Lundquist number is also shown in Fig. 1, and has values from ∼10⁶ to ∼10⁹. These values significantly exceed the critical Lundquist number (10³ ∼ 10⁴) for the onset of plasmoid instability (e.g., Bhattacharjee et al. 2009; Huang & Bhattacharjee 2010; Ni et al. 2013; Leake et al. 2012, 2013).

2.3. Boundary conditions

We used outflow boundaries at the top boundary. The fluid is only allowed to flow out of the computation domain by assuming:

$$\begin{array}{c}\hfill v_{\mathit{yug}}=\{\begin{array}{cc}{v}_{\mathit{yul}}\hfill & \text{if}{v}_{\mathit{yul}}>0\hfill \\ -v_{\mathit{yul}}\hfill & \text{if}{v}_{\mathit{yul}}<0,\hfill \end{array}\end{array}$$

where v_yug and v_yul separately represent the velocities at the two top ghost layers and at the last two layers inside the simulation box in the y-direction. The gradients of velocities in the x-direction and z-direction vanish at the top boundary by assuming $\frac{\partial {\mathit{v}}_x}{\partial \mathit{y}}=0$, $\frac{\partial {\mathit{v}}_z}{\partial \mathit{y}}=0$. Similarly, the gradients of the thermal energy, mass density, and parallel components of magnetic field are set to zero. The perpendicular component of the magnetic field is obtained by divergence-free extrapolation of the magnetic field. We set the dynamic viscosity coefficient $\nu ={10}^4\rho [1+10(1+\tanh \frac{\mathit{y}-1.95L_0}{0.1L_0})]$ to make an enhancement of viscous diffusion around the top boundary.

To drive the reconnection, we mimicked a flux emergence process similar to Ni et al. (2017) by setting the magnetic fields at the bottom boundary as:

$$\begin{array}{cc}\hfill B_{\mathit{xb}}& =\{\begin{array}{cc}-0.85b_0+b_1\frac{(y-y_0)L_0^{1.6}t}{{\left[x^2+({(y-y_0)}^2\right]}^{1.3}{t}_0},& t⩽t_0\\ -0.85b_0+b_1\frac{(y-y_0)L_0^{1.6}}{{\left[x^2+({(y-y_0)}^2\right]}^{1.3}},& t>t_0\end{array}\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill B_{\mathit{yb}}& =\{\begin{array}{cc}-0.35b_0-b_1\frac{x{L}_0^{1.6}t}{{\left[x^2+({(y-y_0)}^2\right]}^{1.3}{t}_0},& t⩽t_0\\ -0.35b_0-b_1\frac{x{L}_0^{1.6}}{{\left[x^2+({(y-y_0)}^2\right]}^{1.3}},& t>t_0\end{array}\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill B_{\mathit{zb}}& =c_0{b}_0\hfill \end{array}$$(23)

where t₀ = 60 s. The emerging flux and the greatest heights it can reach increase monotonically with y₀ and b₁. The divergence-free condition of the magnetic field in the simulation domain is satisfied. We chose y₀ = −3 × 10⁵ m and b₁ = 500 G for all the cases except Cases 5, 5a, and 5b. We set y₀ = −2 × 10⁵ m and b₁ = 600 G in Cases 5, 5a, and 5b so that the emerging flux reaches greater heights. The emerging magnetic field with such a form leads to reconnection magnetic fields in the (x, y) plane around the temperature minimum region of about several hundred Gauss. The small-scale flux emergences with a maximum strength of several hundred Gauss are frequently observed in reconnection events of UV bursts and EBs (e.g., Wang et al. 2020; Chen et al. 2019a), and much stronger magnetic fields of over several thousand Gauss have also been observed in active regions (e.g., Getling & Buchnev 2019; Liu et al. 2020; Leenaarts et al. 2018; Yan et al. 2017, 2020). Therefore, the strength of reconnection magnetic fields in our simulations is reasonable.

The fluid velocity was set to zero, and the density and thermal energy were fixed with the initial values at the bottom boundary. In addition, periodic boundary conditions were used in the horizontal direction.

3. Results and discussions

We tested ten different simulation cases, as summarized in Table 1. As mentioned before, we used a maximum AMR level of 1, 3, 5, and 7 in Cases 1, 2, 3, and 4, respectively. Radiative cooling, heat conduction, and ambipolar diffusion are all included in these cases. The effects of increasing spatial resolution are studied based on these cases.

Then we turned off the radiative cooling, heat conduction, and ambipolar diffusion in Cases 3a, 3b, and 3c, respectively, to investigate the roles of different terms in our simulations; the highest refinement level in these cases is 5. We also ran cases (Case 5, 5a, 5b) with larger emerging flux. In Cases 5 and 5a, the maximum AMR refinement level is 5 and 3 respectively, but excluding ambipolar diffusion and heat conduction. The maximum AMR refinement level in Case 5b is also 3, but ambipolar diffusion is included.

3.1. Simulations with different spatial resolution

We started our simulations at relatively low resolution in order to obtain a general view of the magnetic reconnection processes. The maximum AMR level in Case 1 is 1, giving a minimum grid size of ∼2.86 km. The resolution in Case 1 is the lowest among all the Cases presented in this work, but it is still higher that those in previous 3D simulations in the low solar atmosphere (e.g., Hansteen et al. 2019). Four snapshots of Case 1 are presented in Fig. 2. The original data calculated using the NIRVANA code are transformed to uniform IDL data. First, the new emerging flux rises from the bottom boundary to the computation domain. The emerging flux interacts with the existing overlaying magnetic field, and a current sheet forms along the separatrice. The current sheet is filled with dense low-photosphere plasmas with a density that is larger than the ambient density. The temperature inside the current sheet is not heated above 10 000 K. As the current sheet continues to rise to a higher location, and plasmoids start to appear after it reaches 0.5 Mm above the solar surface. At the same time, the mass density in the current sheet decreases because the plasma is ejected to the exhaust region and then drops back because of gravity. The density and temperature distributions along the current sheet become nonuniform after the plasmoids appear. The cores of these plasmoids fill with dense photosphere plasmas and the plasmas at the edges of the plasmoids are usually more tenuous (see Fig. 3), especially when they are growing to bigger ones. Therefore, these regions with lighter plasmas are heated to higher temperatures (> 20 000 K). Later on, the plasmoids are ejected to the exhaust regions and the current sheet continues to rise. In the end, an oblique current sheet extends from 0 to 1.2 Mm in height and the mass density at the higher part of the current sheet is much lower than that at the lower part of the current sheet. Plasma at greater heights of the current sheet is heated to ∼10⁵ K, while the temperature at lower heights of the current sheet is below 10⁴ K. Such a final stage in Case 1 is similar to the results in Hansteen et al. (2019). In their 3D RMHD simulations, a vertical current sheet extends more than 2 Mm from the photosphere to the low corona, and the cool (< 10⁴ K) and hot (≳10⁵ K) plasmas are located at opposite ends of the current sheet.

Fig. 2. General view of Case 1 with the lowest resolution. The distributions of temperature (a), density (b), total magnetic field strength (c), current density (d), and vertical velocity (e) at t = 226, 375, 569, and 611 s are presented. The white lines in the third column outline the magnetic field lines. The initial background magnetic field is about 500 G.

Based on the results of Case 1, we increased the maximum AMR level step by step to investigate the effects of different computational spatial resolutions. In this work, we selected four cases with maximum AMR levels of 1, 3, 5, and 7, the comparisons among these cases are presented in Figs. 3 and 4. Increasing resolutions in the NIRVANA code shortens the time-step, which causes the runs with higher resolutions to be terminated earlier. However, the stages with plasmoid instability appearing in the reconnection region are well solved in all the cases. Fig. 3 presents the distributions of plasma temperature, density, and the joint probability distribution function (PDF) of temperature and density at a time after plasmoids appear in Cases 1 to 4. Figure 4 shows a zoom in to small scale in order to show the fine structures in one plasmoid. The Level-1 IDL data are used to plot Fig. 3, but the original levels of IDL data are used in Fig. 4.

Fig. 3. Comparisons among Cases 1, 2, 3, and 4 with different maximum AMR levels, of namely 1, 3, 5, and 7, respectively. The distributions of temperature (panels a-d), density (panels e–h), and the joint PDF of temperature and density (panel i–l) at around 376 s are presented. The vertical and horizontal purple dashed lines in panels i–l indicate the positions with the density of 10−4.5 kg m−3 and the temperature of 2 × 104 K.

Fig. 4. Same as Fig. 3 but for a smaller field of view in one plasmoid taken around 352 s.

We find that the locations and morphologies of the current sheets in all four cases are similar. However, the fine structures inside the current sheets are very different. Plasmoid instability appears in Case 1, but only the big plasmoids of first order are identified. These big plasmoids are very smooth, and no turbulent structures appear inside them (see Figs. 3 and 4). More plasmoids appear in Cases 3 and 4, and the turbulent fine structures inside the big plasmoids as shown in Fig. 4 are only clearly observed in these two cases. The numerical results in Case 3 with level 5 AMR are similar to those in Case 4 with level 7 AMR. The density and temperature distributions in the reconnection region in Case 1 are also not uniform, but the differences among different components in Cases 3 and 4 are much larger and more obvious than those in Case 1, and more plasmas are heated to higher temperatures. Before the current sheet emerges above the middle chromosphere, the maximum temperature in the reconnection region in Case 1 is always below 30 000 K. Comparing Figs. 3 and 4, one can find the maximum temperature in the main current sheet reaches above 60 000 K in Case 2, and is about 50 000 K in Case 3 and 4. After the plasmoids are ejected out along the current sheet, they merge with the background plasmas and magnetic fields, both sides of the dome region become very turbulent and many fragment currents appear (see Fig. 5), especially on the right-hand side (see inside the black dotted box in Fig. 5). Figure 5 shows the results for Case 3: a turbulent reconnection outflow region is more clearly identified in the high-resolution runs, and the maximum temperature reaches about 80 000 K in this turbulent region in Case 3.

Fig. 5. Distributions of the current density in z direction (panel a), temperature (panel b) and the synthetic Si IV 1394 Å line emissions at each grid point (panel c) at a time during the later stage of the reconnection process in Case 3. The domain in panel b and panel c corresponds to the region in the dotted box in panel a. The black arrows in panel b represent the total velocity.

In the simulations, both physical diffusion and numerical diffusion contribute to heating. Taking ohmic heating as an example, the heating rate Q can be written as Q = (η + η_num)J² = η_effJ², where η_num is numerical diffusivity and η_eff is the total effective diffusivity. When the grid size becomes smaller, the current density increases and η_num decreases. In Cases 1 and 2, the grids are coarse and the numerical diffusion dominates the heating. When the grid sizes become much smaller in Cases 3 and 4, the numerical diffusion becomes much less and the physical diffusion starts to play a role. Though the current density in Cases 3 and 4 is stronger than that in Case 2, the maximum temperature enhancement is less due to the smaller global η_eff. As shown in Figs. 3 and 4, the results in Case 3 are close to those in Case 4, which indicates that the simulation results in Case 3 and Case 4 are approaching the ones with realistic magnetic diffusion.

3.2. Effects of radiative cooling, heat conduction, and ambipolar diffusion

In order to understand the effects of radiative cooling in our simulations, we turned off radiative cooling and ran the simulation again to obtain Case 3a. Similarly, we dropped heat conduction and ambipolar diffusion in Cases 3b and 3c, respectively. The maximum AMR level in these Cases is 5, and we compared the results from Cases 3a, 3b, and 3c with those from Case 3 to study the three different effects on the magnetic reconnection process. The current sheet in these cases is not emerged to a location higher than 0.7 Mm before they are terminated because of the extremely small time-step. In order to check whether the ambipolar diffusion can effect the magnetic reconnection process at the greater height, we also tried to include this effect in Case 5. However, including the ambipolar diffusion effect in Case 5 leads to termination of the run before plasmoid instability appears. Therefore, we turned off the ambipolar diffusion effect in Case 5. We then performed another two simulations (Case 5a and 5b) using less refinements. The highest AMR level in Cases 5a and 5b is 3. The ambipolar diffusion is turned off in Case 5a and included in Case 5b. These results are presented in Figs. 6–8.

Fig. 6. Comparisons among Cases 3, 3a, 3b, and 3c where different terms are turned off in each. The distributions of temperature (panels a–d) and density (panels e–h) of cases 3, 3a (turning off radiative cooling), 3b (turning off heat conduction), and 3c (turning off ambipolar diffusion) are presented.

Fig. 7. Comparisons among three different cases. The distributions of temperature (panels a–c), total velocity the in x − y plane ($\sqrt{{\mathit{v}}_x^2+{\mathit{v}}_{\mathit{y}}^2}$) (panels d–f), and the joint PDF of temperature and density (panels g–i) are presented. Panels a, d and g are for Case 5, panels b, e and h are for Case 5a, and panels c, f and i are for Case 5b. The maximum AMR refinement levels in Case 5, 5a, and 5b are 5, 3, and 3, respectively. The ambipolar diffusion effect is not included in Case 5 and 5a, but it is included in Case 5b.

Fig. 8. Comparisons of two different energy terms at t = 403 s in Case 5b. Panel a: energy term emD contributed by effective magnetic diffusion, and panel b: energy term eAD contributed by ambipolar diffusion.

Comparing the results of Cases 3 and 3a, one sees that radiative cooling plays an important role in such a magnetic reconnection process. Including radiative cooling makes plasmoids appear earlier (see Fig 6) and the dense cores of the plasmoids colder, which may significantly affect the synthesized spectral lines and images. Previous studies suggested that radiative cooling has little effect on the reconnection processes with a strong guide field (e.g., Uzdensky & McKinney 2011; Ni et al. 2015). However, the guide field in our models is quite strong (500 G) and the radiative cooling still plays an important. The mass density in the current sheet region in our simulations is larger than previous simulations, which is probably one of the reasons for such a discrepancy. Moreover, the reconnection process is driven by the strong flux emergence process in this work, which is different from the previous current sheet studies with small initial perturbations (e.g., Ni et al. 2015). We also find that the upper boundary is unrealistically heated to high temperatures where radiative cooling is excluded. Therefore, including radiative cooling is very important for numerical studies of activities in the low solar atmosphere. In this work, we used the radiative cooling model proposed by Gan & Fang (1990), in which the authors derived their model based on the detailed non-LTE calculations. Using the plasma parameters from the low solar atmosphere model (Vernazza et al. 1981), Gan & Fang (1990) found that the radiative loss calculated from their model is close to that from non-LTE calculations. Therefore, this model can provide a reasonable approximation of the radiative cooling in the photosphere and chromosphere. However, more realistic models and further studies are still necessary. Comparing Figs. 6a and c, one can find that heat conduction does not play a role in the main current sheet region. The maximum temperature ranges from several thousand to tens of thousands of Kelvin, which might mean that heat conduction along the magnetic field line is significant. However, the plasmoids with closed magnetic fields trap the heating within the plasmoids, and the heat conduction loses efficacy.

Previous studies showed that ambipolar diffusion can be very important in heating the chromosphere, and amplifying and transporting the magnetic tension (e.g., Khomenko & Collados 2012; Martínez-Sykora et al. 2017) in the low solar atmosphere. Our simulation results indicate that ambipolar diffusion does not significantly change the reconnection process and local plasma heating in the low chromosphere. Comparing Figs. 6a and d, one finds small differences in the details between Cases 3 and 3c. Fig. 7 also shows that the ambipolar diffusion effect causes the maximum temperature in Case 5b to be slightly higher than that in Case 5a. We calculated the two terms e_mD and e_AD on the right-hand side of the energy equation at t = 403s in Case 5b with ambipolar diffusion, where $e_{\mathrm{mD}}=\mathrm{\nabla }·[\frac{{\eta }_{\mathrm{eff}}}{{\mu }_0}\mathit{B}\times (\mathrm{\nabla }\times \mathit{B})]$ and $e_{\mathrm{AD}}=-\mathrm{\nabla }·[\frac{1}{{\mu }_0}\mathit{B}\times {\mathit{E}}_{\mathrm{AD}}]$. The distributions of the two terms contributed by magnetic diffusion and ambipolar diffusion at this moment are presented in Fig. 8. We find that the maximum values of e_AD and e_mD have about the same order of magnitude. As shown in Figs. 9c and d, the plasma density inside and around the current sheet region is ≳10²⁰ m⁻³. Such a high plasma density causes the maximum ambipolar diffusion to be comparable to the maximum effective diffusion in our simulations. However, the areas with strong e_mD are much larger than the areas with strong e_AD. The heating contributed by compression processes is also very strong. Therefore, the heating contributed by ambipolar diffusion is much less efficient compared to the other processes. As we zoom into the regions with strong e_AD and e_mD, we also find that the locations with strong e_AD are always staggered with the locations with strong e_mD. We should point out that the ionization fraction Y_i is fixed in our simulations. As the realistic ionization fraction Y_i should increase with temperature and the ambipolar diffusion η_AD decreases with increasing Y_i and T, heating caused by ambipolar diffusion is possibly overestimated in our simulations. The ambipolar diffusion effect should become even less efficient than it is shown to be in this work if we consider the time-dependent ionization fraction. Ambipolar diffusion can possibly play important roles in the magnetic reconnection process above the middle chromosphere, where the plasma density is much lower. In recent work by Nóbrega-Siverio et al. (2020), the authors show that the ambipolar diffusion does not significantly change the amount of emerged magnetic flux when nonequilibrium ionization and recombination, and molecule formation of hydrogen are included, but it can efficiently heat the shock structures above the middle chromosphere.

Fig. 9. Comparisons of the main current sheet region in Cases 3 (top panels) and 5 (bottom panels) at the times when the plasmoid instability is well developed. The distributions of temperature (panels a and b), density (panels c and d), vertical velocity Vy (panels d and f) are shown, along with the synthetic Si IV 1394 Å line emission at each grid point (panels g and h). The black arrows in panels e and f represent the total velocity. The background magnetic fields are 500 G and 600 G in Cases 3 and 5, respectively. The amount of emerging flux in Case 5 is greater than that in Case 3. The heat conduction and ambipolar diffusion effects are not included in Case 5. Panels g and h are shown in arbitrary units.

3.3. Synthetic Si IV line profiles

We calculated the emissions and line profiles of the Si IV 1394 Å line of Case 3 and 5 based on the method described in Peter et al. (2006) using the atomic data package CHIANTI (version 9.0;Dere et al. 2019). Level 5 IDL data were used for calculations. First, we calculated the emissivity of the Si IV 1394 Å line at each grid point. The emissivity ε depends on density and temperature and can be written as:

$$\begin{array}{c}\hfill \epsilon =G(n_{\mathrm{e}},T)n_{\mathrm{e}}^2\end{array}$$(24)

where n_e is the electron number density and G(n_e,T) is the contribution function which can be calculated from CHIANTI. Assuming that the line profile at each grid point has a thermal width of ${\mathit{w}}_{\mathrm{th}}=\sqrt{2k_{\mathrm{B}}T/m_i}$, we calculated the spectral line profiles at each grid point with the unit of Doppler velocity:

$$\begin{array}{c}\hfill I\left(v\right)=\frac{\epsilon }{\sqrt{\pi }{w}_{\mathrm{th}}}exp[-\frac{{(v-v_0)}^2}{w_{\mathrm{th}}^2}]\end{array}$$(25)

where v₀ is the component of velocity along the line of sight. Finally, we integrated the spectral line profiles over the whole region shown in Fig. 9 to obtain the total spectrum shown in Fig. 10.

Fig. 10. Synthetic Si IV 1394 Å line profiles in Cases 3 and Case 5 taken along different lines of sight (left panel and right panel, respectively). The Si IV line profiles are integral over the whole region shown in the Fig. 9. Black, purple, blue, cyan, green, and red curves represent the line profiles taken along different lines of sight. The angles between the lines of sight and the vertical direction (y-direction) are 0°, 15°, 30°, 45°, 60°, and 75°, respectively. The dashed arrow represents the x-axis of the computation domain, and the solid arrows with different colors indicate the directions of lines of sight taken for calculating the Si IV line profiles shown with the same colors.

In Case 3, we note that the hydrogen density within the current sheet is as large as 10²¹ m⁻³, and the emissions are underestimated because the Si IV has near-Saha-Boltzmann opacity and its formation temperature drops to 10⁴ ∼ 2 × 10⁴ K at such a high density (Rutten 2016). However, there are still obvious Si IV emissions in the current sheet when using such a calculation method based on the optically thin assumption. The newly formed turbulent region also has strong emissions (see Fig. 5), and this region has a diameter of about 0.4 Mm. The stronger initial background magnetic field and emerged magnetic flux have been applied in the model to simulate Case 5. Comparing the results shown in Fig. 9, we find that the main current sheet with multiple plasmoids emerges to a higher location in Case 5, the plasmas are heated to much higher temperatures, and the maximum velocity is also higher in Case 5. The emission in the reconnection region in Case 5 is much stronger than that in Case 3 (see panels g and h of Fig. 9).

When the line of sight (LOS) is perpendicular to the direction of the reconnection outflow, the synthesized Si IV spectral profile is narrow and close to a single Gaussian profile (see the black solid lines in Fig. 10). As the LOS deviates from the direction perpendicular to the reconnection outflow, the doppler effect appears and the spectral line profile broadens and is no longer a single Gaussian profile. The multiple turbulent plasmoids along the LOS are the cause of the non-Gaussian line profile. However, the line width in Case 3 is only slightly broadened when the LOS passes through multiple plasmoids (see Fig. 10a). We find that the line width of the synthesized Si IV spectral line profile strongly depends on the maximum velocity along the LOS, which depends on both the local Alfvén velocity near the reconnection inflow region and the angle between the LOS and the reconnection outflow. Assuming that the strength of the reconnection magnetic field ranges from 200 to 2000 G, the corresponding plasma density is then between 4 × 10¹⁴ ∼ 1.6 × 10²¹ m⁻³ in order to obtain a local Alfvén velocity of 100 km s⁻¹. Therefore, the reconnection region for generating the UV emissions with broad Si IV line width over 100 km s⁻¹ should generally be located above the solar TMR. Our simulations in Case 5 with reconnection magnetic fields of ∼600 G show that the high-speed plasma that accounts for the broad Si IV line is located around 0.85 Mm above the solar surface (see Figs. 9b, d, f, and h). When the LOS is parallel to the direction of the maximum reconnection outflow velocity, the blue wing of the line profiles extends to about 100 km s⁻¹. If the reconnection region is extended to a higher position or the reconnection magnetic field is stronger than that in Case 5, a broader line width can be generated.

UV bursts with small line widths have also been reported Hou et al. (2016), and their formation mechanisms are still unclear. Our simulations show that the narrow Si IV line profiles (see Fig. 10a) can easily form in the reconnection region around and above the solar TMR, as long as the reconnection magnetic fields are strong enough to cause the plasmas to be heated above 20 000 K. Compared to the UV bursts with broad line widths, the UV bursts with narrow line widths can be formed in a reconnection region at a lower height or with weaker reconnection magnetic fields (Fig. 10a), or the LOS direction and the reconnection outflow direction are almost perpendicular to each other (Fig. 10b).

4. Summaries, conclusions, and discussions

In this work, magnetic reconnection between the emerged and background magnetic fields is simulated to study the UV bursts found to be co-spatial and co-temporal with EBs. The strength of the reconnection magnetic fields with opposite directions is about several hundred Gauss after the current sheet reaches 0.5 Mm above the solar surface. The stratified atmosphere is close to the C7 model (Avrett & Loeser 2008). The initial temperature varies with height, and a temperature minimum region is included. Also, the plasma density changes more than five orders of magnitude from the bottom to 2 Mm above the solar surface. The partial ionization effect including ambipolar diffusion is considered in simulations. An empirical radiative cooling model based on observations is applied, and an anisotropic heat conduction is also included in our simulations. Simulations were performed by using different ARM levels, resulting in different resolutions in the current sheet region in different cases.

Our simulations can explain aspects of the formation of UV bursts. The hot tenuous plasmas and the cold dense plasmas carried upward from the photosphere are mixed up during plasmoid instabilities inside the reconnection region. The ejected plasmoids are then further merged into the background magnetic fields and plasmas, which causes the formation of the new turbulent region with many fragment currents. The temperatures in the plasmoids and in the new turbulent region can range from several thousand to about one-hundred thousand Kelvin, which could explain the different emission and absorption lines observed in a UV burst. The multiple turbulent structures along the LOS give rise to the non-Gaussian line profile of Si IV. However, the line width of the synthesized Si IV spectral line profile strongly depends on the maximum velocity along the LOS. The roundish shape of the plasmoids and the new turbulent region also closely agrees with the shape of UV bursts in observations (e.g., Peter et al. 2014; Tian et al. 2016; Chen et al. 2019b), which is different from the narrow elongated shape in the recent 3D simulations (Hansteen et al. 2019). Our simulations can also explain the association between UV bursts and EBs. The hot plasmas (> 20 000 K) inside the current sheet appear during the later stage of the reconnection process, which explains the fact that UV bursts are usually seen after their EB counterparts (Ortiz et al. 2020). The hot turbulent structures with strong Si IV emissions mainly appear above and around the solar TMR, and the cold plasma extends to much lower locations near the bottom of the current sheet at the solar surface (see Figs. 5 and 9). These phenomena are consistent with the observations that UV bursts tend to be observed at greater heights than corresponding EBs, and that EBs have partial overlap with their corresponding UV bursts in space.

Previous 2.5D simulations also showed the formation of plasmoids and the well-sythesized Si IV spectral line in a stratified solar atmosphere with flux emergence (e.g., Nóbrega-Siverio et al. 2016; Rouppe van der Voort et al. 2017). However, we should point out that plasmoids are formed above the upper chromosphere and the density in UV bursts is ≲10¹⁸ m⁻³ in those previous simulations (Nóbrega-Siverio et al. 2016; Rouppe van der Voort et al. 2017). The surface flows at the bottom boundary can also self-consistently drive the formation of the reconnection region of UV bursts (Peter et al. 2019). Instead of using a realistic density stratification, cases with different plasma densities were tested to study the effects of plasma β on generating UV bursts (Peter et al. 2019). The reconnection magnetic field in the simulations in Peter et al. (2019) is about 50 G, and the plasma with a density > 10¹⁹ m⁻³ is not heated above 20 000 K in their simulations. Observations suggest that the density of UV bursts is ≳10¹⁹ m⁻³ (Young et al. 2018); it could be even higher for those connecting with EBs (e.g., Tian et al. 2016). When the plasma β is the same, the temperature increase is much less in the plasma with a higher density in the lower atmosphere (Ni & Lukin 2018). Therefore, it is more difficult for plasmas to be heated to high temperatures above 20 000 K in the lower atmosphere. In this work, plasmoids are formed below the middle chromosphere, the density in UV bursts is about two to three orders of magnitude higher and reconnecting magnetic fields are also stronger than those in the previous models. Furthermore, the cool cores of the plasmoids and the multi-thermal turbulent fine structures inside the plasmoids and the reconnection outflow regions were not clearly shown in those previous simulations.

In the 3D simulations of Hansteen et al. (2019), the part of the long current sheet with Si IV emission extends downwardly to a position around the middle chromosphere. However, the plasmoid-like structures are not shown in their work. The hot and cool parts are basically located at opposite ends of a long current sheet (Hansteen et al. 2019). Ortiz et al. (2020) further proposed that the coexisting EBs and UV bursts could be part of the same reconnection system, even though they are located far apart on the vertical axis, which is different from the model in the present work. Though our results also show that the cold EB part extends to a lower height, the cold and hot plasmas could be concentrated in one plasmoid or the turbulent outflow region at about the same height. Priest et al. (2018) and Syntelis et al. (2019), Syntelis & Priest (2020) recently studied magnetic reconnection driven by photospheric flux cancelation. These authors pointed out that the hot and cool outflows produced without time difference and spatial offset are also found in a reconnection process without plasmoid-like structures.

A more realistic radiative-transfer process than that applied here has been well included in some of the previous simulations (e.g., Nóbrega-Siverio et al. 2016; Rouppe van der Voort et al. 2017; Hansteen et al. 2019) used to study UV bursts, but we have considered partial ionization effects and included the more realistic magnetic diffusion.

The simulations in this work have coupled the macroscopic scale of several Mm and the microscopic scale down to the ion inertial length (∼45 m), and both the observational features and the hidden small-scale physics are revealed. The results indicate that high-resolution simulations including realistic diffusivity in the stratified low solar atmosphere are essential for better understanding the magnetic reconnection mechanisms and plasma heating during the formation process of UV bursts. Obviously, low resolutions inevitably result in large numerical diffusivity and low Lundquist number. Large numerical diffusivity might result in unrealistic plasma heating in the reconnection region (as shown in Cases 1 and 2). Plasmoid instability will not appear if the Lundquist number is too low (e.g., Bhattacharjee et al. 2009; Ni et al. 2013; Leake et al. 2012), and the mixed hot tenuous and cold dense plasmas in the plasmoids and the turbulent reconnection outflow region will be smoothed out. This work qualitatively reveals the multi-thermal turbulent structures in the reconnection region of a UV burst. However, nonequilibrium ionization-recombination and more realistic radiative transfer in future simulations are needed to find the exact temperatures of the plasmas in these events. The obvious flux cancelations at the solar surface usually appear during the formation process of EBs and UV bursts (e.g., Peter et al. 2014; Tian et al. 2016; Chen et al. 2019b), which should also be included and discussed in future more realistic 3D simulations.

Acknowledgments

We thank Professor Rony Keppens for helpful discussions. This research is supported by the Strategic Priority Research Program of CAS with grants XDA17040507 and QYZDJ-SSWSLH012; the NSFC Grants 11973083; the Youth Innovation Promotion Association CAS 2017; the Applied Basic Research of Yunnan Province in China Grant 2018FB009; the Yunnan Ten-Thousand Talents Plan-Young top talents; the project of the Group for Innovation of Yunnan Province grant 2018HC023; the YunnanTen-Thousand Talents Plan-Yunling Scholar Project; the Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (nsfc2015-460, nsfc2015-463, the second phase); the Computational Solar Physics Laboratory of Yunnan Observatories; the Max Planck Partner Group program. Y.C. is supported by the China Scholarship Council for his stay in Germany. Y.J. and H.T. are supported by NSFC grants 11825301, 11790304(11790300), and the Max Planck Partner Group program.

References

All Tables

All Figures

	Fig. 1. Initial state of the simulations. Panel a: temperature (black), density (blue), and ionization fraction (red) of the atmosphere. Panel b: initial Lundquist number varying with height. The corresponding dashed lines represent the exact values in C7 model.
In the text

	Fig. 2. General view of Case 1 with the lowest resolution. The distributions of temperature (a), density (b), total magnetic field strength (c), current density (d), and vertical velocity (e) at t = 226, 375, 569, and 611 s are presented. The white lines in the third column outline the magnetic field lines. The initial background magnetic field is about 500 G.
In the text

Fig. 3. Comparisons among Cases 1, 2, 3, and 4 with different maximum AMR levels, of namely 1, 3, 5, and 7, respectively. The distributions of temperature (panels a-d), density (panels e–h), and the joint PDF of temperature and density (panel i–l) at around 376 s are presented. The vertical and horizontal purple dashed lines in panels i–l indicate the positions with the density of 10−4.5 kg m−3 and the temperature of 2 × 104 K.

In the text

	Fig. 4. Same as Fig. 3 but for a smaller field of view in one plasmoid taken around 352 s.
In the text

	Fig. 5. Distributions of the current density in z direction (panel a), temperature (panel b) and the synthetic Si IV 1394 Å line emissions at each grid point (panel c) at a time during the later stage of the reconnection process in Case 3. The domain in panel b and panel c corresponds to the region in the dotted box in panel a. The black arrows in panel b represent the total velocity.
In the text

	Fig. 6. Comparisons among Cases 3, 3a, 3b, and 3c where different terms are turned off in each. The distributions of temperature (panels a–d) and density (panels e–h) of cases 3, 3a (turning off radiative cooling), 3b (turning off heat conduction), and 3c (turning off ambipolar diffusion) are presented.
In the text

Fig. 7. Comparisons among three different cases. The distributions of temperature (panels a–c), total velocity the in x − y plane ($\sqrt{{\mathit{v}}_x^2+{\mathit{v}}_{\mathit{y}}^2}$) (panels d–f), and the joint PDF of temperature and density (panels g–i) are presented. Panels a, d and g are for Case 5, panels b, e and h are for Case 5a, and panels c, f and i are for Case 5b. The maximum AMR refinement levels in Case 5, 5a, and 5b are 5, 3, and 3, respectively. The ambipolar diffusion effect is not included in Case 5 and 5a, but it is included in Case 5b.

In the text

	Fig. 8. Comparisons of two different energy terms at t = 403 s in Case 5b. Panel a: energy term emD contributed by effective magnetic diffusion, and panel b: energy term eAD contributed by ambipolar diffusion.
In the text

Fig. 9. Comparisons of the main current sheet region in Cases 3 (top panels) and 5 (bottom panels) at the times when the plasmoid instability is well developed. The distributions of temperature (panels a and b), density (panels c and d), vertical velocity Vy (panels d and f) are shown, along with the synthetic Si IV 1394 Å line emission at each grid point (panels g and h). The black arrows in panels e and f represent the total velocity. The background magnetic fields are 500 G and 600 G in Cases 3 and 5, respectively. The amount of emerging flux in Case 5 is greater than that in Case 3. The heat conduction and ambipolar diffusion effects are not included in Case 5. Panels g and h are shown in arbitrary units.

In the text

Fig. 10. Synthetic Si IV 1394 Å line profiles in Cases 3 and Case 5 taken along different lines of sight (left panel and right panel, respectively). The Si IV line profiles are integral over the whole region shown in the Fig. 9. Black, purple, blue, cyan, green, and red curves represent the line profiles taken along different lines of sight. The angles between the lines of sight and the vertical direction (y-direction) are 0°, 15°, 30°, 45°, 60°, and 75°, respectively. The dashed arrow represents the x-axis of the computation domain, and the solid arrows with different colors indicate the directions of lines of sight taken for calculating the Si IV line profiles shown with the same colors.

In the text