Magnetohydrodynamic waves in braided magnetic fields

TL;DR

This study uses MHD simulations to analyze how complex braided magnetic fields influence Alfvén wave phase mixing and energy dissipation, suggesting enhanced heating mechanisms in the solar corona due to magnetic complexity.

Contribution

It demonstrates that magnetic field complexity significantly affects wave phase mixing, leading to increased energy dissipation over larger areas compared to classical models.

Findings

Wave energy is deposited over larger cross-sections in complex fields.

Small gradients in wave drivers can produce large gradients in plasma.

Wave dissipation rate increases with magnetic field complexity.

Abstract

We consider a series of MHD simulations in which a small amplitude, transverse velocity perturbation is introduced into a complex magnetic field. We analysed the deformation of the wave fronts as the perturbation propagates through the braided magnetic structures and explore the nature of Alfv\'enic wave phase mixing in this regime. Spatial gradients in the local Alfv\'en speed and variations in the length of magnetic field lines ensure that small scales form throughout the propagating wave front due to phase mixing. Additionally, the presence of complex, intricate current sheets associated with the background field locally modifies the polarisation of the wave front. The combination of these two effects enhances the rate of viscous dissipation, particularly in more complex field configurations. Unlike in classical phase mixing configurations, the greater spatial extent of Alfv\'en speed gradients ensures that wave energy is deposited over a larger cross-section of the magnetic structure. Further, the complexity of the background magnetic field ensures that small gradients in a wave driver can map to large gradients within the coronal plasma. The phase mixing of MHD waves in a complex magnetic field will progress throughout the braided volume. As a result, in a non-ideal regime wave energy will be dissipated over a greater cross-section than in classical phase mixing models. The formation rate of small spatial scales in a propagating wave front is a function of the complexity of the background magnetic field. As such, if the coronal field is sufficiently complex it remains plausible that phase mixing induced wave heating can contribute to maintaining observed temperatures. Furthermore, the weak compressibility of the transverse wave and the observed phase mixing pattern may provide seismological information about the nature of the background plasma.

Tables

Table 1. Table 1 : Mean, minimum and maximum values of important plasma and magnetic field parameters for the simulations discussed within this study. The quantities relate to the state of the plasma following the numerical relaxation and prior to the introduction of the wave front.

Simulation	Temperature (MK)			Field Strength (G)			Density (${\rho }_0$)			Plasma-$\beta$
	Min.	Mean	Max.	Min.	Mean	Max.	Min.	Mean	Max.	Min.	Mean	Max.
Uniform	1.9	1.9	1.9	10.0	10.0	10.0	1.0	1.0	1.0	0.13	0.13	0.13
s1	1.9	2.0	2.4	9.7	10.0	10.4	0.9	1.0	1.1	0.12	0.14	0.16
s2	1.9	2.1	8.7	9.5	10.0	10.8	0.3	1.0	1.3	0.13	0.14	0.23
s3	1.9	2.2	8.3	9.3	10.0	11.0	0.4	1.0	1.3	0.12	0.15	0.24
s5	1.9	2.4	8.3	8.9	10.0	11.7	0.4	1.0	1.4	0.12	0.16	0.27
t5 (High Res.)	1.9	2.4	10.3	9.6	10.0	10.7	0.3	1.0	1.4	0.12	0.17	0.24

Equations

$$\frac{\text{D}\rho }{\text{D}t}=-\rho \text{@vec}\nabla \cdot \text{@vec}v,$$

$$\rho \frac{\text{D}\text{@vec}v}{\text{D}t}=\text{@vec}j\times \text{@vec}B-\text{@vec}\nabla P+\text{@vec}{F}_{\text{visc.}},$$

$$\rho \frac{\text{D}ϵ}{\text{D}t}=\eta |\text{@vec}j{|}^2-P(\text{@vec}\nabla \cdot \text{@vec}v)+Q_{\text{visc.}},$$

$$\frac{\text{D}\text{@vec}B}{\text{D}t}=\left(\text{@vec}B\cdot \text{@vec}\nabla \right)\text{@vec}v-\left(\text{@vec}\nabla \cdot \text{@vec}v\right)\text{@vec}B-\text{@vec}\nabla \times \left(\eta \text{@vec}\nabla \times \text{@vec}B\right),$$

$$\Omega (x,y)={\int }_{z_{\text{min}}}^{z_{\text{max}}}\frac{\mathrm{d}s}{v_A},$$

$$I(z)={\int }_A|\nabla v_A|\mathrm{d}x\mathrm{d}y,$$

$$v_{y}=\cases{v}_{0}\sin\left(\frac{2\pi t}{\tau}\right)\hskip 42.67912pt&\text{if }t\leq\tau,\\ 0\hskip 42.67912pt\text{if }t>\tau.{}$$

$$K_x={\int }_V\rho v_x^2\mathrm{d}V,$$

$$v_{y}=\cases{v}_{0}\sin\left(\frac{2\pi t}{\tau}\right)e^{-\left({}^{r}/_{l}\right)^{2}}\hskip 42.67912pt&\text{if }t\leq\tau,\\ 0\hskip 42.67912pt\text{if }t>\tau,{}$$

Full text

*Methods. *We performed a series of three dimensional MHD simulations in which a small amplitude, transverse velocity perturbation is introduced into a complex magnetic field. We analysed the deformation of the wave fronts as the perturbation propagates through the braided magnetic structures and explore the nature of Alfvénic wave phase mixing in this regime. We considered the effects of viscous dissipation in a weakly non-ideal plasma and evaluate the effects of field complexity on wave energy dissipation.

*Results. *Spatial gradients in the local Alfvén speed and variations in the length of magnetic field lines ensure that small scales form throughout the propagating wave front due to phase mixing. Additionally, the presence of complex, intricate current sheets associated with the background field locally modifies the polarisation of the wave front. The combination of these two effects enhances the rate of viscous dissipation, particularly in more complex field configurations. Unlike in classical phase mixing configurations, the greater spatial extent of Alfvén speed gradients ensures that wave energy is deposited over a larger cross-section of the magnetic structure. Further, the complexity of the background magnetic field ensures that small gradients in a wave driver can map to large gradients within the coronal plasma.

*Conclusions. *The phase mixing of transverse MHD waves in a complex magnetic field will progress throughout the braided volume. As a result, in a non-ideal regime wave energy will be dissipated over a greater cross-section than in classical phase mixing models. The formation rate of small spatial scales in a propagating wave front is a function of the complexity of the background magnetic field. As such, if the coronal field is sufficiently complex it remains plausible that phase mixing induced wave heating can contribute to maintaining the observed temperatures. Furthermore, the weak compressibility of the transverse wave and the observed phase mixing pattern may provide seismological information about the nature of the background plasma.

1 Introduction

Over recent years, high spatial and temporal resolution observations have allowed many authors to identify the existence of a wide variety of magnetohydrodynamic (MHD) waves throughout the solar corona (e.g. Aschwanden et al., 1999, 2002; Verwichte et al., 2010; Threlfall et al., 2013; Duckenfield et al., 2018). Of particular interest to the current study is the evidence of waves propagating along coronal structures (e.g. De Moortel et al., 2000; McEwan & De Moortel, 2006; Okamoto et al., 2007; Thurgood et al., 2014; Morton et al., 2015). Whilst estimates of the energy associated with these waves are not well constrained (see, for example, Tomczyk et al., 2007; McIntosh et al., 2011), it is hypothesised that they may contribute to the heating of the coronal plasma and/or the acceleration of the fast solar wind. We refer interested readers to reviews by Erdélyi & Ballai (2007); Parnell & De Moortel (2012); Arregui (2015), and references therein.

Whilst the shuffling of magnetic foot points by photospheric motions may be the source of propagating MHD waves (e.g. Cranmer & van Ballegooijen, 2005; Engvold, 2008; Wang et al., 2009; Hillier et al., 2013), flows at the solar surface may also cause the slow braiding of the coronal magnetic field (Parker, 1972). This may lead to the formation of coronal current sheets and the heating of plasma through Ohmic dissipation and magnetic reconnection, even in a weakly non-ideal regime (e.g. Peter et al., 2004; Klimchuk, 2006; Wilmot-Smith et al., 2011; Reale et al., 2016; O’Hara & De Moortel, 2016). It is hypothesised that the coronal magnetic field may exist in a permanently stressed state and, as such, the background structures that are perturbed by MHD waves may have a complex and highly structured topology (e.g. Longcope & Sudan, 1994; Longbottom et al., 1998; Wilmot-Smith, 2015; Pontin et al., 2016). Indeed, whilst the small scale nature of the coronal field is not well-constrained, many observational studies have highlighted the presence of large, stressed magnetic structures and posited their association with impulsive energy release (e.g. Srivastava et al., 2010; Yan et al., 2014; Joshi et al., 2015; Lim et al., 2016).

Many previous authors have modelled the propagation of waves through a variety of coronal-like plasmas containing a range of magnetic field configurations. These include considerations of wave behaviour around magnetic null points (e.g. McLaughlin, 2016; McLaughlin et al., 2019; Prokopyszyn et al., 2019), wave interactions with large scale fields (e.g. Ofman & Thompson, 2002; Afanasyev & Zhukov, 2018), and oscillations in randomly structured plasmas (e.g. Pascoe et al., 2011; Yuan et al., 2015, 2016; Magyar et al., 2017). Unlike in uniform media, the fast and Alfvén modes are not decoupled in inhomogeneous plasmas and waves with mixed properties are able to propagate (e.g. Goossens et al., 2009, 2013). For example, in cylindrical geometries, there is an extensive body of literature investigating the kink mode and its resonant coupling to azimuthal Alfvén waves (Ionson, 1978; Goossens & Ruderman, 2011). In a more general setting, Lazzaro et al. (2000), present a method for investigating the propagation of MHD waves through finely structured plasmas.

The propagation of waves through a plasma with a transverse gradient in the Alfvén speed can lead to phase mixing (Heyvaerts & Priest, 1983) as wave fronts on adjacent magnetic field lines travel at different velocities. As such, large spatial gradients in the perturbed magnetic and velocity fields can form and will enhance the rate of wave heating dissipation even in the weakly dissipative coronal plasma (e.g. Pagano & De Moortel, 2017). However, in Cargill et al. (2016), the authors highlighted concerns of the suitability of the classical phase mixing regime for heating the corona. In particular, they argued that the typically assumed density profile cannot be self-consistently sustained without an additional background heating term. Additionally, they show that the rate of wave energy dissipation is insufficient unless transport coefficients well in excess of expected coronal values are implemented.

In the classical phase mixing regime, the dissipation time scale is expected to be proportional to $R_{t}^{{}^{1}/_{3}}$, where $R_{t}$ is a combination of the fluid and magnetic Reynolds numbers (Heyvaerts & Priest, 1983). However, an analytic investigation of Alfvén waves in more complex magnetic field structures presented in Similon & Sudan (1989), demonstrates that in some geometries, the dissipation time may be proportional to $\log{R_{M}}$. This may provide a significant enhancement above the classical heating rate and suggests phase mixing may be an important wave energy dissipation mechanism if the coronal field is sufficiently complex.

In this article, we present the results of numerical models investigating the propagation of MHD waves through a braided magnetic field. We investigate the effects of the field complexity, viscous and Ohmic dissipation, and numerical resolution. In Section 2, we outline our numerical model and describe the initial conditions of the simulations. In Section 3 we present our results and finally, in Section 4, we present a discussion of the implications of this study.

2 Numerical method

2.1 Initial conditions

For the numerical simulations presented within this article, we seek initial conditions consisting of complex magnetic field topologies. In order to obtain such fields, we have used the output of a 3D MHD model presented by Reid et al. (2018). Whilst the results of the simulations presented in the aforementioned paper are not directly related to the current work, we present a brief description of the model in order to describe the initial conditions of our experiments.

In Reid et al. (2018), an initially uniform, three dimensional, straight magnetic field is stressed by continuous, counter-rotational motions imposed on both foot points of three cylindrical threads. The configuration of the velocity driver is shown in Fig. 1. These motions generate columns of twisted field, which, at least initially, remain distinct. The central thread is twisted at a faster rate than the other two and hence, reaches the threshold for kink instability sooner. The onset of the instability generates complex current structures and destabilises the remaining two magnetic threads.

Ultimately, the constant injection of Poynting flux by the continuous rotational driving generates distributed current sheets and a complex field structure throughout the domain. In a non-ideal regime, viscous and Ohmic heating increases the plasma temperature during the course of the simulation. We direct the interested reader to the analysis presented in Reid et al. (2018) for further details.

In the present study, we used four simulation times from the model described above in order to obtain fields with varying field complexity and different quantities of magnetic twist. In each case we relaxed the field structure towards a numerical equilibrium under the action of a large viscosity. This causes further plasma heating but allows the wave dynamics (see Results) to be analysed separately from the velocities associated with the previous rotational driving. The numerical relaxation continues until the velocities are small in comparison to the amplitude of the wave driver (see below).

As in Reid et al. (2018), for the numerical simulations presented within this paper, we have implemented the Lagrangian-remap code, Lare3D (Arber et al., 2001). We advance the full, 3D MHD equations in normalised form given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where all variables have their usual meanings. The term $\@vec{F}_{\text{visc.}}$ is the combination of viscous forces and $Q_{\text{visc.}}$ is the associated heating. These dissipative terms can be separated into contributions from a real viscosity, $\nu$, and two small shock viscosity terms which are included within the following simulations to ensure numerical stability. Unless otherwise stated $\nu$ is set to zero (see section 3.2 for discussion of non-zero viscosity).

Our numerical domain has dimensions of 30 Mm $\times$ 30 Mm $\times$ 100 Mm (in $x$, $y$ and $z$, respectively) and in the majority of the simulations presented below, we implemented a numerical grid of 128 $\times$ 128 $\times$ 512. We also considered a higher resolution case with 256 $\times$ 256 $\times$ 1024 grid cells. The equilibrium obtained during the relaxation phase is resolution dependent and narrower current sheets are present when the more refined grid is used (see below for more details).

The model presented in Reid et al. (2018) is driven continuously and the field complexity slowly increases during the course of the simulations. In order to attain various degrees of field complexity, we use the simulation states at 100, 200, 300 and 500 Alfvén times as the initial conditions for the numerical relaxation. Hereafter, we refer to these simulations as s1, s2, s3 and s5 (t5 for the high resolution simulation), respectively. In each case, the rotational driving is ceased in order to allow the field to relax.

In Figs. 2 and 3, we show the configuration of the magnetic field for various cases following the numerical relaxation. In Fig. 2, magnetic field lines are traced from the foot points of each cylinder for the least (s1; left-hand panel) and most (s5; right-hand panel) twisted simulations. The colour of each field line identifies the cylinder containing the lower foot point of the field line. In Fig. 3, we display the projections of magnetic field lines onto the lower boundary of the domain for each of the four post-relaxation states. In this case, the colour of the field line is simply for visualisation purposes. In the least braided simulation (upper left panel), we see that the central and right-hand flux tube have both become unstable and have merged to form a larger structure. However, the left-hand flux tube has remained distinct and we simply observe the helical nature of the constituent field lines. For the field configurations that have been stressed further (see lower row of Fig. 3 and right-hand panel of Fig. 2), the three cylinders are all inter-connected and no separate structures remain.

In all cases, the resultant plasma is highly inhomogeneous with spatial gradients in the Alfvén speed, gas and magnetic pressures, density and temperature. The mean, minimum and maximum values of plasma parameters in the relaxed states are presented in Table 1. Due to the plasma heating described in Reid et al. (2018), and because of the viscous heating that occurs during the relaxation phase, the more complex field simulations exhibit an increase in the mean temperature. We note that the very high temperatures observed in the more complex field simulations (s2-s5 and t5) are typically confined to a few grid points. These locations have been heated by the most energetic current sheets and are not representative of the conditions throughout the remainder of the domain. Despite the inhomogeneity present within much of the numerical grid, for large values of $|y|$, the field remains relatively simple (see Fig. 3) as the plasma was not exposed to rotational driving in these locations (see Fig. 1).

In each case $\beta<1$, and the relaxed field is approximately force-free. Since $\@vec{j}\times\@vec{B}\approx\@vec{0}$, the currents are dominated by the component parallel to $\@vec{B}$. As the relaxed, approximately force-free field is predominantly aligned with the $z$ axis, we have $j_{x}\sim j_{y}\ll j_{z}$. This property of the background currents is important for explaining the modification in wave polarisation discussed in Section 3.1. An isosurface of the initial current associated with the most complex field (s5) is displayed in the left-hand panel of Fig. 4.

In the right-hand panel of Fig. 4, we display an isosurface of the Alfvén speed for the same simulation. The viscous and Ohmic heating generated by the rotational driving is largest close to the centre of the domain. As such, the temperature increases most rapidly here. Therefore, in order to maintain pressure balance, density is expelled (predominantly along magnetic field lines) from the centre of the domain and thus the Alfvén speed is higher here.

Despite this, the Alfvén travel time is not necessarily shorter in the centre of the $x$-$y$-plane, as here, field lines are typically more twisted, and hence, longer. In Fig. 5, we display the travel time for a wave propagating at the local Alfvén speed along magnetic field lines traced from the upper $z$ boundary. In particular we calculate

[TABLE]

where $s$ is an infinitesimal field line element. The spatial scales in the travel time, $\Omega$, provide an estimate for the amount of phase mixing that can be expected as a wave front propagates through the complex field structures.

In Fig. 6, we show an additional measure of the small scales within the background field. For each of the four initial conditions (s1, s2, s3 and s5), at each height, $z$, we integrate the magnitude of gradients in the Alfvén speed over the $x$-$y$-plane. In each case, we caclulate

[TABLE]

and normalise by the maximum of $I(z)$ in the s5 simulation. This quantity acts as a proxy for the rate of deformation (phase mixing) in the wave front as it propagates along magnetic field lines. In Fig. 6, we see that this measure predicts more significant phase mixing in the simulations with more complex field structures.

In addition to the initial states described above, for the purposes of comparison, we will also consider the behaviour of waves in a uniform domain with a magnetic field aligned with the vertical ($z$) axis. In this case, the initial conditions correspond to the plasma prior to any rotational driving is applied (see Reid et al., 2018, or Table 1 for details of plasma parameters).

2.2 Boundary conditions

Following the numerical relaxation, a time-dependent velocity driver was imposed on the lower $z$ boundary in order to introduce a transverse wave into the domain. Hereafter, $t=0$ refers to the start time of this driving in each simulation. In subsequent sections, we conduct an analysis on different wave drivers (detailed below), however, in the most simple case, the imposed velocity has the form $\@vec{v}=(0,v_{y}(t),0)$, where

[TABLE]

Here, $v_{0}$ is the amplitude of the driver and $\tau$ is the period of the driver. For the purpose of these simulations, we set $v_{0}\approx 20\text{ km s}^{-1}$ and $\tau\approx 30$ s. We note that the velocity amplitude is much smaller than the local Alfvén speed ($560\text{ km s}^{-1}<v_{A}<1200\text{ km s}^{-1}$ with a mean of approximately 700 $\text{km s}^{-1}$) and is similar in magnitude to wave amplitudes observed in the solar atmosphere (e.g. McIntosh et al., 2011; Morton et al., 2015). Additionally, we observe that the selected wave period is much shorter than the frequently studied 3-5 minute period which could be excited by p-modes at the photospheric surface (see Tomczyk & McIntosh, 2009; Morton et al., 2016, for example). Such a short period was selected in order to ensure that the height of our domain (100 Mm) will contain several wavelengths of the propagating pulse. Furthermore, we highlight that, in the current form, only a single pulse is introduced into the domain. Numerical simulations implementing continuous wave driving will be the focus of subsequent work.

In Fig. 7, we display the energy injected by the wave driver described by equation 7 (solid curves) for different initial conditions (s1-s5 and the uniform field case). Since there are no flows through the boundary, this corresponds to the time integral of the Poynting flux. In these simulations, the mean Poynting flux of energy into the domain (during the 30 s wave driving period) is around 50 $\text{W m}^{-2}$. Even if all of the injected energy is dissipated, this is approximately an order of magnitude below the expected energy requirements in the quiet Sun (Withbroe & Noyes, 1977). Indeed, in each case, the total energy increase is much less than 0.1% of the initial magnetic and thermal energies.

The small difference between each of the solid curves arises due to the differences in the background field complexities and the weak non-linearity of the wave driver. In particular, for the same velocity driver, a greater Poynting flux is injected for a more complex magnetic field. The dashed red curve corresponds to a modified driver which is described in Section 3.5.

In all simulations, the $x$ and $y$ boundaries are periodic and field lines are able to cross through the corresponding faces (see, for example, red and purple field lines in the right-hand panel of Fig. 2). For the driven $z$ boundaries, the horizontal velocity components are described above, flows through the boundary are not permitted ($v_{z}=0$) and all remaining variables have zero gradients across the boundaries.

3 Results

3.1 Wave dynamics

In this section, we analyse the behaviour of the propagating wave in the simulation with the most complex field, the s5 simulation. The wave front introduced by the velocity driver reaches the upper $z$ boundary after approximately $T_{e}=145$ s. We select this time, $T_{e}$, as the end point of the simulation. The mean propagation speed coupled with the selected driver frequency, ensures that the height of the numerical domain (100 Mm) is approximately 5 times the wavelength.

Since the wave front propagates through a non-uniform plasma, it cannot be identified as a pure Alfvén wave. In particular, it has mixed properties, is mildly compressible and exhibits increasingly complex behaviour as the simulation progresses. In Fig. 8, we display isosurfaces of the magnitude of $v_{y}$ at $t=0.3\,T_{e}$ (left-hand panel) and $t=0.8\,T_{e}$ (right-hand panel). Initially, the wave front is spatially uniform and polarised in the $y$ direction. However, the complex nature of the plasma generates fine spatial variation and transfers energy to the other velocity components. We now describe the two key processes that are associated with the formation of these small spatial scales.

Firstly, the non-uniform Alfvén speed contributes to phase mixing as the wave front propagates at different speeds on neighbouring field lines. Consequently, as time progresses, small length scales develop across the magnetic field. Additionally, the length of field lines vary within the braided region. Therefore, since the wave propagates along magnetic structures, this enhances the rate of phase mixing and deformation of the wave packet.

In the classic case of phase mixing (as described in Heyvaerts & Priest, 1983), magnetic field lines have a constant Alfvén speed along their length, ensuring that a wave front will always become increasingly out of phase as time progresses. However, in the current case, the Alfvén speed varies along field lines and thus phase mixing does not proceed as outlined in Heyvaerts & Priest (1983). Indeed, it is possible for a phase mixed Alfvén wave to become more in phase if it encounters an appropriate Alfvén speed profile. This continuous phase mixing and reverse phase mixing ensures that the wave front does not become as distorted as in the classic case (see the right-hand panel of Fig. 5 in Pascoe et al., 2010, for example). However, in the current paradigm, the spatial distribution of gradients in the Alfvén speed ensures that the phase mixing is spread across a large cross-sectional area of the wave front.

Secondly, the rotational nature of the driver imposed in Reid et al. (2018), ensures that, at $t=0$, in each non-uniform simulation the currents are predominantly aligned in the $z$ direction (see previous section). Meanwhile, the velocity driver induces a perturbation of the $y$ component of the magnetic field. The interaction between this perturbation and the $z$ component of the background currents generates a Lorentz force acting parallel to the $x$ axis. This transfers energy from the initially driven component of the velocity ($v_{y}$) to the other horizontal component ($v_{x}$). This has the effect of locally modifying the polarisation of the wave front. There is also a small transfer of energy to $v_{z}$ associated with horizontal background currents. However in the phase-mixing region the mean of $\lvert v_{z}\rvert$ is approximately $2\%$ of the mean of $\lvert v_{y}\rvert$, whereas the mean of $\lvert v_{x}\rvert$ is approximately $20\%$ of the mean of $\lvert v_{y}\rvert$.

In Fig. 9, we show isosurfaces of the magnitude of $v_{x}$ at $t=0.3\,T_{e}$ and $t=0.8\,T_{e}$. We notice that at both times the spatial extent of the energy conversion is limited to the regions of most complex field (compare to the lower right-hand panel of Fig 3, for example) as these are the locations containing the most significant background currents. Further, as time progresses, energy is continuously transferred to the $x$ component of the velocity, resulting in the increase in volume within the isosurface observed between the two panels in Fig. 9.

The combination of these two, highly localised, effects results in the formation of transverse gradients in the perturbed velocity field. The associated small scales are shown in greater detail in Figs. 10 and 11. In Fig. 10, we show the profile of $v_{y}$ in a horizontal cut through the wave packet at $t=0.8\,T_{e}$. We note that there is spatial variation throughout the entire cross-section but the perturbation is more homogeneous at large $|y|$, where the magnetic field is less complex. The solid and dashed lines shown in Fig. 10 correspond to the locations of the respective lines depicted in Fig. 11. Here, we clearly see the large transverse gradients that form as a result of the complexity of the background plasma. In Fig. 11, we again highlight the relative homogeneity of the dashed line (region of less complex field) in comparison to the solid line (region of more complex field).

The formation of these small spatial scales in the velocity field can be tracked using the vorticity, ${\@vec{\omega}}=\nabla\times\@vec{v}$. The imposed wave driver injects vorticity as it introduces vertical gradients in the driven component of the velocity $\left(\dfrac{\partial v_{y}}{\partial z}\right)$, and hence contributes to an increase in the $\omega_{x}$ term. Meanwhile, the formation of small scales caused by the background field complexity contributes to the formation of transverse gradients of $v_{y}$ and to a lesser extent $v_{x}$. This results in a gradual increase in the magnitude of $\omega_{z}$.

We display this behaviour by showing the change in the volume integrated magnitude of vorticity and of each component in Fig. 12 as the simulation progresses. We see an initial increase in $\omega_{x}$ (red line) as the driver introduces a wave pulse into the domain. However, at $t\approx 29$ s, the driving motion is suspended and the magnitude of $\omega_{x}$ component slowly decreases as energy is transferred from the $y$ component of the velocity (see above). Meanwhile, there is a steady increase in the $z$ component of vorticity (green line) as phase mixing progresses. This increase is much greater than the decrease observed in $\omega_{x}$, resulting in a net increase in total vorticity (black line) during the simulation. This will produce an enhancement in the rate of viscous dissipation as we shall discuss later. We note that there is also a small increase in the volume integral of $\lvert\omega_{y}\rvert$ (blue line) as time progresses. This is predominantly associated with the localised formation of vertical gradients in $v_{x}$ (see Fig. 9).

In Fig. 13, we display an isosurface of $\lvert\@vec{\omega}\rvert$ at $t=0.3\,T_{e}$ (left-hand panel) and $t=0.8\,T_{e}$ (right-hand panel) for the s5 simulation. At both times, we see that the formation of large vorticities is confined to regions of complex field and even within this sub-section of the domain, the vorticity has a fragmented nature. This suggests that any viscous heating associated with this wave will not heat uniformally. Despite this, in contrast with many classical phase mixing simulations, small spatial scales form over a larger proportion of the cross-sectional area of the magnetic structure and not simply in a small region in the boundary of a coronal loop. Thus, in a non-ideal regime, viscous heating will deposit wave energy throughout the complex field. We will consider the rate of wave energy dissipation in a subsequent section.

At this stage, it is important to note that the wave dynamics are also associated with the development of small scales in the magnetic field. As such, for comparable values of resistivity, $\eta$, and viscosity, $\nu$, Ohmic and viscous heating are expected to be equally effective as mechanisms for dissipating wave energy. However, since the background field is already finely structured, it is difficult to separate the wave-induced small scales in the magnetic field from the pre-existing gradients. Therefore, for simplicity, we have restricted our analysis to the evolution of small scales in the perturbed velocity field.

3.2 Field complexity

In order to quantify the effects of varying the field complexity, we now contrast the observed wave dynamics in the s1, s2, s3 and s5 simulations. In Figs. 14 and 15, we show isosurfaces of $\lvert v_{y}\rvert$ and $\lvert v_{x}\rvert$, respectively at $t=0.8\,T_{e}$ for the s2 (left-hand panels) and s3 (right-hand panels) simulations. These correspond to Figs. 8(b) and 9(b) which show the corresponding s5 results and were discussed above. Similar plots for the s1 simulation show an almost perfectly planar wave in the $v_{y}$ component, and very little transfer of energy to the $v_{x}$ component. From this collection of plots we see that the rate of phase mixing and, thus deformation of the propagating wave front, is sensitive to the complexity of the magnetic field and background Alfvén speed profile.

In particular, we observe that the isosurfaces of $\lvert v_{y}\rvert$ exhibit greater variation in the more complex field simulations. Not only are the gradients larger (see Fig. 17 and associated discussion for more detail), but they also occur over a greater cross-section of the propagating wave front. This is a result of the expansion of the braided field region between the simpler and more complex cases (see Fig. 3 for example). Additionally, the transfer of energy to the $v_{x}$ component is more efficient in the s3 and s5 simulations than in the simpler cases. Again, the increased cross-sectional area of stressed field translates to the formation of a significant $v_{x}$ component, and hence modification of wave polarisation, over a larger proportion of the domain.

In Fig. 16, we aim to quantify the rate of energy transfer from the $y$ component of the velocity to the $x$ component. We calculate

[TABLE]

the $x$ component of the kinetic energy integrated over the volume of the domain, $V$. We show the time evolution of this quantity for each of the four complex field configurations. Meanwhile, the uniform field does not induce any transfer of energy between the transverse velocity components. In Fig. 16, we have normalised by the total kinetic energy injected into the domain by the wave driver. As such, we see that only a relatively small (of the order 4% for the most complex field simulation) proportion of the energy is transfered between the $v_{y}$ and $v_{x}$ components. However, this is a highly localised process and in some regions of the wave front, the magnitude of $v_{x}$ is comparable to the size of $v_{y}$. The observational significance of this result will be examined in subsequent work.

We discussed the cause of this energy transfer previously and will simply note that the larger the background currents, the more readily the wave polarisation is modified. Additionally, we note that each simulation will induce small perturbations in the field aligned component due to the non-linear ponderomotive force, which will contribute to the $x$ component of the kinetic energy in regions of non-zero $B_{x}$. Since the strength of the perturbed magnetic field is small in comparison to the background field, we expect this to be a relatively weak effect.

In Fig. 17, we compare the phase mixing gradients that form during the course of the s1 (blue) and s5 (red) simulations. We show $\dfrac{dv_{y}}{dx}$ and $\dfrac{dv_{y}}{dy}$ along horizontal lines through the wave front in the left-hand and right-hand panels, respectively. The term shown in the left-hand panel does contribute to the growth in the magnitude of the vorticity (discussed above). The right-hand panel, on the other hand, depicts a term that is not directly tracked by the vorticity. Instead, this contributes to the compression of the plasma and hence may be associated with an observational signature that is absent from a pure, incompressible Alfvén wave. Whilst detailed analysis of this effect remains beyond the scope of this work, we note that the plasma compression observed in these simulations is highly inhomogeneous and may provide seismological information about the nature of the background magnetic field.

In both cases, it is clear that the larger gradients (both positive and negative) form in the more complex field simulation (red curves). This is indicative of the enhanced phase mixing observed in the s5 simulation. Further, we note that the magnitude of the largest gradients in the left-hand panel are greater than those in the right-hand panel (see different scales for $y$-axes). This is because the compression of plasma and magnetic field associated with the right-hand panel generates a total pressure force which limits the formation of very small scales in the direction of the velocity perturbation.

Additionally, we highlight the relatively small gradients in the right-hand panel for large $\lvert y\rvert$. These are particularly apparent in the s1 simulation (blue curve) and are a result of the relatively simply field and Alfvén speed profile present in this region of the domain. As such we see very little phase mixing and plasma compression at large values of $\lvert y\rvert$. Thus, by comparing the (potentially observable) plasma compression profiles, we can deduce information about the relative complexity of the field at different values of $\lvert y\rvert$. The observational implications of this effect will be explored in future work.

3.3 Dissipation

In this section, we seek to quantify the significance of the previously discussed formation of small spatial scales in relation to the dissipation of wave energy. In particular, we consider the effects of non-zero transport coefficients on the evolution of the simulations. As a non-zero resistivity, $\eta$, will cause the dissipation of energy in the non-potential background field (in addition to any wave heating), we predominantly focus on the effects of viscosity, $\nu$. With this approach, we do not expect significant energy to be extracted from the background conditions.

We repeat the previously described simulations with values of viscosity, $\nu=10^{-3},10^{-4}$ and $\nu=10^{-5}$. These correspond to Reynolds numbers of $10^{3},10^{4}$ and $10^{5}$, respectively, and act in addition to the aforementioned shock viscosities (see Section 2) which also provide dissipative effects.

In Fig. 18, we display the volume integral of the magnitude of the vorticity (also see Fig. 12) for the s5 simulation with various values of $\nu$. For comparison, we also include the results of a simulation with $\eta=10^{-4},\nu=0$ (green curve). We see that the formation of small scales in the velocity field is suppressed by increasing the frictional effects of viscosity (larger values of $\nu$). The similarity between the red and black dashed curves confirms that at $\nu=10^{-5}$, numerical dissipation and the effects of the shock viscosities have greater significance for the wave dynamics. As such, the maximum Reynolds number obtained with the current numerical resolution in these simulations is of the order $10^{4}-10^{5}$.

Counterintuitively, increasing the resistivity enhances the total vorticity observed within the domain (compare green and red curves). This is not a result of the wave dynamics but instead is caused by the diffusion and Ohmic dissipation of the background magnetic field. This process generates small, localised flows and contributes to the increase of the vorticity. The phase mixing of the wave front is still suppressed in this case, however the associated decrease in small scale formation is dominated by the effects caused by the non-zero $\eta$ acting on the background field.

In Fig. 19, we display the cumulative volume integrated viscous heating for the s1, s2, s3, s5 and uniform simulations for $\nu=10^{-3}$ (left-hand panel) and $\nu=10^{-4}$ (right-hand panel). This heating is associated with these values of $\nu$ and the shock viscosities which are not changed between simulations. In the Lare3d code, numerical dissipation does not contribute to the heating of plasma. It is important to note that due to the nature of the wave driver (low amplitude and only a single pulse), very little energy is injected into the domain. As such, we cannot expect significant plasma heating to occur in this case (see Sect. 4 for further discussion). However, continuous driving may provide sufficient energy in wave heating models (e.g. Karampelas & Van Doorsselaere, 2018) and this possibility will be investigated in future studies.

In the left-hand panel of Fig. 19, we see that the viscous disspation of wave energy is a little more than twice as efficient in the s5 simulation (red) line than in the uniform field case (black line). Indeed, we see that increasing the field complexity gradually increases the total viscous heating. Despite this, given that we expect very little heating from a planar Alfvén wave in the high Reynolds number corona, the increase in energy dissipation due to phase mixing seems insufficient to generate significant wave heating. However, in the right-hand panel we see that the relative increase of viscous heating between the uniform and most braided field cases is enhanced at lower values of $\nu$. This is because the formation of small scales due to phase mixing is more efficient at lower viscosities for the s5 simulation (no phase mixing occurs in the uniform field case).

By comparing the two panels of Fig. 19, we see that the viscous heating does not simply scale with the size of $\nu$. In particular, the energy dissipation depicted in the right-hand panel is not an order of magnitude smaller than in the left-hand panel. There are two main reasons for this. Firstly, the formation of small scales in the velocity field is more efficient in the higher Reynolds number case (right-hand panel) and thus the viscosity is acting on larger spatial gradients. Secondly, the magnitudes of the shock viscosities are unchanged between the two sets of simulations and thus the effects of these terms are more significant when the real viscosity is reduced.

3.4 Resolution considerations

In the model presented within Reid et al. (2018), the widths of the current sheets that form are a function of the spatial resolution implemented in the numerical code. In particular, higher resolution simulations allow the formation of narrower current sheets and smaller scale inhomogeneties in the plasma as the numerical dissipation is reduced. In this section we consider a new initial state, t5, which is obtained using a domain with $256\times 256\times 1024$ grid cells. In this case, the same uniform field is driven with identical rotational motions for 500 Alfvén times. The plasma is then allowed to relax numerically under the effects of a high viscosity as with all previously discussed simulations. For comparison with the lower resolution experiments, important characteristics of the initial state are included in Table 1. Once a numerical equilibrium is obtained, an identical, single pulse sinusoidal wave is intoduced at the $z=-50$ Mm (lower) boundary.

In the higher resolution simulation, the fundamental wave dynamics remain unchanged, however, the increased intricacy of the initial conditions induce enhanced phase mixing and allow even smaller spatial scales to form. In Fig. 20, we quantify this effect by plotting the time evolution of the volume integral of the magnitude of the vorticity for the s5 (solid red), t5 (dashed red) and the uniform field (black) simulations. As in Fig. 12, we observe an initial rise in the vorticity in all cases as a result of the gradients introduced by the wave driver. Subsequently, the progression of phase mixing in the complex field simulations (but not in the uniform case) enhances the volume integrated vorticity.

The vorticities in the uniform field case represent the vertical gradients of $v_{y}$ associated with the Alfvén wave. As such, they are inversely proportional to the wavelength of the perturbation. The red curves are a combination of these gradients and the small scales associated with phase mixing. The increase observed in the dashed red line in comparison to the solid red line is simply due to the enhanced phase mixing observed in the high resolution simulation.

The increased field line and Alfvén speed complexity causes a greater deformation in the propagating wave front. Additionally, since the effects of numerical dissipation are reduced in the t5 simulation, larger gradients in the velocity field (and hence greater vorticities) are accessible. As a result of the formation of smaller scales in the t5 experiment (in comparison to the s5 case), we expect an additional increase in the efficiency of wave energy dissipation compared to the uniform field case in non-ideal simulations.

3.5 Non-uniform wave driver

In this section, we consider the effects of imposing a transverse wave driver that is not spatially uniform across the lower boundary of the domain. As waves propagate along magnetic structures, the complexity of the field ensures that small gradients at the foot points of field lines can map to large gradients within the numerical domain. It is reasonable to expect that two field lines which are well-separated within the lower solar atmosphere may be excited by different wave drivers. If these field lines converge towards each other at higher altitudes, then as the wave packets propagate into the corona, the transverse gradients will be enhanced.

We consider a modified driver of the form

[TABLE]

where $r^{2}=x^{2}+y^{2}$. We set $l=6$ Mm to ensure that the imposed velocity is approximately [math] close to the $x$ and $y$ boundaries of the domain. The form of the imposed wave driver is shown in Fig. 21 and the total Poynting flux through the lower boundary of the domain is shown in Fig. 7 (dashed red line). It is imposed on two initial conditions; the s5 field and, for comparison, the straight field case.

In Fig. 22, we show the transverse velocity profile in a horiztontal cross-section through the wave front at $t=0.8\,T_{e}$. We show the magnitude of the horizontal velocity perturbation (contour plot) and the direction of these flows (overplotted vectors). We note that although the wave driver has the same maximum amplitude as previously, the flows are smaller in this case. This is because the perturbed plasma has to displace surrounding fluid that is not moved by the imposed velocity driver. This transfers wave energy from the centre of the domain to the surrounding plasma. Additionally, there is some evidence of a return flow being generated at large $\lvert x\rvert$ as the non-driven plasma moves into the space vacated by the perturbed fluid.

The small scales observed in Fig. 22 are once again indicative of the phase mixing that progresses as the wave propagates along magnetic field lines. However, since the wave energy injected by the driver is smaller in this case (compare the dashed and solid red lines in Fig. 7), the total vorticity is also reduced. Therefore, in order to compare the growth in small scales due to phase mixing for the two wave drivers, we normalise the vorticity by the maximum value obtained in the corresponding straight field simulations.

In Fig. 23, we show the magnitude of the vorticity, integrated over the computational domain for the two forms of the driver. We show the uniform driver acting on the straight field and s5 initial conditions (red and black lines, respectively). For both of these curves, we normalise by the maximum of the red line. Additionally, we show the non-uniform driver on these two initial conditions (dashed and solid blue lines, respectively). For these two curves, we have normalised by the maximum of the dashed blue line.

In Fig. 23, we see that once the reduction in the injected vorticity (for the non-uniform driver) is accounted for, the integral of $\lvert\omega\rvert$ increases by a greater amount in the non-uniform cases. This is largely due to the phenomenon described above by which complex magnetic fields can ensure small gradients in the wave driver can map to large gradients in the numerical domain. As such, in a non-ideal regime, wave energy dissipation will be significantly more efficient for a wave driver that has (possibly small) spatial gradients.

4 Discussion and conclusions

In this article we have presented the results of MHD simulations of transverse wave packets propagating through complex magnetic field configurations. A single sinusoidal pulse was driven on the lower boundary of the domain for initial equilibria containing varying field complexity.

The presence of an approximately helical component in the magnetic field induces a localised change in the polarisation of the wave front, which is sensitive to the complexity of the initial currents. Additionally, the inhomogeneous Alfvén speed profile associated with the background plasma and magnetic field generates a complex phase mixing pattern. Enhancing the field complexity increases the average rate of phase mixing and the formation of small scales throughout the wave front.

As a result of these effects, large spatial gradients are induced in both the perturbed velocity and magnetic fields. However, the pre-existing field complexity ensures it is relatively difficult to identify the formation of additional small scales in the magnetic field as the wave propagates through the numerical domain. On the other hand, the background flows are very small (in comparison to the driver amplitude) and thus we have focussed our analysis on the evolution of the velocity field. In particular, we highlight the evolution of the vorticity as a useful measure for tracking the generation of small scales.

The single pulse examined in the current study contains very little energy and thus, even in non-ideal regimes, does not lead to significant plasma heating. The Poynting flux could certainly be enhanced by adopting a larger amplitude driver or by increasing the background field strength. Alternatively, a continuous driver will inject substantially more energy into the domain and the potential excitation of field line resonances can further increase the energy input. Since the field exhibits a wide range of natural field line frequencies, it is relatively easy to excite resonances regardless of the driving frequency. The additional energy input and the potential for even more complex dynamics ensures wave heating is more significant in a continuously driven case. A more thorough investigation of such simulations will be considered in subsequent work.

We have discussed the enhanced efficiency of kinetic energy dissipation in non-ideal regimes within the complex field structures. In a relatively low Reynolds number plasma, the increase in heating is limited in comparison to the dissipation associated with an Alfvén wave propagating through a uniform field. However, at higher Reynolds numbers, the enhancement is greater (compared to an analagous uniform field case) due to the increased effects of phase mixing for reduced viscosities. Additionally, the implementation of a non-uniform driver generates further small scales due to the convergence of initially distant magnetic field lines within the corona. As such, in the presence of continuous wave driving, this model may allow significant heat to be deposited within the solar atmosphere. At present, this hypothesis remains untested and will be the focus of future research.

Many classical phase mixing models rely on a prescribed density profile which are not self-consistently generated by the heating process. Indeed, Cargill et al. (2016) argued that even with transport coefficients augmented by many orders of magnitude above expected coronal values, the heating profile established by phase mixing will not be able to sustain the dense loops observed within the solar atmosphere. Since, in this model, the spatial distribution of energy dissipation is not confined to regions of strong density gradient, it seems that the arguments presented by Cargill et al. (2016) do not preclude phase mixing as a coronal heating mechanism in this regime. It is important to highlight that, in these simulations, wave energy is dissipated across the entire cross-section of the braided magnetic structure and not simply in narrow layers within the flux tubes.

Since the transverse wave propagates through a highly inhomogeneous domain, it cannot be identied as a pure Alfvén wave. Indeed, in this regime the wave is weakly compressible and, as such, given sufficient telescope sensitivity, may be directly observable in the solar corona. Subsequent work will seek to identify the necessary wave amplitude for such plasma compression to be detected with current observational capabilities.

In any case, the wave front becomes highly deformed and transfers energy to smaller scales which will be beyond the spatial resolving power of observing instruments. As such, it remains unclear whether the localised phase mixing pattern and modification of the wave polarisation have implications for estimating the coronal wave energy budget. In particular, it may be difficult to accurately estimate the energy associated with observations of similar waves in the solar atmosphere.

The compressibility of the propagating wave front is not homogeneous and is sensitive to the nature of the background field. For example, it is associated with greater density variation in regions of greater field complexity. Therefore, it is plausible that observations of the weak plasma compression as a similar wave propagates through coronal field will indicate the complexity of the background Alfvén speed and/or the nature of magnetic field lines. As such synthetic observables derived from the models presented within this article are expected to contain information about the initial conditions. It remains to be seen whether any seismological techniques can be employed on the wave dynamics to infer properties of the background medium. If so, observations of propagating transverse coronal waves may be used to infer the complexity of the coronal magnetic field.

Acknowledgements. The authors would like to thank the referee for their helpful comments. The research leading to these results has received funding from the UK Science and Technology Facilities Council (consolidated grants ST/N000609/1 and ST/S000402/1) and the European Union Horizon 2020 research and innovation programme (grant agreement No. 647214). IDM acknowledges support from the Research Council of Norway through its Centres of Excellence scheme, project number 262622 and JR acknowledges the support of the Carnegie Trust for the Universities of Scotland.