Jet formation in solar atmosphere due to magnetic reconnection

Using numerical simulations, we show that jets with features of type II spicules and cold coronal jets corresponding to temperatures $10^{4}$ K can be formed due to magnetic reconnection in a scenario in presence of magnetic resistivity. For this we model the low chromosphere-corona region using the C7 equilibrium solar atmosphere model and assuming Resistive MHD rules the dynamics of the plasma. The magnetic filed configurations we analyze correspond to two neighboring loops with opposite polarity. The separation of the loops' feet determines the thickness of a current sheet that triggers a magnetic reconnection process, and the further formation of a high speed and sharp structure. We analyze the cases where the magnetic filed strength of the two loops is equal and different. In the first case, with a symmetric configuration the spicules raise vertically whereas in an asymmetric configuration the structure shows an inclination. With a number of simulations carried out under a 2.5D approach, we explore various properties of excited jets, namely, the morphology, inclination and velocity. The parameter space involves magnetic field strength between 20 and 40 G, and the resistivity is assumed to be uniform with a constant value of the order $10^{-2}\Omega\cdot m$


INTRODUCTION
Magnetic reconnection is a topological reconfiguration of the magnetic field caused by changes in the connectivity of its field lines (Priest 1984;Priest et al. 2000). It is also a mechanism of conversion of magnetic energy into thermal and kinetic energy of plasma when two antiparallel magnetic fields encounter and reconnect with each other. Magnetic reconnection can occur in the chromosphere, photosphere and even in the convection zone. In particular, the chromosphere has a very dynamical environment where magnetic features such as H α upward flow events (Chae et al. 1998) and erupting mini-filaments (Wang et al. 2000) happen. The dynamics of the chromosphere at the limb region is dominated by spicules (Beckers et al. 1968) and related flows such as mottles and fibrils on the disk (Hansteen et al. 2006;De Pontieu et al. 2007a). Spicular structures are also visible at the limb in many spectral lines at the transition region temperatures (Mariska 1992;Wilhelm 2000), and some observations suggest that coronal dynamics are linked to spicule-like jets Tsiropoula et al. 2012;Cheung et al. 2015;Skogsrud et al. 2015;Tavabi et al. 2015a;Narang et al. 2016). With the large improvement in spatiotemporal stability and resolution given by the Hinode satellite (Kosugi et al. 2007), and with the Swedish 1 m Solar Telescope (SST) (Scharmer et al. 2008), two classes of spicules were defined in terms of their different dynamics and timescales (De Pontieu et al. 2007c).
The so-called type I spicules have lifetimes of 3 to 10 minutes, achieve speeds of 10-30 km s −1 , and reach heights of 2-9 Mm (Beckers et al. 1968;Suematsu et al. 1995), and typically involve upward motion followed by downward motion. In , the authors studied in detail the propagating shocks using simplified onedimensional models. In (Hansteen et al. 2006;De Pontieu et al. 2007a;Heggland et al. 2007), the authors study the propagation of shocks moving upwards passing through the upper chromosphere and transition region toward the corona. They also describe how the spicule-driving shocks can be gener-ated by a variety of processes, such as collapsing granules, pmodes and dissipation of magnetic energy in the photosphere and lower chromosphere. In (Matsumoto & Shibata 2010) the authors state that spicules can be driven by resonant Alfvén waves generated in the photosphere and confined in a cavity between the photosphere and the transition region. Another studies in this direction, like (Murawski & Zaqarashvili 2010;Murawski et al. 2011), where they use ideal MHD and perturb the velocity field in order to stimulate the formation of type I spicules and macro-spicules. Furthermore in (Scullion et al. 2011), the authors simulate the formation of wave-driven type I spicules phenomena in 3D trough a Transition Region Quake (TRQ) and the transmission of acoustic waves from the lower chromosphere to the corona.
Type II spicules are observed in Ca II and Hα, these spicules have lifetimes typically less than 100s in contrast with type I spicules that have lifetimes of 3 to 10 min, are more violent, with upward velocities of order 50-100 km s −1 and reach greater heights. They usually exhibit only upward motion (De Pontieu et al. 2007b), followed by a fast fading in chromospheric lines without observed downfall. Spicules of type II seen in the Ca II band of Hinode fade within timescales of the order of a few tenths of seconds (De Pontieu et al. 2007a). The type II spicules observed on the solar disk are dubbed "Rapid Blueshifted Events" (RBEs) (Langangen et al. 2008;Rouppe van der Voort et al. 2009). These show strong Doppler blue shifted lines in the region from the middle to the upper chromosphere. The RBEs are linked with asymmetries in the transition region and coronal spectral line profiles ). In addition the lifetime of RBEs suggests that they are heated with at least transition region temperatures (De Pontieu et al. 2007c;Rouppe van der Voort et al. 2009). Type II spicules also show transverse motions with amplitudes of 10-30 km s −1 and periods of 100-500 s (Tomczyk et al. 2007;McIntosh et al. 2011;Zaqarashvili & Erdélyi 2009), which are interpreted as upward or downward propagating Alfvenic waves (Okamoto & De Pontieu 2011;Tavabi et al. 2015b), or MHD kink mode waves (He et al. 2009;McLaughlin et al. 2012;Kuridze et al. 2012).
As mentioned above, there are several theoretical and observational results about type II spicules, but there is little consensus about the origin of type II spicules and the source of their transverse oscillations. Some possibilities discussed suggest that type II spicules are due to the magnetic reconnection process (Isobe et al. 2008;De Pontieu et al. 2007c;Archontis et al. 2010), oscillatory reconnection process (McLaughlin et al. 2012), strong Lorentz force (Martínez-Sykora et al. 2011) or propagation of p-modes (de Wijn et al. 2009). Moreover, type II spicules could be warps in 2D sheet like structures (Judge et al. 2011). A more recent study suggests another mechanism, for instance in (Sterling& Moore 2016), the authors suggested that solar spicules result from the eruptions of small-scale chromospheric filaments.
The limited resolution in observations and the complexity of the chromosphere make difficult the interpretation of the structures, and even question the existence of type II spicules as a particular class, for instance in (Zhang et al. 2012). In consequence, these difficulties spoil the potential importance of magnetic reconnection as a transcendent mechanism in the solar surface. Nevertheless, there is evidence that magnetic reconnection is a good explanation of chromospheric anemone jets (Singh et al. 2012), which are observed to be much smaller and much more frequent than surges (Shibata et al. 2007). A statistical study performed by  showed that the chromospheric anemone jets have typical lengths of 1.0-4.0 Mm, widths of 100-400 km, and cusp size of 700-2000 km. Their lifetimes is about 100-500 s and their velocity is about 5-20 km s −1 . Other types of coronal jets can be generated by magnetic reconnection, for example in (Yokoyama & Shibata 1995, 1996 the authors using twodimensional numerical simulations study the jet formation, or in (Nishizuka et al. 2008), it is shown that emerging magnetic flux reconnects with an open ambient magnetic field and such reconnection produces the acceleration of material and thus a jet structure. The reconnection seems to trigger the jet formation in a horizontally magnetized atmosphere, with the flux emergence as a mechanism (Archontis et al. 2005;Galsgaard et al. 2007).
Another approach uses a process that produces a magnetic reconnection using numerical dissipation of the ideal MHD equations, and the atmosphere model is limited to have a constant density and pressure profiles and assumes there is no gravity (Pariat et al. 2009(Pariat et al. , 2010(Pariat et al. , 2015Rachmeler et al. 2010). In this paper we show that magnetic reconnection can be responsible for the formation of jets with some characteristics of Type II spicules and cold coronal jets (Nishizuka et al. 2008), for that i) we solve the system of equations of the Resistive MHD subject the solar gravitational field, ii) we assume a completely ionized solar atmosphere consistent with the C7 model. The resulting magnetic reconnection accelerates the plasma upwards by itself and produces the jet.
The paper is organized as follows, in Section 2 we describe the resistive MHD equations, the model of solar atmosphere, the magnetic field configuration used the numerical simulations and the numerical methods we use. In Section 3, we present the results of the numerical simulations for various experiments. Section 4 contains final comments and conclusions.

The system of Resistive MHD equations
The minimum system of equations allowing the formation of magnetic reconnection is the resistive MHD. In this paper we follow (Jiang et al. 2012) to write the dimensionless Extended Generalized Lagrange Multiplier (EGLM) resistive MHD equations that include gravity as follows: where ρ is the mass density, v is the velocity vector field, B is the magnetic vector field, E is the total energy density, where γ = 5/3, the plasma pressure p is described by the equation of state of an ideal gas; g is the gravitational field, J is the current density, η is the magnetic resistivity tensor and ψ is a scalar potential that helps damping out the violation of the constraint. Here c h is the wave speed and c p is the damping rate of the wave of the characteristic mode associated to ψ. In this study we consider an uniform and constant magnetic resistivity, because the solar chromosphere is fully collisional and anomalous or space dependent resistivity -which is the result of various collisionless processes-may not be expected . We normalize the equations with the quantities given in Table 1, which are typical scales in the solar corona.
In the EGLM-MHD formulation, equation (7) is the magnetic field divergence free constraint. As suggested in (Dedner et al. 2002), the expressions for c h and c p are where ∆t is the time step, ∆x, ∆y and ∆z are the spatial resolutions, c c f l < 1 is the Courant factor, c d is a problem dependent coefficient between 0 and 1, this constant determines the damping rate of divergence errors. The parameters c h and c p are not independent of the grid resolution and the numerical scheme used, for that reason one should adjust their values. In our simulations we use c p = √ c r c h , with c r = 0.18 and c h = 0.01. For our analysis we use a 2.5D model, but we solve the 3D resistive MHD equations with high resolution in the xz plane and with a 4 cells along the y direction, so that the speeds c h and c p only depend on ∆x and ∆z. All the state variables depend on x, y and z (González-Avilés& Guzmán 2015; González-Avilés et al. 2015).
Temperature in red and mass density in green as a function of height for the C7 equilibrium solar atmosphere model. Notice the steep jump in Temperature, which makes the C7 model a realistic one.

Numerical methods
We solve numerically the resistive EGLM-MHD equations given by the system of equations (1)-(7) on a single uniform cell centered grid, using the method of lines with a third order total variation diminishing Runge-Kutta time integrator (Shu & Osher 1989). In order to use the method of lines, the RHS of resistive MHD equations are discretized using a finite volume approximation with High Resolution Shock Capturing methods (LeVeque 1992). For this, we first reconstruct the variables at cell interfaces using the Minmod limiter. The numerical fluxes are calculated using the Harten-Lax-van-Leer-Contact (HLLC) approximate Riemann solver (Li 2005).

Model of the solar atmosphere
We choose the numerical domain to cover part of the photosphere, chromosphere and corona. We consider the atmosphere in hydrostatic equilibrium and study the evolution on a finite xz domain, where x is a horizontal coordinate and z labels height. The temperature field is assumed to obey the semiempirical C7 model of the chromosphere (Avrett & Loeser 2008) and is distributed to obtain optimum agreement between calculated and observed continuum intensities, line intensities, and line profiles of the SUMER (Curdt et al. 1999) atlas of the extreme ultraviolet spectrum. The profiles of T (z) and ρ(z) are shown in Fig. 1, where the expected gradients at the transition region can be seen.

The magnetic field
The magnetic field in the model is chosen as a superposition of two neighboring loops. Following (Priest 1982;Del Zanna et al. 2005) we construct a loop with the vector potential where B 01 is the photospheric field magnitude at the foot points x = ±L/2 and k = π/L. Here L is the distance between the two foot points of the loop and k defines the nodes of the potential. In this model the components of the magnetic field can be represented as In order to superpose two loops we use a modified version of (8): where l 0 defines the location of the foot points for each loop, B 01 and B 02 are the magnetic field strengths of the left and right loop respectively. With the parameter l 0 it is possible to control the separation between the two magnetic loops, which in turn will influence the thickness of a current sheet. The case B 01 = B 02 describes two neighboring loop configurations with the same magnetic field strength, whereas B 01 = B 02 describes two nearby loops with different magnetic field strength. A schematic picture of these two configuration is shown in Fig.  2. 3. RESULTS OF NUMERICAL SIMULATIONS Putting altogether, we run a number of simulations with the magnetic field configuration given by equation (11) in a scenario with constant resistivity η = 5 × 10 −2 Ω · m across the whole domain, which is a realistic value estimated for a fully ionized solar atmosphere (Priest 2014). We experiment with various values of the magnetic field strength and separation of the loops as indicated in Table 2.
We fix k = π/L with L = 8 Mm and the simulations were carried out in a domain x ∈ [−4, 4], y ∈ [0, 1], z ∈ [0, 10] in units of Mm, covered with 300×4×375 grid cells. In all the numerical simulations we use out flux boundary conditions, which in our approach translate into copy boundary conditions applied to the conservative variables (Toro 2009).

Symmetric configurations
For the first case show the results for the case B 01 = B 02 for three values of magnetic strengths 20, 30, 40 G and l 0 = 2.5, 3.0, 3.5 Mm. Representative simulations for this case are shown in Fig. 3 for l 0 = 3.5 Mm. For Run #1 correspond to the typical formation of a jet, with a special feature at the top with a bulb related to a Kelvin-Helmholtz type of instability which is contained due to the presence of the magnetic field. This jet reaches a height of 9 Mm with a maximum speed v z,max ≈ 34 km/s at time t = 180 s. After this time the jet starts falling down until it disperses away by time t ≈ 400 s. In Fig. 4 we show T inside , ρ inside and v z,max as a function of time estimated inside the jet, and the method used to estimate this properties of the jet. The temperature during the evolution is a useful scalar at determining the location of the jet, because it shows a minimum precisely located where the head of the jet is and is used to estimate the maximum height achieved.
We also show in Fig. 3 the y component of current density J y , which is the most significant component and shows the formation of an elongated structure representing the reconnection process. At the middle panel of Fig. 3 we show the results of Run #3 (the magnetic field strength is 20 G). Due to smaller magnetic field strength the resulting jet reaches a height of 5.5 Mm with a maximum vertical velocity v z,max ≈ 32.5 km/s at time t = 180 s. Finally at the lower panel of Fig. 3 we show Run #6, in this case the effect of a weaker magnetic field is seen in the appearance of a small jet that reaches a height of 3 Mm with a maximum speed v z,max ≈ 21.32 km/s at time t = 180 s. These results indicate that the height of the jet is stronger for bigger values of B 01 = B 02 , and the sharpness of the jet is also more clear for magnetic fields. We present the snapshots of all the cases at t = 180 s for comparison, the maximum height and velocity are different in each case.
The separation of the loops is also important due to the dependence of the current sheet parameters on it. In the case of configurations with a larger separation, the plasma is accelerated rapidly, which produces diffusion and consequently the jet does not form. However in the case of closer loops, the plasma is accelerated slowly, which allows the formation of a jet later on. According to our results, for the parameters we analyzed, the most effective separation between the loops to trigger a jet is l 0 = 3.5 Mm.

Non-symmetric configurations
In this case we show the results for a more realistic magnetic configuration of the numerical simulations corresponding to the case of non-symmetric magnetic field loops, i.e. when the magnetic strengths of the left and right loops are different (B 01 = B 02 ) for the combinations of magnetic field strengths 20, 30 and 40 G and l 0 = 2.5, 3.0, 3.5 Mm. In order to illustrate the effect of the asymmetry in the formation of jets, we show the results for Runs # 7, # 9 and # 13 in Fig. 5. At the top panel we present the result for Run # 7, that shows the inclination of the jet toward the loop with the weak magnetic field. Similar to the previous case, the top part of the jet exhibits a bulb which appears due to Kelvin-Helmholtz instability (Kuridze et al. 2016). This jet reaches a height of 8.5 Mm and a maximum vertical velocity v z,max ≈ 24.65 km/s at t = 180 s. Like in the symmetric case, the jet starts to weaken and finally vanishes. At the middle panel of Fig. 5 we show Run # 9, in this case the jet shows a more significant inclination to the direction of the weaker magnetic field. In this case the jet reaches a height of 8.0 Mm and a maximum veloc- Table 2 Maximum velocity, temperature and density of jets in terms of the magnetic field strength of the two loops. We include symmetric (B 01 = B 02 ) and non-symmetric configurations (B 01 = B 02 ).
Run # B 01 (G) B 02 (G) l 0 (Mm) vz,max (km/s) T inside (K) ρ inside (kg/m 3 ) ity v z ≈ 44.15 km/s at time t = 180 s. The inclination is also shown in the y component of the current density J y . At the bottom we show Run # 13, in this case the jet is small due to the weaker magnetic field of the loops. The maximum height of the jet is about 5 Mm and its maximum vertical speed is v z ≈ 32.14 km/s at t = 180s. The inclined propagation of the jet produces a distortion of the magnetic field lines, that can be seen by comparing the lines at initial time with the lines at the time of the snapshot. The velocity vector field shows the plasma moving in the direction of the weaker magnetic field. The y component of the current density J y shows the formation of an elongated current sheet, directly related to the elongated shape of the jet and similar shape like those observed by Hinode for instance in Fig. 1 of (Tavabi et al. 2015b). We do not show the density of the plasma in these Figures, however its shape is pretty much that of the Temperature profile.
We summarize the results with combinations of the magnetic field configurations presented in Table 2 for symmetric and non-symmetric magnetic loops. In the same Table we also show the resulting values for the maximum velocity of the plasma along the vertical direction v z,max , plasma temperature estimate T inside and density estimate ρ inside measured along the jet.
In the case of two non-symmetric magnetic loops the inclination depends of the magnetic field strength between the two loops as shown in Fig. 5. In order to see this dependence more clearly we show the inclination angle of the jet as a function of the ratio B 02 /B 01 for different separation parameters in Fig. 6 at the time when the jet reaches the maximum height for each case.

CONCLUSIONS
In this paper we present the numerical solution of the equations of the resistive MHD submitted to the solar constant gravitational field, and simulate the formation of narrow jet structures on the interface low-chromosphere and corona. For this we use a magnetic field configuration of two superposed loops in a way that a current sheet is formed that allows the magnetic reconnection process, which in turn accelerates the plasma.
An ingredient of our simulations is that we use a realistic atmospheric model that includes the transition region. The rarefaction of the environment above the transition region helps the acceleration of the plasma. We can summarize our findings in the schematic picture shown in Fig.7. This rarefied atmosphere allows the formation of a bulb at the top of the  Here we illustrate the way we measure the properties of the jet using diagnostics of the case B 01 = B 02 = 40 G and l 0 = 3.5 Mm. On the top-left we show snapshots of the plasma temperature at t = 50, 100, 250 s along the z-axis. The plasma moving upwards during the evolution is cold with respect to the temperature of the corona and therefore the temperature shows a minimum where the front of the jet is developing. We say the Temperature of the jet is that of the minimum and changes with time and position. We define the temperature of the jet T inside as that of the minimum at every time, independently of the location of such minimum in space and the result appears in the top-right panel. With a similar procedure we measure vz of the jet as the maximum along the z-axis (remember that in the case of T it was a minimum) and plot this quantity in time vz,max that we show in the bottom-left panel. Finally we also estimate the jet density ρ inside with a similar procedure and show the result in the bottom-right. jet due to Kelvin-Helmholtz instability (Kuridze et al. 2016), which is contained and stabilized by the magnetic field as found in (Flint et al. 2014;Zaqarashvili et al. 2014).
We consider symmetric and asymmetric magnetic field configurations. In the symmetric case different jet properties where found in terms of separation and magnetic field strength of the loops. The magnetic field used ranges from 20 to 40 G, leading to the conclusion that the stronger the magnetic field the higher and faster the jet. The separation of the loops' feet was also found to be important, because it determines the thickness of the current sheet that later on produces the magnetic reconnection. With our parameters, a separation with l 0 = 4Mm is so large that no jet is formed anymore. The Temperature within the jet structure is of the order of 10 4 K, which is within the observed range of a cool jet (Nishizuka et al. 2008). An illustrative example is that of B 01 = B 02 = 40G and l 0 = 3.5Mm, which shows a height of 7Mm measured from the transition region, and a maximum vertical velocity of v z ≈ 34 km s −1 , parameters similar to those of Type II spicules (De Pontieu et al. 2007c). The evolution of these structures indicate that they last about 200 s, which is a lifetime similar to type II spicules.
In the case of asymmetric magnetic field configurations we also simulated the formation of jets with similar properties of Temperature, velocity and height. These jets show a considerable inclination toward the loop with the weaker magnetic field. We found that the inclination of the jet depends on the magnetic field ratio of the two loops. According to the results of this paper, a good model for the formation of realistic jets mimicking type II spicules is to have two magnetic loops close together with opposite polarity. This produces a current sheet at chromospheric level capable to trigger magnetic reconnection. A key ingredient in the process is the inclusion of magnetic resistivity, which is a mechanism consistent with the magnetic reconnection process. Another feature of these jets is that they are based at the level of the transition region, which is characterized by a sharp gradient in density and temperature.   Figure 7. We show the two configurations of the magnetic coronal loops with foot points at the photospheric level. Initially two magnetic loops close together have opposite polarity which produces magnetic reconnection. In the case of two symmetric magnetic loops (left) the jet appears at the middle of the configuration and moves straight upward until it reaches the solar corona and diffuses later on. In the case of two asymmetric loops (right) the jets appears at the middle of the configuration, but in this case the loop with the stronger magnetic field pushes the jet toward the weak magnetic field and the jet reaches the solar corona. In both cases there is an elongated current sheet represented with the density current J in the perpendicular direction to the jet.