© ESO 2022

1. Introduction

Turbulent convection is crucial for a variety of solar processes, including the formation and decay of magnetic concentrations, magnetic flux cancellation, and so on (Nordlund et al. 2009; Stein 2012). Vortices are a major component of any turbulent flow and the Sun’s highly turbulent intergranular plasma is no exception. Though vortex flows were first seen on the solar surface, their detection is not restricted to the near-surface levels: they have also been observed in the solar chromosphere (Wedemeyer-Böhm & Rouppe van der Voort 2009). Solar observations and realistic 3D magnetoconvection simulations have both revealed vortices existing across a wide range of spatial and temporal scales (Kitiashvili et al. 2012; Yadav et al. 2020, 2021a). A distinction has been made between small-scale vortex flows (Bonet et al. 2008, 2010; Park et al. 2016), large-scale vortex flows (Brandt et al. 1988; Attie et al. 2009; Requerey et al. 2018), and solar tornadoes (Wedemeyer-Böhm et al. 2012), depending on their locations and spatial extents. Ever smaller vortices are being detected with improved observational and simulation capabilities. Turbulence is thought to create vortices of varying sizes, although the exact mechanism is unknown and remains an unresolved question in solar physics (Battaglia et al. 2021).

In the solar atmosphere, magnetic fields serve to guide and influence vorticity propagation. Vorticity generated by torsional motions at the surface propagates to the upper layers of the solar atmosphere in the form of torsional Alfvén waves (Liu et al. 2019). These upwardly propagating torsional MHD waves carry a lot of energy that could be dissipated in the solar atmosphere’s upper layers, making vortex flows a possibly important source of solar atmospheric heating. Various mode transformations and mode couplings occur as these torsional Alfvén waves propagate through the solar atmosphere. Despite the fact that they are difficult to detect in observations due to their incompressible nature, some recent observations have confirmed their presence (Jess et al. 2009; De Pontieu et al. 2014; Morton et al. 2015; Srivastava et al. 2017). Since the Alfvén waves are highly incompressible, in order to heat the plasma, they require a mechanism to convert their magnetic energy into thermal energy by mode conversion, phase mixing, or a non-linear turbulent cascade. Furthermore, because Alfvén waves have magnetic tension as their restoring force, they are less affected by large thermodynamic gradients in the transition zone than the other MHD wave modes. Unlike acoustic waves, they have no propagation cutoff (Musielak et al. 2007).

Magnetic null points, such as vortices, are a common occurrence in the solar atmosphere, particularly in the solar corona. Magnetic nulls are locations where the magnetic field strength and, hence, the Alfvén speed is effectively zero. Because they are truly local, small-scale features, they are still impossible to observe with existing instruments. Numerical calculations and extrapolations of photospheric magnetograms, on the other hand, suggest that null points will be found throughout the solar chromosphere and corona (Longcope & Parnell 2009). They can be classified as positive or negative nulls, based on the sign of eigenvalues in the magnetic field’s Jacobian matrix. For example, a null point is called positive if one eigenvalue is negative and the other two are positive (or have positive real parts). The real and imaginary parts of the eigenvectors corresponding to these two eigenvalues span a plane such that all field lines are originating from the null point. These field lines form a surface called the fan surface of the null (Parnell et al. 1996; Longcope 2005). Multiple null points can also develop in certain magnetic configurations: as a result of a local bifurcation, a null point may split into two separate nulls or two nulls may coalesce to form a new null (Brown & Priest 1999).

In particular, MHD waves have been intensively investigated due to their potential as a diagnostic tool for coronal seismology and their contribution to solar atmospheric heating (e.g., see reviews in Parnell & De Moortel 2012; Arregui 2015; Nakariakov & Kolotkov 2020; Srivastava et al. 2021). Despite the existence of various wave-based heating theories, reconnection-based heating mechanisms have often been preferred since the wave energy flux is usually found to be insufficient to overcome radiative losses (Fossum & Carlsson 2005). Magnetic reconnection and heating due to waves are not independent mechanisms, however, since magnetic reconnection can be a source of waves (Verwichte et al. 2004; Jess et al. 2008; Luna et al. 2008; Li & Zhang 2012; Provornikova et al. 2018), and waves can produce instabilities needed to trigger reconnection processes (Isobe & Tripathi 2006; Lee et al. 2014; McLaughlin et al. 2009). Moreover, recent advances in observational instrumentation and simulation capabilities imply that wave-based heating could be a viable process for maintaining the high temperatures of the solar atmosphere (Abbasvand et al. 2020, 2021; Yadav et al. 2021b). These new results have drawn renewed interest in wave-heating mechanisms, both observationally and theoretically. Combinations of MHD wave aspects, vortices, and null point ingredients have only been studied to a certain extent (McLaughlin et al. 2011, 2019; Thurgood & McLaughlin 2013; Shelyag et al. 2013; Mumford et al. 2015; Tarr et al. 2017; Tarr & Linton 2019; Liu et al. 2019; Tziotziou et al. 2020; Battaglia et al. 2021). Because vortex flows, resulting torsional Alfvén waves, and localized null points are so common, their co-occurrence is quite likely, opening up new challenges in the study of wave and energy propagation.

Furthermore, MHD wave propagation in the solar atmosphere is strongly influenced by the magnetic field topology and it undergoes numerous wave transformations as it encounters obstacles such as the Alfvén-acoustic equipartition layer, the transition region, and magnetic nulls, among others. There have been several modeling efforts to investigate wave conversion and wave coupling as MHD waves travel through the different layers of the solar atmosphere with (or without) a magnetic null (Bogdan et al. 2003; McLaughlin & Hood 2004; Khomenko & Collados 2006; Gruszecki et al. 2011; Fedun et al. 2011; Vigeesh et al. 2012; Pascoe et al. 2013; Santamaria et al. 2015; Tarr et al. 2017; Tarr & Linton 2019; Sabri et al. 2020). In particular, Sabri et al. (2020) investigated the nonlinear effects associated with interaction of an Alfvén wave with a magnetic null point in a 2.5D simulation setup. They found that an incoming Alfvén pulse excites magneto-acoustic waves which in return, induce current density and are responsible for creation of plasma flows that are possible sources of jets in the solar corona.

Many previous studies are limited to 2D, where Alfvén waves are omitted as they have a perturbation normal to the plane of the adopted background magnetic field. Inclusion of a third dimension such as in 2.5D simulations or in full 3D simulations allows Alfvén waves to exist. In a 2.5D simulation, it is relatively easy to isolate the three fundamental MHD waves, as slow and fast modes can be separated by projecting velocities into the directions parallel and perpendicular to the local magnetic field in the plane of the magnetic field (at least in a low-beta plasma), while the Alfvén wave is along the normal to the plane. However, in full 3D setups it becomes really challenging to isolate the three modes if wave identification is desired. At the same time, full 3D MHD simulations have proven to be a very powerful tool to investigate 3D wave propagation and transformations in the solar atmosphere. As a result, numerous wave decomposition methods are in use that allow for the contributions of the three basic MHD waves to be separated, which is valid under certain assumptions (discussed in detail in Sect. 3). Being able to do such a wave mode decomposition has important applications. In direct solar observations, it helps us understand the role of specific wave types in various spectacularly dynamic and intricate phenomena. Depending on the type of the wave, we know that certain dissipation mechanisms (not included in our ideal MHD study here) might be more efficient. For example, in a 2D geometry, partial ionization effects are a more efficient mechanism of dissipation for fast modes than for slow modes (Popescu Braileanu & Keppens 2021; Khomenko & Cally 2019). Moreover, by comparing observed wave properties to analytical or numerical models, coronal seismology helps us to determine the physical properties of the coronal regions through which the waves are propagating (Jess et al. 2015). Even nonlinearly, the ideal MHD description features various shock types that in a complicated way generalize the linear wave modes: MHD shocks form a three-parameter family of solutions, while linear slow, Alfvén, and fast waves can be considered as shocks of infinitesimal strength, although the Alfvén waves more properly relate to rotational (Alfvén) discontinuities (Goedbloed et al. 2019). Therefore, although in any stratified medium, there are no pure slow, Alfvén and fast modes, it still remains meaningful to identify how a signal decomposes in its pure MHD wave types locally. This is done in practice by a local, mutually orthogonal decomposition, allowing to discuss wave transformations and their possible role in plasma heating.

In the present paper, we investigate the propagation of linear MHD waves generated by wave-like torsional motion at the solar surface to the higher layers of the solar atmosphere, including the photosphere, chromosphere, transition region, and the solar corona. Our magnetic field configuration resembles a local parasitic-polarity magnetic field embedded within a homogeneous magnetic field of the opposite polarity. This introduces a 3D null in the coronal regions. The goal of this research is to learn how MHD waves driven by ubiquitous vortices at the solar surface propagate across the stratified solar atmosphere, specifically in the vicinity of an isolated 3D magnetic null. We propose a new approach for identifying the three wave modes and compare it to previously applied methods. We discuss the advantages of our method over existing strategies. We compute the wave energy fluxes and investigate the impact of the presence of a null point on energy flux accumulation on various locations for the three MHD wave types. We also examine the specific locations where strong currents will be generated as these will be the sites for preferential heating when resistivity is included.

The paper is organized as follows. We describe the MHD equations and how they are solved numerically, along with the adopted stratification of the solar atmosphere and the specific magnetic topology in Sect. 2. We discuss the proposed wave decomposition method and its comparison with existing methods in Sect. 3. In Sect. 4, we discuss the results related to wave energy fluxes and current accumulation. We present our conclusions in Sect. 5.

2. Numerical setup

2.1. Governing equations

We performed three-dimensional (3D) numerical simulations of wave propagation through the stratified solar atmosphere. The simulations were carried out with the MPI-parallelised Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Porth et al. 2014; Xia et al. 2018; Keppens et al. 2021). We used a single-fluid MHD module of the code, where a background stratified equilibrium is explicitly removed from the equations, to better simulate the linear to weakly nonlinear wave regime that is of interest here. The well-known governing ideal MHD equations are as follows:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·\left(\rho \mathit{v}\right)=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial (\rho \mathit{v})}{\partial t}+\mathrm{\nabla }·[\rho \mathit{v}\mathit{v}+(p+\frac{B^2}{2})\mathbb{I}-\mathit{BB}]=\rho \mathit{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial }{\partial t}(e_{\mathrm{int}}+\frac{\rho v^2}{2}+\frac{B^2}{2})+\mathrm{\nabla }·[\mathit{v}(e_{\mathrm{int}}+\frac{\rho v^2}{2}+B^2+p)\hfill \\ \hfill & -\mathit{B}\mathit{B}·\mathit{v}]=\rho \mathit{v}·\mathit{g},\hfill \end{array}$$(3)

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

where ρ, p, v, and B represent the mass density, thermal pressure, fluid velocity, and magnetic field vector, respectively. The internal energy density is

$$\begin{array}{c}\hfill e_{\mathrm{int}}=p/(\gamma -1),\end{array}$$(5)

where γ = 5/3 is the adiabatic index. We actually split the MHD variables B, ρ and p into split-off, prescribed, time-independent parts, and their deviations, such that ρ = ρ₀ + ρ₁, p = p₀ + p₁, B = B₀ + B₁, where we subtract the magneto-hydrostatic (time-independent background) equilibrium (where possibly a non-vashing current density distribution J₀ = ∇ × B₀ enters) such that:

$$\begin{array}{c}\hfill \frac{\partial {\rho }_0}{\partial t}=0,\frac{\partial p_0}{\partial t}=0,\frac{\partial {\mathit{B}}_0}{\partial t}=\mathbf{0},-\mathrm{\nabla }{p}_0+{\rho }_0\mathit{g}+{\mathit{J}}_0\times {\mathit{B}}_0=\mathbf{0},\end{array}$$

Then Eqs. (1)–(4) become:

$$\begin{array}{cc}& \frac{\partial {\rho }_1}{\partial t}+\mathrm{\nabla }·\left(\rho \mathit{v}\right)=0,\hfill \end{array}$$(6)

$$\begin{array}{cc}& \frac{\partial (\rho \mathit{v})}{\partial t}+\mathrm{\nabla }·[\rho \mathit{vv}+(p_1+\frac{B_1^2}{2}+{\mathit{B}}_0·{\mathit{B}}_1)\mathbb{I}\hfill \\ \hfill & -({\mathit{B}}_0{\mathit{B}}_1+{\mathit{B}}_1{\mathit{B}}_0+{\mathit{B}}_1{\mathit{B}}_1)]={\rho }_1\mathit{g},\hfill \end{array}$$(7)

$$\begin{array}{cc}& \frac{\partial }{\partial t}(e_{\mathrm{int}1}+\frac{\rho v^2}{2}+\frac{B_1^2}{2})+\mathrm{\nabla }·[\mathit{v}(e_{\mathrm{int}}+\frac{\rho v^2}{2}+p_1)-(\mathit{v}\times \mathit{B})\times {\mathit{B}}_1]\hfill \\ \hfill & ={\rho }_1\mathit{v}·\mathit{g}-p_0\mathrm{\nabla }·\mathit{v}-(\mathrm{\nabla }\times {\mathit{B}}_0)·(\mathit{v}\times {\mathit{B}}_1),\hfill \end{array}$$(8)

$$\begin{array}{cc}& \frac{\partial {\mathit{B}}_1}{\partial t}+\mathrm{\nabla }·(\mathit{v}\mathit{B}-\mathit{B}\mathit{v})=\mathbf{0}.\hfill \end{array}$$(9)

2.2. Initial condition

The computational domain is a 3D box in Cartesian geometry (x, y, z) such that z corresponds to the vertical direction, that is, it is normal to the solar surface. Our simulation domain extends from −30 ≤ x ≤ 30 Mm, −30 ≤ y ≤ 30 Mm and 0 ≤ z ≤ 10 Mm and is resolved by 200 × 200 × 300 grid points.

We used the temperature profile from VAL-C model of the solar atmosphere in the photosphere and the chromosphere (Vernazza et al. 1981), and joined it smoothly to a constant temperature of 1 MK for the solar corona. This slightly smoothes the transition region variation (which is not well represented anyway as non-adiabatic effects such as thermal conduction are ignored here), as seen in Fig. 1. The background state is a full magnetohydrostatic equilibrium, namely,

$$\begin{array}{c}\hfill -\mathrm{\nabla }{p}_0+{\mathit{J}}_0\times {\mathit{B}}_0+{\rho }_0\mathit{g}=0,\end{array}$$(10)

Fig. 1. Height variation of temperature (green), density (blue), and pressure (red). The left vertical axis displays the temperature and the pressure, and the right vertical axis displays the mass density in CGS units.

where p₀ is the background kinetic pressure, g is the solar gravitational acceleration, and ρ₀ is the background density. In fact, for simplicity, we use a potential magnetic field, B₀, so that the current term, J₀, and the background Lorentz force vanish. The resulting pressure and density are then just one-dimensional profiles, and the height variation of pressure and density, along with our input temperature profile are shown in Fig. 1. The initial potential configuration of the atmospheric magnetic field is similar to that of Candelaresi et al. (2016). It includes an isolated 3D magnetic null point located at the apex of a fan surface. This magnetic field topology represents a parasitic polarity embedded in a dominant opposite polarity magnetic field. The analytic expressions for the three components of the magnetic field vector in dimensionless form are:

$$\begin{array}{cc}& B_x=-B_0\frac{3x(z-d)}{{[x^2+y^2+{(z-d)}^2]}^{5/2}},\hfill \end{array}$$(11)

$$\begin{array}{cc}& B_y=-B_0\frac{3y(z-d)}{{[x^2+y^2+{(z-d)}^2]}^{5/2}},\hfill \end{array}$$(12)

$$\begin{array}{cc}& B_z=B_0+B_0\frac{[x^2+y^2-2{(z-d)}^2]}{{[x^2+y^2+{(z-d)}^2]}^{5/2}},\hfill \end{array}$$(13)

where B₀ and d are the parameters that decide the magnetic field topology and the location of the magnetic null point. We used B₀ ≃ 11.2 G and d = −6 Mm for the present simulations. The map of the background magnetic field strength (in real units) across a vertical cut passing through the magnetic null is displayed in Fig. 2. To get a 3D view of our magnetic topology, we numerically construct a few constant magnetic flux tubes shown by their field-line boundaries, that is, we show “flux surfaces” in Fig. 3. They enclose the same amount of magnetic flux over their whole length. The equilibrium magnetic field is rotationally symmetric about the vertical (z in our case) axis. The propagation of MHD waves and their energy fluxes will be studied in further detail using the three inner flux surfaces depicted in red, green, and blue colors.

Fig. 2. Spatial map of the magnetic field strength at the y = 0 plane. The red lines represent magnetic field lines, and the yellow lines are contours for β = 1: this happens in the chromosphere and near the null point.

Fig. 3. 3D visualization of various magnetic flux surfaces.

Height variations of Alfvén speed (green curves) and sound speed (red curve) are displayed in Fig. 4. Since the sound speed depends only on the local plasma temperature (from pressure and mass density), it is independent of the horizontal location and is only height-dependent. In contrast, the Alfvén speed includes the local magnetic field that has 3D spatial structuring and hence depends both on the height and horizontal position. For illustration, we show the height variation of Alfvén speed through the central (solid green curve) and corner (dashed green curve) position of the domain.

Fig. 4. Height variation of Alfvén and sound speeds in the equilibrium atmosphere. By construction, the sound speed is only height-dependent. On the other hand, the Alfvén speed profile changes in a two-dimensional fashion as it depends on the magnetic field strength; we show both the central and corner variation.

2.3. Boundary conditions and wave driver

To avoid the reflection of waves from the upper boundary, we employed a uniform absorbing layer of a thickness of 1 Mm at the top boundary, such that the time-dependent variables are damped exponentially in this layer:

$$\begin{array}{c}\hfill u^{\prime }=u\exp (-a\frac{z-z_{\mathrm{db}}}{z_{\mathrm{tb}}-z_{\mathrm{db}}}),\end{array}$$(14)

where z_db and z_tb are the bottom and the top of the absorbing layer respectively, and a is the damping coefficient. This kind of absorbing layer has been successfully used in previous numerical studies related to wave propagation in the solar atmosphere (Felipe et al. 2010; Popescu Braileanu & Keppens 2021). Lateral boundaries are kept periodic for all variables. The bottom boundary includes the wave driver in the horizontal velocity components, all other time dependent variables viz. p₁, ρ₁, B₁, and v_z are kept at zero.

The wave driver represents an oscillating spiral vortex at a circular region of a radius 1 Mm and has the following analytical form:

$$\begin{array}{cc}& v_x=-v_0\sin (\omega t)\sin (\varphi +{\varphi }_0)\sqrt{x^2+y^2},\hfill \end{array}$$(15)

$$\begin{array}{cc}& v_y=v_0\sin (\omega t)\cos (\varphi +{\varphi }_0)\sqrt{x^2+y^2},\hfill \end{array}$$(16)

where v₀ = 1 km s⁻¹ is the amplitude, ω is the driver frequency, $\varphi ={\tan }^{-1}\left(\frac{\mathit{y}}{x}\right)$ is the azimuth angle in the bottom plane, and ${\varphi }_0=\frac{\pi }{10}$ is the expansion parameter of the spiral vortex flow. This velocity driver at the bottom boundary perturbs the equilibrium to excite MHD waves.

We used a velocity driver with a frequency ω ≃ 26.2 mHz, which equates to a time period of 240 s. Since the driver frequency is much greater than the Brünt–Väisälä frequency (∼2.2 mHz) in our corona, dispersion due to gravity will not be strong enough to change the nature of pure MHD waves. Figure 5 displays the spatial profile of the velocity driver in the bottom plane at the peak amplitude. We used a velocity driver with a small amplitude such that despite its propagation through a gravitationally stratified atmosphere, perturbations remain small compared to the background values. We do this to focus on the linear to weakly nonlinear regime in the present study and to avoid nonlinear effects, such as wave steepening and nonlinear coupling.

Fig. 5. Spatial profile of the velocity driver applied in the ghost cells of the simulation domain at peak amplitude. Green lines are streamlines of the velocity vector field, over-plotted on the map of velocity magnitude.

We note that our driver does not use a prechosen (vertical) spatial wavelength directly, but fixes the temporal variation instead. In general, it is only when the wavelength is much smaller than the scale on which background quantities vary that we can approximate the background plasma to be locally uniform and MHD waves to be identifiable in their pure linear slow, Alfvén and fast form. In our case, because of gravitational stratification, the wavelength of propagating waves may increase and become comparable to the local density scale height of the solar corona. However, our simulation domain extends to a distance of about one fifth of the scale height that is approximately 50 Mm for coronal plasma at 1 MK. Therefore, background plasma parameters remain more or less uniform in our coronal part, as can also be seen in Fig. 1.

3. Wave mode decomposition

3.1. MHD wave modes

In a compressible, homogeneous and uniformly magnetized plasma of infinite extent, MHD theory supports three fundamental linear MHD waves with regard to slow magneto-acoustic waves, fast magneto-acoustic waves, and Alfvén waves. In an inhomogeneous setting, one can still identify the corresponding slow, Alfvén and fast subspectrum, which organize the entire collection of MHD eigenmodes (Goedbloed et al. 2019). The three wave modes become coupled in the presence of a gravitational stratification and a non-uniform magnetic field. In particular, magneto-acoustic waves traveling through the solar atmosphere get coupled and undergo linear mixing in the layer where the sound and the Alfvén speed are comparable (Schunker & Cally 2006; Cally 2006). These waves, on the other hand, are mostly decoupled in places where either magnetic or thermal pressure dominate. In such situations, one of the waves is acoustic, while the other is magnetic. Decomposing the plasma motions into locally slow magneto-acoustic, fast magneto-acoustic, and Alfvén waves is essential to examine the numerous wave transformations occurring in the solar environment and to quantify wave energy fluxes that relate to the three fundamental MHD waves.

A great deal of intuitive work related to MHD wave properties deals with their known behavior in uniform media. Magneto-acoustic waves have a parallel component of plasma motion and are driven by total pressure (gas pressure and magnetic pressure) and magnetic tension. They do not propagate parallel vorticity and have compressible plasma motions. Alfvén waves are driven purely by magnetic tension and exhibit transverse plasma movements. They propagate vorticity and have incompressible plasma motions. Invoking the compression due to parallel and transverse velocity motions and the parallel vorticity, Cally (2017) constructed proxies that can be used to separate out the various wave types. However, in the gravitationally stratified solar atmosphere and in a non-uniform plasma, these waves are no longer pure modes and can propagate both compression and vorticity (Goossens 2003; Roberts 2019).

In much of the solar chromosphere and corona, the regions we are dealing with in the present paper, magnetic pressure dominates over thermal pressure and the plasma beta is much less than unity. In such low-beta plasma, where the Alfvén speed is much larger than the sound speed, magneto-acoustic oscillations can meaningfully be split into slow and fast components, where slow magneto-acoustic waves are predominantly guided along the magnetic field vector. From Fig. 4, one notices that this is true over most of our domain, except in the lowest atmospheric regions and near the null.

3.2. An overview of MHD mode decompositions

A few approaches for isolating individual, pure MHD wave modes in low-beta plasma regions have been proposed and used in the literature, although most of them are only suitable to 2D or 2.5D systems. In 2D simulations, it is relatively easy to separate out the wave modes, as Alfvén waves do not exist while slow and fast mode can be separated by taking projections of the velocity vector along and normal to the background magnetic field (Rosenthal et al. 2002; Bogdan et al. 2003; Khomenko et al. 2009). In 2.5D simulations, there is an invariant direction such that the magnetic pressure gradient vanishes and the only driving force in the invariant direction is magnetic tension. Therefore, plasma motions along this invariant direction are associated with Alfvén waves. Fast and slow waves can be distinguished by taking projections along and normal to the background field in the same way as in 2D simulations (Thurgood & McLaughlin 2012, 2013). If the magnetic field is in the x − y plane, the velocity perturbations corresponding to the three MHD waves can be calculated by the projections of the velocity vector v as shown below:

$$\begin{array}{cc}& v_{\mathrm{alf}}=\mathit{v}·\widehat{\mathit{z}},\hfill \end{array}$$(17)

$$\begin{array}{cc}& v_{\mathrm{slow}}=\mathit{v}·\widehat{\mathit{s}},\hfill \end{array}$$(18)

$$\begin{array}{cc}& v_{\mathrm{fast}}=\mathit{v}·(\widehat{\mathit{z}}\times \widehat{\mathit{s}}),\hfill \end{array}$$(19)

where v is a linear velocity field on top of a magnetohydrostatic equilibrium and $\widehat{\mathit{s}}$ is the unit vector along the (in-plane) magnetic field vector.

In full 3D simulations, where the magnetic field is not confined to a 2D geometry, the isolation of MHD waves becomes non-trivial as there are an infinite number of vectors normal to the magnetic field vector. However, under certain assumptions, it is possible to derive analytical expressions for characteristic directions for velocity fluctuations associated with the different MHD modes. We discuss a few methods that are widely used in the literature below.

One possible way to separate the three different MHD waves is by taking velocity projections onto the following three characteristic directions (Khomenko & Cally 2011, 2012):

$$\begin{array}{cc}& \widehat{\mathit{s}}=[\cos \varphi \sin \theta ,\sin \varphi \sin \theta ,\cos \theta ],\hfill \end{array}$$(20)

$$\begin{array}{cc}& \widehat{\mathit{\varphi }}=[-\cos \varphi {\sin }^2\theta \sin \varphi ,1-{\sin }^2\theta {\sin }^2\varphi ,-\cos \theta \sin \theta \sin \varphi ],\hfill \end{array}$$(21)

$$\begin{array}{cc}& \widehat{\mathit{n}}=[-\cos \theta ,0,\cos \varphi \sin \theta ],\hfill \end{array}$$(22)

where a plane parallel atmosphere is permeated by the magnetic field such that θ is the inclination angle, measured from the vertical and ϕ is the azimuth angle, measured from the x − z plane. Here, the first projection ($\widehat{\mathit{s}}$) corresponds to the slow magneto-acoustic wave, the projection along $\widehat{\mathit{\varphi }}$ corresponds to the Alfvén wave, according to the asymptotic polarization direction derived by Cally & Goossens (2008), and the projection along $\widehat{\mathit{n}}$ selects the fast magneto-acoustic wave. This specific decomposition method is derived assuming that the wave propagation vector is confined to the x − z plane such that the y-component in $\widehat{\mathit{n}}$ is always zero. However, as an approximation, this decomposition has been also applied to fully three dimensional setups where the wave propagation vector is not constrained to lie in the 2D plane (Felipe et al. 2010).

Shelyag et al. (2016) instead decomposed the different MHD wave modes in a low-beta solar atmosphere by projecting the velocity vector in the three characteristic directions defined by the following unit vectors:

$$\begin{array}{cc}& \widehat{\mathit{s}}=[\cos \varphi \sin \theta ,\sin \varphi \sin \theta ,\cos \theta ],\hfill \end{array}$$(23)

$$\begin{array}{cc}& \widehat{\mathit{\varphi }}=[\sin \varphi ,-\cos \varphi ,0],\hfill \end{array}$$(24)

$$\begin{array}{cc}& \widehat{\mathit{n}}=[\cos \varphi \cos \theta ,\sin \varphi \cos \theta ,-\sin \theta ].\hfill \end{array}$$(25)

This decomposition can only be used when the flux tube is symmetric with regard to the z axis that defines the direction of gravity throughout the domain. Indeed, we see here how the Alfvén velocity fluctuations are strictly in the x − y plane. However, in general cases when the flux tubes are inclined or twisted, this method loses merit as velocity perturbations corresponding to Alfvén waves can also have components along the z direction.

Riedl et al. (2019) employed a method similar to the one used by Shelyag et al. (2016) to decompose velocity components. However, these authors calculated the direction vectors numerically instead of using analytic expressions. For example, to calculate the direction vector for Alfvén waves, they found the iso-contour of the magnetic field for all horizontal layers and then calculated a linear fit at each pixel. This method works perfectly when the flux tube profile is known beforehand and is symmetric about the z axis. However, for curved magnetic flux tubes or in general for any complex 3D time-dependent magnetic field, the azimuthal direction vector calculated using isocontours of the magnetic field at a fixed geometric height would not be fully consistent.

Instead, Mumford et al. (2015) invoked the concept of “flux tubes” and “flux surfaces” to split the velocity perturbation vector into three waves. A flux tube is a well defined theoretical construct, mostly used in analytical models where we have, for instance, cylindrical (or translational) geometry that contains a constant amount of magnetic flux at each cross-section. A flux surface is defined as any smooth surface with a local normal $\widehat{\mathit{n}}$ such that locally, $\mathit{B}·\widehat{\mathit{n}}=0$; namely, a surface buildup by individual magnetic field lines. To distinguish between the three MHD modes in such geometry, they took projections of the local velocity vector along three directions: along the magnetic field vector (for the slow magneto-acoustic wave); normal to the flux tube (for the fast magneto-acoustic wave); and along the azimuthal vector that is perpendicular to the magnetic field and parallel to the flux surface (Alfvén wave). To construct the magnetic flux surface numerically in their simulation domain, they selected seed points on the circumference of a circular area and traced magnetic field lines in the 3D space. To calculate the normal vector for their axially symmetric flux surface, they made use of neighboring field lines and computed normal vectors by constructing polygons. Once the normal vector, $\widehat{\mathit{n}}$, is known, the azimuthal unit vector can readily be calculated as $\widehat{\mathit{\varphi }}=\widehat{\mathit{s}}\times \widehat{\mathit{n}}$. This procedure is shown in Fig. 6 in a cartoon representation. This triplet of mutually orthogonal directions has also been used for detailed MHD spectroscopy of inhomogeneous 2.5D configurations, namely, fusion devices (Goedbloed et al. 2019) and 2.5D flux ropes in stratified atmospheres (Blokland & Keppens 2011a,b).

Fig. 6. Cartoon representation of wave decomposition methods. Left: wave decomposition approach utilized by Mumford et al. (2015). Right: our proposed wave decomposition method based on the field-line geometry. The unit direction vectors for motions associated with Alfvén, fast magneto-acoustic, and slow magneto-acoustic waves are represented by red, blue, and green vectors, respectively, while the fluid velocity vector is represented by a grey vector.

The direction vectors calculated by the aforementioned method depend on the chosen flux surface, therefore, we refer to this method as the “flux-surface method”. Mumford et al. (2015) constructed axially symmetric magnetic flux surfaces that are appropriate given their axisymmetric magnetic field configuration. However, in reality, a magnetic configuration can be significantly more complex, for instance, in fully 3D magneto-convection simulations as well as in actual solar observations. Observations of solar active regions, for example, show that the overall emerging flux is the result of several small flux bundles emerging concurrently or in succession, rather than a single cohesive, monolithic flux bundle (Emonet & Moreno-Insertis 1998). Therefore, the choice of flux surface can greatly influence these direction vectors and, hence, the velocity components corresponding to MHD waves.

All of the methods discussed above work well in special cases when a clear flux tube geometry is present or when we have prior knowledge of the nature of the wave propagation vector. However, in general cases, for instance, in realistic simulations or in real solar observations, the methods discussed above face some limitations. A method for decomposing waves that can be used in more realistic scenarios is based on the compressive versus rotational nature of the velocity vectors (Cally 2017; Khomenko et al. 2018; Khomenko & Cally 2019). In this method, using the (assumed linear) velocity from the simulations, the following quantities are calculated:

$$\begin{array}{cc}& f_{\mathrm{alf}}=\widehat{\mathit{s}}·(\mathrm{\nabla }\times \mathit{v}),\hfill \end{array}$$(26)

$$\begin{array}{cc}& f_{\mathrm{slow}}=\widehat{\mathit{s}}·\mathrm{\nabla }(\mathit{v}·\widehat{\mathit{s}}),\hfill \end{array}$$(27)

$$\begin{array}{cc}& f_{\mathrm{fast}}=\mathrm{\nabla }·(\mathit{v}-(\mathit{v}·\widehat{\mathit{s}})\widehat{\mathit{s}}).\hfill \end{array}$$(28)

Here, the quantities computed are not velocities but relate to their derivatives. Therefore, these proxies cannot directly be used to calculate the wave energy fluxes.

3.3. Field line-based mode decomposition

To overcome the limitations of the methods discussed above, we suggest a new technique that is inspired by widely known Frenet-Serret formulae. In the proposed approach, individual magnetic field lines are used for wave identification rather than constructing constant magnetic flux surfaces.

In this approach, just as in the other decomposition strategies discussed above, the unit vector in the direction of the background magnetic field vector corresponds to the slow magneto-acoustic waves ($\widehat{\mathit{s}}$). Then, we calculate the unit vector corresponding to the Alfvén wave and utilize the mutual orthogonality for computing the third unit vector for the fast magneto-acoustic wave. To compute a uniquely defined unit vector corresponding to the Alfvén wave, we first constructed a local plane at any point on the curved magnetic field-line by connecting it to two neighboring points on each side of the field-line and then we calculated the unit vector normal to this plane and associated it to the Alfvén wave ($\widehat{\mathit{\varphi }}$). Since the solar atmosphere is gravitationally stratified, the magnetic field expands with height and, hence, the magnetic field lines always possess substantial curvature. Finally, the unit vector corresponding to the fast magneto-acoustic wave is calculated by $\widehat{\mathit{n}}=\widehat{\mathit{\varphi }}\times \widehat{\mathit{s}}$. To be consistent throughout the domain, we implement a check such that our unit vector for Alfvén wave does not flip direction after crossing the point of inflection. Moreover, since we wish to compare our results with the method used by Mumford et al. (2015), we additionally require that the unit vector for fast magneto-acoustic wave is always orientated outwards of the flux surface. A cartoon representation of this proposed method is illustrated in Fig. 6. The advantage of this method is that whenever we have field lines available (either numerically or analytically), it is always possible to define this frame on every point of every field line, and we never need to resort to flux surfaces or tubes. Moreover, the projections are well defined for all curves and can directly lead to a quantification of wave energy fluxes along field lines.

4. Results and discussion

4.1. Equivalence of flux surface and field-line based decompositions

To compare and validate the suggested approach of wave decomposition, we select three flux surfaces (red, green, and blue in Fig. 3) and then isolate the velocity perturbations corresponding to the three MHD waves using both the proposed “field-line based method” and the “flux-surface method” used by Mumford et al. (2015) for a fixed time instance. The choice of flux surfaces is based on their spatial topology with respect to the magnetic null point. The flux surface shown in red stays quite far from the spine and fan surface of the null configuration, while the surface in green reaches the vicinity of the inner spine and the fan surface. The flux surface in blue samples parts of the outer spine and fan surface and it is closest to the magnetic null among the three selected surfaces. The three selected magnetic flux surfaces are shown in row a of Fig. 7: we note that these three flux surfaces are shown on the same scale, but their actual physical extent and relative configurations are different (as seen in Fig. 3).

Fig. 7. Comparison and validation of the proposed wave mode decomposition method with the previously used method: selected magnetic flux surfaces are shown in row a; velocity component associated with slow magneto-acoustic wave is shown in row b; velocity component associated with fast magneto-acoustic wave and Alfvén wave using our new method is shown in rows c and e, respectively; velocity component associated with fast magneto-acoustic wave and Alfvén wave using the method used by Mumford et al. (2015) is displayed in rows d and f, respectively. Over-plotted yellow and turquoise solid line curves represent the Alfvén-acoustic equipartition layer and the transition region, respectively.

In our proposed field-line based wave decomposition approach, we first need to trace the magnetic field lines in our 3D simulation domain. In fact, we selected three circular seed areas at the bottom of our simulation domain (z = 0) with varying radii and placed 500 seed points on each of them to trace multiple field lines. Due to the axisymmetric nature of the background magnetic field, tracing the associated magnetic field lines through the simulation domain provides three constant magnetic flux surfaces as shown in row a of Fig. 7. In general, seed points can be randomly placed, however, here they are placed on a circular region to provide constant magnetic flux surfaces, allowing us to compare directly to earlier approaches that relied on calculating normal vectors of flux surfaces. To calculate the velocity perturbation associated with slow magneto-acoustic waves, we calculate the corresponding tangent unit direction vectors. The component of local velocity vector in this direction constitutes the velocity perturbation associated with slow magneto-acoustic waves (v_slow). This is shown in row b of Fig. 7.

To isolate the other two modes in our field-line based method, we use the fact that slow and fast magneto-acoustic waves have perturbations in the plane formed by joining three consecutive points on the field-line while the Alfvén wave has velocity perturbations in the direction normal to this plane in which the background field is lying. Once the unit direction vectors have been computed, we compute the velocity perturbation associated with various waves by taking projections of velocity vector onto these directions. Our simulations run for ∼48 min in total and reach a periodically repeating state after around 15 min, where the system is time-varying, with a quasi-steady background state and oscillations around its mean value. In this quasi-steady state, energy input to the domain through the applied wave driver is balanced by the energy dissipation through damping at the upper absorbing layer (mimicking an upper open boundary, to some extent) or through numerical dissipation. To compare and validate our proposed approach, we chose data at a selected time (∼36 min). To obtain the velocity perturbation associated with the Alfvén wave (v_alf) and fast magneto-acoustic wave (v_fast), velocity components along the azimuthal unit vector ($\widehat{\mathit{\varphi }}$) and normal vector ($\widehat{\mathit{n}}$) are calculated and displayed in rows c and e of Fig. 7, respectively.

For comparison and validation of our proposed field-line based wave identification method, we also calculated the velocity components for three MHD waves using the flux-surface method. In this method, the direction vector for slow magneto-acoustic waves is calculated in the same way namely, along the local magnetic field vector and is, hence, identical to the one shown in panel b. To calculate the direction vector corresponding to fast magneto-acoustic and Alfvén waves that lie in the plane normal to the background magnetic field vector, a small plane is constructed at each point on the magnetic flux surface, by joining it to the two adjacent points of the neighboring field line. Then, the normal vector to this plane is calculated for the fast magneto-acoustic wave. For consistency, this normal vector is always oriented away from the central axis of the background magnetic field. The unit vector corresponding to the Alfvén wave is then computed using the orthogonality condition, that is, $\widehat{\mathit{\varphi }}=\widehat{\mathit{s}}\times \widehat{\mathit{n}}$. Velocity perturbations along these directions are computed and shown in rows d and f of Fig. 7, which represent velocity perturbations corresponding to Alfvén wave and fast magneto-acoustic waves, respectively. We note that the selected neighboring magnetic field-line must be part of the constant magnetic flux surface, which was easy to construct here due to axisymmetry.

If we compare panels in rows c and e that are calculated by field-line based method, with rows d and f that are calculated by flux-surface method, we see that both methods produce similar results. Small differences between them are purely numerical and their accuracy relies on different factors, for instance, on the number of seed points chosen to construct the surface in the flux-surface method, or in our proposed method on the step-size used for spatial integration while tracing a field line. Thus, we suggest that both our proposed method and the flux-surface method works equally well in decomposing linear MHD waves in the solar atmosphere, whenever we have a situation characterized by axisymmetric magnetic fields. It is important to note that the wave decomposition used by Shelyag et al. (2016) would also work equally well here as the magnetic field is symmetric about the z axis. However, it would give inconsistent results if the magnetic field would possess finite twist, instead of expanding only in the vertical direction. In axisymmetric twisted flux tube scenarios, the flux-surface method and the one proposed here would work. However, in realistic magneto-convection simulations and in actual solar observations, magnetic fields are almost never fully axisymmetric as magnetic flux elements are continuously jostled by the convective flows. In these circumstances, it is impossible to create a unique flux surface at any given location, however, a magnetic field-line can be traced in the domain in a unique way. Therefore, the proposed field-line based method has merits over other methods to investigate wave propagation and wave mode transformations taking place in more realistic models and in actual solar observations as it only relies on the individual magnetic field-line geometry to identify the three MHD waves.

4.2. Analysis of wave fluxes and effect of the null point

To quantify the energy content in the various MHD modes, we take a spatial average of the energy fluxes over all the field lines constituting the three selected surfaces and then additionally average over three wave periods. In ideal linear MHD, the wave energy density is a conserved quantity and is related with the wave energy flux of the following form (Bray & Loughhead 1974; Leroy 1985; Bogdan et al. 2003):

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{waves}}=p_1\mathit{v}+{\mathit{B}}_1\times (\mathit{v}\times {\mathit{B}}_0),\end{array}$$(29)

where the quantities with subscript ‘1’ represent the perturbation part that evolve during the simulation while the quantities with subscript ‘0’ correspond to the steady background. The first term on the right hand side of the above equation represents the acoustic flux, while the second term represents the Poynting or magnetic flux.

In a low beta plasma, found in the majority of our simulation domain and in the (upper) solar chromosphere and corona throughout, the slow magneto-acoustic waves are predominantly of acoustic nature while fast magneto-acoustic waves and Alfvén waves are of a magnetic nature. Therefore, comparing the acoustic and magnetic wave energy fluxes separates the slow wave energy flux from the combined energy flux associated with fast and Alfvén waves. However, this does not provide any insight about the energy flux distribution among fast magneto-acoustic waves and Alfvén waves. Moreover, it is important to separate out the energy flux contribution from different waves in order to estimate the available energy fluxes in these waves and to investigate energy transfer among various wave modes in the process of mode transformation. For example, in the case of standing waves, the total acoustic and magnetic flux might be zero but individual wave energy fluxes will always give positive quantities with no cancellation (Felipe 2012).

Therefore, to get an estimate of the contribution of individual MHD modes, we calculate the energy flux contained in the three different MHD waves using the following relations:

$$\begin{array}{cc}& F_{\mathrm{slow}}=\rho v_{\mathrm{slow}}^2{c}_{\mathrm{s}},\hfill \end{array}$$(30)

$$\begin{array}{cc}& F_{\mathrm{alf}}=\rho v_{\mathrm{alf}}^2{v}_{\mathrm{A}},\hfill \end{array}$$(31)

$$\begin{array}{cc}& F_{\mathrm{fast}}=\rho v_{\mathrm{fast}}^2\sqrt{c_{\mathrm{s}}^2+v_{\mathrm{A}}^2}.\hfill \end{array}$$(32)

These expressions just combine the ram pressure in the decomposed wave modes (identified from the MHD mode decomposition) with the relevant group speed of the mode in question: in a low-beta plasma, this group speed is the sound speed for the slow, the Alfvén speed for the Alfvén mode, and the combined sound-Alfvén speed for the fast mode.

In a few recent works, the total wave energy flux vector from Eq. (29) is decomposed into three components by taking projections in the same way as done for velocity components associated with the different MHD wave modes (Mumford et al. 2015; Leonard et al. 2018). However, the wave energy flux component computed in this way is not always a correct representation and should be interpreted with caution. For slow magneto-acoustic waves, wave energy can be computed as the projection of the acoustic component of the total wave energy flux vector along the magnetic field vector, as it travels predominantly along the magnetic field in a low beta plasma. However, the parallel component of the total wave energy flux vector also gets a contribution from the magnetic component associated with Alfvén and fast magneto-acoustic waves and therefore does not necessarily quantify the purely slow wave energy flux. In order to be able to separate the net wave energy flux associated with the different MHD waves, knowledge of the wave propagation vector is a pre-requisite, which is not known because fast magneto-acoustic waves travel (nearly) isotropically, originating from various possible sources. Therefore, in the present work, we restrict ourselves to calculating the “available” energy flux using the expressions given in Eqs. (30)–(32).

The three rows in Fig. 8 shows the energy flux comparison for the three selected surfaces as can be distinguished by the background color. The panels in the left column compare the “available” wave energy flux corresponding to the different MHD waves (normalized to the sum of the wave energy fluxes at that location) while the panels at right compare magnitudes of acoustic and magnetic parts of the wave energy flux vector, F_wave (normalized to the sum of acoustic and magnetic flux at that location). The color coding for various curves is the same for the three panels in each column (according to the inserted legends in row a). Another point to notice is that the length of the horizontal axis corresponds to the total length of the field-line constituting the selected surface, which, hence, varies from row to row. The right end of the horizontal axis always corresponds to the end of the field-line that is farthest from the axis while lying in the z = 0 plane. The left end of the axis, that is, s = 0 corresponds to z = 0 for a red and green surface, while for a blue surface, it corresponds to z = 10 Mm, that is, the top of the simulation domain. The vertical lines represent the locations of the Alfvén-acoustic equipartition layer (yellow) and the transition region (turquoise). There are two vertical yellow and two vertical turquoise lines for red and green panels as the field lines return back to the surface and encounter the transition region and Alfvén-acoustic equipartition layer twice. While for the blue surface, there is only one vertical yellow and one vertical turquoise line, as the field lines for this surface do not return back to the surface and end on the top of the simulation box. Here, we must note that although the blue surface crosses very close to the null, it does not cross the equipartition layer surrounding the null. Since the wave decomposition method proposed here is not fully applicable for the high beta regions, we restrict our discussion of wave energy fluxes to the regions confined within the yellow vertical lines for red and green surfaces, while for the blue surface, we consider the region on the left of the yellow line.

Fig. 8. Energy flux comparison for the three selected surfaces also indicated as background colors: red (first row), green (second row), blue (third row). Left column: energy fluxes calculated using Eqs. (30)–(32): slow (red line), fast (green line), and Alfvén (blue line). Right column: acoustic (orange line) and magnetic (brown line) fluxes, calculated as the first and second terms in Eq. (29), respectively.

Row a shows that the energy flux is distributed among magneto-acoustic waves and Alfvén waves as the applied driver is not necessarily in the direction normal to the magnetic field and contains both magnetic and acoustic flux. Since the fast wave energy flux is expected to get focused at the null point, there is not much fast wave energy flux on this surface as it is far from the null point (more discussion on selective fast wave energy accumulation follows).

Row b displays the results for the green surface which reaches closer to the null point and the fan surface than the red surface. Here in the left panel, we see that the energy flux is dominated by the Alfvén wave and consistent with that, the right panel is dominated by magnetic flux. The accumulation of Alfvén waves at the fan surface has been previously shown by both analytical and numerical investigations where twist propagates towards a null point along the spine axis (Galsgaard et al. 2003; Thurgood & McLaughlin 2013). These studies, on the other hand, were constrained to a uniform background density, and they solved the ideal MHD equations with the β = 0 approximation. Similar results are found for the blue surface shown in row c of Fig. 8. Just like the green surface, a large part of this surface is close to the fan surface, and hence is dominated by Alfvén wave energy flux in the left panel and by magnetic flux in the right panel. Interestingly, we see a sharp increase in the slow and fast wave energy fluxes on the left panel at around s = 5 Mm, and an increase in the acoustic flux in the right panel at the location of closest proximity to the null point.

Since the Alfvén speed drops down to nearly zero near the null point, Alfvén waves cease to propagate and spread over the fan surface as magnetic field lines diverge. This explains the sharp dip in Alfvén energy flux when the blue surface is in close proximity for the null point. Moreover, a fast wave undergoes refraction due to an inhomogeneous Alfvén speed and gets focused near the null point. This focusing results in an increase of fast wave energy flux at the null point clearly seen here for the blue surface. This fast wave accumulation near the null point can also be seen in the rows e and f of Fig. 7; especially in the panels corresponding to the blue surface, as it is closest to the null point. Increase in the slow wave energy flux can be a possible result of mode conversion taking place at the Alfvén-acoustic equipartition layer surrounding the null point (McLaughlin & Hood 2004). Thus, the results obtained here are consistent with previous analytical and numerical investigations, showing the wave accumulation at predictable locations of the magnetic topology (McLaughlin & Hood 2006; Thurgood & McLaughlin 2013). We note, hence, that our quantification of the wave energy fluxes by the two means shown in Fig. 8 provide complimentary info on the true wave transformations occurring throughout the domain.

4.3. Current density accumulation

Though we did not include resistivity in our model, to gain a notion of the locations where preferential heating might take place in such a system, we calculate the current density throughout the domain. Figure 9 displays the magnitude of the current density vector for a slice (y = 0) that is calculated by taking a time-average over three wave periods. We find strong current buildup at the spine and also at the fan surface. This is consistent with the wave energy analysis, as Alfvén waves spread along the fan surface, and create gradients of the magnetic field that result in strong currents. Similar results were obtained by Thurgood & McLaughlin (2013) in which they investigated behavior of Alfvén waves near a null point in a comparatively simpler model by ignoring gravitational stratification and acoustic effects. This increase in current is extremely important, because it suggests that finite resistivity will eventually dissipate the wave energy and contribute to plasma heating at these locations in the solar corona.

Fig. 9. Spatial map of magnitude of current (time averaged over three wave periods) shown for a y = 0 slice through the domain.

Many of the previous models showed preferential current accumulation at the null point (McLaughlin & Hood 2004, 2006) that is not present in our simulations. The reason for this lack of current accumulation at the null point in our case is related to the nature of our applied wave driver. Our driver is predominately of Alfvénic nature and has very little wave power in the fast mode, as can be seen in the left columns of Fig. 8. In the earlier studies, their wave driver is in the form of fast waves, and wave energy gets focused at the null point and leads to current accumulation close to the null. These findings imply that in the actual solar case, where waves are stochastically excited due to turbulent convection, the nature of the excited waves cannot be predicted in advance and, thus, current accumulation and plasma heating can occur in different characteristic regions of a 3D null topology, such as the null, spine, or fan plane.

5. Conclusions

We performed three-dimensional ideal MHD simulations of waves travelling through a gravitationally stratified solar atmosphere extending from the solar surface to the solar corona. The magnetic field configuration was chosen to be representative of a uni-polar magnetic field region with a parasitic magnetic patch of opposite polarity at its center. This creates a null point in the corona, and such magnetic configurations are believed to be common in the lower solar atmosphere. To study the propagation of waves excited by ubiquitous vortex flows found at the solar surface, we employed a rotational velocity driver at the bottom boundary of our simulation domain and place it in the center of the parasitic polarity region. Most of our simulation domain contains low-beta plasma where the slow magneto-acoustic wave is predominantly acoustic and the fast magneto-acoustic wave is predominantly magnetic in nature. However, we do have the sound-Alfvén equipartition height in our chromosphere, and the coronal null introduces another equipartition surface centered on this null. We introduced and applied an MHD wave decomposition method to see how wave energy fluxes traverse this fairly realistic stratified magnetic topology.

In a fully three-dimensional setup, the magnetic component of the total wave energy flux vector also gets a contribution from Alfvén waves and, hence, it is not always easy to separate out the wave energy fluxes for all three different MHD waves. To tackle this issue, various wave identification methods have been used in the literature but they have certain limitations with respect to the magnetic field configuration or the wave propagation vector. Here we present another wave identification method that can avoid these limitations, so that different wave modes can be isolated in more general scenarios. Our proposed method makes use of the local Frenet frame vectors on a magnetic field-line to calculate the components of velocity perturbations corresponding to the three different MHD waves. We selected an axially symmetric magnetic field configuration, so that we can compare our proposed method with an existing method that requires us to construct or acquire knowledge of constant and nested magnetic flux surfaces.

For comparing our method and the flux-surface method, we selected three constant flux surfaces in our domain and compared the velocity components corresponding to the different MHD modes. We showed that both methods give identical results, as expected in this setting. However, the proposed method has advantages over the flux-surface method as it can be applied to more realistic scenarios where the magnetic field configuration does not have an axis of symmetry. We must note that the proposed field-line based method is only applicable for identifying pure linear MHD waves in low-beta plasma regions. However, in the proximity of the null point, magnetic gradients are steep leading to strong variations in the phase speeds of the waves, thus modifying the propagation of MHD waves. Moreover, plasma beta is more than unity near the null point and thus the proposed wave decomposition does not hold true. We further calculated the time-averaged (over three full time period of the input driver) wave energy fluxes associated with slow, fast, and Alfvén waves after the simulations reached a stationary state. We found that fast-wave energy flux and Alfvén wave energy flux accumulate at predictable locations in the null configuration, that is, near the null for fast waves and near the fan surface for Alfvén waves. Thus, we show that even in a more realistic setup, the obtained results are consistent with previous analytical and simulation studies that used several approximations in terms of background parameters or ignored acoustic effects. We found current buildup at the spine and fan surface of a null topology. In the presence of finite resistivity, excessive currents can lead to preferential heating at specific locations and may thus have important implications in triggering solar jets (Shen 2021).

In the present simulations, we have a realistic solar atmosphere in the sense that excited waves pass through various mode transformation and reflecting layers, viz. the acoustic cutoff layer, two Alfvén-acoustic equipartition layers (one near the surface and one surrounding the null) and also the transition region. Therefore, several interesting mode conversion, reflection and transmission processes take place at these layers. The amplitude of our wave driver was intentionally kept small in order to investigate the wave propagation and transformations in a linear regime. We compared the available energy fluxes in various MHD waves in this work, however, we will further investigate the mode transformations in detail in a follow-up study and compare them with a nonlinear case. Also, in the nonlinear case for the present setup, additional wave excitation sources might play an important role due to nonlinear effects. We also plan to extend the presented study by including different profiles of wave drivers in terms of their spatial localization, frequency, and polarization. Since the proposed method of wave decomposition can be applied to magnetic fields obtained by extrapolations of photospheric magnetograms, we plan to perform wave studies using a variety of wave drivers on extrapolated magnetic fields that host a magnetic null. Moreover, in the present simulations, we excited waves with a fixed frequency, while it is well known that waves excited by turbulent convection and granular buffeting cover a broadband spectrum. Therefore, including multi-frequency or stochastic drivers in the present setup can also be an interesting follow-up study. We also intend to apply the proposed wave decomposition technique in realistic magnetoconvection simulations. As shown by Yadav et al. (2021a), slow magneto-acoustic waves carry enough wave energy flux to heat the solar chromosphere. However, these waves get dissipated in the chromosphere due to shock formation and do not contribute much to the higher layers. When it comes to magnetic waves, these propagate to higher layers and may heat the solar corona. By employing the wave identification method proposed here, we intend to calculate energy flux associated with fast magneto-acoustic and Alfvén waves in such realistic simulations and estimate their energy content.

Acknowledgments

N.Y. and R.K. are supported by Internal funds KU Leuven, project C14/19/089 TRACESpace. R.K. further received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 833251 PROMINENT ERC-ADG 2018), a joint FWO-NSFC grant G0E9619N and FWO grant G0B4521N. B.P. is supported by the FWO grant 1232122N.

References

All Figures

	Fig. 1. Height variation of temperature (green), density (blue), and pressure (red). The left vertical axis displays the temperature and the pressure, and the right vertical axis displays the mass density in CGS units.
In the text

	Fig. 2. Spatial map of the magnetic field strength at the y = 0 plane. The red lines represent magnetic field lines, and the yellow lines are contours for β = 1: this happens in the chromosphere and near the null point.
In the text

	Fig. 3. 3D visualization of various magnetic flux surfaces.
In the text

	Fig. 4. Height variation of Alfvén and sound speeds in the equilibrium atmosphere. By construction, the sound speed is only height-dependent. On the other hand, the Alfvén speed profile changes in a two-dimensional fashion as it depends on the magnetic field strength; we show both the central and corner variation.
In the text

	Fig. 5. Spatial profile of the velocity driver applied in the ghost cells of the simulation domain at peak amplitude. Green lines are streamlines of the velocity vector field, over-plotted on the map of velocity magnitude.
In the text

Fig. 6. Cartoon representation of wave decomposition methods. Left: wave decomposition approach utilized by Mumford et al. (2015). Right: our proposed wave decomposition method based on the field-line geometry. The unit direction vectors for motions associated with Alfvén, fast magneto-acoustic, and slow magneto-acoustic waves are represented by red, blue, and green vectors, respectively, while the fluid velocity vector is represented by a grey vector.

In the text

Fig. 7. Comparison and validation of the proposed wave mode decomposition method with the previously used method: selected magnetic flux surfaces are shown in row a; velocity component associated with slow magneto-acoustic wave is shown in row b; velocity component associated with fast magneto-acoustic wave and Alfvén wave using our new method is shown in rows c and e, respectively; velocity component associated with fast magneto-acoustic wave and Alfvén wave using the method used by Mumford et al. (2015) is displayed in rows d and f, respectively. Over-plotted yellow and turquoise solid line curves represent the Alfvén-acoustic equipartition layer and the transition region, respectively.

In the text

Fig. 8. Energy flux comparison for the three selected surfaces also indicated as background colors: red (first row), green (second row), blue (third row). Left column: energy fluxes calculated using Eqs. (30)–(32): slow (red line), fast (green line), and Alfvén (blue line). Right column: acoustic (orange line) and magnetic (brown line) fluxes, calculated as the first and second terms in Eq. (29), respectively.

In the text

	Fig. 9. Spatial map of magnitude of current (time averaged over three wave periods) shown for a y = 0 slice through the domain.
In the text