Spicule jets in the solar atmosphere modeled with resistive MHD and thermal conduction

Using numerical simulations, we study the effects of magnetic resistivity and thermal conductivity in the dynamics and properties of solar jets with characteristics of Type II spicules and cool coronal jets. The dynamic evolution of the jets is governed by the resistive MHD equations with thermal conduction along the magnetic field lines on a 2.5D slice. The magnetic field configuration consists of two symmetric neighboring loops with opposite polarity, used to support reconnection and followed by the plasma jet formation. In total 10 simulations were carried out with different values of resistivity and thermal conductivity, that produce jets with different morphological and thermal properties we quantify. We find that an increase in magnetic resistivity does not produce significant effects on the morphology, velocity and temperature of the jets. However, thermal conductivity affects both temperature and morphology of the jets. In particular, thermal conductivity causes jets to reach greater heights and increases the temperature of the jet-apex. Also, heat flux maps indicate the jet-apex and corona interchange energy more efficiently than the body of the jet. These results could potentially open a new avenue for plasma diagnostics in the Sun's atmosphere.


INTRODUCTION
Solar spicules are small-scale, jet-like plasma features observed ubiquitously in the solar chromosphere (Beckers 1972;Sterling 2000;De Pontieu et al. 2004, 2011. Spicules may play an important role in energy and material supply to the upper layers of the solar atmosphere (De Pontieu et al. 2011;Samanta et al. 2019).
Similarly, there are obseva-tions that highlight the evidence of magnetic reconnection in the generation of jets (Borrero et al. 2013;Martínez Pillet et al. 2011). In particular, Type II spicules are collimated jets that reach maximum heights of 3-9 Mm and even longer in coronal holes, have a typical lifetime of 50-150 s (De Pontieu et al. 2007a;Pereira et al. 2012), although these features can survive up to 500-800 s, while appearing and dissapearing in multiple chromospheric passbands, due to thermal evolution (Pereira et al. 2014). Type II spicules show apparent upward motions with speeds of order 30-110 km s −1 and temperatures of order 10 4 K (Sterling 2000;Sterling et al. 2010). At the end of their life they usually exhibit rapid fading in chromospheric lines (De Pontieu et al. 2007b, 2017a. Numerical modelling is an important way for styling and analyzing various dynamical plasma processes in the solar atmosphere. In particular, it is a powerful tool for better understanding of transient phenomena such as jets. For instance Takasao et al. (2013), study the acceleration mechanism of chromospheric jets associated with emerging fluxes using 2D MHD simulations. Similarly, there are more sophisticated models of jet formation in 2D and 3D, which have been performed by (Isobe et al. 2006;Pariat et al. 2009;Archontis et al. 2010;Jiang et al. 2012a), that take into account physical effects such as magnetic resistivity and thermal conductivity. In particular, models including magnetic resistivity and thermal conductivity, are close to realistically describe the conditions in the solar atmosphere, as examples of this type of simulations we have (Botha et al. 2011), where the authors use 3D MHD simulations to show that thermal conduction plays an essential role in the kink instability of coronal loops and cannot be ignored. Another example can be found in Fang et al. (2014), where the authors study the formation of coronal jets through the numerical simulation of the emergence of a twisted flux rope and found that field-aligned thermal conduction efficiently distributes the energy release, which is essential for comparing with synthetic emission. Apart from the ingredients of resistivity and thermal conductivity, there are other more sophisticated numerical simulations of Type II spicule formation that include the effect of radiation, partial ionization and ambipolar diffusion (Martínez-Sykora et al. 2009, 2011, 2017aDe Pontieu et al. 2017a).
In this paper we continue the analysis carried out by González-Avilés et al. (2017), by including the thermal conductivity flux term in the resistive MHD equations, in addition to that we consider a wider range of resistivity values. In particular, we use a set of 4 realistic values of magnetic resistivity and thermal conductivity to analyze their effects on the morphology, maximum height, vertical velocity, thickness, temperature of the jet-apex and lifetime of the jets modelled. Aside of state of the art contribution in simulations, like partial ionization, radiation and ambipolar difussion such as in De Pontieu et al.  Yang et al. (2013), we show that thermal conductivity can modify the temperature, maximum height and width of the jets with some characteristics of Type II spicules and cool coronal jets. The carrying out the study in 2.5D allows us to control computing time in a flexible way and in the use of various parameter combinations to analyze the numerical simulations.
The paper is organized as follows. In Section 2 we describe the resistive MHD equations with thermal conduction, the numerical methods we use, the model of solar atmosphere and the magnetic field configuration.
In Section 3 we describe the parameters analyzed and the results of the numerical simulations for various experiments. Section 4 contains conclusions and final comments.

MODEL AND METHODS
2.1. The system of resistive MHD equations with thermal conduction The model we consider to drive the plasma dynamics and the jet formation is the resistive MHD with thermal conduction and in particular we choose the Extended Generalized Lagrange Multiplier (EGLM) (Jiang et al. 2012b) to be effective at keeping the evolution of jets under control (González-Avilés et al. 2017. The system of equations we use for the evolution of the plasma is given in Jiang et al. (2012b), whose dimensionless version reads as follows: where ρ is the mass density, v is the velocity field, B is the magnetic field, E is the total energy density, the plasma pressure p is described by the equation of state of an ideal gas p = (γ − 1)ρe, where e is the internal energy and γ its adiabatic index, g is the gravitational field at the solar surface, J is the current density, η is the magnetic resistivity and ψ is a scalar potential that helps damping out the violation of the constraint ∇ · B = 0.
Here c h is a wave speed and c p is the damping rate of the wave of the characteristic mode associated to ψ. In our simulations we use c p = √ c r c h , with c r =0.18 and c h =0.001, that have shown useful for 2.5D and 3D simulations.
The contribution of thermal conduction is included in the equation for the energy (3), through the heat flux vector that allows the heat propagation along the magnetic field lines (see e.g. Jiang et al. 2012a) where κ is the thermal conductivity of the plasma, T its temperature and B is the magnetic field.

Numerical methods
The evolution of jets is analyzed with a 2.5D approach, which consists in the restriction of the dynamics under consideration to have slab symmetry along one of the two horizontal directions of the spatial domain. We use Cartesian coordinates x and z describing the twodimensional domain, and y is the direction of the slab symmetry.
We solve numerically the resistive EGLM-MHD equations with thermal conduction given by the system of equations (1)-(5) using the Newtonian CAFE code (see e.g. González-Avilés et al. 2015; González-Avilés & Guzmán 2018), on a uniform cell centered grid, using the method of lines with a third order total variation diminishing Runge-Kutta time integrator described in Shu & Osher (1989). In order to use the method of lines, the right hand side of equations (1)-(5) are discretized using a finite volume approximation with High Resolution Shock Capturing methods, e.g., LeVeque (1992). For this, we first reconstruct the variables at cell interfaces using the minmod limiter, and numerical fluxes are calculated using the Harten-Laxvan Leer-Contact (HLLC) approximate Riemann solver (see e.g. Li 2005).

Model of the solar atmosphere and magnetic field configuration
At initial time of simulation, we assume the solar atmosphere to be in hydrostatic equilibrium, with the temperature field consistent with the semiempirical C7 model of the chromosphere transition region (Avret & Loeser 2008). The density and temperature profiles used to start the simulations are shown at the top of Figure 1. For the equation of state we assume an adiabatic index γ = 5/3 and the gravitational field is set to g = −gẑ with g = 274m s −2 in the equations of momentum and energy, more details can be found in González-Avilés et al. (2017).
The magnetic field configuration is a superposition of two neighboring loops and is constructed from a potential. Based on Priest (1982); Del Zanna et al. (2005) and González-Avilés et al. (2017), the magnetic field potential for two symmetric loops, that decreases with height exponentially, is given by where l 0 parametrizes the position of the foot points for each of the loops and B 0 is the magnetic field strength.
In this paper, we use the parameters B 0 = 40 G and l 0 = 3.5 Mm, because they are physically sound and produce jets with the some properties of Type II spicules successfully (González-Avilés et al. 2017). In this manner we concentrate on the properties of the ejected structure. A schematic picture of the magnetic configuration is shown to the bottom of Figure 1.
Further details of the set up are that we fix k = π/L with L = 8 Mm, where L is the distance between the two footpoints of the loop. Then we execute simulations in the 2D domain x ∈ [−4, 4] × z ∈ [0, 10] in units of Mm, covered with 300×375 grid cells in x and z direction correspondingly. Since we are using a three dimensional code, we cover the additional direction y with four cells. The boundary conditions used are outflow at all faces of the domain.

Parameters analyzed
We study the effects of magnetic resistivity and thermal conductivity on jet formation process with characteristics of Type II spicules. The following parameters were analyzed in details for different runs: the maximum height, width of jet-apex, average temperature of the jet-apex at maximum height and time when the jet reaches the maximum height obtained with the various combinations of the resistivity η and thermal conductivity κ.
With the various combinations of these parameter values we define 10 cases specified in Table 1, the combinations were chosen taking into account the level of realism of the values. For ease, in what follows, we will use the code unit values of resistivity and thermal conductivity to define and describe the various scenarios in our analysis.

Results of numerical simulations
Of all the runs for simulations and summarized in Table 1, we select 3 illustrative ones, e.g. Run #1: η = 3.97 × 10 −8 , κ = 0, Run #6: η = 3.97 × 10 −8 , κ = 3275 and Run #10: η = 8 × 10 −8 , κ = 32755, which are highlighted in light gray. In Figure 2 we show snapshots of the temperature, vertical velocity with the vector field distribution and the magnitude of the heat flux |q| given by equation (6) with magnetic field lines at time t = 210 s for the three illustrative Runs. For example, in Figures 2(a), (b) and (c) we show the results for the values η = 3.97 × 10 −8 and κ = 0, which are practically the same snapshots corresponding to the top of Figure 3 of the paper González-Avilés et al. (2017), these snapshots will be useful for comparison with cases where thermal conductivity is included. Figure 2(c) shows the magnetic field lines with |q| = 0, corresponding to κ = 0, for completeness.
In Figures 2(d), (e) and (f) we show the results for Run #6, with the combination of parameters η = 3.97×10 −8 and κ = 3275. In this case according to the temperature map, we can see that the jet is wider in its lower part and thinner in the upper part just below its apex compared to the Run #1, and reaches a height of 7.5 Mm, which represents 0.2 Mm larger than the jet of Run #1. In addition, the vertical velocity is higher at the sides of the jet compared to that of Run #1. The magnitude of the heat flux |q| is high at the top of the jet-apex, and at the sides beneath the jet-apex, while inside the jet structure the value tend to zero. The reason is that q is a basically a projection of ∇T along the field lines, then it is maximum when ∇T is large has an important component parallel to field lines. At the top of jet-apex the gradient of temperature is nearly radial (i.e perpendicular to the ball shape of the jet-apex) as seen in Fig.  2(a) and magnetic field lines have a radial component; another case is that heat flux is nearly zero aside the jet from z between 2 to 7 Mm because the gradient of T is nearly horizontal whereas the magnetic field lines are nearly vertical and therefore effect of thermal conductivity is minimal. In Figures 2(g), (h) and (i) we show the results for Run #10 with values η = 8 × 10 −8 and κ = 32755, which is the combination that includes the highest value of thermal conductivity. In this case we can also see that the bottom of the jet is about 0.1 Mm wider than the bottom of the jet of Run #1 and jet-apex is 0.28 Mm thinner compared to the Run #1, and 0.10 Mm thinner compared to the Run #6. This jet reaches a height of about 7.7 Mm, which is 0.4 Mm higher compared to Run #1 and 0.2 Mm higher compared to Run #6, and in the same way it is seen that the vertical speed is greater aside the jet. In Figure 2(i), |q| is higher at the at the top of the jet, similar to Figure 2(f), however in this case |q| reaches higher values, which correspond with bigger κ used in this numerical run. The jet during development reaches maximum speeds of the order 100 km s −1 at early times, between 0 and 50 s, while at later times after 60s the speed decreases to values up to 15-30 km s −1 . Therefore in the first stage of the jet's evolution the speeds are comparable to those of Type II spicules, whereas at later times the speeds are smaller than the lower limit of the observed velocities. Similarly, the vertical speed of our jets have similarities with the observed velocities of the Rapid Redshifted and Blueshifted Excursions (RREs, RBEs), which are in the range of 50-150 km s −1 (Langangen et al. 2008). In addition, the vector field shows the appearance of vorticity near the top of the jets. According to the results reported in Table 1, the width of the jet-apex at the maximum height vary in the range 0.8-1.1 Mm, which is four times greater than the width of 0.25 Mm that has been observed in RREs and RBEs spicule features (Kuridze et al. 2015), however the widths of the RREs and RBEs are of the entire observed structure. In fact, if we estimate the width of the jet structure from our simulations, we obtain that the width varies in the range 0.2-0.6 Mm from the bottom of the jet to below its apex, these widths are close to the observed values. Regarding the spicules observed at the limb, it has been estimated cross-sectional widths in the range 0.27-0.36 Mm (Sharma et al. 2018), which are again smaller than the width of the jet-apex, but they are close to the widths estimated in the entire jet structure of our simulations. Similarly to the results obtained in González-Avilés et al. (2017), in this paper we find that jets show a special feature at the apex with a bulb possibly related to the formation of a Kelvin-Helmholtz (KH) type of instability. However as it was shown in González-Avilés et al. (2017), this instability is suppressed by the magnetic field.
To see more clearly the differences in morphology for the cases shown in Figure 2, in Figure 3 we show a zoom of snapshots of temperature, vertical velocity v z together with the velocity field, magnitude of the heat flux |q| with magnetic field lines and the y-component of the vorticity (∇×v) y with the velocity field at the time when the jets are at the same height. For example, in Figures  3(a), (b) and (c) we can see that the jets have different morphology when are at the same height, in particular the jet with the highest value of thermal conductivity is smaller and thinner. The cold material develops a horizontal structure connected to the jet-apex that is more notorious for higher κ, which can also be seen in Figure 2. In Figures 3(d), (e) and (f) we show that the vertical component of velocity is higher on the side of the jet for the cases when thermal conductivity is higher. In Figure 3(g), we show the case of Run #1 when |q| is zero. The magnetic field lines have the shape of the jet shown in Figure 2(a). In Figures 3(h) and (i) we show that |q| is high near to the top of the jet-apex for Runs #6 and #10, where the magnetic field lines have significant horizontal component. We can also see that the magnetic field lines follow the jet structure in Run #1, whereas for the other two cases the magnetic field lines tend to flatten near the jet-apex. This is consistent with the fact that β < 1 outside the jet and β > 1 inside, where hydrodynamical effects dominate. Finally, in Figures 3(j), (k) and (l) we can notice that at the sides of the jet-apex, the flux develops vorticity, which apparently does not change with the increase in thermal conductivity. The appearance of the vorticity could be due to the interaction of the cold gas with the hotter plasma in the corona.
In order to quantify the influence of resistivity and thermal conductivity, we average the temperature and vertical velocity along the jet and show it as functions of time. We also calculate an average of the jet-apex temperature at time t = 210 s. For instance, in Figure  4(a), we show the temperature along the jet as a function of time for different values of η and κ = 0. In this case, we notice that temperature does not show significant variations in all cases. For example, in Figure 4(b) we have found that vertical velocity along the jet is less sensitive to the increase of resistivity. In Figure 4(c), we show the temperature along the jet as a function of time for different values of κ and η = 3.97×10 −8 . In this case, we can see that the temperature increases slightly when the value of the thermal conductivity is higher. In Figure 4(d), we show that the vertical component of velocity along the jet is higher about time t ∼ 100 s when the value of κ is higher, for all other times the velocity remains without significant variations. In Figures 4(b) and 4(d), we can see that the vertical velocity along the jet shows a bump about the time t ∼ 100 s, this is due to the rapid acceleration of the plasma produced by magnetic reconnection, which is triggered by the magnetic loops close together with opposite polarity and the inclusion of resistivity as indicated in González-Avilés et al. (2017). In Figure 4(e), we show the average of jet-apex temperature as a function of η for κ = 0 at time t = 210 s. For the values of η used, the temperature increases by 2.1%. Finally, in Figure 4(f), we show the average of jetapex temperature as a function of κ and η = 3.97×10 −8 . In this case we can see that temperature of the jet-apex increases by ∼ 8.1% with respect to the temperature obtained with the smallest value of conductivity.

CONCLUSIONS AND FINAL COMMENTS
In this paper we perform numerical simulations to explore the effects of magnetic resistivity and thermal conductivity in the dynamics and parameters of jet structures with some characteristics of Type II spicules. This paper is a continuation of the work presented in González-Avilés et al. (2017), where we show that jets with features of Type II spicules and cool coronal jets can be formed as a result of magnetic reconnection in a scenario with magnetic resistivity.
According to the analysis carried out in this paper, we found that the inclusion of thermal conductivity along the magnetic field lines affects the morphology of the jets, in particular the combination κ = 32755 and η = 8×10 −8 makes the jet structure 0.1 Mm wider at its lower part and 0.28 Mm thinner at its apex, measured with respect to jets where only resistivity is considered. In addition, the increase in thermal conductivity makes the jets reach maximum heights of about 7.7 Mm, such heights are 0.4 Mm larger compared to jets with only resistivity. According to the results of the Runs #5 to #8 shown in Table 1, we can see that an increase of the order of 10 times the value of thermal conductivity causes the jet to reach a maximum height 0.2 Mm larger, it also makes the jet-apex 0.28 Mm thinner, this means that the inclusion of thermal conductivity makes the jet have a width closer to the value of observations. One of the most important results is that the jet-apex heat up and reach temperatures of the order of 73000 K for the combination η = 3.97 × 10 −8 and κ = 32755 in comparison with the combination η = 3.97 × 10 −8 and κ = 0, in which the jet-apex temperature is of the order 58000 K. These results are similar to those obtained in Kuźma et al. (2017), however in this paper the formation and evolution of solar spicules use numerical simulations triggered with a vertical velocity pulse that is launched from the upper chromosphere. Instead In panels a, b and c we show the results for Run #1 (η = 3.97 × 10 −8 and κ = 0) at time t = 210 s. In the panels d, e and f we show the results for the Run #6 (η = 3.97 × 10 −8 and κ = 3275) at time t = 210 s. Finally, in the panels g, h and i, we show the results for Run #10 (η = 8 × 10 −8 and κ = 32755) at time t = 210 s. Note that as the value of thermal conductivity increases the jet gets thinner closer to the observations. In this figure we only show temperature maps, but the mass density has the same morphology. Heat flux is particularly consistent, notice it is nearly zero at the body of the jet where ∇T is perpendicular to the field lines, and becomes important at the top of the jet, where ∇T is more parallel to the shape of the jet-apex and field lines nearly horizontal. It is also seen that heat transfer produces a stream of cold material from the jet-apex to the sides. We use different density of magnetic field lines in each case, to capture some of the differences, specially at the jet-apex. Figure 3. From left to right we show a zoom of the region around the jets for Runs #1, Run #6 and Run #10. We show the temperature in Kelvin (panels a, b and c), vertical component of velocity vz in km s −1 with the velocity vector field shown as black arrows (panels d, e and f), magnitude of the heat flux |q| in W m −2 as well as magnetic field lines (panels g, h and i) and y-component of vorticity (∇ × v)y in s −1 (panels j, k and l). Note that we only show the temperature, but the mass density has the same morphology. In all the figures the snapshot is taken when the jet is at the same height.   in our analysis, we do not perturb the solar atmosphere with any pulse at initial time, rather it is the magnetic reconnection that accelerates the plasma and forms the jet. Another difference is that we use a range of resistivity and thermal conductivity values consistent with a fully ionized solar atmosphere. We have also found, that the increase in the resistivity does not affect the morphology of the jets, and although it slightly increases the jet-apex temperature it was found that it does not modify the behavior of the temperature along the whole jet structure.
Note that although the model of two magnetic loops close together with opposite polarity is very simple and approximate, this configuration can support the formation of jets mimicking some properties of Type II spicules and cool coronal jets. Furthermore, the inclusion of resistivity and thermal conductivity is consistent with the physical properties found in the solar atmosphere. In particular, the resistivity supports development of the magnetic reconnection process, and thermal conductivity helps the heat to propagate more efficiently along the magnetic field lines.
These results are important to understanding of nature of spicules. The reason is that two different values of thermal conductivity produce spicules with different temperature and maximum height. From other side, the problem could be degenerate, because two jets with the same temperature and height could be obtained with and without thermal conductivity at the price of modifying for example the magnetic field, or the temperature model of the chromosphere-corona interface. Eventually the addition of ingredients to a model-simulation will have to face the observational restrictions that will in turn refine the parameter values of simulations and the degeneracy of the problem can be an important prob-lem in itself. Finally, it is important to mention the role of heat transfer prior jet propagation, in this case the time scale of physical process related to the thermal conductivity is much smaller than the time needed for jet development and reach the maximum height. At the beginning of simulation runs, the temperature gradient has it's maximum value at the interface, and at the same time the field lines are nearly parallel to this gradient near the foot points, therefore heat transfer influences on initial numerical background in terms of temperature changes even before the jet propagates up. This is important that deserves special attention.