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

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 into contact (see, e.g., Shibata et al. 2007). The magnetic reconnection acts as a mechanism of conversion of the magnetic 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 1968), which are ubiquitous, highly dynamic jets of plasma (Secchi 1878; De Pontieu et al. 2007a; Tsiropoula et al. 2012). Improvements in the resolutions of the observations by the Hinode satellite and the Swedish 1 m Solar Telescope on La Palma (Scharmer et al. 2003; 2008 Kosugi et al. 2007) have suggested the existence of two classes of spicules.

The first type of spicule is so-called Type I, which reaches maximum heights of 4–8 Mm, maximum ascending velocities of 15–40 km s⁻¹, has lifetimes of 3–6.5 minutes (Pereira et al. 2012), and shows upward and downward motions (Beckers 1968; Suematsu et al. 1995). Type I spicules are probably the counterparts 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 mangnetoacoustic 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 spicule (Type II) reaches maximum heights of 3–9 Mm (longer in coronal holes) and has shorter lifetimes of 50–150 s, compared to Type I spicules (De Pontieu et al. 2007b; Pereira et al. 2012). Type II spicules show apparent upward motions with speeds of the order of 30–110 km s⁻¹. At the end of their lives they usually exhibit rapid fading in chromospheric lines (De Pontieu et al. 2007a). 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 10 minutes (Pereira et al. 2014; Skogsrud 2015). Also, in Zhang et al. (2012) the authors established 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 (De Pontieu et al. 2009, 2011). Other 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 of 10–30 km s⁻¹ and periods of 100–500 s (De Pontieu et al. 2007a; Tomczyk et al. 2007; Zaqarashvili & Erdélyi 2009; McIntosh et al. 2011; Sharma et al. 2017), suggesting generation of upward, downward, and standing Alfvén waves (Okamoto & De Pontieu 2011; Tavabi et al. 2015), and the generation of MHD kink mode waves or Alfvén waves due to magnetic reconnection (Nishizuka et al. 2008; He et al. 2009; Kuridze et al. 2012; McLaughlin et al. 2012). Suematsu et al. (2008) suggested that some spicules show multi-thread structure as a result of possible rotation. Other possible motions that Type II spicules show are the torsional motions suggested in Beckers (1972) and established using high-resolution spectroscopy at the limb (De Pontieu et al. 2012). According to the latter, Type II spicules show torsional motions with 25–30 km s⁻¹ speeds.

There are observational results and theoretical models for the Type II spicules, but our understanding of their physical origins remains limited. Some possibilities are that Type II spicules are due to magnetic reconnection (De Pontieu et al. 2007a; Isobe et al. 2008; Archontis et al. 2010; González-Avilés et al. 2017), 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 the formation of a jet with the characteristics of a Type II spicule. We assume that (i) a completely ionized solar atmosphere that is governed by the resistive MHD equations is subject to a constant gravitational field, and (ii) we can 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. In Section 4, we present our final comments and conclusions.

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:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\cdot \left(\left(p+\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2}\right){\boldsymbol{I}}+\rho {\boldsymbol{vv}}-{\boldsymbol{BB}}\right)\\ &&\qquad =\,-({\rm{\nabla }}\cdot {\boldsymbol{B}}){\boldsymbol{B}}+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\cdot \left({\boldsymbol{v}}\left(E+\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2}+p\right)-{\boldsymbol{B}}({\boldsymbol{B}}\cdot {\boldsymbol{v}})\right)\\ &&\quad =\,-{\boldsymbol{B}}\cdot ({\rm{\nabla }}\psi )-{\rm{\nabla }}\cdot ((\eta {\boldsymbol{J}})\times {\boldsymbol{B}})+\rho g\cdot {\boldsymbol{v}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+{\rm{\nabla }}\cdot (B{\boldsymbol{v}}-{\boldsymbol{vB}}+\psi {\boldsymbol{I}})=-{\rm{\nabla }}\times (\eta {\boldsymbol{J}}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \psi }{\partial t}+{c}_{h}^{2}{\rm{\nabla }}\cdot {\boldsymbol{B}}=-\displaystyle \frac{{c}_{h}^{2}}{{c}_{p}^{2}}\psi ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\boldsymbol{J}}={\rm{\nabla }}\times {\boldsymbol{B}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&E=\displaystyle \frac{p}{(\gamma -1)}+\displaystyle \frac{\rho {{\boldsymbol{v}}}^{2}}{2}+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{2},\end{eqnarray} \tag{ 7 }$$

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.

Table 1.  Normalization Units

Variable	Quantity	Unit	Value
x, y, z	Length	l0	106 m
ρ	Density	ρ0	10−12 kg m−3
B	Magnetic field	B0	11.21 G
v	Velocity	${v}_{0}={B}_{0}/\sqrt{{\mu }_{0}{\rho }_{0}}$	106 m s−1
t	Time	t0 = l0/v0	1 s
η	Resistivity	η0 = l0μ0v0	1.256 × 106 m2 s−1 N A−2

Download table as:  ASCIITypeset image

In the EGLM-MHD formulation, Equation (5) is the magnetic field, divergence-free constraint. As suggested in Dedner et al. (2002), the expressions for c_h and c_p are

$$\begin{eqnarray*}&&{c}_{h}=\displaystyle \frac{{c}_{{cfl}}}{{\rm{\Delta }}t}\min ({\rm{\Delta }}x,{\rm{\Delta }}y,{\rm{\Delta }}z),\,\,\,{c}_{p}=\sqrt{\left|\displaystyle \frac{-{\rm{\Delta }}t}{\mathrm{ln}{c}_{d}}\right|{c}_{h}^{2}},\end{eqnarray*}$$

where Δ t is the time step, Δx, Δy, and Δz are the spatial resolutions, c_cfl < 1 is the Courant factor, and c_d is a problem-dependent coefficient between 0 and 1, a constant that 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, and for that reason one should adjust their values. In our simulations we use ${c}_{p}=\sqrt{{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 can be many orders of magnitude smaller than the 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⁻²) with the entropy density equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial S}{\partial t}+{\rm{\nabla }}\cdot (S{\boldsymbol{v}})=(\gamma -1){\rho }^{1-\gamma }\eta {{\boldsymbol{J}}}^{2},\end{eqnarray} \tag{ 8 }$$

where $S=\tfrac{p}{{\rho }^{\gamma -1}}$ is the entropy density and ${{\boldsymbol{J}}}^{2}={J}_{x}^{2}+{J}_{y}^{2}+{J}_{z}^{2}$. 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).

2.2. 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 the MURaM code (Vögler et al. 2005; Shelyag et al. 2012). 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 an unsigned vertical magnetic field strength of 200 G. This field configuration was inserted into a well-developed 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 a strength of ∼1.5 kG (Shelyag et al. 2012).

The potential field extrapolation is based on a vector-potential, 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. We do this because 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.1 Mm are shown in the top panel of Figure 1. 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.

Zoom In Zoom Out Reset image size

2.3. Numerical Methods

2.4. Model of the Solar Atmosphere

We choose the numerical domain to cover part of the interconnected solar photosphere, chromosphere, and corona (see the top left panels of Figures 1 and 2). The atmosphere is initially assumed to be in hydrostatic equilibrium. The temperature field is considered to obey the semi-empirical 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 steep gradients.

Zoom In Zoom Out Reset image size

We carried out a numerical simulation within a specific domain with magnetic fields constructed with the MURaM code, which contained a region with 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 values of the variables set to their initial condition values at a ghost boundary that is 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 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, which later on traduces into the burst of material upward. 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 with 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 that accelerates the plasma. Between t = 30 s and t = 45 s, the jet continues to develop and moves upward. The most representative time of the jet formation is t = 60 s; at this time we can see a structure with a similar morphology to a Type II spicule, which reaches a height of about z ≈ 7 Mm measured from the transition region (Tavabi et al. 2015), and a vertical velocity of about v_z ≈ 130 km s⁻¹, as shown in Figure 6. These characteristics are similar to those of a Type II spicule (De Pontieu et al. 2007b). At time t = 90 s, the spicule-like structure reaches the top of the domain located at z = 10 Mm.

Zoom In Zoom Out Reset image size

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 the 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 the transition region (see Figure 4), which is in agreement with the observed heights of the Type II spicules, between 3 and 9 Mm (Pereira et al. 2012; Tavabi et al. 2015). The structure of the spicule obtained at time t = 60 s is similar to that 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.

Zoom In Zoom Out Reset image size

In order to locate regions where magnetic reconnection can take place, we show 2D perspectives of the evolution of $| {\boldsymbol{J}}|$ (A m⁻²) 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 at 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 upward. 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.

Zoom In Zoom Out Reset image size

As it has been reported in a number of observational papers (De Pontieu et al. 2007b; 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⁻¹), the vector velocity field, and temperature contours (K) in Figure 6. At time t = 15 s, the spicule is moving upward with a maximum vertical velocity v_z ∼ 190 km s⁻¹. At time t = 30 s, the spicule continues to move upward with a velocity of the order of v_z ∼ 178 km s⁻¹. At time t = 60 s, the maximum vertical velocity of the spicule is of the order of v_z ∼ 148 km s⁻¹, 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⁻¹ at the top of the domain.

Zoom In Zoom Out Reset image size

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 the pressure gradient $| {\boldsymbol{J}}\times {\boldsymbol{B}}| /| {\rm{\nabla }}p|$. The results of the evolution of the ratio $| {\boldsymbol{J}}\times {\boldsymbol{B}}| /| {\rm{\nabla }}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. Note 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 the spicule move upward. 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.

Zoom In Zoom Out Reset image size

In order to identify more clearly the behavior of the magnetic reconnection process, we estimate the ratio between magnetic energy density ${E}_{\mathrm{mag}}=\tfrac{| {B}^{2}| }{2{\mu }_{0}}$ and kinetic energy density ${E}_{\mathrm{kin}}\,=\tfrac{\rho {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.

In addition, we calculate the velocity, 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 (see left panel of Figure 8). 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 signs, 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Another important diagnostic of Type II spicules is whether they are twisted, rotating, or have torsional flows. Observations on the Doppler shift of various emission lines in the limb suggest that Type II spicules are rotating (De Pontieu et al. 2012; 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. We consider the plane at z = 5 Mm located approximately at the middle of the spicule. We show the magnitude of ω, the 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 vorticity value and a low temperature. This vortex is related to the motion of the spicule structure.

Zoom In Zoom Out Reset image size

Vorticity and Doppler. In the plane defined above we estimate 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 redshift and blueshift. This method is illustrated in Figure 12. We show a zoomed-in view 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 redshift is of the order of 15 km s⁻¹, whereas the blueshift has an amplitude of the order of 25 km s⁻¹. The results of the estimation of the Doppler effect due to tangent motion Δv_Dy are also shown in Figure 12.

Zoom In Zoom Out Reset image size

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 the characteristics of a Type II spicule, specifically the morphology, upward velocity range, and formation timescale. This result provides a simple explanation and is in contrast with that in Martínez-Sykora et al. (2017), where, using 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 is important, the magnetic pressure, which is related to full Lorentz force being 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 with 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 accelerate the structure upward. It is also expected that the pressure gradient at the transition region contributes to accelerating the plasma upward.

This 3D model reveals complexity, since a solar atmosphere containing the transition region, in combination with a magnetic field with a complex topology, more effectively render the complexity of the solar atmosphere.

Our findings also indicate 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 blueshift and redshift. A detailed analysis on the torsional properties of the spicule, generated waves, and rotational and radial displacements will be presented in a separate paper (J. J. González-Avilés et al. 2018, in preparation).

In order to contrast our simulations with similar analyses, we note 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.

This research is partly supported by the following grants: Royal Society-Newton Mobility grants NI160149, CIC-UMSNH 4.9, and CONACyT 258726 (Fondo Sectorial de Investigación para la Educación). The simulations were carried out in the facilities of the CESCE-UNAM, Iceberg HPC Cluster and the Big Mamma cluster at the LIASC-IFM. V.F. and G.V. thank the STFC for their financial support. Visualization of the simulations data was done with the use of the VisIt software package.