Numerical Simulations of Torsional Alfv\'en Waves in Axisymmetric Solar Magnetic Flux Tubes

We investigate numerically Alfv\'en waves propagating along an axisymmetric and non-isothermal solar flux tube embedded in the solar atmosphere. The tube magnetic field is current-free and diverges with height, and the waves are excited by a periodic driver along the tube magnetic field lines. The main results are that the two wave variables, the velocity and magnetic field perturbations in the azimuthal direction, behave differently as a result of gradients of physical parameters along the tube. To explain these differences in the wave behavior, the time evolution of the wave variables and the resulting cutoff period for each wave variable are calculated, and used to determine regions in the solar chromosphere where strong wave reflection may occur.


Introduction
Observations by several recently launched spacecrafts have revealed the ubiquitous presence of oscillations in the solar atmosphere, which can be interpreted as magnetohydrodynamics (MHD) waves (e.g. Nakariakov and Verwichte, 2005) or specifically as Alfvén waves whose signatures were observed in prominences, spicules and X-ray jets by Okamoto (2007), De Pontieu (2007) and Cirtain D.Wójcik e-mail: dwojcik@kft.umcs.lublin.pl (2007), respectively. Moreover, observational evidence for the existence of torsional Alfvén waves in the solar atmosphere was given by Jess (2009), see however Dwivedi and Srivastava (2010) and Okamoto and De Pontieu (2011) who reported the presence of propagating Alfvén waves in solar spicules. Alfvén waves have been of a particular interest because they can carry energy and momentum along magnetic field lines to the solar corona, where the wave energy can heat the corona and the wave momentum may accelerate the solar wind.
There is a large body of literature devoted to Alfvén waves and their propagation in the solar atmosphere (Zhugzhda and Locans, 1982;Hollweg, Jackson and Galloway, 1982;An et al., 1989;Hollweg, 1990;Musielak, 1992;Musielak and Moore, 1995;Kudoh and Shibata, 1999;Hollweg and Isenberg, 2007;Musielak, Routh and Hammer, 2007;Murawski and Musielak, 2010;Routh, Musielak and Hammer, 2010;Webb et al., 2012;Chmielewski et al., 2013;Perera, Musielak and Murawski, 2015). There are two main problems considered in these papers, namely, the propagation conditions for Alfvén waves and the dissipation of energy and momentum carried by these waves. The first problem involves the concept of cutoff frequency whose existence is caused by the presence of strong gradients of physical parameters in the solar atmosphere, and it has been explored for linear Alfvén waves by Murawski and Musielak (2010) and Perera, Musielak and Murawski (2015), and for torsional Alfvén waves in thin magnetic flux tubes by Musielak, Routh and Hammer (2007) and Routh, Musielak and Hammer (2010). The second problem deals with the coronal heating (Suzuki and Inutsuka, 2005) and it involves different mechanisms of Alfvén wave dissipation, like phase-mixing (Ofman and Davila, 1995) or nonlinear mode coupling (Ulmschneider, Zähringer and Musielak, 1991), and wave momentum deposition (Dwivedi and Srivastava, 2006;Chmielewski et al., 2013Chmielewski et al., , 2014; however, a realistic modeling of both Alfvén wave propagation and dissipation is difficult to perform (Banerjee, Hasan and Christensen-Dalsgaard, 1998;O'Shea, Banerjee and Doyle, 2005;Bemporad and Abbo, 2012;Chmielewski et al., 2013;Jelínek et al., 2015). Murawski and Musielak (2010) considered impulsively generated Alfvén waves in the one-dimensional solar atmosphere with a smoothed step-wise temperature profile and a vertical magnetic field. Perera, Musielak and Murawski (2015) studied analytically and numerically the case of periodically driven Alfvén waves and their propagation in an isothermal solar atmosphere. Musielak, Routh and Hammer (2007) and Routh, Musielak and Hammer (2010) investigated torsional Alfvén waves propagating in thin magnetic flux tubes embedded in the isothermal and non-isothermal solar atmosphere, respectively. The main aim of this paper is to extend the work of Murawski and Musielak (2010) and Perera, Musielak and Murawski (2015) to an axi-symmetric solar magnetic flux tube embedded in the solar atmosphere with a realistic temperature profile of Avrett and Loeser (2008) and curved magnetic field lines (Low, 1985), perform numerical simulations of the propagation of torsional Alfvén waves, and to compare the obtained numerical results to the analytical ones obtained by Routh, Musielak and Hammer (2010). This paper is organized as follows. The numerical model of the solar atmosphere is described in Section 2. The numerical and analytical results are presented and discussed in Sections 3 and 4, respectively. Our discussion of the analytical and numerical result is given in Section 5, and a short summary of the results and some concluding remarks are presented in Section 6.

Numerical Model
We consider a magnetically structured and gravitationally stratified solar atmosphere, which is described by the following set of ideal MHD equations: where ̺ is mass density, p is gas pressure, V is the plasma velocity, B the magnetic field, and g = (0, −g, 0) is the gravitational acceleration of magnitude 274 m s −2 . The symbol T denotes temperature, m is particle mass, specified by a mean molecular weight of 1.24, k B is the Boltzmann's constant, γ = 1.4 is the adiabatic index, and µ the magnetic permeability of plasma.
We consider the axisymmetric case in which all plasma variables are invariant in the azimuthal direction θ, but the θ-components of the perturbed magnetic field B θ and velocity V θ are not identically zero. We assume that the equilibrium magnetic field is current-free, ∇ × B = 0, and its components along the radial (r), azimuthal (θ), and vertical (y), directions are specified by the following expressions: where a and S are free parameters which denote the vertical location of the magnetic moment and the vertical magnetic field strength, respectively. Equations (5)-(7) comprise a special (potential and axisymmetric) case of the threedimensional model which was developed by Low (1985). We set a = −0.75 Mm choose S such that at r = 0 Mm and y r = 6.0 Mm the magnetic field is B e = B ey = 9.5 Gs. Here, y = y r denotes the reference level. The magnetic field vectors resulting from these equations are drawn in Figure 1. The magnetic field lines are curved and the curvature decays with height y, it grows with the radial distance r, and the magnitude of B e decays with y and r. As the equilibrium magnetic field, specified by Equations (5)-(7) is currentfree, the equilibrium gas pressure and mass density are hydrostatic and they are given as (e.g. : where p 0 = 0.01 Pa is the gas pressure evaluated at y = y r , is the pressure scale-height, and the symbol T h (y) denotes the hydrostatic temperature profile specified in the model of Avrett and Loeser (2008). In the model, the solar photosphere occupies the region 0 < y < 0.5 Mm, the chromosphere resides in 0.5 Mm < y < 2.1 Mm, and the solar corona is represented by higher atmospheric layers, which start from the transition region that is located at y ≈ 2.1 Mm. Note that in equilibrium the gas pressure and mass density do not vary with the radial distance r, while they fall off with the height y (not shown). Figure 2 shows the vertical profile of the Alfvén speed c A = B e / √ µ̺ h along the central axis, r = 0 Mm, with c A = 12 km s −1 at y = 0 Mm. Higher up c A grows with y, reaching maximal value of 2715 km s −1 at y = 2.44 Mm, and then c A falls off with height. This characteristic shape of the Alfvén velocity profile will have a prominent effect on the wave behavior in our model as shown and discussed in Section 3.

Numerical Results
Numerical simulations are performed using the PLUTO code, which is based on a finite-volume/finite-difference method (Mignone et al., 2007(Mignone et al., , 2011, designed to solve a system of conservation laws on structured meshes. For all considered cases we set the simulation box as 0 to 5.12 Mm in the r direction and −0.1 Mm to 30 Mm in the y direction. We impose the boundary conditions by fixing at the top and bottom of the simulation region all plasma quantities to their equilibrium values, while at the right boundary, outflow boundary conditions are implemented. These outflow boundary conditions do no lead to any incoming signal reflection as we use the background magnetic field splitting (Mignone et al., 2007). At r = 0 Mm axisymmetrical boundaries are implemented. Along the r-direction this simulation box is divided into 1024 cells and along the ydirection into 1536 cells for −0.1 < y < 7.58 Mm and into 512 cells in the range of 7.58 < y < 30 Mm. For our problem the Courant-Friedrichs-Lewy number is set to 0.3 and the Harten-Lax-van Leer Discontinuities (HLLD) approximate Riemann solver (Toro, 2009) is adopted. At the bottom boundary the periodic driver is additionally set as where P d is the period of the driver, A V is the amplitude of the driver and w is its spatial width. We set and hold fixed w = 100 km and A V = 5 km s −1 , while allowing P d to vary. Figure 4 illustrates the r-y spatial profiles of V θ (left) and B θ (right) for P d = 35 s which corresponds to an effective amplitude of the driver of about  0.2 km s −1 . It is seen in V θ (r, y) that this driver excites Alfvén waves which penetrate the chromosphere, and at t ≈ 70 s reach the transition region and the solar corona. As a result of the Alfvén speed profile, displayed in Figure 2, the Alfvén waves accelerate up to the altitude where c A attains its maximal value, and higher up the Alfvén waves decelerate. Note that as a consequence of the divergence of the magnetic field with height (Figure 1), the wave front spreads with height and magnetic shells develop in time  and a standing wave pattern is present in the low photosphere, below y ≈ 0.5 Mm. A qualitatively similar scenario can be seen in the low atmospheric layers in profiles of B θ (r, y) (Figure 4, right). However, the perturbations in B θ decay with altitude, remaining largest in low atmospheric regions. This conclusion confirms the former findings of Murawski and Musielak (2010). The patterns of standing waves and magnetic shells are easily seen below the transition region. Figure 3 displays the vertical profiles of V θ (left) and B θ (right) at r = 0.1 Mm. The penetration of Alfvén waves into the solar corona is seen in the V θ profiles (left). The presented results clearly show that the wave variables behave differently, namely the gradients in physical parameters have stronger effects on B θ than on V θ because the former decays rapidly with height while the latter reaches the solar corona; the fact that the Alfvén wave variables behave differently is a well-known phenomenon (e.g. Hollweg, Jackson and Galloway, 1982;Musielak, 1992;Musielak and Moore, 1995;Musielak, Routh and Hammer, 2007;Murawski and Musielak, 2010;Routh, Musielak and Hammer, 2010;Perera, Musielak and Murawski, 2015). The presence of gradients is the main physical reason for such a different behavior of V θ and B θ , and it also leads to wave reflection, which occurs below y ≈ 1.75 Mm and decays with time; the waves that undergo reflection become trapped and this trapping dominates in the upper parts of the solar chromosphere. Similar to V θ (r = 0.1 Mm, y, t), the profiles of B θ reveal the transient phase, which occurs for t 150 s but later on the oscillations in B θ reach a quasi-stationary stage (bottom).
The above numerical results have important physical implications, namely, they show that Alfvén waves do lose their identity because their wave velocity and magnetic field variables behave differently. Since the wave velocity variable reaches the corona and the wave magnetic field variable cannot, the Alfvén wave in the solar corona does not obey any more its equipartition of energy between the two wave variables. This means that the physical nature of the Alfvén wave in the solar corona is different than in the solar chromosphere. The problem is out of the scope of this paper but it does require future studies.

Cutoff Periods for Torsional Waves in Thin Flux Tubes
To compare the results of our numerical simulations to the analytically determined conditions for the propagation of torsional Alfvén waves in a solar magnetic flux tube, we follow Musielak, Routh and Hammer (2007) and Routh, Musielak and Hammer (2010). With the wave variables v = v θ (r, y, t)θ and b = b θ (r, y, t)θ, we follow Musielak, Routh and Hammer (2007) and write the θ-components of the linearized induction equation of motion and induction as and According to Musielak, Routh and Hammer (2007), the thin flux tube approximation gives the following relationship: which allows writing Equations (11) and (12) as and ∂b θ ∂t Assuming the tube being in temperature equilibrium with its surroundings, Routh, Musielak and Hammer (2010) used horizontal pressure balance and obtained the following wave equations for torsional Alfvén waves: and where H(y) = c 2 s (y)/γg is the pressure scale height with c s (y) being the sound speed, and c A (y) = B ey (y)/ ρ h (y) is the Alfvén velocity.
As demonstrated by Routh, Musielak and Hammer (2010), the critical frequencies are given by and Ω 2 cr,b (y) = 1 2 1 2 and the resulting turning-point frequencies become and where t ac is actual wave travel time expressed by with y b being an atmospheric height at which the wave is initially generated. Finally, the cutoff frequency for torsional Alfvén waves propagating in a thin and non-isothermal magnetic flux tube embedded in the solar atmosphere is Ω cut,τ (y) = max [Ω tp,τ,v (y), Ω tp,τ,b (y)] .
Having obtained the turning-point frequencies and the cutoff frequency, we may define the turning-point periods as and the cutoff period is Figure 5 Cutoff period for Alfvén waves vs height y calculated from Equations (23) and (25).

Discussion of Analytical and Numerical Results
Figure 5 displays the cutoff period, P cut , vs height y. At y = 0, which is the bottom of the photosphere, P cut = 10 s. Higher up it grows, reaching its maximum value of about 300 s at y ≈ 0.2 Mm, and subsequently it falls off with height to its global minimum of about 0.03 s, which occurs at the solar transition region. In the solar corona P cut increases again with y but its values are lower than 1 s. This general trend of P cut (y) is consistent with that found by Murawski and Musielak (2010), however, the exact values of P cut (y) are different in the two cases. In particular, Murawski and Musielak (2010) found that for the straight vertical equilibrium magnetic field P cut was higher than 40 s, while Figure 5 reveals the corresponding values being a fraction of 1 s. This clearly shows how strongly the cutoff depends on the geometry and strength of the background magnetic field. Our numerical results show that Alfvén waves -or at least their wave velocity perturbation -penetrate the transition region, whereas the magnetic field perturbation does not, as shown in Figures 3 and 4. The analytical results imply the same as they also show that the two wave variables behave differently. However, there are some discrepancies between the analytical and numerical results because the thin flux tube approximation, which is the basis for obtaining the cutoff frequency, does not apply to our model of the solar corona because the tube becomes thick and the condition given by Equation (13) is not satisfied any longer. The agreement between the analytical theory and our numerical model is better in the solar chromosphere where the tube remains thin. Because of the different behavior of the two wave variables, one may look directly at the values of the turning-point frequencies given by Equations (20) and (21) for each wave variable. Since in the photosphere for 0.2 Mm < y < 2.5 Mm, P tp,b = 0, a signal in magnetic field is not able to penetrate this region, which explains why B θ is absent in the solar corona (Figure 3, right); this is an important result, which clearly shows the connections between the turningpoint frequencies and the different behavior of the two wave variables. However, a signal can be seen in this region in V θ (Figure 3, left). Indeed, for P d = 35 s the condition of P < P tp,v is satisfied only up to y ≈ 0.3 Mm, while for y ≈ 0.1 Mm P < P tp,b . For larger values of y, we have P > P tp,bv and the Alfvén waves become evanescent. Figure 6 illustrates spectral power corresponding to P d = 35 s for V θ (left) and B θ (right). A sudden fall-off of power at y ≈ 0.2 Mm corresponds to energy being spent on excitation of a non-zero signal in B θ . Thus energy goes from V θ , in which it was originally generated, into B θ . We point out that for linear Alfvén waves propagating in homogeneous media, there is a perfect equipartition of the wave kinetic energy associated with V θ and the wave magnetic energy associated with B θ ; our results show that this perfect equipartition of wave energy is strongly violated in inhomogeneous media (see our comments at the end of Section 3).
In the solar corona the value of the P d is larger then a local value of P tp,v and P tp,b so the waves are evanescent in this region. However, our numerical results show a tunneling of V θ (left) to the solar corona with the signal in B θ (right) falling off with height more strongly, with and B θ being practically zero above y = 1.5 Mm.   Figure 8). Note that oscillations take place until the zero of V θ and higher up V θ becomes evanescent as it falls off with altitude above its local maximum. The oscillatory behavior of B θ until the last zero and its subsequent strong fall-off are clearly seen.
The turning points for V θ and B θ , which separate oscillatory behavior of the wave variables from non-oscillatory behavior, are illustrated in Figure 9. Note that the turning points are located at lower heights for larger values of P d with the falling-off trend confirming the numerical findings of Perera, Musielak and Murawski (2015), and that the turning points for B θ are located at lower altitudes (see Figure 9). There is strong consistency between our conclusions based on the cutoff and on the concept of turning points, and the reason is that all these two concepts are related to each other (e.g. Musielak, 2006).

Summary and Concluding Remarks
In this paper, we simulate numerically Alfvén waves which are driven by a periodic driver operating in the solar photosphere 100 km below the solar surface. Our findings can be summarized as follows (c) This different behavior of the two wave variables can be explained by the fact that the variables for Alfvén waves have different turning point frequencies, and that these frequencies are used to define the cutoff frequency (or period), which sets up the conditions for the wave propagation in the solar atmosphere.
(d) The presence of the cutoff for Alfvén waves implies that the waves with periods close to the cutoff are being reflected and that interaction between incoming and reflected waves leads to the formation of standing wave patterns.
(e) Our numerical and analytical results demonstrate that the Alfvén wave loses its identity while traveling in the solar atmosphere because the wave velocity perturbation reaches the solar corona but the wave magnetic field perturbation does not; as a result, the wave variables do not satisfy equipartition of energy and thus the Alfvén wave has a different nature in the upper atmospheric layers.
The obtained results extend the previous work of Murawski and Musielak (2010) and Perera, Musielak and Murawski (2015) to describe torsional Alfvén waves propagating along solar magnetic flux tubes. From a mathematical point of view, the previous and current work is different, namely, our study is mathematically more complex. However, from a physical point of view, there are some differences and similarities between Alfvén waves studied in the previous work and the torsional Alfvén waves investigated here. In particular, the previous work by Musielak, Routh and Hammer (2007) and Routh, Musielak and Hammer (2010) was used to obtain the cutoff frequency for torsional Alfvén waves propagating along a thin and non-isothermal solar magnetic flux tube. The approach formally takes into account both B er (r, y) and B ey (r, y), however, only within the thin flux tube approximation. Our analytical results in Section 4 show that the turning-point frequencies are the same as those obtained by Murawski and , who considered regular Alfvén waves propagating along a uniform background magnetic field, however, the cutoff wave periods are different. The results presented here clearly show that the wave variables behave Figure 9 The vertical location of the turning point for B θ (triangles) and V θ (squares) vs driving period P d .
differently and that the atmospheric inhomogeneities have stronger effects on the wave magnetic variable than on the wave velocity variable; again, this result is consistent with that obtained previously by Murawski and Musielak (2010) and Perera, Musielak and Murawski (2015). Moreover, we demonstrate that one of the best ways to account for these differences is to use the turning point frequencies because their values uniquely determine the behavior of each wave variable in the solar atmosphere, and because they are used to define the cutoff frequency (or period) for Alfvén waves. It is also important to emphasize that the different behavior of the Alfvén wave variables in inhomogeneous media leads to violation of the equipartition of the wave kinetic and magnetic energies in Alfvén waves and that this fact makes a clear distinction between Alfvén waves propagating in homogeneous and inhomogeneous media.