I. Jet Formation and Evolution due to 3D Magnetic Reconnection

Using simulated data-driven three-dimensional resistive MHD simulations of the solar atmosphere, we show that magnetic reconnection can be responsible of the formation of jets with characteristic of Type II spicules. For this, we numerically model the photosphere-corona region using the C7 equilibrium atmosphere model. The initial magnetic configuration is a 3D potential magnetic field, extrapolated up to the solar corona region from a dynamic realistic simulation of solar photospheric magnetoconvection model which is mimicking quiet-Sun. In this case we consider a uniform and constant value of the magnetic resistivity of $12.56 ~\Omega~{\rm m}$. We have found that formation of the jets depends on the Lorentz force, which helps to accelerate the plasma upwards. Analyzing various properties of the jet dynamics, we found that the jet structure shows Doppler shift near to regions with high vorticity. The morphology, upward velocity, covering a range up to 100 $\rm km$ $\rm s^{-1}$, and life-time of the estructure, bigger than 100 ${\rm s}$, are similar to those expected for Type II spicules.


INTRODUCTION
Jet-like emissions of plasma in the solar atmosphere have been extensively observed over a range of wavelengths, e.g. X-ray, EUV and Hα, that usually occur in active regions and polar coronal holes. It is believed that many plasma jets are produced directly by magnetic reconnection, when oppositely directed magnetic fields come in contact (see e.g. Shibata et al. 2007). The magnetic reconnection acts as a mechanism of conversion of the magnetic field energy into thermal and kinetic energy of the ejected plasma and can occur from the convection zone to the solar corona. In particular, the observed chromospheric dynamics at the solar limb is dominated by spicules (Beckers et al. 1968), which are ubiquitous, highly dynamic jets of plasma (Secchi 1878;Tsiropoula et al. 2012;De Pontieu et al. 2007b). The improvement in the resolution of the observations by the Hinode satellite and Swedish 1 m Solar Telescope (SST) on La Palma (Kosugi et al. 2007;Scharmer et al. 2003Scharmer et al. , 2008 has suggested the existence of two classes of spicules. The first type of spicules are so-called Type I, which reach maximum heights of 4-8 Mm, maximum ascending velocities of 15-40 km s −1 , have a lifetime of 3-6.5 minutes (Pereira et al. 2012), and show up and downward motions (Beckers et al. 1968;Suematsu et al. 1995). These Type I spicules are probably the counterpart of the dynamic fibrils on the disk. They follow a parabolic (ballistic) path in space and time. In general the dynamics of these spicules is produced by mangneto-acoustic shock wave passing or wave-driving through the chromosphere (Shibata et al. 1982;De Pontieu et al. 2004;Hansteen et al. 2006;Martínez-Sykora et al. 2009;Matsumoto & Shibata 2010;Scullion et al. 2011). The second type of spicules (Type II) reach maximum heights of 3-9 Mm (longer in coronal holes) and have shorter lifetimes of 50-150 s than Type I spicules (De Pontieu et al. 2007a;Pereira et al. 2012). These Type II spicules show apparent upward motions with speeds of order 30-110 km s −1 .
At the end of their life they usually exhibit rapid fading in chromospheric lines (De Pontieu et al. 2007b). However the timescale of both types of spicules depends on the temperature, i.e., Ca II observations show short spicules, whereas Mg II or transition region lines show lifetimes of the order of ten minutes (Pereira et al. 2014;Skogsrud et al. 2015). Also in Zhang et al. (2012), the authors stablished the complexity of differentiating between Type I and Type II, so in general we can say that the Spicules are not well understood. It has been suggested from observations that Type II spicules are continuously accelerated while being heated to at least transition region temperatures . Another observations indicate that some Type II spicules also show an increase or a more complex velocity dependence with height (Sekse et al. 2012).
Apart from the upward motion, Type II spicules show swaying or transverse motions at the limb with velocity amplitudes of the order 10-30 km s −1 and periods of 100-500 s (De Pontieu et al. 2007b;Tomczyk et al. 2007;Zaqarashvili & Erdélyi 2009;McIntosh et al. 2011;Sharma et al. 2017), suggesting generation of upward, downwards and standing Alfvén waves (Okamoto & De Pontieu 2011;Tavabi et al. 2015), the generation of MHD kink mode waves or Alfvén waves due to magnetic reconnection (Nishizuka et al. 2008;He et al. 2009;McLaughlin et al. 2012;Kuridze et al. 2012). Also, Suematsu et al. (2008) suggest that some spicules show multi-thread structure as result of possible rotation. Another possible motions that Type II spicules show are the torsional motions as suggested in (Beckers 1972), and established using high-resolution spectroscopy at the limb . According to the latter, Type II spicules show torsional motions with 25-30 km s −1 speeds.
There are observational results and theoretical models for the Type II spicules, however our understating of their physical origins remains limited. Some possibilities are that Type II spicules are due to magnetic reconnection (Isobe et al. 2008;De Pontieu et al. 2007b;Archontis et al. 2010;González-Avilés et al. 2017a), oscillatory reconnection processes (Heggland et al. 2009;McLaughlin et al. 2012), strong Lorentz force (Martínez-Sykora et al. 2011;Goodman 2012) or propagation of p-modes (de Wijn et al. 2009). More recently, (Martínez-Sykora et al. 2017) showed that spicules occur when magnetic tension is amplified and transported upward through interaction between ions and neutrals or ambipolar diffusion. The tension is impulsively released to drive flows, heat plasma, and generate Alfvénic waves.
In this paper, we show that 3D magnetic reconnection may be responsible for formation of a jet with characteristics of a Type II spicule. For that (i) we assume a completely ionized solar atmosphere which is governed by the resistive MHD equations subject to a constant gravitational field, (ii) we model the solar atmosphere based on the C7 model in combination with a 3D potential magnetic field configuration extrapolated from a realistic photospheric quiet-Sun model.
The system of equations, the magnetic field configuration, the numerical methods and the model of the solar atmosphere are described in detail in Section 2. The results of the numerical simulations are presented in Section 3. Finally in the Section 4, we present the final comments and conclusions.
2. MODEL AND NUMERICAL METHODS 2.1. The system of Resistive MHD equations We solve the dimensionless Extended Generalized Lagrange Multiplier (EGLM) resistive MHD (Jiang et al. 2012) equations that include gravity: where ρ is the mass density, v is the velocity vector field, B is the magnetic vector field, E is the total energy density and γ = 5/3 is the adiabatic index. 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 aims at damping out the violation of the constraint ∇ · B = 0. Here c h is the wave speed and c p is the damping rate of the wave of the characteristic mode associated with ψ. In this study we consider uniform and constant magnetic resistivity for simplicity. The system of Equations (1)-(7) was normalized by the quantities given in Table 1, which are typical scales in the solar atmosphere.
In the EGLM-MHD formulation, Equation (5) is the mag- Velocity netic 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.1.
In this work we solve the 3D resistive MHD equations with resolutions ∆x, ∆y and ∆z. The gas pressure is computed using the thermal energy, which is obtained by subtracting the kinetic and magnetic energy from the total energy, defined by the total energy Equation (7). In the solar corona region, the plasma-β can become very small, and the thermal energy could be many orders of magnitude smaller than magnetic energy. Therein, small discretization errors in the total energy can produce unphysical negative pressure. We fix this problem by replacing the total energy density Equation (3) in low-beta regions (β ≤ 10 −2 ) with the entropy density equation.
where S = p ρ γ−1 is the entropy density and J 2 = J 2 x + J 2 y + J 2 z . In this way, we calculate the pressure directly using the entropy, which, by definition, is a positive quantity. The entropy density equation is used to maintain the positivity of gas pressure in the context of the ideal MHD simulations (Balsara & Spicer 1999;Li 2008;Derings et al. 2016), and is also used in some resistive MHD simulations of the solar corona (Takasao et al. 2015). In the ideal MHD limit, equation (8) is an advection type of equation, whereas in the case of the resistive MHD equations the Ohmic dissipation is added as a source term. This entropy equation is consistent with the second law of thermodynamics in the continuum limit (Derings et al. 2016).

The magnetic field
As an initial magnetic configuration, we use a 3D potential (current-free) magnetic field extrapolated from a simulated quiet-Sun photospheric field. The latter has been obtained from a large-scale, high-resolution self-consistent simulation of solar magnetoconvection in a bipolar photospheric region with MURaM code (Shelyag et al. 2012;Vögler et al. 2005). The original computational box had a size of 480×480×400 pixels with a spatial resolution of 25 km in the horizontal directions and 10 km in the vertical direction. The initial magnetic field was created as a checkerboard (positive-negative) pattern with the unsigned vertical magnetic field strength of 200 G. This field configuration was inserted into a welldeveloped non-magnetic photospheric convection model and evolved for 20 minutes of physical time. During this simulation phase the magnetic field partially cancelled and partially concentrated in the intergranular lanes forming the intergranular magnetic field concentrations with random polarities and with the strength of ∼ 1.5 kG (Shelyag et al. 2012).
The potential field extrapolation is based on a vectorpotential Grad-Rubin-like method as described in Amari et al. (1997). The potential field extrapolation uses open boundary conditions on the side and top of the computational box: the first derivative of the magnetic field component normal to the surface of the box vanishes. We select a 3D domain of 6×6×10 Mm containing a topology of interest to perform our numerical simulations, the reason is that with such structure it is likely that reconnection may happen and lead to jet generation. The magnetic field lines of the 3D configuration and the magnitude of the magnetic field at z = 0.1Mm are shown on the top panel of Figure 1. The all three components of the magnetic field B x , B y , B z in the plane z =0.1 Mm are shown in the bottom panels, where dipolar structures can be observed. In our convention the xy plane is horizontal and z labels height. These plots show the region used to simulate the evolution of the system, which contains magnetic dipoles at around the location (x, y, z)∼(1.4, 2.3, 0.1) Mm.

Numerical methods
The implementation is the same High Resolution Shock Capturing method as used in (González-Avilés et al. 2017a), based on finite volume approximation. However, in the present paper we exploit the full three-dimensional capabilities of the Newtonian CAFE code (González-Avilés et al. 2015). A summary of the specific numerical methods is as follows. We solve numerically the system of Equations (1)-(8) on a uniform cell centered grid, using the method of lines with a third order Runge-Kutta time integrator (RK3) (Shu & Osher 1989). The discretization of the resistive MHD equations above is based on finite volume approximation. We use the MINMOD and MC limiters for the flux reconstruction, and a combination of the HLLE and HLLC approximate flux formulas (Einfeldt 1988;Harten et al. 1983;Li 2005). The combinations of limiters and flux formulas is adaptive and depends on the magnitude of the discontinuities and shocks formed during the evolution, using the maximum dissipative combination MINMOD-HLLE in zones where β < 10 −2 and the least dissipative combination MC-HLLC otherwise.

Model of the solar atmosphere
We choose the numerical domain to cover part of the interconnected solar photosphere, chromosphere and corona (see top left panel of Figure 1 and Figure 2). For this the atmosphere is initially assumed to be in hydrostatic equilibrium. The temperature field is considered to obey the semiempirical C7 model of the chromosphere-transition region (Avrett & Loeser 2008) and is distributed consistently with observed line intensities and profiles from the SUMER atlas of the extreme ultraviolet spectrum (Curdt et al. 1999). The photosphere is extended to the solar corona as described by Fontela et al. (1990) and Griffiths et al. (1999). The temperature T (z) and density ρ(z) as functions of height z are shown in Figure 2, where the transition region is characterized by the steep gradients.

RESULTS OF NUMERICAL SIMULATIONS
We carried out a numerical simulation within a specific domain with magnetic fields constructed with the MURaM code, which contained a region with a high magnetic field strength dipoles. We define the numerical domain to be x ∈ [0, 6], y ∈ [0, 6], z ∈ [0, 10] Mm, covered with 240×240×400 grid cells, i.e., the effective resolution is 25 km in each direction. In the faces of the numerical box we set fixed in time boundary conditions, which keep the value of the variables set to their initial condition value at a ghost boundary three ghost cells out from the six faces of the physical boundary.
Once we set the magnetic field and the atmosphere model described above, we start evolving the plasma according to the Equations (1-8). We do not apply any explicit perturbation to the system, instead, the round-off errors suffice to trigger the instability of the whole system, including the magnetic field and hydrodynamic equilibria, that later on traduces into the burst of material upwards. The reconnection happens and is accompanied by the introduction of a finite magnetic resistivity η = 12.56 Ω m.
We focus on the process of jet formation and track the temperature evolution that helps understanding the dynamics of the system. In Figure 3 we show snapshots of temperature on the plane x = 2.5 Mm and the magnetic field lines in 3D at different times. For instance, at time t = 15 s the jet-like structure starts to develop in the region of magnetic reconnection which accelerates the plasma. Between t = 30 s and t = 45 s, the jet continues to develop and moves upwards. The most representative time of the jet formation is t = 60 s, at this time we can see a structure with a similar morphology of a Type II spicule, which reaches a height of about z ≈7 Mm measured from the transition region (Tavabi et al. 2015) and vertical velocity of about v z ≈ 130 km s −1 as shown in Figure 6, these characteristics are similar to those of a Type II spicule (De Pontieu et al. 2007a). At time t = 90 s, the spicule-like structure reaches the top of the domain located at z = 10 Mm.
We show a 2D perspective of the process with a cut of the 3D domain at the plane x = 0.1 Mm in Figure 4, where various snapshots of the evolution of the temperature (in Kelvin) and the magnetic field lines are shown. For instance, at time t = 15 s the jet starts to develop at the transition region level z ≈ 2.1 Mm where there is a strong current density, which may be an indication of reconnection happening. The location of the exact reconnection process turns out to be crucially different (see e. g. (Pontin 2012)). Between t = 30 s and t = 45 s the jet continues to form. At time t = 60 s a jet with features of a Type II spicule appears with a basis located at z ≈2 Mm and reaches a height of about z ≈ 7 Mm measured from transition region (see Figure 4), which is in agreement with the observed heights of the Type II spicules, between 3-9 Mm (Pereira et al. 2012;Tavabi et al. 2015). The structure of the spicule obtained at time t = 60 s is similar to the obtained in Figure 5 in Martínez-Sykora et al. (2011). At time t = 90 s the spicule reaches the top of the domain and the magnetic field lines tend to be uniform.
In order to locate regions where magnetic reconnection can take place, we show 2D perspectives of the evolution of |J| (A m −2 ) and temperature contours (K) in Figure 5. For instance, at time t = 15 s, which is the time when the spicule starts to develop, we can see regions of strong current density located the transition region and chromosphere. Between t = 30 s and t = 45 s the stronger current density regions are located at the basis of the spicule, which can accelerate the plasma upwards.  At time t = 60 s, when the spicule is well formed, the stronger current locates around (y, z) ∼(2,2) Mm, at the basis of the spicule, which is consistent with the results shown in Figure 4 at t = 60 s. At the next two snapshots t = 75 s and t = 90 s, the regions of stronger current are still located at the bottom of the spicule. This analysis shows that magnetic reconnection mainly happens at the chromosphere and transition region.
As it has been reported in a number of observational papers (De Pontieu et al. 2007a;Anan et al. 2010;Pereira et al. 2012;Zhang et al. 2012), the upward velocity of spicules is important, therefore we monitor this quantity in our simulation. For the analysis, we show 2D maps of the vertical velocity v z (km s −1 ), the vector velocity field and temperature contours (K) in Figure 6. At time t = 15 s, the spicule is moving upwards with a maximum vertical velocity v z ∼ 190 km s −1 . At time t = 30 s, the spicule continues to move upwards with a velocity of the order v z ∼ 178 km s −1 . At time t = 60 s, the maximum vertical velocity of the spicule is of the order v z ∼ 148 km s −1 , which is slightly above the range of observed upward velocities of a Type II spicule. At times t = 75 s and t = 90 s, the velocity reaches a value of v z ∼ 116 km s −1 at the top of the domain.
In order to understand the physics behind the modeled spicule formation, it is important to identify the dominant force(s) acting during the formation and development of the spicule. For this we compare the forces due to the magnetic field and hydrodynamics, thus we calculate the ratio between the magnitude of the Lorentz force and the magnitude of pressure gradient |J × B|/|∇p|. The results of the evolution of the ratio |J × B|/|∇p| and temperature contours (K) are shown in Figure 7 in the cross cut at the plane x =0.1 Mm of the 3D domain. Notice that at time t = 15 s, which is the time when the    spicule starts to form, that Lorentz force dominates. At times t = 30 s and t = 45 s the dominance of the Lorentz force helps to the spicule moving upwards. At time t = 60 s, the Lorentz force is stronger exactly where the spicule forms. This dominance is still clear at times t = 75 s and t = 90 s. This analysis shows that Lorentz force is an important ingredient of the jet formation.
In order to identify more clearly the behavior of the magnetic reconnection process, we calculate the velocity and magnetic field components, the gas pressure and mass density as functions of y along the constant line z = 2.1 Mm, which is the location of the base of the jet (see Figure 4) and where the current density is strong. For instance in Figure 9 we show the velocity and magnetic field components along this line at time t = 45 s, we can see that v x and v y change sign indicating a bidirectional flow, which is characteristic of a current sheet region. We can see at the bottom of Figure 9 that magnetic field components B x , B y and B z also change sign, in particular the vertical magnetic field component B z indicates a current sheet region. We also analyze the behavior of the gas pressure p and mass density ρ at time 45 s in Figure 10, measured along the line z =2.1 Mm. These plots show an increase in density and pressure near the reconnection region.
In addition, we estimate the ratio between magnetic energy density E mag = |B| 2 2µ0 and kinetic energy density E kin = ρv 2 2 at the point A = (0.1, 1.75, 2.1) Mm shown in the right panel of Figure 8, which is located in a region where the reconnection can be triggered. From the estimation we obtain that magnetic energy density is being converted into kinetic energy during the evolution, which is an indication of a reconnection process.
Another important diagnostics of Type II spicules is whether they are twisted, rotate or show torsional flows. Observations on the Doppler shift of various emission lines in the limb suggest that Type II spicules are rotating Sekse et al. 2013;Sharma et al. 2017). Thus we calculate the vorticity ω = ∇ × v and the vector velocity field in order to look for rotational motion in the spicule region. For this we consider the plane at z = 5 Mm located approximately to the middle of the spicule. We show the magnitude of ω, velocity field and temperature contours (K) in Figure 11. By t = 15 s we can see regions where the magnitude of vorticity is high, the vector velocity field starts to circulate and the temperature is low. At time t = 60 s we can see a region with a high value of the vorticity and a low value of the temperature. This vortex is related to the motion of the spicule structure.
Vorticity and Doppler. We estimate in the plane defined above the Doppler effect related to the dynamic of the spicule in a simple way. We specifically estimate this effect in a small region where the vorticity is high, the velocity vector field is circulating and the temperature is low. In order to estimate the Doppler effect we define a center in the region mentioned above where the velocity is v c . Then we chose a set of points to the left and to the right along the x direction from the center (it could have been any other), with velocities v L and v R , respectively. Then we calculate the difference in the y component of these velocities with respect to that of the center, specifically ∆v D = v L,R − v c , where v L,R , which is an estimate of the tangent velocity of the points around the center and therefore a measure of a red and blue shift. This method is illustrated in Figure 12. We show a zoom in of the vortex where the circulation of the vector velocity field is more evident. In this particular case we calculate a plot of ∆v D for the y component of the velocity v y as a function of the distance d c from the center to the right or left, along the blue or red line. The amplitude of the red shift is of the order of 15 km s −1 , whereas the blue shift has an amplitude of the order 25 km s −1 . The results of the estimation of the Doppler effect due to tangent motion ∆v Dy are shown also in Figure 12. 4. CONCLUSIONS In this paper we have presented a 3D numerical simulation on a small region of the solar atmosphere, showing the formation of a jet structure with characteristics of a Type II spicule, specifically the morphology, upward velocity range and timescale formation. This result provides a simple explanation and is in contrast with that in Martínez-Sykora et al. (2017), where out of 2D simulations the formation of spicules is explained in terms of the amplification of the magnetic tension and the interaction between ions and neutrals. In our simulation we show that even if magnetic tension might be important, the magnetic pressure, which is a part of full Lorentz force is important as well, which is consistent with the results obtained in the simulations of vortex tubes (Kitiashvili et al. 2013) and in the formation of solar chromospheric jets (Iijima & Yokoyama 2017). A quantitative distinction between the components of the different forces involved, would require the development of new analysis tools for time-dependent structures.
For this, we solve the equations of the resistive MHD submitted to the solar constant gravitational field. We use a 3D magnetic field configuration extrapolated up to the solar corona region from a simulated quiet-Sun photospheric field. This magnetic field configuration contains bipolar regions with a strong magnetic field strength at the bottom, which helps the development of the magnetic reconnection process from the photospheric level.
A key result of our analysis is that the Lorentz force dominates over the pressure gradient in the region where the spicule takes place and helps accelerating the structure upwards. It is also expected that the pressure gradient at the transition region contributes at accelerating the plasma upwards.
This 3D model, reveals the complexity, since a solar atmosphere containing the transition region in combination with a magnetic field with a complex topology sketch better the complexity of the solar atmosphere.
Our findings include also that the vorticity near the spicule is important. By looking at the velocity field in a specific cross-cut of the spicule we can track the circular displacement of plasma that eventually can be identified with blue-red shifts. A detailed analysis on the torsional properties of the spicule, generated waves, rotational and radial displacements will be presented in a separate paper (González-Avilés et al. 2017b).
In order to contrast our simulations with other similar analyses, we mention that our simulations are limited in the sense that we do not consider thermal conductivity, radiation and partial ionization as in Martínez-Sykora et al. (2017), however our simulation uses a topologically complex magnetic field in full 3D.  their financial support.