© ESO 2019

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 et al. 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_{\text{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 Sect. 2, we outline our numerical model and describe the initial conditions of the simulations. In Sect. 3 we present our results and finally, in Sect. 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.

Fig. 1. Schematic of rotational driving implemented in Reid et al. (2018).

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

$$\begin{array}{cc}\hfill \frac{\mathrm{D}\rho }{\mathrm{D}t}& =-\rho \mathbf{\nabla }·\mathit{v},\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \rho \frac{\mathrm{D}\mathit{v}}{\mathrm{D}t}& =\mathit{j}\times \mathit{B}-\mathbf{\nabla }P+{\mathit{F}}_{\mathrm{visc}.},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \rho \frac{\mathrm{D}ϵ}{\mathrm{D}t}& ={\eta |\mathit{j}|}^2-P(\mathbf{\nabla }·\mathit{v})+Q_{\mathrm{visc}.},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \frac{\mathrm{D}\mathit{B}}{\mathrm{D}t}& =(\mathit{B}·\mathbf{\nabla })\mathit{v}-(\mathbf{\nabla }·\mathit{v})\mathit{B}-\mathbf{\nabla }\times (\eta \mathbf{\nabla }\times \mathit{B}),\hfill \end{array}$$(4)

where all variables have their usual meanings. The term F_visc. is the combination of viscous forces and Q_visc. is the associated heating. These dissipative terms can be separated into contributions from a real viscosity, ν, and two small shock viscosity terms which are included within the following simulations to ensure numerical stability. Unless otherwise stated ν is set to zero (see Sect. 3.2 for discussion of non-zero viscosity).

Our numerical domain has dimensions of 30 Mm × 30 Mm × 100 Mm (in x, y and z, respectively) and in the majority of the simulations presented below, we implemented a numerical grid of 128 × 128 × 512. We also considered a higher resolution case with 256 × 256 × 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.

Fig. 2. Configuration of magnetic field lines for the least (s1; left-hand panel) and most (s5; right-hand panel) complex initial conditions.

Fig. 3. Projections of magnetic field lines onto x–y-plane for the four field configurations. Upper left: s1, upper right: s2, lower left: s3 and lower right: s5.

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 β <  1, and the relaxed field is approximately force-free. Since j × B ≈ 0, the currents are dominated by the component parallel to B. As the relaxed, approximately force-free field is predominantly aligned with the z axis, we have j_x ∼ j_y ≪ j_z. This property of the background currents is important for explaining the modification in wave polarisation discussed in Sect. 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.

Fig. 4. Initial configuration for the most complex field (s5) simulation. In the right-hand panel, the Alfvén speed is largest within the volume.

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

Fig. 5. Alfvén travel time, Ω (see Eq. (5)), along magnetic field lines traced from the upper z boundary.

$$\begin{array}{c}\hfill \mathrm{\Omega }(x,y)={{\displaystyle \int }}_{z_{\min }}^{z_{\max }}\frac{\mathrm{d}s}{v_{\mathrm{A}}},\end{array}$$(5)

where s is an infinitesimal field line element. The spatial scales in the travel time, Ω, 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

$$\begin{array}{c}\hfill I(z)={{\displaystyle \int }}_{\mathrm{A}}|\mathrm{\nabla }{v}_{\mathrm{A}}|\mathrm{d}x\mathrm{d}y,\end{array}$$(6)

Fig. 6. Measure of the field complexity, I (see Eq. (6)), in each of the four field configurations as a function of z.

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 v = (0, v_y(t),0), where

$$\begin{array}{c}\hfill v_y=\{\begin{array}{cc}{v}_{0}sin\left(\frac{2\pi t}{\tau }\right)\hfill & \mathrm{i}ft\le \tau ,\hfill \\ 0\hfill & \mathrm{i}ft>\tau .\hfill \end{array}\end{array}$$(7)

Here, v₀ is the amplitude of the driver and τ is the period of the driver. For the purpose of these simulations, we set v₀ ≈ 20 km s⁻¹ and τ ≈ 30 s. We note that the velocity amplitude is much smaller than the local Alfvén speed (560 km s⁻¹ < v_A <  1200 km s⁻¹ with a mean of approximately 700 km s⁻¹) 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 min 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 Eq. (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 W m⁻². 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.

Fig. 7. Energy injected into the domain by the wave driver. We display the time integral of the Poynting flux through the z = −50 Mm (lower) boundary for different initial conditions with the drivers described by Eqs. (7) and (9) (solid and dashed lines, respectively). Here we have normalised by the maximum of the s5 (solid red) curve.

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 Sect. 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.

Fig. 8. Isosurfaces of the magnitude of vy corresponding to |vy| ≈ 0.8 v0.

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 |v_z| is approximately 2% of the mean of |v_y|, whereas the mean of |v_x| is approximately 20% of the mean of |v_y|.

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.

Fig. 9. Isosurfaces of the magnitude of vx corresponding to |vx| ≈ 0.25 v0.

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).

Fig. 10. y component of the velocity field in a horizontal cut through the wave front at t = 0.8 Te.

Fig. 11. y component of the velocity field along the two horizontal lines shown in Fig. 10.

The formation of these small spatial scales in the velocity field can be tracked using the vorticity, ω = ∇ × v. The imposed wave driver injects vorticity as it introduces vertical gradients in the driven component of the velocity $\left(\frac{\partial v_y}{\partial z}\right)$, and hence contributes to an increase in the ω_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 ω_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 ω_x (red line) as the driver introduces a wave pulse into the domain. However, at t ≈ 29 s, the driving motion is suspended and the magnitude of ω_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 ω_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 |ω_y| (blue line) as time progresses. This is predominantly associated with the localised formation of vertical gradients in v_x (see Fig. 9).

Fig. 12. Volume integrated |ω| as a function of time for simulation, s5. Additionally, we display the volume integral of the magnitude of each component and the dashed black line corresponds to the time t = τ when the wave driving ceases.

In Fig. 13, we display an isosurface of |ω| 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.

Fig. 13. Isosurfaces of the magnitude of ω in the s5 simulation.

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, η, and viscosity, ν, 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 |v_y| and |v_x|, respectively at t = 0.8 T_e for the s2 (left-hand panels) and s3 (right-hand panels) simulations. These correspond to Figs. 8 and 9 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.

Fig. 14. Isosurfaces of the magnitude of vy corresponding to |vy| ≈ 0.8 v0 at t = 0.8 Te.

Fig. 15. Isosurfaces of the magnitude of vx corresponding to |vx| ≈ 0.25 v0 at t = 0.8 Te.

In particular, we observe that the isosurfaces of |v_y| 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

$$\begin{array}{c}\hfill K_x={{\displaystyle \int }}_V\rho v_x^2\mathrm{d}V,\end{array}$$(8)

Fig. 16. Time evolution of the volume integrated x component of the kinetic energy, Kx (see Eq. (8)), for each of the four complex field simulations.

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 $\frac{\mathrm{d}{v}_y}{\mathrm{d}x}$ and $\frac{\mathrm{d}{v}_y}{\mathrm{d}y}$ 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 |y|. 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 |y|. Thus, by comparing the (potentially observable) plasma compression profiles, we can deduce information about the relative complexity of the field at different values of |y|. 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, η, 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, ν. 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, ν = 10⁻³, 10⁻⁴ and ν = 10⁻⁵. These correspond to Reynolds numbers of 10³, 10⁴ and 10⁵, respectively, and act in addition to the aforementioned shock viscosities (see Sect. 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 ν. For comparison, we also include the results of a simulation with η = 10⁻⁴, ν = 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 ν). The similarity between the red and black dashed curves confirms that at ν = 10⁻⁵, 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⁴ − 10⁵.

Fig. 17. Transverse gradients of vy for the s1 (blue) and s5 (red) simulations. We show $\frac{\mathrm{d}{v}_y}{\mathrm{d}x}$ along y = 0 Mm (left-hand panel) and $\frac{\mathrm{d}{v}_y}{\mathrm{d}y}$ along x = 0 Mm (right-hand panel) at t = 0.8 Te through the location of the wave front. In each case we have normalised by the maximum transverse gradient of vy in the horizontal plane containing these lines.

Fig. 18. Time evolution of |ω| integrated over the numerical domain for the s5 simulation with different transport coefficients. We note that the early stages of the simulation are omitted here and we have normalised the integral using the maximum value attained in the ideal s5 simulation (as in Fig. 12).

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 η 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 ν = 10⁻³ (left-hand panel) and ν = 10⁻⁴ (right-hand panel). This heating is associated with these values of ν 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.

Fig. 19. Cumulative volume integrated viscous heating for the braided and uniform field simulations. We show the heating for ν = 10−3 in the left-hand panel and for ν = 10−4 in the right-hand panel. In each case we have normalised by the total heating obseverd in the s5 simulation with ν = 10−3.

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 ν. 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 ν. 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 × 256 × 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.

Fig. 20. Volume integrated |ω| as a function of time for the s5, t5 and uniform field simulations. In each case, we have normalised by the maximum value of the integral for the s5 simulation (see also Figs. 12 and 18).

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

$$\begin{array}{c}\hfill v_y=\{\begin{array}{cc}{v}_{0}sin\left(\frac{2\pi t}{\tau }\right)e^{-{\left({}^r{/}_l\right)}^2}\hfill & \mathrm{i}ft\le \tau ,\hfill \\ 0\hfill & \mathrm{i}ft>\tau ,\hfill \end{array}\end{array}$$(9)

where r² = x² + y². We set l = 6 Mm to ensure that the imposed velocity is approximately 0 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.

Fig. 21. Non-uniform driver profile imposed at the lower z boundary.

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 |x| as the non-driven plasma moves into the space vacated by the perturbed fluid.

Fig. 22. Transverse velocity profile in a cross-section of the wave packet at t = 0.8 Te for the non-uniform driver.

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.

Fig. 23. Volume integrated |ω| for the two forms of the wave driver in the s5 and straight field simulations. For the uniform driver cases, we have normalised by the maximum of the red curve. For the non-uniform driver cases, we have normalised by the maximum of the dashed blue curve.

In Fig. 23, we see that once the reduction in the injected vorticity (for the non-uniform driver) is accounted for, the integral of |ω| 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.

Acknowledgments

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.

References

All Tables

All Figures

	Fig. 1. Schematic of rotational driving implemented in Reid et al. (2018).
In the text

	Fig. 2. Configuration of magnetic field lines for the least (s1; left-hand panel) and most (s5; right-hand panel) complex initial conditions.
In the text

	Fig. 3. Projections of magnetic field lines onto x–y-plane for the four field configurations. Upper left: s1, upper right: s2, lower left: s3 and lower right: s5.
In the text

	Fig. 4. Initial configuration for the most complex field (s5) simulation. In the right-hand panel, the Alfvén speed is largest within the volume.
In the text

	Fig. 5. Alfvén travel time, Ω (see Eq. (5)), along magnetic field lines traced from the upper z boundary.
In the text

	Fig. 6. Measure of the field complexity, I (see Eq. (6)), in each of the four field configurations as a function of z.
In the text

	Fig. 7. Energy injected into the domain by the wave driver. We display the time integral of the Poynting flux through the z = −50 Mm (lower) boundary for different initial conditions with the drivers described by Eqs. (7) and (9) (solid and dashed lines, respectively). Here we have normalised by the maximum of the s5 (solid red) curve.
In the text

	Fig. 8. Isosurfaces of the magnitude of vy corresponding to |vy| ≈ 0.8 v0.
In the text

	Fig. 9. Isosurfaces of the magnitude of vx corresponding to |vx| ≈ 0.25 v0.
In the text

	Fig. 10. y component of the velocity field in a horizontal cut through the wave front at t = 0.8 Te.
In the text

	Fig. 11. y component of the velocity field along the two horizontal lines shown in Fig. 10.
In the text

	Fig. 12. Volume integrated |ω| as a function of time for simulation, s5. Additionally, we display the volume integral of the magnitude of each component and the dashed black line corresponds to the time t = τ when the wave driving ceases.
In the text

	Fig. 13. Isosurfaces of the magnitude of ω in the s5 simulation.
In the text

	Fig. 14. Isosurfaces of the magnitude of vy corresponding to |vy| ≈ 0.8 v0 at t = 0.8 Te.
In the text

	Fig. 15. Isosurfaces of the magnitude of vx corresponding to |vx| ≈ 0.25 v0 at t = 0.8 Te.
In the text

	Fig. 16. Time evolution of the volume integrated x component of the kinetic energy, Kx (see Eq. (8)), for each of the four complex field simulations.
In the text

Fig. 17. Transverse gradients of vy for the s1 (blue) and s5 (red) simulations. We show $\frac{\mathrm{d}{v}_y}{\mathrm{d}x}$ along y = 0 Mm (left-hand panel) and $\frac{\mathrm{d}{v}_y}{\mathrm{d}y}$ along x = 0 Mm (right-hand panel) at t = 0.8 Te through the location of the wave front. In each case we have normalised by the maximum transverse gradient of vy in the horizontal plane containing these lines.

In the text

	Fig. 18. Time evolution of |ω| integrated over the numerical domain for the s5 simulation with different transport coefficients. We note that the early stages of the simulation are omitted here and we have normalised the integral using the maximum value attained in the ideal s5 simulation (as in Fig. 12).
In the text

	Fig. 19. Cumulative volume integrated viscous heating for the braided and uniform field simulations. We show the heating for ν = 10−3 in the left-hand panel and for ν = 10−4 in the right-hand panel. In each case we have normalised by the total heating obseverd in the s5 simulation with ν = 10−3.
In the text

	Fig. 20. Volume integrated |ω| as a function of time for the s5, t5 and uniform field simulations. In each case, we have normalised by the maximum value of the integral for the s5 simulation (see also Figs. 12 and 18).
In the text

	Fig. 21. Non-uniform driver profile imposed at the lower z boundary.
In the text

	Fig. 22. Transverse velocity profile in a cross-section of the wave packet at t = 0.8 Te for the non-uniform driver.
In the text

	Fig. 23. Volume integrated |ω| for the two forms of the wave driver in the s5 and straight field simulations. For the uniform driver cases, we have normalised by the maximum of the red curve. For the non-uniform driver cases, we have normalised by the maximum of the dashed blue curve.
In the text