EDP Sciences
Free Access
Issue
A&A
Volume 597, January 2017
Article Number A133
Number of page(s) 8
Section The Sun
DOI https://doi.org/10.1051/0004-6361/201628747
Published online 18 January 2017

© ESO, 2017

1. Introduction

Spicules are thin, cool, and dense structures that are observed in the solar limb (Beckers 1968, 1972; Suematsu 1998; Sterling 2000; Zaqarashvili & Erdélyi 2009). They are seen to emerge from the chromospheric background at an altitude of about 2000 km above the solar surface where they reveal a speed of 25 km s-1, reach a maximum level, and then either disappear or sink down to the chromosphere. A typical lifetime of spicules is within the range of 5–15 min with an average value of ~7 min (Pasachoff et al. 2009). Spicules seem to consist of double-thread structures (Tanaka 1974; Dara et al. 1998; Suematsu et al. 2008), and they reveal bidirectional flows (Tsiropoula et al. 1994; Tziotziou et al. 2003, 2004; Pasachoff et al. 2009). Typical electron temperature and electron density in spicules are (15–17) × 103 K and 2 × 10113.5 × 1010 cm -3 at altitudes of 4–10 Mm above the solar surface (Beckers 1968) with a diameter estimated as 660 ± 200 km (Pasachoff et al. 2009). As a result, spicules are much cooler and denser than ambient coronal plasma. High-resolution observations by the Solar Optical Telescope onboard Hinode have revealed another type of spicules with many features different from those of classical limb spicules, and they are referred to as type II spicules (De Pontieu et al. 2007). The type II spicules are distinguished by (a) smaller diameters (200 km) in the Ca II H line and a significantly shorter height of 4 Mm; (b) a lifetime of 10–150 s; (c) the evolution, which shows an upflow and then disappears; and finally by (d) much higher speeds of 50–100 km s-1 (e.g., De Pontieu et al. 2007, 2009; Rouppe van der Voort et al. 2009; Kuridze et al. 2015). Macrospicules were recently discussed by Nóbrega-Siverio et al. (2016), who considered whether the relevance of the entropy sources in the surges, such as the optically thin losses, can be applied to similar phenomena as macrospicules.

In spite of various theoretical models that have been brought forward to explain the spicule ejection in the lower solar atmosphere, many recent numerical methods have been developed to simulate the solar spicules or macrospicules with an energy input at their base in the photosphere such as a gas pressure pulse or an Alfvén wave that steepens into a shock wave (Sterling 2000, and references therein). Hansteen et al. (2006) and De Pontieu et al. (2007) simulated the formation of dynamic fibrils that are due to slow magneto-acoustic shocks through two-dimensional (2D) numerical simulations. They suggested that these shocks are formed when acoustic waves generated by convective flows and global p-modes in the lower lying photosphere leak upward into the magnetized chromosphere. Heggland et al. (2007) used the initial periodic piston to drive the upward propagating shocks in 1D simulation, and Martinez-Sykora et al. (2009) considered the emergence of new magnetic flux, but the drivers of spicules originate from collapsing granules, energy release in the photosphere, or in the lower chromosphere. However, these simulations were unable to mimic the double structures and bidirectional flows in spicules. On the other hand, Murawski & Zaqarashvili (2010) performed 2D numerical simulations of magnetohydrodynamic (MHD) equations and showed that the 2D rebound shock model of Hollweg (1982) may explain both the double structures and bidirectional flows. They used a single initial velocity pulse, which led to the formation of consecutive shocks as a result of the nonlinear wake in the stratified atmosphere. However, they considered a simple model of atmospheric temperature that was approximated by a smoothed step function for the temperature profile.

Understanding exact drivers of spicules requires further investigation, and more than one mechanism may trigger their evolution depending on the local plasma and magnetic field conditions. The goal of this paper is to contribute to the above-mentioned studies by performing simulations of the generation and evolution of the spicules in the solar atmosphere and compare our results for adiabatic and non-adiabatic MHD equations. The method we chose to trigger a spicule is a localized vertical velocity pulse launched from the upper chromosphere. This method is similar to the calculations performed by Shibata (1982), Sterling et al. (1993), Murawski & Zaqarashvili (2010), and Guerreiro et al. (2013). This approach differs from the models that attempt to model spicules with a disturbance in the photosphere (e.g., Suematsu et al. 1982; Hollweg 1982). Guerreiro et al. (2013) studied the mid-chromospheric energy inputs of earlier simulations by adding additional physics to the radiative loss term and including hydrogen ionization and recombination. They concluded that it would be difficult to produce spicules through those previously suggested mechanisms (specifically that of Sterling et al. 1993). Our simulations do not include the detailed energy losses of Guerrero et al. However, the energy input is different from that assumed by Guerrero et al. and Sterling et al., who adopted localized increase in the heating rate. Since the form of input energy is different, it is not clear whether the losses Guerrero et al. took into account would have any effect on the velocity of pulse-driven spicules. This will be tested in the future studies.

This paper is organized as follows. A numerical model is presented in Sect. 2, and the corresponding numerical results are shown in Sect. 3. Our paper is concluded by a summary of the numerical results in Sect. 4.

2. Physical model of the solar atmosphere

2.1. MHD equations

We consider a gravitationally stratified and magnetically confined solar plasma that is governed by the following set of non-adiabatic MHD equations: where ϱ is the mass density, p the gas pressure, V represents the plasma velocity, B is the magnetic field, T the temperature, q the anisotropic thermal conduction flux, L(ϱ,T) radiatively thin cooling terms (Mignone et al. 2007), and H(ϱe,Te) denotes the external heating therm that balances L and ∇·q at the equilibrium that is specified in Sect. 2.2. This term depends only on the equilibrium plasma quantities and therefore it does not vary in time. The symbol kB denotes the Boltzmann constant, γ = 5/3 is the adiabatic index, m is the particle mass that is specified by a mean molecular weight of 0.6, and g = (0,−g,0) is the gravitational acceleration. The value of g is equal to 274 m s-2.

thumbnail Fig. 1

Hydrostatic solar atmospheric temperature vs. height y.

Open with DEXTER

thumbnail Fig. 2

Vertical profile of plasma β (left) and sound speed cs (right) at the plasma equilibrium.

Open with DEXTER

2.2. Equilibrium solar atmosphere

In a static solar atmosphere all plasma quantities are time-invariant, which means that fe/∂t = 0, where fe denotes a plasma quantity and the subscript e corresponds to the equilibrium. Then, from Eqs. (1)–(4) it follows that for a still (Ve = 0) medium the Lorentz force must be balanced by the gravity force and the gas pressure gradient, (5)and the heating term must compensate for the radiative losses and thermal conduction, (6)This model of the solar atmosphere corresponds to a quiet Sun.

thumbnail Fig. 3

Magnetic field lines at the plasma equilibrium.

Open with DEXTER

2.2.1. Force-free magnetic field of the hydrostatic atmosphere

A hydrostatic atmosphere corresponds to the force-free ((∇ × Be) × Be = 0) magnetic field. We additionally assume a current-free (∇ × Be = 0) magnetic field whose radial Ber, azimuthal Beθ, and vertical Bey components are given as (7)where a and S are free parameters corresponding to the vertical location of the singularity in the magnetic field and the magnetic field strength, respectively. We set a = −1 Mm and S in such way that at the reference point (r = 0,y = yr = 6) Mm the magnitude of magnetic field Be = 9.5 Gauss. The corresponding magnetic field lines are displayed in Fig. 3. We note that the magnetic lines diverge with height and Be is vertical along the symmetry axis, r = 0 Mm.

For a force-free magnetic field it follows from Eq. (5) that the gas pressure gradient has to be balanced by the gravity force, (8)The subscript h corresponds to a hydrostatic quantity. With the use of the ideal gas law and the vertical y-component of Eq. (8), we express the hydrostatic gas pressure and mass density as (9)where (10)is the pressure scale height, and p0 denotes the gas pressure at the reference level, y = yr.

For simplicity reasons we assume that Th varies with height y only, and it specifies a hydrostatic atmosphere that is determined by the semi-empirical model of Avrett & Loeser (2008) that is extrapolated into the solar corona (Fig. 1). In this model, the temperature attains a value of about 7 × 103 K at the top of the chromosphere, y ≈ 2.0 Mm. At the transition region, which is located at y ≃ 2.1 Mm, Th exhibits an abrupt jump (Fig. 1) and grows to about 1.0 × 106 K in the solar corona at y = 10 Mm. Higher up in the solar corona, the temperature increases very slowly, tending to its asymptotic value of about 2 MK at y = 40 Mm. The temperature profile uniquely determines the equilibrium mass density and gas pressure profiles, which decrease with height (not shown).

We specify the plasma β as the ratio of gas to magnetic pressures, (11)The vertical profile of plasma β is illustrated in Fig. 2 (left panel). We note that for the coronal plasma, the value of plasma β is lower than 1 within the displayed region. The vertical profile of the sound speed (12)is displayed in the right panel of Fig. 2. Below the transition region, cs ≈ 10 km s-1. Higher up cs grows, first fast right above the transition region, and higher up slowly, reaching a value of about 100 km s-1 at y ≈ 20 Mm.

thumbnail Fig. 4

Temporal evolution of log (ϱ(r,y)) at t = 20 s, t = 50 s, t = 110 s, t = 175 s, t = 220 s, t = 300 s (from top left to bottom right), for adiabatic MHD equations and for the case of AV = 40 km s-1 and y0 = 1.75 Mm. Arrows represent velocity vectors in the ry plane, [Vr,Vy].

Open with DEXTER

2.2.2. Perturbation

Initially, at t = 0 s we perturb the model equilibrium by the initial pulse in the y-component of velocity, which is expressed as follows: (13)where y0 is the vertical position of the initial pulse, w is its width, and AV its amplitude. We set and hold fixed w = 0.25 Mm, while allowing other parameters to vary. For our studies, the initial position of y0 varies between 1.5 Mm and 1.75 Mm, and the amplitude AV varies between 30 km s-1 and 50 km s-1. The detailed studies were performed for the case of AV = 40 km s-1 and y0 = 1.75 Mm. The value of AV may be associated with reconnection of magnetic field lines.

3. Numerical simulations of MHD equations

To solve Eqs. (1)–(3) numerically, we used the PLUTO code (Mignone et al. 2007, 2012). In our problem, we set the Courant-Friedrichs-Levy number equal to 0.3 and chose piecewise TVD linear interpolation in a second-order Runge-Kutta method, which leads to second-order accuracy in space and time. Additionally, we adopted the Harten-Lax-van Leer discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano 2005).

Our simulation box in (r,y) was set as (0.0,5.12) Mm ×(1.0,40.0) Mm, where y = 0 denotes the bottom of the photosphere. For our study we used the uniform grid within the region (0.0 ≤ r ≤ 5.12) Mm ×(1.0 ≤ y ≤ 11.24) Mm, which is covered by 1024 × 2048 grid points. This leads to a resolution of 5 km in the lower regions of the simulation box. Above this region, namely within the box (0.0 ≤ r ≤ 5.12) Mm ×(11.24 ≤ y ≤ 40.0) Mm, we implemented a stretched grid along the y-direction; this box was divided into 648 cells whose size grows with y. Such a stretched grid plays the role of a sponge as it absorbs incoming signal and allows us to avoid significant reflections from upper boundary. We imposed open boundary conditions at r = 5 Mm, but at the bottom and top, we fixed all plasma quantities to their equilibrium values. The left boundary (r = 0) was set as axisymmetric.

The heating source term H in Eq. (3) was implemented as follows. The code computed the residuum R0 using the initial condition and then it subtracted R0 also at later times during the update over time step Δt as (14)where Un is a plasma vector state at time t = nΔt, n = 1,2,... This automatically mimics the heating term without any need to explicitly write it.

thumbnail Fig. 5

Temporal evolution of ϱ(r = 0,y) (top) and Vy(r = 0,y) (bottom) for the case of adiabatic (left) and non-adiabatic (right) MHD equations and for AV = 40 km s-1 and y0 = 1.75 Mm.

Open with DEXTER

To fully understand the influence of non-adiabatic effects on simulated spicules, we performed first simulations of the adiabatic case, that is, without the thermal conduction and cooling terms. Figure 4 shows the spatial profiles of logarithm of ϱ(r,y) at six instants of time. We note that the system is axisymmetric along r = 0 Mm. The initial pulse in Vy is launched from r = 0 Mm, y = 1.75 Mm, which is located about 0.35 Mm below the transition region. The amplitude of the initial pulse is AV = 40 km s-1. The shock front that results from the initially launched pulse arrives at the transition region and triggers the plasma jet (top panels), which reaches its maximum height of 4.3 Mm (top right panel). At later moments in time, the chromospheric plasma injected into the corona begins to fall toward the transition region (bottom left panel). The falling plasma then triggers the second pulse (bottom middle panel), which results in the second spicule (bottom right panel).

thumbnail Fig. 6

Temporal evolution of ϱ(r = 0,y = 5) Mm for the case of AV = 40 km s-1 and y0 = 1.75 Mm.

Open with DEXTER

Figure 5 shows the temporal evolution of ϱ(r = 0,y) (top panels) and Vy(r = 0,y) (bottom panels) for adiabatic (left panels) and non-adiabatic (right panels) cases. The rise time of the spicule to its maximum height in the adiabatic case is 110 s, which is shorter by about 15 s than in the non-adiabatic case (125 s), and the chromospheric plasma reaches a height of 4.3 Mm. This value is lower than in non-adiabatic case (4.7 Mm). In the adiabatic case we spot three waves (Fig. 5, top left panel): the leading wave is a shock wave, which travels with velocity 170 km s-1; this shock is followed by the contact wave, and the rarefaction wave. We note that all plasma quantities are continuous across a contact wave; the only exception is the mass density, which is discontinuous. As a result, all three waves are seen on the mass density profiles, while Vy(r = 0,y) profiles reveals only two of these waves. As in the non-adiabatic case, in later moments, the system becomes thermally unstable. The right panels show the initial phase of the system evolution. The cross-section of Fig. 5 for the adiabatic case for a fixed value of y = 3.5 Mm is presented in Fig. 6. The corresponding three waves (shock, contact, and rarefaction) are represented by steep gradients at t ≈ 30 s, t ≈ 55 s, and t ≈ 75 s. As the formation mechanism of the spicule is a shock wave, we assume that it is the crest-shock-type jet (Shibata 1982). The mass density at the top of the spicule is ~102 higher than in the ambient coronal plasma (Fig. 6, t ≈ 75 s), and it matches the predicted value for the crest-shock-type jet.

thumbnail Fig. 7

Radial profile of ϱ(r,y = 3.1 Mm) at t = 35 s (solid line), t = 80 s (dashed line), and t = 155 s (dotted line) in the case of adiabatic MHD equations. The case of AV = 40 km s-1 and y0 = 1.75 Mm.

Open with DEXTER

thumbnail Fig. 8

Maximum height of spicule ymax vs. amplitude of the initial pulse AV for y0 = 1.75 Mm. The case of adiabatic (+) and non-adiabatic (+?) MHD equations.

Open with DEXTER

Figure 7 shows the effect caused by the rarefaction wave. The horizontal plane is located at y = 3.5 Mm. The plasma inside the spicule-like structure becomes rarefied. The time at which this rarefaction occurs is clearly shown: the dense plasma at t = 35 s (solid line) experiences rarefaction within the central part of the spicule at t = 80 s (dashed line), which keeps its trend up to t = 155 s (dotted line). At this time, plasma is so rarefied that it may be missed by the detectors. This may mimic the disappearance of II type spicules.

Figure 8 illustrates spicule height ymax vs. amplitude of the initial pulse AV with a fixed position of the initial pulse y0 = 1.75 Mm. Hence we infer that a larger amplitude pulse results in chromospheric plasma reaching a higher maximum altitude, which is intuitively expected as a larger amplitude corresponds to more energy being launched initially. This more energetic pulse leads to higher jets. Figure 9 shows a similar trend in the arise time of the spicule for a growing amplitude of the pulse. The comparison of results for adiabatic (+) and non-adibatic (*) MHD equations reveals that the non-adiabatic effects have a minor influence on both spicule height and rise time of the spicule. The non-adiabatic terms result in a growth of the maximum height by about ~10–15%, while the rise time of the spicule grows by about ~15–20%.

Figure 10 illustrates the influence of the vertical position of the initial pulse on the maximum height of the spicule. For a pulse launched closer to the transition region, the mass density of the upward-pushed plasma is lower. Thus such a pulse exhibits a lower momentum, and we can expect it to penetrate lower coronal regions. We also expect that the rise time of the spicule becomes shorter with a higher value of y0. Indeed, Fig. 11 illustrates this declining trend. The obtained values of spicule height lie within the range of the observed height of the spicules.

thumbnail Fig. 9

Rise time of spicule tr vs. amplitude of the initial pulse AV for y0 = 1.75 Mm. The case of adiabatic (+) and non-adiabatic (+?) MHD equations.

Open with DEXTER

thumbnail Fig. 10

Maximum height of spicule ymax vs. vertical position of the initial pulse y0 for AV = 40 km s-1 and in the case of adiabatic MHD equations.

Open with DEXTER

thumbnail Fig. 11

Rise time of spicule tr vs. vertical position of the initial pulse y0 for AV = 40 km s-1 and in the case of adiabatic MHD equations.

Open with DEXTER

thumbnail Fig. 12

Spicule velocity Vy vs. time for y0 = 1.75 Mm, AV = 40 km s-1 and adiabatic MHD equations.

Open with DEXTER

Figure 12 displays the dependence of spicule velocity on time. The location of the top of the spicule is determined by the vertical position of the contact wave. The upward movement of the transition region plasma starts at t ≈ 10 s, and at this time the spicule reaches its speed of about 45 km s-1. At later moments in time, this spicule decelerates until at t ≈ 100 s the spicule speed attains zero, which occurs at the highest location of the spicule. At the lowest location of the spicule, which is reached at t = 200 s, the vertical velocity of the decreasing plasma is about 30 km s-1. After the first rise and fall, the transition region subsequently experiences oscillations.

thumbnail Fig. 13

Time-distance plot of ϱ(r = 0,y,t) vs. time for y0 = 1.75 Mm, AV = 30 km s-1 and adiabatic MHD equations.

Open with DEXTER

thumbnail Fig. 14

Acoustic cut-off wave period Pac vs. altitude y.

Open with DEXTER

The transition region exhibits oscillations that are due to the rebound shocks (Hollweg 1982). It is interesting that these oscillations occur with a wave period of about 250 s (Fig. 13). This wave period is close to the acoustic cut-off period, (15)which for y = y0 = 1.75 Mm attains a value of about 259 s (Fig. 14).

4. Summary and conclusions

We performed numerical simulations of a spicule by launching initially (at t = 0 s) a localized pulse in the y-component of a plasma velocity in the upper chromosphere. The initial magnetic field configuration was force-free and the initial stratification remained in static equilibrium. We simulated both ideal and non-ideal MHD equations to reveal non-adiabatic effects such as thermal conduction and radiation. Numerical simulations showed that an upward-propagating signal quickly steepened into a shock that propagated into the corona along the magnetic field lines. This shock was followed by the cold and dense chromospheric plasma jet, which exhibited properties of a contact wave and reached a certain height (typically 4–5 Mm) and then returned to the chromosphere. The mean up-flow speed was 20–25 km s-1. The obtained values match those given for the spicules observed by Beckers (1968, 1972). However, subsequent shocks come from the chromosphere with a periodicity of almost the chromospheric acoustic cut-off period, which is the result of the nonlinear wake behind the pulse that propagates in the stratified atmosphere (Kuridze et al. 2009). These shocks again lift up the chromospheric plasma into the corona and cause the quasi-periodic appearance of plasma jets. This is consistent with the rebound shock model suggested by Hollweg (1982). Numerical solutions show that the rarefaction wave that follows the shock-wave results in a reduction of the mass density inside the structure. Both ideal and non-ideal simulations give similar results for the maximum height, upward speed, rise time, and periodicity. Therefore, non-adiabatic effects do not significantly affect the dynamics of jets, but slightly increase the maximum heights and rise time.

In conclusion, our numerical simulations of the solar spicules approximately mimicked the observed properties of the spicules. Since we did not attempt to synthesize observables, an exact comparison with observations cannot be achieved here.

Acknowledgments

This work (B.K. and K.M.) was financially supported by the project from the Polish National Foundation (NCN) Grant No. 2014/15/B/ST9/00106. The work of T.V.Z. was supported by the Austrian “Fonds zur Förderung der Wissenschaftlichen Forschung” (FWF) under projects P26181-N27 and P25640-N27, by Austrian Scientific-Technology Collaboration (WTZ) Grant PL15/2013 and by FP7-PEOPLE-2010-IRSES-269299 project- SOLSPANET. The PLUTO code used in this work was developed by Andrea Mignone. Numerical simulations were performed on the Lunar/Solaris cluster at Institute of Mathematics of M. Curie-Skłodowska University. The visualizations of the simulation variables have been carried out using the IDL (Interactive Data Language) and VisIt software packages.

References

All Figures

thumbnail Fig. 1

Hydrostatic solar atmospheric temperature vs. height y.

Open with DEXTER
In the text
thumbnail Fig. 2

Vertical profile of plasma β (left) and sound speed cs (right) at the plasma equilibrium.

Open with DEXTER
In the text
thumbnail Fig. 3

Magnetic field lines at the plasma equilibrium.

Open with DEXTER
In the text
thumbnail Fig. 4

Temporal evolution of log (ϱ(r,y)) at t = 20 s, t = 50 s, t = 110 s, t = 175 s, t = 220 s, t = 300 s (from top left to bottom right), for adiabatic MHD equations and for the case of AV = 40 km s-1 and y0 = 1.75 Mm. Arrows represent velocity vectors in the ry plane, [Vr,Vy].

Open with DEXTER
In the text
thumbnail Fig. 5

Temporal evolution of ϱ(r = 0,y) (top) and Vy(r = 0,y) (bottom) for the case of adiabatic (left) and non-adiabatic (right) MHD equations and for AV = 40 km s-1 and y0 = 1.75 Mm.

Open with DEXTER
In the text
thumbnail Fig. 6

Temporal evolution of ϱ(r = 0,y = 5) Mm for the case of AV = 40 km s-1 and y0 = 1.75 Mm.

Open with DEXTER
In the text
thumbnail Fig. 7

Radial profile of ϱ(r,y = 3.1 Mm) at t = 35 s (solid line), t = 80 s (dashed line), and t = 155 s (dotted line) in the case of adiabatic MHD equations. The case of AV = 40 km s-1 and y0 = 1.75 Mm.

Open with DEXTER
In the text
thumbnail Fig. 8

Maximum height of spicule ymax vs. amplitude of the initial pulse AV for y0 = 1.75 Mm. The case of adiabatic (+) and non-adiabatic (+?) MHD equations.

Open with DEXTER
In the text
thumbnail Fig. 9

Rise time of spicule tr vs. amplitude of the initial pulse AV for y0 = 1.75 Mm. The case of adiabatic (+) and non-adiabatic (+?) MHD equations.

Open with DEXTER
In the text
thumbnail Fig. 10

Maximum height of spicule ymax vs. vertical position of the initial pulse y0 for AV = 40 km s-1 and in the case of adiabatic MHD equations.

Open with DEXTER
In the text
thumbnail Fig. 11

Rise time of spicule tr vs. vertical position of the initial pulse y0 for AV = 40 km s-1 and in the case of adiabatic MHD equations.

Open with DEXTER
In the text
thumbnail Fig. 12

Spicule velocity Vy vs. time for y0 = 1.75 Mm, AV = 40 km s-1 and adiabatic MHD equations.

Open with DEXTER
In the text
thumbnail Fig. 13

Time-distance plot of ϱ(r = 0,y,t) vs. time for y0 = 1.75 Mm, AV = 30 km s-1 and adiabatic MHD equations.

Open with DEXTER
In the text
thumbnail Fig. 14

Acoustic cut-off wave period Pac vs. altitude y.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.