ACCELERATION OF TYPE 2 SPICULES IN THE SOLAR CHROMOSPHERE. II. VISCOUS BRAKING AND UPPER BOUNDS ON CORONAL ENERGY INPUT

In a review of small scale chromospheric structures, Tsiropoula et al. (2012) point out that the existence of spicules in the chromosphere, eventually identified as mainly vertical supersonic flows, has been inferred from observations for more than 130 yr. As spatial and temporal resolution has increased, the observed maximum flow speeds of spicules have increased, and the observed minimum length scales over which these flows are localized orthogonal to their velocity have decreased. This trend is similar to that of observing, and using semi-empirical models to infer the existence of, increasingly stronger magnetic field strengths and total magnetic energy on increasingly smaller horizontal scales in the photosphere (de Wijn et al. 2009; Sánchez Almeida & Martínez González 2011; Stenflo 2012, 2013). These trends can be connected in that Lorentz forces on a given scale can drive flow on that scale. There is no sign that the smallest scales of magnetic fields or supersonic flows in any region of the solar atmosphere are resolved.

These trends suggest that type 2 spicules, observed on the limb and on the disk (as rapid blueshifted excursions—RBEs) (De Pontieu et al. 2007a, 2007b, 2007c; Langangen et al. 2008; Rouppe van der Voort et al. 2009; Sekse et al. 2012; De Pontieu et al. 2012; de Wijn 2012) are the sub-set of chromospheric flows with the highest speeds (∼ 50–150 km  s⁻¹), and shortest durations (∼ 10–150 s) that can be detected with the highest available spatial and temporal resolution of ∼100 km and 5–8 s (van Noort & Rouppe van der Voort 2006; De Pontieu et al. 2007c; Tavabi et al. 2011; Pereira et al. 2012).

This paper uses a magnetohydrodynamic (MHD) model that improves upon the model in Goodman (2012, henceforth G12) to determine conditions under which flows can be accelerated to vertical speeds comparable to type 2 spicule speeds in the chromosphere. The model in G12 assumes a pre-spicule state without flow, with a constant vertical magnetic field, and neglects the effects of viscosity. The model presented here includes a pre-spicule state with 2D velocity and magnetic fields, and the viscous force due to H i–H i collisions.

Assume cylindrical coordinates (R, θ, z). Let t be time. All quantities are independent of θ, and z is height above some to be specified point in the chromosphere. Assume a weakly ionized H plasma with constant temperature T. Then $p=\rho k_B T/m_p= \rho V_s^2$, where p, ρ, V_s, m_p, and k_B are the pressure, mass density, sound speed, proton mass, and Boltzmann's constant. Assume p = p₁(R, t)exp (− z/L) and ρ = ρ₁(R, t)exp (− z/L). Here L = k_BT/m_pg is the pressure scale height, where g = 2.74 × 10⁴ cm s⁻². Then the vertical pressure gradient and gravitational forces cancel.

Let V be the center of mass (bulk flow) velocity. Assume the magnetic field B = b(R, t)exp (− z/2L). Include the viscous force due to H i–H i collisions in the momentum equation. For the momentum and mass conservation equations an exact solution for the height dependence of V is that it is independent of height. This form of the solution is assumed here, so V = V(R, t). The assumed form of the height dependence of B, V, and ρ implies the flow is accelerated by height independent Lorentz and viscous forces per unit mass acting on a vertical column of gas. The current density J = j(R, t)exp (− z/2L) = c∇ × B/4π. It is assumed that V_θ = 0.

Under the preceding assumptions, the momentum and mass conservation equations are

$$\begin{equation} \dot{V}_R + V_R V_R^\prime + V_s^2 \frac{\rho _1^\prime }{\rho _1}= \frac{j_\theta b_z-j_z b_\theta }{c \rho _1} + \frac{4 \nu }{3 \rho _1} \left(V_R^{\prime \prime } + \frac{V_R^\prime }{R} - \frac{V_R}{R^2}\right). \end{equation} \tag{ 1 }$$

$$\begin{equation} j_z b_R - j_R b_z=0. \end{equation} \tag{ 2 }$$

$$\begin{equation} \dot{V}_z + V_R V_z^\prime = \frac{j_R b_\theta -j_\theta b_R}{c \rho _1} + \frac{\nu }{\rho _1 R} (R V_z^\prime)^\prime. \end{equation} \tag{ 3 }$$

$$\begin{equation} \dot{\rho }_1 + \rho _1 \frac{(R V_R)^\prime }{R} +V_R \rho _1^\prime = \frac{\rho _1 V_z}{L}. \end{equation} \tag{ 4 }$$

Here the dot and prime denote ∂/∂t and ∂/∂R, and ν = 5(m_pk_BT)^1/2/(16π^1/2d²) is the viscosity of HI, where d is an estimate of the atomic diameter (Mihalas & Mihalas 1984). Here d = 2a₀, where a₀ is the Bohr radius.

The |V_R| is now restricted to being sufficiently small to simplify solving the model. This restriction, together with the restriction made above that V_θ = 0, is consistent with the observation that the primary acceleration is in the vertical direction. These restrictions prevent the model from reproducing horizontal and torsional flows, with speeds up to ∼30 km s⁻¹, exhibited by some type 2 spicules (Sekse et al. 2012; De Pontieu et al. 2012).

Constraints are imposed on V_R to justify omitting the terms involving V_R in the momentum Equations (1) and (3). Inspection of those equations indicates the constraints |V_R| ≪ |V_z|, and |V_R| < V_s should be adequate for the solutions presented in Section 6. As shown in Section 3, the constraint on V_R places an upper bound on R₀, which is the characteristic radius within which J is confined. Imposing these constraints reduces Equations (1) and (3) to

$$\begin{equation} V_s^2 \rho _1^\prime = \frac{j_\theta b_z-j_z b_\theta }{c}, \end{equation} \tag{ 5 }$$

$$\begin{equation} \dot{V}_z = \frac{j_R b_\theta -j_\theta b_R}{c \rho _1} + \frac{\nu }{\rho _1 R} (R V_z^\prime)^\prime. \end{equation} \tag{ 6 }$$

In order that the solution for V_z be finite at R = 0 it is necessary that the boundary condition $V_z^\prime (0,t)=0$ be imposed. One more boundary condition is necessary to uniquely determine the solution for V_z. This is chosen as V_z(R = 100 km, t) = 0 since the observed diameter of type 2 spicules is ≲ 200 km. For the solutions in Section 6, the plasma acceleration, and resistive, viscous, and compressive heating occur almost entirely within a radius ≪100 km from the origin.

2.1. Accelerating and Compressive Lorentz Forces

If the Lorentz force plays an important role in accelerating spicules then it must be localized on the space- and timescales observed to characterize spicules. The vertical Lorentz force is (J_RB_θ − J_θB_R)/c. Here a functional form for j_z is assumed, given by Equation (7) below, that is localized on time and radial scales t₀ and R₀ to model the temporal and radial localization of type 2 spicules. As shown below, the assumed form for j_z together with the equations ∇ · J = 0, and ∇ × B = 4πJ/c determine j_R and b_θ, and hence the contribution J_RB_θ/c to the vertical Lorentz force. For the two spicule solutions presented in Section 6, this is the accelerating Lorentz force in that it drives all of the upward acceleration of the plasma. It does this against a much weaker downward component, −J_θB_R/c, of the vertical Lorentz force present in the background (BG) state. The values of t₀ and R₀ are chosen so the solutions exhibit a plasma acceleration time, and radial localization of the accelerated plasma consistent with observations. The observations of type 2 spicules referenced in Section 1 show they have durations ∼10–150 s, and characteristic diameters ≲ 200 km. In Section 6, the choice t₀ = 33.3 s yields an acceleration time of ∼100 s to maximum vertical speeds of ∼66 and 77 km s⁻¹, and the corresponding choices of R₀ = 5 and 1.25 km yield flows mainly confined within diameters ∼80 and 50 km.

j_z is assumed to have the following form.

$$\begin{equation} j_z = j_0 \exp (-\bar{t})(1-\exp (-\bar{t})) \exp (-\bar{R}^2). \end{equation} \tag{ 7 }$$

Here $\bar{R} = R/R_0$ and $\bar{t} = t/t_0$.

The specification of J_z, T, and V_θ causes the complete MHD model to be over-determined. Consequently, three equations must be omitted from the complete MHD model. The remaining equations determine the solution self-consistently. Here the mass and momentum equations, Faraday's law, the R component of a multi-fluid Ohm's law, the ideal gas equation of state, and an NLTE Saha equation are used to determine the solution to the model. The energy equation, and the θ and z components of the Ohm's are omitted from the model, except that the complete Ohm's law is used to estimate the Joule heating rate.

It follows from ∇ · J = 0 that

$$\begin{equation} j_r = \frac{j_0 R_0}{4 L} \exp (-\bar{t})(1-\exp (-\bar{t})) \frac{(1-\exp (-\bar{R}^2))}{\bar{R}}. \end{equation} \tag{ 8 }$$

From ∇ × B = 4πJ/c it follows that

$$\begin{equation} b_\theta =b_{\theta 0} \exp (-\bar{t})(1-\exp (-\bar{t})) \frac{(1-\exp (-\bar{R}^2))}{\bar{R}}, \end{equation} \tag{ 9 }$$

where b_θ0 ≡ 2πj₀R₀/c.

With respect to $\bar{t}$, the maximum values of j_z, j_r, and b_θ are reached at $\bar{t}=\ln (2) \sim 0.693147$, for which $\exp (-\bar{t})(1-\exp (-\bar{t})) = 1/4$. With respect to $\bar{R}$, the maximum values of j_r and b_θ are reached at $\bar{R}= 1.1209$, for which $(1-\exp (-\bar{R}^2))/\bar{R} = 0.6382$. The maximum value of b_θ, and hence of B_θ is then 0.6382 b_θ0/4 = 0.1596 b_θ0. Then specifying b_θ0 determines the maximum magnitude of B_θ, and determines j₀R₀. Figure 1 shows B_θ/B_{θ, max} over the effective time interval $0 \le \bar{t} \le 4$ of Lorentz force driven acceleration, where B_{θ, max} is the maximum value of B_θ.

Zoom In Zoom Out Reset image size

2.2. Background (BG) State

This is a state in which B_θ = 0. It is determined by B_R, B_z, and J_θ in that once these quantities are known they determine ρ, V_z, and V_R for this state. Equation (2), and ∇ · B = 0 imply

$$\begin{equation} B_R(R,z,t)= \frac{B_0(t) R_0}{4 L \bar{R}} (1-\exp (-\bar{R}^2)) \exp (-z/2L) \end{equation} \tag{ 10 }$$

$$\begin{equation} B_z(R,z,t) = B_0(t) \exp (-\bar{R}^2) \exp (-z/2L), \end{equation} \tag{ 11 }$$

where B₀(t) is an arbitrary function of time, specified in Section 3. For the spicule solutions presented in Section 6, $B_\theta \gtrsim 10 (B_R^2 + B_z^2)^{1/2}$ during the main phase of the acceleration. Such a highly twisted field may be susceptible to the kink instability in the real chromosphere. The specification of the time dependence of j_z and B₀ forces the model solutions to be stable. A model driven only by initial and boundary conditions is required to determine the stability of such field configurations in the chromosphere.

From Equations (10) and (11) it follows that

$$\begin{eqnarray} J_\theta &=& \frac{- B_0(t) c \bar{R}}{2 \pi R_0} \left[\left(\frac{R_0}{4 L \bar{R}}\right)^2 (1-\exp (-\bar{R}^2))\right.\nonumber\\ &&\left. {-} \exp (-\bar{R}^2)\right] \exp (-z/2L). \end{eqnarray} \tag{ 12 }$$

The simplifying assumptions in Section 2 that V_θ = 0 and ∂/∂θ = 0 cause the θ component of the momentum equation to reduce to Equation (2). This causes the BG state to become un-coupled from the spicule in that the BG state influences the spicule acceleration process, but this process does not affect the BG state. In reality, there is bi-directional coupling between the pre-spicule state of the atmosphere, and the subsequent spicule acceleration process. The strength of this coupling is not known. A model that allows for bi-directional coupling, for example one that does not require that V_θ = 0, or that V_R and R₀ be small in the sense defined in Sections 2 and 3, is needed to estimate the importance of this coupling. A potentially important effect omitted due to the un-coupling of the spicule from the BG plasma is the distortion of the BG magnetic field due to partial flux freezing, which increases with temperature and degree of ionization.

The model presented here assumes the existence of a current that accelerates plasma through the Lorentz force. The source of this current is not specified. It is proposed here that the source is the emergence of current carrying magnetic flux through the photosphere in inter-granular lanes on scales ≲ 10² km. These magnetic structures may rise into the chromosphere with unbalanced Lorentz forces that accelerate plasma as part of the process of relaxation of B toward a force free, minimum energy state. Magnetic flux emergence, driven by convection, occurs continuously on inter-granular scales, and as observations have improved, more flux has been observed on increasingly smaller scales (Lites et al. 1996, 2008; Orozco Suárez et al. 2007; Lites 2009). These emerging magnetic structures have or develop a vertical, quasi-cylindrical geometry, similar to the geometry of type 2 spicules. The proposition that magnetic flux emergence drives type 2 spicules is discussed in more detail in Section 7.

Driving by magnetic reconnection (e.g., Rouppe van der Voort et al. 2009) is unlikely since type 2 spicules are often observed in regions that appear to be unipolar (McIntosh & De Pontieu 2009; Martínez-Sykora et al. 2011), and since in reconnection an upward jet is expected to be accompanied by a downward jet, whereas there do not appear to be observations of type 2 spicules as multi-jet phenomena. The model presented here includes a radial Lorentz force that compresses the plasma, increasing its density by factors up to ∼10², but the model excludes the possibility of this compression causing a vertical pressure gradient that contributes to spicule acceleration. It is possible that such a pressure gradient is generated in a more realistic model. MHD simulations support this possibility, showing the development of type 2 spicule like jets in regions with intense electric currents associated with magnetic flux emergence on granulation spatial scales (Martínez-Sykora et al. 2009, 2011).

Solving Equation (5) for ρ₁ gives

$$\begin{eqnarray} \rho _1 &=& \rho _{1 \infty } \left[1+ A(\bar{R})((b_{\theta 0}\exp (-\bar{t})(1-\exp (-\bar{t})))^2\right. \nonumber\\ &&\left. + (\alpha B_0)^2) - \frac{B_0(t)^2}{8 \pi V_s^2 \rho _{1 \infty }} \exp (-2 \bar{R}^2) \right]. \end{eqnarray} \tag{ 13 }$$

Here α = R₀/4L,

$$\begin{equation} A(\bar{R})=\frac{I_0(\bar{R})}{4 \pi V_s^2 \rho _{1 \infty }}, \end{equation} \tag{ 14 }$$

$$\begin{equation} I_0(\bar{R})= \int _{\bar{R}^2}^\infty \frac{(1-\exp (-y))\exp (-y)}{y} \; dy, \end{equation} \tag{ 15 }$$

and ρ_1∞ ≡ ρ₁(R = ∞) is the time independent BG density far from the acceleration region.

It is assumed the BG state evolves on the characteristic granulation turnover timescale t_bg = 600 s. Then this state evolves slowly compared with the spicule acceleration timescale of ∼100 s for the solutions in Section 6. Choose $B_0(t)= b_0 \exp (-t/t_{bg}) = b_0 \exp (-\bar{t} t_0/t_{bg})$, where b₀ is constant. The choice of b₀ is constrained by the requirement that ρ > 0, and is discussed in Section 4.

With ρ₁ known, V_z is determined by numerically integrating Equation (6). This is done using the Matlab function pdefun that is adaptive in time, and uses a fixed, specified spatial grid. Here the spatial grid is specified by the range 0 ⩽ R ⩽ 100 km, and grid spacing 0.005R₀. From Equation (4),

$$\begin{equation} V_R=\frac{R_0}{\rho _1(\bar{R},t) \bar{R}} \int _0^{\bar{R}} x \left(\frac{\rho _1(x,t) V_z(x,t)}{L} - \dot{\rho }_1(x,t) \right) \; dx. \end{equation} \tag{ 16 }$$

It is found numerically that V_R increases with R₀, roughly linearly for the solutions presented in Section 6.¹ This is where the constraint on R₀ enters the solution. For given values of the other input parameters, R₀ must be chosen sufficiently small so that |V_R| ≪ |V_z|, and |V_R| < V_s.

The electron and H i number densities n_e and n_H are needed to compute terms in the Ohm's law. Here n_e is computed using the NLTE, statistical equilibrium Saha equation derived in Goodman & Judge (2012, Equation (11)), and n_H is estimated as ρ/m_p.

The Ohm's law for the partially ionized plasma is assumed to be (Mitchner & Kruger 1973)

$$\begin{equation} {\bf E} + \frac{ {\bf V} \times {\bf B}}{c} = \eta \left({\bf J} + M_e \frac{{\bf J} \times {\bf B}}{B} + \Gamma {\bf J}_\perp \right). \end{equation} \tag{ 17 }$$

Here η is the Spitzer resistivity, J_⊥ is the component of J⊥B, the magnetization factor Γ = M_eM_p(ρ_n/ρ)², where M_e, M_p, and ρ_n are the electron and proton magnetizations,² and the neutral density. In the chromosphere ρ_n ∼ ρ. Explicitly,

$$\begin{equation} \eta = m_e^{1/2} \left(\frac{4 (2 \pi)^{1/2} e^2 \ln (\lambda)}{3 (k_B T)^{3/2}} + \frac{\sigma \rho (k_B T)^{1/2}}{n_e m_p e^2}\right), \end{equation} \tag{ 18 }$$

$$\begin{equation} \lambda = \frac{3(k_B T)^{3/2}}{2 e^3(\pi n_e)^{1/2}}, \end{equation} \tag{ 19 }$$

$$\begin{equation} \eta \Gamma = \frac{B^2 m_p^{1/2}}{c^2 \sigma \rho n_e (k_B T)^{1/2}}, \; {\rm and} \end{equation} \tag{ 20 }$$

$$\begin{equation} \frac{\eta M_e}{B} = \frac{1}{n_e e c}. \end{equation} \tag{ 21 }$$

Here m_e, e, and σ(= 5 × 10⁻¹⁵ cm²) are the electron mass and charge magnitude, and the charged-neutral particle scattering cross section (Osterbrock 1961). The MHD Joule heating rate Q_J = J · (E + (V × B)/c), where E is the electric field. From the Ohm's law, $Q_J = \eta (J^2 + \Gamma J_\perp ^2)$. The Pedersen current dissipation rate $Q_P \equiv \eta \Gamma J_\perp ^2$, where ηΓ is the magnetization induced resistivity. The second term on the right hand side of the Ohm's law Equation (17) is the Hall electric field. It is not dissipative since it is ⊥J. It is used in the model to compute E_R, which is used in Section 6.1.3 to compute the z component of the Poynting flux.

The compressive heating rate per unit volume Q_comp = −p∇ · V.

The viscous heating rate per unit volume may be written as³

$$\begin{eqnarray} Q_{{\rm vis}} &=& \nu \left({V_z^\prime }^2 + \frac{4}{3} \left(\frac{Q_{{\rm comp}}}{\rho _1 V_s^2}\right)^2 + \frac{4 V_R}{R}\right. \nonumber\\ &&\left. \times \left(\frac{V_R}{R} - \frac{V_z}{L} + \frac{(\dot{\rho _1} + \rho _1^\prime V_R)}{\rho _1} \right) \right). \end{eqnarray} \tag{ 22 }$$

Faraday's law reduces to

$$\begin{equation} E_\theta (R,t) = - \frac{2 L}{c} \dot{B_R}, \end{equation} \tag{ 23 }$$

$$\begin{equation} E_z(R,t) = E_z(0,t) + \int _0^R \left(\frac{\dot{B_\theta }}{c} - \frac{E_R}{2 L} \right) d\alpha, \end{equation} \tag{ 24 }$$

where E_z(0, t) = η(0, t)J_z(0, t), and E_R is assumed to be given by the R component of the Ohm's law Equation (17).

The solutions in this section are for spicules accelerated in the upper chromosphere. The inputs to the model are R₀, t₀, t_bg, b_θ0, b₀, ρ_1∞, L, and T. Two solutions are presented. For both solutions, t₀ = 33.3 s, t_bg = 600 s, and T = 8000 K. Then L = 241 km, V_s = 8.13 km s⁻¹, and ν = 2.1 × 10⁻³ poise.

The choice of b₀ is constrained by the requirement ρ₁(R, t) > 0. Equation (13) shows that ρ₁ is a minimum in the BG state, for which b_θ0 = 0. For the solutions considered here, R₀ = 1.25 km or 5 km. This determines the possible values of α(= R₀/4L). It can then be shown from Equation (13) with b_θ0 = 0 that for these values of α, $d \rho _1/d\bar{R} > 0$ for $0 \le \bar{R}< \bar{R}_\ast$, where $\bar{R}_\ast \sim 2.92\hbox{--}3.82$, and that $d \rho _1/d\bar{R} \le 0$ for $\bar{R} \ge \bar{R}_\ast$, where equality holds only at $\bar{R} = \bar{R}_\ast$. However, for $\bar{R} = \bar{R}_\ast$, ρ₁ is already essentially equal to ρ_1∞. It follows that the minimum value of the BG density occurs at $\bar{R}=0$. Then, noting that the maximum value of B₀ is b₀, a necessary and sufficient condition for ρ₁ to be positive is

$$\begin{equation} b_0 < 2 V_s \left(\frac{2 \pi \rho _{1 \infty }}{1 - 2 \alpha ^2 I_0(0)}\right)^{1/2}, \end{equation} \tag{ 25 }$$

where I₀(0) = 0.6931. Since α² ∼ 10⁻⁵–10⁻³ for the values of R₀ and L used here, the denominator in Equation (25) is essentially unity. Then Equation (25) reduces to b₀ < 2V_s(2πρ_1∞)^1/2 ∼ 1.6662(n_∞/10¹¹ cm⁻³)^1/2 G, where n_∞ = ρ_1∞/m_p. For upper chromospheric densities this constrains b₀ to be no more than a few Gauss.

For all solutions n_∞ = 2.575 × 10¹¹ cm⁻³. The constraint on b₀ is then b₀ < 2.674 G. Choose b₀ = 2.5 G for all solutions. This implies the density at $\bar{R}=\bar{t}=0$ is 3.24 × 10¹⁰ cm⁻³.

It remains to specify b_θ0 and R₀ to determine the solutions. All plots are at the reference height z = 0. The reference height is the height at which n(R, z, t) = n_∞ at t = 0 and R = ∞. This follows from Equation (13) and ρ(R, z, t) = ρ₁(R, t)exp (− z/L).

The total magnetic field strengths for the spicule solutions presented here and in G12 during the main phase of the acceleration are ∼10–25 G. This is within the range inferred from observations of spicules, typically at heights ∼1400–5000 km above the photosphere (López Ariste & Casini 2005; Trujillo Bueno 2005; Ramelli et al. 2006; Centeno et al. 2010).

6.1. Solution 1

For this solution b_θ0 = 78.35 G and R₀ = 5 km. Then B_θ ⩽ 12.5 G.

6.1.1. BG State

This state is unremarkable for the given input parameters. The magnetic field is essentially vertical. V_z and V_R are <0, with maximum magnitudes 2.2 km s⁻¹ and 34 m s⁻¹. There is a density depression localized in the region $\bar{R} \lesssim 0.7$, with the maximum depression at R = 0. The density at R = 0 is ∼3 × 10¹⁰ cm⁻³. It increases to 1.13 × 10¹¹ cm⁻³ at $\bar{t}=4$, and continues to increase toward n_∞ as the BG magnetic field decays.

The total thermal energy E_J generated by Joule heating in the BG state is estimated as the integral of Q_J over the volume V defined by 0 ⩽ z ⩽ ∞, 0 ⩽ R ⩽ 100 km, and 0 ⩽ θ ⩽ 2π, and over the time interval T defined by 0 ⩽ t ⩽ 4t₀. Then E_J = 9.14 × 10¹⁸ erg. Since J⊥B, E_J is entirely due to Pedersen current dissipation.

The MHD Poynting theorem may be written as (8π)⁻¹∂B²/∂t + ∇ · S = −Q_J − R_ke. Here R_ke ≡ V · (J × B)/c is the rate per unit volume at which electromagnetic energy is transformed into bulk flow kinetic energy by the action of the Lorentz force, S is the Poynting flux, and Q_J + R_ke = J · E. The total amount of electromagnetic energy transformed into bulk flow kinetic energy is denoted by W, and is defined as the integral of R_ke over V and T. Then W = −3.84 × 10¹⁷ erg. The negative value of W indicates that bulk flow kinetic energy is converted into electromagnetic energy. The net amount of electromagnetic energy converted into particle energy within V during T is then E_J + W ∼ 8.76 × 10¹⁸ erg.

Similarly, the thermal energy generated by viscous and compressive heating is estimated as the integrals of Q_vis and Q_comp over V and T, and denoted E_vis and E_comp. For the BG state (E_vis, E_comp) = (1.52 × 10¹⁸, 2.01 × 10¹⁹) erg.

6.1.2. Spicule Acceleration

Set b_θ0 = 78.35 G. Then B_θ ⩽ 12.5 G. Again generate the solution for $0 \le \bar{t} \le 4$.

Figure 2 shows the number density n = ρ/m_p, V_z, and V_R. During acceleration the plasma is compressed to a density ∼100 times its pre-acceleration value. The compression is confined within a diameter ∼20 km, and is a maximum at t ∼ 0.6932 t₀ = 23 s. The maximum of V_z = 65.6 km s⁻¹. If the viscous force is neglected, the maximum value of V_z is 124.1 km s⁻¹. The constraint |V_R| ≪ V_z is satisfied almost everywhere, in particular where almost all acceleration occurs and where V_z is a maximum. The maximum radial speed is 4.4 km s⁻¹, and the radial flow changes from an inflow to an outflow during the acceleration. The flow is localized within a diameter ∼80 km.

Zoom In Zoom Out Reset image size

Figure 3 shows Q_J, Q_vis, and Q_comp, and a panel showing their volume integrals together with the magnetic energy, which is the volume integral of B²/8π. The plots of the heating rates suggest that different heating and cooling processes are important at different times during the acceleration process. Mainly for $\bar{t} \lesssim 0.8$, and some intervals in the range $\bar{R} \lesssim 1$, Q_J ≳ Q_vis, but otherwise Q_vis tends to significantly dominate Q_J. The plots of the volume integrated Joule and viscous heating rates show that Joule heating dominates viscous heating for $\bar{t} \lesssim 0.5$, but not by much. Overall, compressive heating and cooling dominate the heating effect of Q_J and Q_vis. Q_comp is relatively large and positive for $\bar{t} \lesssim 0.5$, and negative and relatively large in magnitude at later times, suggesting a strong cooling effect due to rarefaction (∇ · V > 0⇒dρ/dt < 0). Characteristic mean chromospheric heating rates, assuming emission over a height range ∼10³ km, are ∼(4 × 10⁶–2 × 10⁷) erg cm⁻² s⁻¹/10³ km = 0.04–0.2 erg cm⁻³ s⁻¹ (Withbroe & Noyes 1977; Vernazza et al. 1981; Anderson & Athay 1989). Then, locally and for tens of seconds, the estimated values of Q_J, Q_vis, and especially the magnitude of Q_comp are comparable to, or far exceed the mean chromospheric heating rate.

Zoom In Zoom Out Reset image size

Determining E_J, W, E_vis, and E_comp as in Section 6.1.1 gives (E_J, W, E_vis, E_comp) = (2.48 × 10²⁰, 4.23 × 10²¹, 2.98 × 10²¹, −8.16 × 10²¹) erg. About 95% of the magnetic energy of E_J + W = 4.48 × 10²¹ erg that is expended in accelerating the plasma, and heating it by resistive dissipation is converted into bulk flow kinetic energy. The value of E_J + W is of the same order of magnitude as the larger values of the magnetic energy (∼ 10²¹ erg) in Figure 3. This is consistent with the fact that the magnetic field provides the energy for the acceleration and heating process. About 98% of E_J is due to Pedersen current dissipation. Viscous heating exceeds Joule heating by a factor ∼12. The net compressive heating rate is relatively large and negative, indicating a net cooling effect during the acceleration process.

The energy estimates presented here are crude because they are not constrained by an energy equation. The question of the degree to which these estimates are realistic can be addressed in a meaningful way only by using models that include an energy equation. That equation couples Q_J, Q_comp, and Q_vis to one another, and to radiative cooling, and diffusive and convective thermal energy flow. Solving an energy equation self-consistently as part of the model is also the only way to obtain a meaningful estimate of the variation of T in space and time.

6.1.3. Energy and Mass Fluxes into the Corona

Ji et al. (2012) report observations of impulsive coronal heating associated with upward mass flows originating in photospheric regions of strong magnetic field in intergranular lanes. The flows from the photosphere into the corona are observed to occur along magnetic loops with diameters ∼100 km, comparable to the reported diameters of type 2 spicules. Some of these flows are observed as RBEs. The observations of Ji et al. are consistent with the suggestion that type 2 spicules are an important source of mass and energy for the corona. Semi-empirical analyses of De Pontieu et al. (2011) and Klimchuk (2012) respectively support and refute the possibility that type 2 spicules make a significant contribution of mass and energy to the corona.

The question of how much energy and mass type 2 spicules contribute to the corona is addressed here by computing the horizontal area and time averaged vertical mass flux, and electromagnetic, convective thermal, and bulk kinetic energy fluxes for the solutions presented here. These fluxes are evaluated at the height above the reference height (z = 0) at which the density is 10¹¹ cm⁻³, defined here as the top of the chromosphere. They are upper limits to the mass and energy fluxes into the corona.

Let R_* = 100 km, and t_* = 4t₀. The average mass flux over the area $\pi R_\ast ^2$, and time interval 0 ⩽ t ⩽ t_* is

$$\begin{equation} F_M = \frac{2}{R_\ast ^2 t_\ast } \int _0^{t_\ast } dt \int _0^{R_\ast } dx \; \rho V_z x. \end{equation} \tag{ 26 }$$

The average convective thermal energy flux F_TE = 3k_BTF_M/(2m_p).

The average bulk kinetic energy flux is

$$\begin{equation} F_{KE} = \frac{1}{R_\ast ^2 t_\ast } \int _0^{t_\ast } dt \int _0^{R_\ast } dx \; \rho V_z^3 x. \end{equation} \tag{ 27 }$$

The average electromagnetic energy flux F_EM is computed using the z component S_z = c(E_RB_θ − E_θB_R)/(4π) of the Poynting flux. Equations (17) and (23) are used to compute E_R and E_θ. Then

$$\begin{equation} F_{EM} = \frac{2}{R_\ast ^2 t_\ast } \int _0^{t_\ast } {dt} \int _0^{R_\ast } dx \; S_z x. \end{equation} \tag{ 28 }$$

Given the estimated total energy flux F_t = F_TE + F_KE + F_EM into the corona due to a single spicule, and assuming the spicules occur with a frequency f such that they balance a coronal energy loss flux F_c over the entire surface area A_s of the Sun, the type 2 spicule occurrence frequency for a given solution is estimated as

$$\begin{equation} f=\left(\frac{F_c A_s 100 \; {\rm s}}{t_\ast F_t \pi R_\ast ^2}\right) \left(100 \; {\rm s}\right)^{-1}. \end{equation} \tag{ 29 }$$

Here f is expressed in number per 100 s because this is a characteristic observed spicule lifetime. Choose F_c = 10⁶ erg cm⁻² s⁻¹.

For the solution in Section 6.1.2, (F_M, F_TE/F_t, F_KE/F_t, F_EM/F_t, F_t) = (1.33 × 10⁻⁷ g  cm⁻² s⁻¹, 14.9%, 70.3%, 14.8%, 8.88 × 10⁵ erg cm⁻² s⁻¹), and f = 1.64 × 10⁸ (100 s)⁻¹. A characteristic solar wind mass flux F_sw at the base of the corona is 10⁻¹¹ g cm⁻² s⁻¹ (Withbroe & Noyes 1977). Then F_M ∼ 10⁴F_sw, F_KE dominates the energy flux, and F_t is comparable to the average coronal energy loss during a time t_*. The value of F_M implies that if all of F_t contributes to balancing coronal energy loss, then essentially all the spicule mass injected into the corona must return to the chromosphere, but this mass must be in the corona long enough to transfer all of its energy to coronal plasma. The value of f is ∼10 times larger than that inferred from the semi-empirical estimate of Judge & Carlsson (2010) that there are ∼2 × 10⁷ type 2 spicules distributed over the surface of the Sun at any given time.

The value F_c = 10⁶ erg cm⁻² s⁻¹ used above is a characteristic value for the coronal energy loss. The value of this loss varies over the surface of the Sun. Withbroe & Noyes (1977) give values of F_c = (3 × 10⁵, 8 × 10⁵, 10⁷) erg cm⁻² s⁻¹ for quiet Sun, coronal hole, and active regions. Using these values in Equation (29) gives values for f ∼ 2.5, 6.6, and 82 times larger than the value inferred from Judge & Carlsson (2010). Given observational and model uncertainties, this suggests type 2 spicules may make a significant contribution to coronal energy input in quiet Sun regions, possibly also in coronal holes, but not in active regions.

6.2. Solution 2

The H i viscous force can have a strong braking effect on type 2 spicule acceleration if it is driven by Lorentz forces associated with current densities localized within diameters ∼5–20 km, corresponding to 4R₀ for the values of R₀ used here. If the actual characteristic radii within which the current density is localized are much larger, then the effect of viscous braking is expected to be significantly smaller. This follows from a dimensionless form of the momentum Equation (3), which shows that the Lorentz force is ∝R₀, while the viscous force is $\propto R_e^{-1}$, where $R_e = \rho _c R_0^2/t_0 \nu$ is a Reynold's number, and ρ_c is a characteristic density. Then if all other input parameters are held constant, the ratio of the Lorentz force to the viscous force is $\propto R_0^3$.

The maximum vertical speeds of ∼66 and 77 km s⁻¹ for the solutions presented here are at the lower end of type 2 spicule speeds. Significantly higher speeds might be possible in models not restricted to using values of R₀ sufficiently small so the effects of V_R in the momentum equation can be neglected.

Joule and compressive heating dominate during the early stage of acceleration. Viscous heating dominates at later times when velocity gradients are largest, but the cooling rate due to rarefaction, corresponding to Q_comp < 0, more than cancels this heating rate. The total thermal energy generated by viscous dissipation during the acceleration process is an order of magnitude larger than that due to Joule dissipation, but the total cooling due to rarefaction can be comparable to or larger in magnitude than this energy. These are crude estimates since the model does not include an energy equation.

If all the upward spicule energy flux in the upper chromosphere contributes to coronal energy input, the energy injected by type 2 spicules into the corona may provide a significant fraction of the coronal energy input in quiet regions, possibly also in coronal holes, but not in active regions.

Magnetic flux emergence through the photosphere in inter-granular lanes is a likely source of non-force free current systems that accelerate type 2 spicules through the associated Lorentz force. This conclusion is based on the following complementary sets of observations, and on the MHD simulations cited in Section 2.2. (1) Type 2 spicules are localized on scales ≲ 10² km perpendicular to their direction, consistent with inter-granular lane spatial scales. (2) Magnetic flux continually emerges through the photosphere in inter-granular lanes, with quiet Sun field strengths up to several hundred Gauss in internetwork, ∼10³ G in network, and with the emergence first appearing as a region of mixed polarity, nearly horizontal flux (Lites et al. 1996, 2008; Orozco Suárez et al. 2007; Lites 2009). The emerging magnetic structures become more vertical and quasi-cylindrical as they rise, similar to the geometry of type 2 spicules. (3) On the orders of magnitude larger scales of emerging sunspots, pores, and active regions, B emerges carrying currents generated below the photosphere (Leka et al. 1996; Lites et al. 1998; Burnette et al. 2004; Sun et al. 2012). Flaring phenomena in these regions include acceleration of plasma up to speeds ∼10–10² times greater than those of type 2 spicules. The field emerges as mixed polarity, nearly horizontal flux ropes with bipolar geometries, and field strengths ∼200–600 G. The topology becomes more complex as the field rises into the chromosphere. Lites et al. (1998) emphasize that sections of these flux ropes can attain kilo-Gauss strength after they become nearly vertical as they flow out of the flux emergence region and rise into the chromosphere. This vertical, quasi-cylindrical geometry is similar to type 2 spicule geometry on smaller scales. By analogy, the observations of flux emergence with plasma acceleration on these relatively larger, better resolved scales support the proposition that emerging non-potential inter-granular flux is a source of current systems that drive type 2 spicules. At present there do not appear to be any other type 2 spicule acceleration mechanisms that are as strongly supported by observations.

The chromosphere is a highly resistive medium with respect to the Pedersen current density J_P (∼J_⊥ in the chromosphere) due to the combination of weak ionization and strong magnetization, which distinguishes the chromosphere from the weakly ionized, weakly magnetized photosphere, and the strongly ionized, strongly magnetized corona (Goodman 2000, 2004). Here "highly resistive" means the Pedersen resistivity η_P(∼ ηΓ) is orders of magnitude larger than η. This raises the following questions. (1) If emerging magnetic flux is considered as a driver of type 2 spicules, are the associated currents strong enough to cause type 2 spicule acceleration? (2) If the answer is yes, is the J_⊥ associated with emerging flux sufficiently dissipated by Pedersen current dissipation below the height at which spicules form so that the Lorentz force cannot be a major component of the force that accelerates type 2 spicules? At present these questions cannot be answered for the following reasons: (1) The formation heights of type 2 spicules are not known, other than that they form somewhere in the chromosphere. (2) An approximate, meaningful answer to these questions can be obtained from simulations of instances of the flux emergence process that include the essential physics, and use spatial and temporal resolutions sufficient to accurately compute J_⊥, and the associated Lorentz force and Q_P. Such simulations do not yet exist. The 3D simulations of flux emergence by Arber et al. (2007) are a step toward developing such a simulation, and suggest that as flux rises into the chromosphere, there is strong dissipation of J_P, consistent with general results in Goodman (2000, 2004), and the expectation that the atmosphere becomes increasingly force free with increasing height.

The author acknowledges support from grants ATM-0650443 and ATM-0848040 from the Solar-Terrestrial Research Program of the National Science Foundation to the West Virginia High Technology Consortium Foundation. This research made use of NASA's Astrophysics Data System (ADS). The author thanks both referees for thorough and constructive reviews.