THE NONLINEAR EVOLUTION OF A TWIST IN A MAGNETIC SHOCKTUBE

Magnetic flux tubes are ubiquitous within the solar atmosphere. They are known as the building blocks that may form small or large scale structures. Flux tubes on different scales often exhibit mass flow. Examples include siphon flows in coronal loops (Orlando et al. 1995a, 1995b), counterstreaming (Lin et al. 2003), and field-aligned flows within filament channels (Lin et al. 2005), spicules (Hollweg et al. 1982; De Pontieu et al. 2004; Zaqarashvili & Erdélyi 2009; Scullion et al. 2011), Evershed flows within sunspots (Montesinos & Thomas 1997; Plaza et al. 1997). In addition, non-thermal broadenings have been observed within the chromosphere (Beckers 1968, 1972; Pishkalo 1994; Jess et al. 2009, 2015). These non-thermal broadenings are usually interpreted as Alfvén waves and have been considered as a possible mechanism responsible for Type-I spicule formation (Hollweg et al. 1982; Sterling & Hollweg 1988; Hollweg 1992; Kudoh & Shibata 1999; James et al. 2003; Erdélyi & James 2004) though it is possible that other mechanisms such as leakage of p-mode oscillations along inclined field lines are responsible (De Pontieu et al. 2004). For a review, see Sterling (2000).

Using Hinode De Pontieu et al. (2007) proposed the existence of Type-II spicules, which are more explosive than the classical, Type-I variety. De Pontieu et al. (2012) used the Swedish 1 m Solar Telescope to establish and quantify three types of motions within Type-II spicules: (i) field aligned flows of 50–100 km s⁻¹ (ii) swaying back and forth with speeds of 15–20 km s⁻¹ and (iii) torsional motions between 25 and 30 km s⁻¹. Type-II spicules are often associated with regions of rapid blueshift and redshift. The launch of IRIS (Interface Region Imaging Spectrograph) has led to observations of small-scale twists that propagate along these structures. The typical amplitudes are 10–30 km s⁻¹ and the duration is less than a minute (De Pontieu et al. 2014).

Recently, high-resolution observations from Crisp Imaging Spectro-Polarimeter (Wedemeyer-Böhm et al. 2012) detected a number of small-scale chromospheric swirls. It was also shown numerically by (non)convection simulations (Fedun et al. 2011; Wedemeyer-Böhm et al. 2012) that granular-driven magnetic vortices in the photosphere (Brandt et al. 1988; Bonet et al. 2008, 2010) are capable of driving such motions up to the transition region and further into the lower corona, which provides an alternative method for transmitting energy.

Li et al. (2012) witness flow channels within a solar tornado that rise and fall along curved trajectories and suggest that the motions are along helical magnetic field lines. They also claim that flows and density waves may propagate along the magnetic field lines within the solar tornado. In intergranular lanes, similar motions have been seen in the form of magnetic whirlpools of downflowing plasma (Bonet et al. 2008, 2010). However, it is also worth mentioning that Panasenco et al. (2014) show rapid swirling motions and solar tornadoes may be an illusion caused by 2D projection effects at the solar limb.

Work by Vranjes (2014) has shown that partial ionization of the chromosphere can lead to Alfvén waves being severely damped between the photosphere and corona. Studies by James et al. (2003) and Erdélyi & James (2004) have shown that partial ionization of the atmosphere may aid Alfvén waves with the heating and formation of chromospheric structures such as spicules.

Magnetohydrodynamics (MHD) instabilities in the solar atmosphere play an important role in the amplification of perturbations, which may drastically alter the system. Taroyan (2008) discovered the possibility of a new MHD instability which is based on the interaction between an incompressible Alfvénic perturbation and a compressible plasma flow. The perturbations amplify due to over-reflection—a concept first investigated by Acheson (1976) in regards to the Kelvin–Helmholtz instability.

The Alfvén instability was investigated for siphon flows within asymmetric coronal loops and found that linear torsional perturbations driven at the footpoints are exponentially amplified (Taroyan 2009). Taroyan (2011) considered the Alfvén instability for smooth plasma flow within open and expanding flux tubes in the presence of gravity. The analysis was extended by Taroyan (2015) to include non-isothermal flows, a body force term in the momentum equation, and a full energy equation.

The study of shock waves can be dated back to the 17th century with Sir Isaac Newton's Principia in 1687. Poisson's memoirs in 1808 laid the foundations for the first shock wave theory and later work by De Laplace (1816) and subsequently Rankine (1859) improved Newton's theory by assuming that sound is an adiabatic process. In Hugoniot's (1889) memoirs, the equation which relates the internal energy to the kinetic energy of a thermodynamic system is introduced. Fast forward a couple of centuries and shocks are a topical study in astrophysics and the Sun. Orlando et al. (1995b) investigate the physical conditions within coronal loops that allow the formation of isothermal and adiabatic shocks for shocked critical solutions and shocked supersonic solutions when a siphon flow is present. Shocks appear in abundance in the solar system with shock waves in the solar atmosphere (Hollweg 1992; Ryutova & Tarbell 2003; Centeno et al. 2009), CMEs (Magdalenić et al. 2014; Fainshtein & Egorov 2015; Lugaz et al. 2015), interplanetary shocks within the solar wind (Grygorov et al. 2014) and at the solar wind termination shock (Baranov et al. 1996).

High resolution observations have shown that the Evershed flow is mainly concentrated into the dark penumbral filaments where the magnetic field is nearly horizontal (Title et al. 1992). The flow speed can reach up to 20 km s⁻¹ or more in individual filaments. In some filaments the flow is observed to slow abruptly at a location of enhanced emission, which is strongly suggestive of a standing shock in supersonic siphon flow (Maltby 1975; Georgakilas et al. 2003; Borrero et al. 2005). However, none of these studies have focused on the interaction between a standing shock and an Alfvén wave, until now.

The present paper is the first, pure numerical study that confirms the existence of the Alfvén instability. The study extends the investigation of the phenomenon to the nonlinear regime. The interaction of a small Alfvénic twist with a classical gas-dynamic shock; the amplification of the twist and consequences of nonlinearity; and the role the Alfvén wave plays in energy transfer are analyzed.

In the next section the model set-up is presented. This includes the geometry used, choice of boundary conditions, replication of an Alfvén wave, setting up a stationary shock, and normalization. In Section 3, the conditions required for the Alfvén instability and thus the amplification of an Alfvén wave are discussed. The model results are presented in Section 4 and discussed in Sections 5 and 6.

The model presented in this paper is a 1.5-D axisymmetric magnetic shocktube, which was introduced by Hollweg et al. (1982) and subsequently employed by Sterling & Hollweg (1988), Kudoh & Shibata (1999), Matsumoto & Shibata (2010), and others. It is discussed by Hollweg (1981) that a single field-line is modeled which resides close to but not on the axis of symmetry, such that for typical cylindrical coordinates, $r\ne 0$ at any point. The equations used in these models may describe torsional and shear Alfvén waves in the nonlinear regime depending on whether the chosen geometry is cylindrical or Cartesian (Priest 2014, Section 4.3).

A uniform background magnetic field in the longitudinal direction, B_z which does not vary in space or time is defined. The magnetic field is highly inclined such that the field-line is virtually horizontal.

In the forthcoming analysis, θ denotes the azimuthal direction and as it is assumed that none of the quantities vary in the azimuthal direction, it follows that $\frac{\partial }{\partial \theta }=0$. The model is generated using the MHD equations in VAC (Versatile Advection Code; see, e.g., Tóth 1997). The numerical scheme used to resolve the system is CD4 (4th order central differencing) with a TVDLF (Total Variance Diminishing Lax–Friedrich) predictor step and minmod limiter. The shocktube geometry can be seen in Figure 1.

Zoom In Zoom Out Reset image size

2.1. Normalization

2.2. Model Equations

For the 1.5-D problem, the MHD equations implemented in VAC that describe the model are shown below:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+\displaystyle \frac{\partial }{\partial z}({v}_{z}\rho )=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(\rho {v}_{z})+\displaystyle \frac{\partial }{\partial z}({v}_{z}\rho {v}_{z})=-\displaystyle \frac{\partial }{\partial z}({p}_{\mathrm{tot}}),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}(\rho {v}_{\theta })+\displaystyle \frac{\partial }{\partial z}({v}_{z}\rho {v}_{\theta })=\displaystyle \frac{{B}_{z}}{{\mu }_{0}}\displaystyle \frac{\partial }{\partial z}({B}_{\theta }),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial e}{\partial t}+\displaystyle \frac{\partial }{\partial z}({v}_{z}e) & = & \displaystyle \frac{{B}_{z}}{{\mu }_{0}}\displaystyle \frac{\partial }{\partial z}({v}_{z}{B}_{z}+{B}_{\theta }{v}_{\theta })\\ & & -\displaystyle \frac{\partial }{\partial z}({v}_{z}{p}_{\mathrm{tot}}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {B}_{\theta }}{\partial t}+\displaystyle \frac{\partial }{\partial z}({v}_{z}{B}_{\theta })=\displaystyle \frac{{B}_{z}}{{\mu }_{0}}\displaystyle \frac{\partial }{\partial z}({v}_{\theta }),\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{p}_{\mathrm{tot}}=\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2}+p,\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&p=(\gamma -1)\left(e-\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}-\displaystyle \frac{1}{2}{{\boldsymbol{B}}}^{2}\right).\end{eqnarray} \tag{ 7 }$$

For brevity, ${\boldsymbol{B}}={B}_{z}{{\boldsymbol{e}}}_{z}+{B}_{\theta }{{\boldsymbol{e}}}_{\theta }$, and ${\boldsymbol{v}}={v}_{z}{{\boldsymbol{e}}}_{z}+{v}_{\theta }{{\boldsymbol{e}}}_{\theta }$. The plasma density is denoted by ρ, whereas the longitudinal and transverse velocities are presented as v_z and v_θ. The internal energy is denoted by e, while the thermal, and total pressures are p and p_tot. B_z and B_θ are the longitudinal and transverse magnetic field components. In this first study, the gravitational force is not included in the momentum equation, which implies a highly inclined longitudinal magnetic field. Such magnetic field configurations are observed in the photosphere and in the chromosphere. Examples include sunspot penumbrae and plage regions in the chromosphere.

The last terms of Equations (2) and (4) describe the nonlinear coupling between the z and θ variables. Due to the small initial twist the perturbation of the magnetic field and corresponding derivatives remain small. Therefore, the initial system behavior is linear. With amplification of the perturbations, the p_tot terms become more significant and lead to nonlinear behavior by producing momentum perturbations in the z components.

The expression for energy conservation, which is of the same form as that used by Hollweg (1992) is obtained through combining Equations (4), (6), and (7):

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}\left[\displaystyle \frac{p}{\gamma -1}+\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\displaystyle \frac{1}{2}{B}_{\theta }^{2}\right]\\ &&\quad +\displaystyle \frac{\partial }{\partial z}\left({v}_{z}\left[\displaystyle \frac{\gamma p}{\gamma -1}+\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}\right]\right)\\ &&\quad +\displaystyle \frac{\partial }{\partial z}({v}_{z}{B}_{\theta }^{2}-{v}_{\theta }{B}_{\theta }{B}_{z})=0.\end{eqnarray} \tag{ 8 }$$

The terms within the temporal derivative of Equation (8) give the thermal, kinetic, and magnetic energy densities. The two terms within the first spatial derivative give the convection of enthalpy and kinetic energy whereas the two terms in the second spatial derivative give the Poynting flux in the inertial frame.

2.3. Shock Equilibrium

In order to set up the required stationary gas-dynamic shock in the center of the shocktube, a longitudinal background flow, which is supersonic yet sub-Alfvénic, is introduced. The plasma either side of the gas-dynamic shock is described by the Rankine–Hugoniot jump conditions (9)–(12). Subsequently, the background flow becomes subsonic once it passes through the shock interface whereas the density, pressure, and temperature of the plasma increase:

$$\begin{eqnarray}&&{v}_{z2}={v}_{z1}\displaystyle \frac{{\rho }_{1}}{{\rho }_{2}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rho }_{2}={\rho }_{1}\displaystyle \frac{{v}_{z1}}{{v}_{z2}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{p}_{2}={p}_{1}-{\rho }_{2}{v}_{z2}^{2}+{\rho }_{1}{v}_{z1}^{2},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{e}_{2}-{e}_{1}=\displaystyle \frac{1}{2}({v}_{z1}^{2}-{v}_{z2}^{2})+\displaystyle \frac{{p}_{1}{\rho }_{2}-{p}_{2}{\rho }_{1}}{{\rho }_{1}{\rho }_{2}},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{v}_{z1}={M}_{c1}{c}_{S},\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&{c}_{S}=\sqrt{\displaystyle \frac{\gamma {P}_{1}}{{\rho }_{1}}}=\sqrt{\displaystyle \frac{\gamma {RT}}{\xi }}.\end{eqnarray} \tag{ 14 }$$

Subscript 1 denotes the plasma upstream of the shock and subscript 2 denotes the post-shock plasma. M_c1 and c_S are the sonic Mach number of the plasma flow, and sound speed in region 1. R is the molar gas constant, T is the plasma temperature, and ξ is the molar mass. From Priest (2014), Equations (9)–(11) can be expressed in terms of the sonic Mach number of region 1:

$$\begin{eqnarray}&&{v}_{z2}=\displaystyle \frac{(\gamma -1){M}_{c1}^{2}+2}{(\gamma +1){M}_{c1}^{2}}{v}_{z1},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{\rho }_{2}=\displaystyle \frac{(\gamma +1){M}_{c1}^{2}}{(\gamma -1){M}_{c1}^{2}+2}{\rho }_{1},\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{p}_{2}=\displaystyle \frac{2\gamma {M}_{c1}^{2}-(\gamma -1)}{\gamma +1}{p}_{1}.\end{eqnarray} \tag{ 17 }$$

Analytically, the shock location is represented by a sharp discontinuity in the density, pressure, and velocity variables. The computational modeling of such a structure is difficult due to numerical instabilities appearing at the shock region. This difficulty in replicating an analytical shock is true regardless of the number of grid-points incorporated as there will always be an associated error. These errors lead to the Gibbs phenomenon at the shock with a trailing wake of oscillations downstream of the discontinuity. In order to eliminate this and introduce stability throughout the simulation, smoothing is introduced using a smooth cubic function. The smooth profile adjusts itself into a stable, static shock, where the Rankine–Hugoniot jump conditions are completely satisfied either side of this region (Figure 2). The resulting profile is then taken as the initial equilibrium from which the study begins.

Zoom In Zoom Out Reset image size

2.4. Boundary Conditions

There are different approaches as far as photosphere boundary conditions are concerned. In studies of thin flux tubes it is customary to introduce a driven boundary at the photospheric level. For example, a driven photospheric boundary was introduced by Hollweg (1981) to study loop resonances and by Hollweg (1986), Van Ballegooijen et al. (2011) to study Alfvén wave turbulence via nonlinear interactions. On the other hand, in their study of wave propagation and conversion in 30–40 Mm wide sunspot slabs, Khomenko et al. (2009) adopted a model which extends from 10 Mm below the photosphere to the low chromosphere and has a gradual decrease in the density.

Two different types of wave excitation were implemented for the simulation. First, consider a single sine-squared driver, F which is implemented as a source term in Equation (3) with the following form:

$$\begin{eqnarray}\rho F & = & A{\mathrm{sin}}^{2}\left(\pi \displaystyle \frac{t-{t}_{\mathrm{min}}}{{t}_{\mathrm{max}}-{t}_{\mathrm{min}}}\right)\\ & & {\mathrm{sin}}^{2}\left(\pi \displaystyle \frac{z-{z}_{\mathrm{min}}}{{z}_{\mathrm{max}}-{z}_{\mathrm{min}}}\right),\end{eqnarray} \tag{ 18 }$$

and in Equation (4) as ρF · v_θ. Note that, t_min < t < t_max, and z_min < z < z_max. The start time for the simulation is denoted as t_min and is chosen for any time after the equilibrium stabilizes. The source term is active until time, $t={t}_{\mathrm{max}}\ ={t}_{\mathrm{min}}+0.2\ L/{c}_{S1}$. A is the amplitude of the driver and z_min, z_max denote the region in which the driver is active.

It is assumed that the ghost cells for ρ, ρv_z, e, B_z, and B_θ are all set as symmetric (reflects from close-by mesh cells). This corresponds to a zero derivative boundary condition, i.e., a Neumann boundary condition. This means that the reflected quantity is "squashed" by the dense photosphere until the photosphere "pushes" that quantity back into the atmosphere due to Newton's 3rd law of motion.

The boundary type for ρv_θ is set as asymmetric, which is the same as symmetric with the addition that the reflection is multiplied by −1. The multiplication by −1 means that the twist velocity changes sign and direction as the pulse interacts with the boundary. This corresponds to a Dirichlet boundary and allows the twist to reflect back into the direction it has just propagated from rather than allowing the twist to "skip" off the surface and twist in the same direction continually. This choice of driver and boundary conditions generates a smooth Alfvénic pulse, which upon backward propagation reflects at the photospheric boundary.

Second, a driven boundary for ρv_θ was tested to drive the Alfvénic perturbations from the photosphere. The driven boundary is prescribed as:

$$\begin{eqnarray}&&\rho {v}_{\theta }{| }_{z=0}={At}\mathrm{exp}\left(-\displaystyle \frac{t}{T}\right),\end{eqnarray} \tag{ 19 }$$

where the decay time of the wave is T = 1/1000 L/c_S1 and t > t_min. The remaining quantities are set to continuous at the left boundary. Additionally, it is worth mentioning that we also tested a driven boundary for $e{| }_{z=0}={e}_{1}+{e}_{k\theta }$ when driving the $\rho {v}_{\theta }{| }_{z=0}$ boundary. However, no discernible difference could be made between this incorporation and solely driving the $\rho {v}_{\theta }{| }_{z=0}$ boundary.

Convection in general is a complex compressible fluid motion with the energy spectrum corresponding to the turbulence. Cranmer & Van Ballegooijen (2005) took the power spectrum of transverse magnetic bright point motions in the photosphere as the lower boundary condition for their model of Alfvénic fluctuations. Here, we are interested in the basic physics of the twist evolution. In order to have a clear picture of the process, a single, decaying pulse of the form (19) is introduced at the footpoint. As the results are largely similar, the principle of Occam's razor may be applied in that the boundary condition (19) with fewest assumptions is best.

That is, the ghost cells at the base of the magnetic shocktube are set as continuous for all variables, except ρv_θ, which is driven through Equation (19) to excite Alfvén waves at the photospheric level. This is an inhomogeneous Dirichlet-like boundary condition for ρv_θ, which results in the boundary becoming reflective. At the other side of the shocktube, the boundary type for the ghost cells is set to continuous. This allows the propagation of waves and flows to pass through the end of the simulated region without any reflections.

Alfvénic pulses generated in region 1 will become over-reflected and amplified at the shock interface if

$$\begin{eqnarray}&&{v}_{z1}\gt {v}_{z2}+{c}_{A2},\end{eqnarray} \tag{ 20 }$$

where c_A denotes the Alfvén speed in the corresponding regions (Taroyan 2008) and is described by

$$\begin{eqnarray}&&{c}_{A}=\displaystyle \frac{B}{\sqrt{\mu \rho }},\end{eqnarray} \tag{ 21 }$$

where μ is the permeability of free space.

In the case of a sub-Alfvénic flow, v_z1 < c_A1 which also satisfies condition (20), the shock tube is unstable with respect to Alfvénic perturbations as they bounce back and forth between the shock interface and the photospheric driver upstream of the shock. The shock front reflects and transmits Alfvénic perturbations upstream and downstream with amplified amplitudes. It acts as an amplifier leading to the Alfvén instability. The flow velocity jump across the shock interface is expressed through the Rankine–Hugoniot relation (15).

Condition (20) contains equilibrium quantities both upstream and downstream of the shock that are dependent. Using the Rankine–Hugoniot relation for velocity (15) it is possible to rewrite the over-reflection condition (20) in terms of the independent sonic and Alfvénic Mach numbers in the upstream region 1. First, note that the Alfvén speeds on either side of the shock are related through the equation:

$$\begin{eqnarray}&&\displaystyle \frac{{v}_{z1}}{{c}_{A1}^{2}}=\displaystyle \frac{{v}_{z2}}{{c}_{A2}^{2}},\end{eqnarray} \tag{ 22 }$$

which is a consequence of the continuity condition (10). Substituting the expressions (15) and (22) into the over-reflection condition (20) yields

$$\begin{eqnarray}&&{v}_{z1}-\displaystyle \frac{(\gamma -1){M}_{c1}^{2}+2}{(\gamma +1){M}_{c1}^{2}}{v}_{z1}\gt \sqrt{\displaystyle \frac{{v}_{z2}}{{v}_{z1}}}{c}_{A1}.\end{eqnarray} \tag{ 23 }$$

Using the Rankine–Hugoniot relation (15) and the inequality (23), the following Alfvén instability criterion is derived in terms of the sonic Mach number, M_c1, and the Alfvénic Mach number, M_A1 = v_z1/c_A1:

$$\begin{eqnarray}&&\displaystyle \frac{\sqrt{(\gamma +1){M}_{c1}^{2}(\gamma -1){M}_{c1}^{2}+2(\gamma +1){M}_{c1}^{2}}}{2({M}_{c1}^{2}-1)}\\ &&\quad \lt {M}_{A1}\lt 1.\end{eqnarray} \tag{ 24 }$$

This condition (24) shows that the Alfvén instability is not present when M_A1 ≳ 1 because the flow becomes super-Alfvénic. For γ ≳ 1 the instability condition (24) is reduced to

$$\begin{eqnarray}&&\displaystyle \frac{{M}_{c1}}{{M}_{c1}^{2}-1}\lt {M}_{A1}\lt 1.\end{eqnarray} \tag{ 25 }$$

The above derived inequality (25) shows that the range of flow speeds for which the system is unstable becomes broader with increasing sonic Mach numbers.

The implemented driver (19), which mimics the torsional component of the photospheric motions, such as granular buffeting, generates a small-amplitude torsional Alfvén wave. This twist propagates along the shocktube at a constant speed of ${c}_{A1}+{v}_{z1}$ until it interacts with the gas-dynamic shock.

At this point, if the flow is too low to satisfy the instability criterion (20), the pulse is partially reflected and transmitted with damped amplitudes. The reflected wave propagates against the supersonic plasma flow at a constant speed of ${c}_{A1}-{v}_{z1}$ until it reaches the photosphere. The wave bounces back and forth between the dense photosphere and static shock until the pulse vanishes after a few iterations of this process.

However, if the Alfvén instability is present, the Alfvén wave is amplified through over-reflection and transmission at the shock position. The pulse propagates backwards until it returns to the photosphere where it is reflected, and the process repeats (Supplementary Movie 1).

The amplification of the Alfvénic perturbation leads to the formation of secondary waves (Figures 3 and 4). This occurs due to nonlinear coupling of the transverse variables with the longitudinal variables in Equations (2) and (5). The variable p_tot (6) translates twists in the magnetic field to form perturbations in the longitudinal variables. Initially, the Alfvénic perturbation is small, so the effects of this nonlinear coupling between the two directions is negligible—i.e., this corresponds to the linear stage of the amplification process. The secondary waves generated by the initial twist correspond to fast- and slow-magnetoacoustic waves propagating along the magnetic field lines with phase speeds of c_A and c_S in the flow reference frame. This is discussed in detail in the book by Priest (2014).

During the growth of the Alfvén wave via over-reflection, the nonlinear coupling becomes more important, as the magnetoacoustic waves form into propagating shocks. These are described by Hollweg (1992) as an increase in $| {B}_{\theta }|$ across the shock for a fast shock and a decrease in $| {B}_{\theta }|$ across the shock for a slow one. The fast-shock corresponds to the leading edge of the Alfvénic pulse and travels at ${c}_{A}+{v}_{z}$ upon forward propagation, and ${c}_{A}-{v}_{z}$ on backward propagation. The slow-shock travels at the ${c}_{S}+{v}_{z}$ on forward propagation and corresponds to the trailing edge of the Alfvénic pulse.

However, unlike the fast-shock, the reflected slow-shock does not propagate against the plasma flow once it reaches the gas-dynamic shock, nor does it propagate through the shock interface and out of the low-chromosphere. Instead, it remains at the gas-dynamic shock location as it is too slow to propagate against the upcoming flow, yet too fast for the post-shock flow to overcome its backward propagation.

Note that unlike previous studies (Hollweg et al. 1982; Hollweg 1992; and others), the formation of slow and fast shocks does not require the presence of gravity. Instead, the formation of these shocks is caused by the instability mechanism.

Each subsequent slow-shock that is generated by the photospheric reflection of the Alfvénic pulse is trapped in the same manner, with each one merging with the gas-dynamic shock and contributing to the formation of a standing shock in v_θ and B_θ. This evolution of the slow-mode waves with the gas-dynamic shock can be seen in Figures 5 and 6 and Supplementary Movie (1) as initially there is no discernible difference between the pre-shock and post-shock values—and later it is clear that a standing shock is created in the transverse variables. For details on the confluence of shocks, see Whitham (1975; pp. 52–53 and 110–112).

In addition to the generation of shocks, the nonlinear coupling also affects the sound speed (Figure 7), and Alfvén speed (Figure 8) of the shocktube plasma. It is clear that a change in sound speed also means a change in the plasma temperature (see Equation (14)). Therefore the system evolves into a configuration of ever-increasing complexity until the Alfvénic pulse ceases to amplify and twist the shocktube due to the propagation of the gas-dynamic shock. Figure 8 also shows that as the damping time of the driver becomes long, the amplification process and the nonlinear evolution occurs over a shorter period of time. However, qualitatively the response of the system to the twist remains similar.

Zoom In Zoom Out Reset image size

Once the gas-dynamic shock propagates out of the modeled magnetic flux tube, the Alfvén wave disappears, and a new equilibrium is formed. The resulting equilibrium is shown in Figure 8 and Supplementary Movies (1) and (2) after t = 4.5 L/c_S1. The magnetic field has become highly twisted from the entrapment and amplification of the Alfvén wave between the gas-dynamic shock and photosphere. The highly twisted magnetic field leads to the Alfvén speed doubling when the nonlinear coupling is at its strongest, and remains raised thereafter (Figure 8). The plasma flow remains supersonic, though its speed has been reduced by the Alfvén wave amplification through over-reflection at the shock interface. The density and pressure have been raised to ρ = 1.14ρ₁ and p = 1.33p₁. This is a consequence of the nonlinear coupling, which has altered the plasma properties, i.e., heating of the flux tube, increased magnetic field strength, Alfvén and sound speeds increasing and a decrease in flow speed upstream of the shock.

4.1. Gas-dynamic Shock Propagation

4.2. The Conversion of Energy and Generation of a Global Twist

4.3. Energy Flux and Energy Density

Following Hollweg (1992), Equation (8) can be split into longitudinal and transverse components to yield the energy fluxes for z and θ:

$$\begin{eqnarray}&&{F}_{z}=\displaystyle \frac{1}{2}\rho {v}_{z}^{3}+\displaystyle \frac{{v}_{z}\gamma p}{\gamma -1},\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&{F}_{\theta }=\displaystyle \frac{1}{2}\rho {v}_{\theta }^{2}{v}_{z}+\displaystyle \frac{1}{\mu }({B}_{\theta }^{2}{v}_{z}-{B}_{z}{v}_{\theta }{B}_{\theta }),\end{eqnarray} \tag{ 27 }$$

where μ is the permeability of free space.

It is worth noting that the integral of the internal energy may be expressed as

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{\displaystyle \int }_{{z}_{1}}^{{z}_{2}}e\;{dz}=F({z}_{1})-F({z}_{2}).\end{eqnarray} \tag{ 28 }$$

This shows that the difference in the energy flux between two spatial locations is equal to the temporal derivative of the energy content in the region between those two locations. The energy content in a region increases when the right-hand side of the above equation is positive. Conversely, a negative flux difference corresponds to a decrease in the energy content.

Figure 10 shows that during the linear stage of the evolution the torsional flux upstream and downstream of the gas-dynamic shock are approximately equal and there is no change in the energy content. However, as the system becomes nonlinear, the θ-flux downstream of the gas-dynamic shock does not amplify as much as the flux upstream of the gas-dynamic shock (Figure 10). The explanation for this is related to the propagation of the fast shock. That is, the nonlinear coupling term in Equation (2) translates the magnetic perturbations into z-momentum. Thus, some of the θ-flux generated by the Alfvénic perturbations at z = 0.3 L is not seen by the time the perturbations reach z = 0.99 L. The decrease in the θ-flux from left to right is represented by the difference between the black and blue curves in Figure 8. According to Equation (28) this leads to an increase in the corresponding energy content, namely, the magnetic energy and the θ component of the kinetic energy. Indeed, Figure 9 shows an increasing trend for the θ components of the energy density. Similarly, the increase in the z-flux from left to right represented by the difference between green and red curves in Figure 8. According to Equation (9), the result is a decrease in the corresponding energy content represented by the thermal energy and the z-component of the kinetic energy. Again, this is evidenced by Figure 9. The slight increase in the kinetic energy density at z = 0.99 L is explained above. These results indicate that the instability mechanism extracts the z components of energy and largely converts those into θ components.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.4. The Effect of Flow Speed on the System Evolution

The effect the initial flow speed v_z has on the evolution with respect to the Alfvénic instability discussed in Section 3 is also investigated. For an initial Alfvén speed of 3.89 c_S1, the range for the flow-speed that will yield the Alfvénic instability is 2.25 c_S1 ≤ v_z1 ≤ c_A1. All parameters are kept constant, with the same values as used before, but with v_z1 varying from 2.25 c_S1 to 3.75 c_S1 in increments of 0.25 c_S1. Figure 11 shows that as the initial flow-speed approaches the initial Alfvén speed, the time required for the Alfvénic pulse to amplify and "push" the gas-dynamic shock out of the low-chromosphere decreases. The speed of this process more than doubles when v_z1 = 3.75 c_S1 compared to when the initial speed is taken as the lower velocity limit, v_z1 = 2.25 c_S1. It can also be seen in Figure 11 that not only does the time required for this process to occur decrease but the amplitude of the Alfvénic perturbations increase with increasing flow-speed. These results are interesting as one could initially assume that as the flow-speed approaches the Alfvén speed, it may inhibit the backward propagation of the reflected waves/shocks due to ${c}_{A}-{v}_{z}$ approaching 0. When the initial flow-speed approaches the Alfvén speed, the strength of the over-reflection of the pulse at the gas-dynamic location must sufficiently increase the Alfvén speed upon reflection so that the backward propagation is not inhibited. It may be trivial to mention, but it is worth remembering that the fast- and slow-shock forward propagation will only increase as the flow-speed increases, which is another contributing factor into the decreasing timespan for the instability to amplify the Alfvén waves and "push" the gas-dynamic shock out of the shocktube.

Zoom In Zoom Out Reset image size

It has been shown that a small twist may amplify through over-reflection at a gas-dynamic shock. This occurs for a supersonic yet sub-Alfvénic flow. The Alfvén wave amplifies through harnessing kinetic energy of the supersonic flow upon reflection at the shock location and converts it into magnetic twists. Due to this strong reflection at the shock, and the inability for perturbations to propagate into the photosphere, the Alfvénic pulse becomes trapped between the photosphere and gas-dynamic shock. Here, it continually amplifies and broadens until a nonlinear system evolves due to the instability. The magnetic twists reach a maximum velocity of v_θ = 15.4 km s⁻¹ with the magnetic field eventually reaching a global twist of B_θ = 1.04 B_z. During the linear stage of the evolution, the twists propagate downstream quasi-periodically with each pulse lasting for about 10 s. However, as the amplitudes increase and nonlinearity sets in, the periodicity deteriorates and the pulses last longer. The propagating nonlinear twists are a signature of fast shocks accompanied by fluctuations in density, temperature and flow speed. The entire process of evolution from an undisturbed shocktube to a globally twisted tube takes about 50 minutes. This is within the lifetimes of Evershed flows (Rimmele 1994).

The nonlinear coupling means that the Alfvén pulse converts the kinetic energy of the plasma flow into magnetic energy upstream of the stationary shock. Downstream of the shock, these twists are somewhat converted back into the flow as kinetic energy. The Alfvén pulse thus acts as a means of energy transport from the upstream to the downstream plasma.

The conversion of kinetic energy from the background flow to magnetic energy in the upstream plasma leads to an imbalance in the total pressure gradient as a slow shock is set-up. At first, the gas-dynamic shock remains quasi-stationary as the Rankine–Hugoniot jump conditions are disturbed then restored. However, as the over-reflection of the Alfvénic pulse continues, the time taken for the jump conditions to restore increases until eventually they cannot be restored before the leading edge of the pulse has returned to the gas-dynamic shock. At this point, the gas-dynamic shock continually propagates until it exits the flux tube.

Before the shock is lost, the Alfvénic perturbations have sufficiently converted kinetic energy of the supersonic flow into magnetic twists such that the Alfvén instability is lost. The perturbations are now damped upon reflection and eventually cease all together.

The magnetic field, unlike the torsional velocity, continually twists in the same direction during this study. This is regardless of the propagation direction of the torsional perturbations. Similar twisting motions have been seen as magnetic swirls in solar tornadoes (Li et al. 2012; Su et al. 2014). Wedemeyer-Böhm & Steiner (2014) showed numerically that such spiral motions are caused by the rotation of magnetic field structures, which are driven by photospheric vortices. However, previous work by Shelyag et al. (2013) demonstrates that these motions are due to torsional Alfvén waves.

The twisting motions are not only confined to the region between the photosphere and gas-dynamic shock. It is seen in Figure 6 that the portion of the Alfvénic pulse, which is transmitted through the gas-dynamic shock, induces a global magnetic twist. This means that the mechanism could permanently twist a larger flux tube without the gas-dynamic shock propagating along the entirety of the flux tube. The final state of the shocktube presented is similar to the equilibrium assumed by Zaqarashvili et al. (2010) and Díaz et al. (2011). They consider a magnetic flux tube with a helical magnetic field and mass flow along the z-axis. This means the magnetic field structure is "pulled" along the flux tube without any deformation. These twisted flux tubes are prone to the kink instability, which may lead to explosive events such as solar flares and CMEs.

The amplification of the Alfvén wave results in velocities ≈10–15 km s⁻¹, which is within the lower limits of the ranges observed in the chromosphere (De Pontieu et al. 2007, 2014). The observed twists last for less than a minute and propagate at the local Alfvén speed. However, the torsional amplitude of the waves presented in this model is somewhat lower than the reported values in Type-II spicules, which have torsional velocities of 25–30 km s⁻¹ (De Pontieu et al. 2012). That said, the field-aligned flow of Type-II spicules observed by De Pontieu et al. (2012) is between 50 and 100 km s⁻¹, which is 2–3 times the initial background flow employed in this study. Also, the amplitudes are expected to increase with increasing altitude when gravitational stratification is taken into account.

The θ-energy flux density observed during this study is 4–5x larger than that seen in spicular studies (Hollweg et al. 1982; Kudoh & Shibata 1999; Matsumoto & Shibata 2010). However, this is not necessarily a negative attribute as these studies either meet the lower end of the energy requirement to heat the quiet corona (Kudoh & Shibata 1999) or fail to meet even these minimum requirements for coronal heating (Hollweg 1992). This increased flux density is due to the background flow, the energy of which is extracted and transferred downstream by amplified twists. That said, the values employed are all consistent with the lower chromosphere.

Incorporating gravity and a magnetic flux tube that expands as a function of height, such that rB_z ≈ constant, where r is the flux tube radius, leads to increased torsional velocities as an Alfvén wave propagates forwards (Hollweg 1992; James et al. 2003; Erdélyi & James 2004). These additions to the model will only serve to further increase the amplitude of the torsional motions seen, and possibly strengthen the effects of the Alfvénic instability along the shocktube.

Using the 1.5-D approximation for a magnetic shocktube, it has been shown that a gas-dynamic shock traps and amplifies a twist when the flow is supersonic and sub-Alfvénic. Propagating fast shocks are produced due to the nonlinear coupling in Equations (2) and (4) once the Alfvén wave is sufficiently amplified. The amplification and broadening of the Alfvénic pulse leads to magnetic twisting motions that are always in the same direction, much like the swirls seen in solar tornadoes (Li et al. 2012). This continuous twisting in the same direction leads to a global twist along the flux tube where B_θ becomes, and remains, comparable to B_z.

The amplification strengthens the nonlinear effects and destabilizes the gas-dynamic shock. Initially, the shock sways back and forth as the total pressure balance is repeatedly disturbed by the impelling signals.

Eventually, the pressure balance can no longer be restored and the gas-dynamic shock propagates indefinitely. The twisting motions cease to amplify at this point as the Alfvénic pulse becomes damped in the absence of the Alfvén instability. The end result is a highly twisted magnetic flux tube with a sub-Alfvénic mass flow which pulls the magnetic field structure along the z-axis. This is much like the flux tubes used as a starting point for kink instability studies by Zaqarashvili et al. (2010) and Díaz et al. (2011).

The torsional energy flux generated by this mechanism is comparable to spicular models that provide sufficient energy to heat the corona (Kudoh & Shibata 1999; Matsumoto & Shibata 2010). However, adding gravitational stratification, losses and expansion of the flux tube will affect the Alfvén instability and torsional flux generated by the model and will be addressed in a future study.

The authors would like to thank the referee for their helpful comments and T.W. would also like to thank the STFC for their financial support.