Varying Driver Velocity Fields in Photospheric MHD Wave Simulations

Torsional motions are ubiquitous in the solar atmosphere. In this work, we perform 3D numerical simulations which mimic a vortex-type photospheric driver with a Gaussian spatial profile. This driver is implemented to excite MHD waves in an axially symmetric, 3D magnetic flux tube embedded in a realistic solar atmosphere. The Gaussian width of the driver is varied and the resulting perturbations are compared. Velocity vectors were decomposed into parallel, perpendicular and azimuthal components with respect to pre-defined magnetic flux surfaces. These components correspond broadly to the fast, slow and Alfven modes, respectively. From these velocities the corresponding wave energy fluxes are calculated, allowing us to estimate the contribution of each mode to the energy flux. For the narrowest driver ($0.15$ Mm) the parallel component accounts for $\sim 55 - 65\%$ of the flux. This contribution increases smoothly with driver width up to nearly $90\%$ for the widest driver ($0.35$ Mm). The relative importance of the perpendicular and azimuthal components decrease at similar rates. The azimuthal energy flux varied between $\sim 35\%$ for the narrowest driver and $<10\%$ for the widest one. Similarly, the perpendicular flux was $\sim 25 - 10\%$. We also demonstrate that the fast mode corresponds to the sausage wave in our simulations. Our results therefore show that the fast sausage wave is easily excited by this driver and that it carries the majority of the energy transported. For this vortex-type driver the Alfven wave does not contribute a significant amount of energy.


INTRODUCTION
Magnetohydrodynamic (MHD) waves are ubiquitous in the solar atmosphere and it is considered likely by many that they contribute to solar atmospheric heating by transporting energy from the photosphere up through the lower solar atmosphere and into the low corona. There have been numerous observations in various magnetic structures of each of the MHD wave modes -fast, slow and Alfvén. The fast mode, in particular, is frequently seen having been observed in sunspots (e.g. Dorotovič et al. 2014), pores (e.g. Morton et al. 2012;Dorotovič et al. 2014;Freij et al. 2014) and other magnetic structures in the chromosphere (Morton et al. 2012). Dorotovič et al. (2014) also observed the slow mode. Alfvén waves have been observed in a group of bright points by Jess et al. (2009), andMcIntosh et al. (2011) claim to have detected them in the corona. For reviews of the wide range and variety of wave observations see e.g. Nakariakov and Verwichte (2005) Torsional motions have great potential to excite Alfvén waves, the favourite candidate for energy transport in solar MHD (see e.g. Mathioudakis et al. (2013) for a review of Alfvén wave observations and theory). Therefore torsional motions have been searched for and have been successfully observed at e.g. intergranular lanes in the form of resolution-limited small-scale vortices (e.g. Bonet et al. 2008;Wedemeyer-Böhm and Rouppe van der Voort 2009;Bonet et al. 2010). It is widely accepted that these vortices form due to turbulent convection (e.g. Shelyag et al. 2011;Wedemeyer-Böhm et al. 2012;Kitiashvili et al. 2013).
Given the ubiquity of these vortex motions in the photosphere, it is important to understand how the waves they excite contribute to the heating of the lower solar atmosphere. To this end, several three-dimensional simulations have been performed by (e.g. Fedun Mumford et al. 2015). These studies implemented torsional motions at the base of a realistic magnetic flux tube and analysed the resulting perturbations. In each case it was found that such a driver excites fast and slow magnetoacoustic waves and the Alfvén wave, and that in all but one case the sausage and kink modes were both present. Vigeesh et al. (2012) and Mumford et al. (2015) also quantified the energy flux of waves produced by torsional motions and found that the azimuthal components of the waves made a greater contribution to the flux than the perpendicular or parallel components.
The present work is a continuation of the work of Mumford et al. (2015), which investigates the effects of varying driver parameters on the wave motions stimulated by those drivers in the low solar atmosphere. In this case, we now implement a spiral velocity driver and investigate how varying that driver's width scales the wave energy transport from the driver into the lower solar atmosphere. Since a range of vortex sizes are observed, we wish to investigate whether this variation causes different waves, as this information will have implications for atmospheric heating. We also outline a new way to unambiguously demonstrate the presence of sausage and kink modes (whether slow or fast, depending on the equilibrium conditions) in our simulations by calculating the displacement of the magnetic flux surface from its original position.
For this study we use the Sheffield Advanced Code (SAC; Shelyag et al. 2008), which is built on the basis of the Versatile Advection Code (VAC; Tóth 1996). SAC separates variables into background and perturbed components, allowing the simulation of highly gravitationally stratified media such as the solar atmosphere.
This paper is structured as follows: in Section 2 we describe the background atmosphere and the properties of the photospheric drivers employed in the simulations; in Section 3 we describe the simulation parameters and the analysis method; in Section 4 we present the results of the simulations and the analysis; in Section 5 we discuss these results and present our conclusions.

BACKGROUND ATMOSPHERE AND PHOTOSPHERIC DRIVERS
Here, we use a three dimensional background atmosphere ( Figure 1) based on the VAL IIIC model (Vernazza et al. 1981), and implement an axisymmetric magnetic flux tube modelled based on the self-similar approach (Schluter and Temesvary 1958;Deinzer 1965;Schüssler and Rempel 2005;Fedun et al. 2011a;Gent et al. 2013). The full-width at half-maximum of the magnetic flux tube (FWHM) was approximately 90 km in the photosphere. The footpoint of the magnetic flux tube was centred at x, y, z = (1.0, 1.0, 0.0) Mm, and had a magnetic field strength of 143.6 mT (see Figure 1). For more details on this initial configuration, see Mumford et al. (2015) -who use the same background atmosphereand references therein.
In each simulation, perturbations to the background atmosphere are driven by introducing a velocity field in the horizontal plane close to the footpoint of the flux tube. These drivers are intended to mimic different kinds of velocity fields that may be found in the photosphere, as a result of granu-  Figure 2. Normalised horizontal profile of the velocity field generated by the implemented spiral driver. Colour-coding indicates the magnitude of the velocity and the cyan arrows indicate the direction. This velocity profile is multiplied by the magnitude of the driver (see text) and varies with both height and time.
lation. This paper uses the same logarithmic spiral velocity fields as Mumford et al. (2015), to study the excitation of torsional waves in the atmosphere. The logarithmic spiral shape is based on observations by Bonet et al. (2008) and others of vortex flows in intergranular lanes (see Section 1). The logarithmic spiral driver is centred on the point x, y, z = (1.0, 1.0, 0.1) Mm, the spatial extent of which is determined by a Gaussian profile in each direction. The velocity at a given point and at time t is described by: where is the Gaussian profile with width ∆x, ∆y and ∆z in the x-, y-and z-directions, respectively. A is the driver amplitude, P is the driver period, θ = tan −1 (y/x) is the angle around  Table 1. Driver width ∆x = ∆y and corresponding driver amplitude values used to ensure the same input of total kinetic energy to each simulation. The middle column indicates the ratio of the width of the driver to the FWHM of the flux tube.
the flux tube axis and φ = tan −1 (1/0.15) determines the expansion of the spiral. In a series of numerical simulations, (Mumford 2016, Chapter 6) found that the period of such a swirly motion has only a relatively minor effect on the results of this kind of study, with the contribution from the Alfvén mode, for instance, varying by less than 20%. The value of the period for this study was therefore mostly chosen so that a few periods would fit into the run-time of the simulation, and was set to 90 s. The choice of the expansion parameter, 0.15, was also largely arbitrary and was selected to allow a few rotations of the spiral within the driver. We use five values for the horizontal Gaussian width of the driver, ∆x = ∆y, as indicated in Table 1. This range of parameter values was chosen to correspond to the range of major axes found by Sánchez Almeida et al. (2004) for a sample of 126 magnetic bright points (MBPs) observed in intergranular lanes, and are consistent with Bonet's observations that these vortexes have sizes of 0.5 Mm. Each driver had the same vertical width, ∆z = 0.05 Mm.
All drivers were designed to supply the same total amount of energy, ET , into the simulation. To ensure this, the amplitude of the driver was adjusted according to the width of the driver. The exact relation is In the above, V is the volume of the computational domain, ρ is the density, and n is the integer number of periods in the simulation run-time. The actual amplitudes corresponding to the widths implemented in the simulations are listed in Table 1.

SIMULATIONS
The simulation domain ranged from 0.0 Mm to 2.0 Mm in the x-and y-directions, and from 0.0 Mm to 1.6 Mm in the zdirection, with a mesh size of 128, resulting in 128 3 grid cells.
The boundaries of the domain were set to the 'continuous' setting in SAC (i.e.: the gradient of each variable was zero across the boundaries). Each simulation was run for 270 s of simulation time, equal to three full driver periods. This amount of time is approximately equal to the lifetime of vortex flows observed in the photosphere by Bonet et al. (2008). We can therefore be confident that the flux tube would reasonably remain stable within the runtime of the simulation.

Velocity vector decomposition
A flux surface is constructed by selecting seed points on a circle near the top of the domain centred on the flux tube axis. Field lines are traced down through the domain from those seed points to the bottom of the domain using the method of Mumford et al. (2015). These field lines then enclose a constant amount of magnetic flux at any given height and thus describe the surface of a flux tube. We refer to these field lines and flux surfaces by the radius of the circle of seed points, as a fraction of the maximum radius possible in the domain (64 grid cells). These field lines are retraced from these advected seeds at each time-step and new flux surfaces are calculated. This treatment allows us to separate the velocities into components which are locally parallel to the direction of the magnetic field, perpendicular to the magnetic flux surface and azimuthal around the flux tube. These components correspond broadly to the fast and slow MHD waves and the Alfvén wave, respectively.
Of course, with this interpretation of wave modes, one has to bear in mind the local value of the plasma-β, where β is the ratio of kinetic to magnetic pressure. Given the fairly weak magnetic field of our flux tube, plasma-β is large (> 1) everywhere in the simulation domain except for a small region at the top of the domain close to the flux tube axis. Therefore the slow mode propagates mainly along magnetic field lines with the local Alfvén speed, vA. The fast mode is allowed to propagate in any direction. Along field lines the fast mode travels with the sound speed cs, while in the direction perpendicular to the magnetic field it propagates with the phase speed vp = c 2 s + v 2 A .

Velocity and flux calculation
Following Mumford et al. (2015), we interpolate the decomposed velocity vectors onto a single field line in order to study how waves propagate along this field line throughout the simulation. In addition to the velocity components, we define the wave energy flux using the following equation (Leroy 1985;Bogdan et al. 2003;Mumford et al. 2015): wherep k is the kinetic pressure perturbation, Here, v is the velocity, B is the magnetic field, e is the total energy density, µ0 is the permeability of free space and γ is the adiabatic index of the plasma. Background and perturbed components of quantities are indicated by a subscript b and a tilde, respectively. Once calculated, the energy flux vector can be decomposed into its parallel, perpendicular and azimuthal components in the same way as the velocity vector.

Flux surface displacement
We calculate the distance at a number of angular positions around the axis by finding the intersection of the flux surfaces with lines through the axis at those angles. These distances are calculated for a number of heights in the domain and   The sausage wave, which distorts the flux tube in the same direction (i.e.: towards or away from the axis) at all angles, will manifest as the displacement of any two points having the same sign. Conversely, the kink mode moves the flux tube towards the axis on one side and away from it on the opposite side, resulting in negative and positive displacement of points on those sides, respectively.
Torsional motions could be detected by inspecting azimuthal velocity, v θ , using the same method, but doing so is not trivial due to technical limitations -this will be addressed in a later work. This part of the analysis is therefore intended only to determine the presence of sausage and kink motions.

Wave mode identification
The displacement of flux surfaces with respect to their original positions was calculated as described above for twenty magnetic flux surfaces throughout the domain. The radii of the seed point circles for these flux surfaces were equally spaced between r = 0.047 and r = 0.984 at intervals of 3 grid cells, and each seed point was initially located 10 grid cells below the top of the domain at z = 1.475 Mm. Figure 3 shows the displacement of these flux surfaces at three different heights in the domain. The plots in the left-hand column of this figure show the displacement of each point from its original position. In these plots the radial and azimuthal axes indicate the position of each point with respect to the flux tube axis and the colour scale shows the displacement of the flux surface at that point from its original position. Motion away from and towards the axis of the flux tube are shown in red and blue, respectively. The plots in the right-hand column show this information only for points at θ = 136.8 • and at θ = 316.8 • , that is, points on opposite sides of the flux tube. These angles are chosen to align with the direction of displacement of the flux tube to most clearly show the wave motion. To aid the analysis of these plots the magnitude of velocity is also plotted in Figure 4 for a vertical slice through the centre of the domain.
Near the bottom of the domain (Figure 3a, b), we see clear motion towards the axis on one side and away from it on the other, indicating a kink wave. However, near the top (Figure 3e, f), the motion on either side of the axis is in phase, either towards the axis on both sides or away from it on both sides. This demonstrates the presence of a sausage wave. In the middle of the domain (Figure 3c, d), both waves are visible in different parts of the domain. In this case, the kink wave is dominant close to the axis ( 230 km) and the sausage mode becomes more dominant further away.
This interpretation is consistent with the velocity magnitudes plotted in Figure 4, which shows two distinct wave fronts. One of these propagates almost isotropically and has reached the top of the domain, indicating that it is a fast mode. The other wave front has not reached as great a height and is more closely confined to the magnetic field, and must therefore be the slow mode. Since we can see in Figure 3 that only the sausage wave is visible at the top of the domain, this must correspond to the fast mode. Similarly, the kink wave must correspond to the slow mode, since both are only seen close to the flux tube and in the lower half or so of the domain.
From this analysis we identify fast sausage waves and slow kink waves in our simulations without ambiguity. This identification provides useful context to the rest of our analysis.

Velocity components
We select a field line at r = 0.469. The changes in the velocities along this field line with time are plotted in the time-distance diagrams shown in Figure 5 for the narrowest and widest drivers used in the simulations. We also calculate the values of the fast speed (v f ), slow speed (also called the tube speed, vt), sound speed (cs) and Alfvén speed (vA) along this field line (Equations 5) and plot these for comparison. These values are calculated thus: where p is the total (background plus perturbed) kinetic pressure. Figure 5 shows that for the narrowest driver, 0.15 Mm (1.67 FWHM), the azimuthal velocity component is most dominant with an absolute value of ∼ 30 ms −1 , compared to ∼ 16 ms −1 for both the parallel and perpendicular components. For the widest driver, meanwhile, the absolute value of the azimuthal perturbation is lower (∼ 16 ms −1 ), whereas the parallel and perpendicular perturbations have increased (to ∼ 60 ms −1 and ∼ 30 ms −1 , respectively). We see then that the azimuthal component is most dominant for a narrow driver, whereas a wider driver produces perturbations in which the parallel component is greatest.

Flux contribution
Time-distance diagrams for the contribution of each wave flux component are shown in Figure 6, again for the narrowest and widest drivers. This contribution is expressed as the square of each flux component as a fraction of the total square flux, so that the sum of the three component contributions is equal to unity.
In the case of the widest driver, 0.35 Mm (3.89 FWHM) we can see that the majority of the wave flux is contained in the parallel component. This is particularly clear in the region above the height reached by slow magnetoacoustic and Alfvén waves, where some contribution comes from the perpendicular component but almost none can be seen in the azimuthal component. The case for the narrowest driver is similar at these heights. Below this, however, the relative importance of the perpendicular and azimuthal components is much greater compared to the widest driver. Figure 7 plots the percentage square wave flux, averaged over the full simulation run-time for each flux component and for each driver width. This gives a broad indication of how dominant each component is over the simulation as a whole. This comparison is plotted for each of three flux surfaces at different distances from the flux tube axis. The parallel compenent varied between ∼ 50% and ∼ 90%, the perpendicular component between ∼ 20% and ∼ 10%, and the azimuthal component ∼ 35% and ∼ 10%. As spiral width increases, the influence of the parallel component increases, while the roles of the other components decrease. Close to the axis of the flux tube, the contribution from the parallel component reaches almost 90% for the widest driver, compared to around 65% on the flux surface furthest from the axis. The contributions of the perpendicular and azimuthal components are below 40% for all flux surface radii and all driver widths, and both are almost universally much closer in value to each other than to the parallel component, apart from medium-and large distances from the flux tube axis for low-width drivers.

CONCLUSIONS
The aim of this study was to investigate how the width of a photospheric spiral velocity driver affected the excited MHD waves in the lower solar atmosphere. To achieve this, velocity profiles with a range of different widths between 0.15 Mm (1.67 FWHM) and 0.35 Mm (3.89 FWHM) were implemented to excite perturbations in a localised magnetic flux tube similar to one that might be found above a magnetic bright point (MBP). The resulting perturbations were decomposed into parallel, perpendicular and azimuthal components and projected on to flux surfaces, and the corresponding wave energy fluxes calculated. The relative contributions to the wave energy flux from these components were compared and evaluated.
First, these simulations do not include the transition region, where in reality waves might be reflected. Whether or not such reflection takes place, and to what extent, will clearly have an impact on the amount of energy transmitted through the transition region into the corona. Additionally, given the rapid expansion of the flux tube, it is likely that flux tubes which expand more gradually would display slightly different behaviour. Slower expansion would lead to greater magnetic field strength near the top of the domain, assuming that other variables were kept the same. This would change the plasma-β and the fast, slow and Alfven speeds, all of which would affect the propagation of waves and thus may have implications for the amount of energy transfered to the corona. The perpendicular component was found to have a minimal contribution for each spiral width, particularly on flux surfaces further from the centre of the domain. Its contribution was greatest for the innermost flux surface and the narrowest driver, but even here it is quite small (< 30%). The azimuthal component behaves similarly, in that its contribution decreases for wider drivers, but unlike the perpendicular component, its contribution is greatest close to the centre of the domain. Both components vary least on the largest flux surface.
The parallel component, on the other hand, has a significant flux contribution (> 50%) for all drivers and all flux surfaces. This contribution is greatest close to the flux tube axis and increases with driver width, reaching ∼ 90% for the widest driver.
The effective excitation of the parallel wave component by these drivers is an important result, since we have shown it indicates the presence of a fast sausage mode. It has been shown that this mode is ubiquitous in the quiet Sun and may carry enough energy to meet heating requirements in the chromosphere and low corona (Morton et al. 2012). Our results present a mechanism by which such waves could be excited by photospheric spiral velocity swirls consistent with observations.