Exact Nonlinear Decomposition of Ideal-MHD Waves Using Eigenenergies

The core idea behind the proposed decomposition method is the use of the ideal-MHD eigensystem to construct local transformations of the solutions onto the characteristic curves, where the fundamental MHD modes are discernible. We stress the word “solutions” as throughout this paper it is assumed that the MHD solutions are fully known a priori, whether by direct numerical simulation or other means, and here we merely use those solutions as input to analyse the MHD modes that exist therein. Once the local transformation is known, we may break down the rate of change of the total energy density and ultimately tease out the energy contribution from each fundamental mode. For the sake of brevity, throughout this paper we will use the word energy to mean energy density.

The full set of ideal-MHD equations describing a polytropic plasma in the absence of gravity can be written in SI units by separating the temporal and spatial derivatives as follows,

where the quantities have their usual meaning, i.e., $\rho$ is the mass density, $\boldsymbol{v}$ is the plasma velocity, $\boldsymbol{B}$ is the magnetic field, $\boldsymbol{j}=\mbox{\boldmath$\nabla$}{\mathbf{\times}}{\boldsymbol{B}}/\mu_{0}$ is the current density, $\mu_{0}$ is the permeability of vacuum, $\gamma$ is the heat capacity ratio (taken to be $5/3$ throughout this study), and $p=\kappa\rho^{\gamma}$, with $\kappa$ assumed to be constant (indicating an isentropic process). Note that the precondition on $\boldsymbol{B}$ to remain divergence-free is built into the induction equation (1c), and hence if $\mbox{\boldmath$\nabla$}{\mathbf{\cdot}}{\boldsymbol{B}}=0$ initially, then it shall so remain throughout.

The system of equations above is hyperbolic, and hence can be recast into the quasi-linear form as follows (Roe & Balsara, 1996)

in which $\boldsymbol{P}=\left(\rho,v_{x},v_{y},v_{z},B_{x},B_{y},B_{z},p\right)^{T}$ is the primitive plasma state vector (the solution vector mentioned above which is assumed to be known at all times), $M_{q}$ are $8\times 8$ time- and space-dependent flux matrices in the direction $q$ (for which, from here on out, we assume $q\in\left(x,y,z\right)$, unless stated otherwise), and summation is implied over the repeated indices. Without imposing any simplifying assumption, the flux matrices in 3D Cartesian coordinates are found to be

Given the above-mentioned hyperbolicity of the equations, the eigensystem associated with these matrices is analytical everywhere, and may be used to extract the required transformation matrices.¹¹1Note that the MHD equations are not strictly hyperbolic, as the eigensystem becomes degenerate in certain limits. However, the system remains well behaved in those limiting cases for judicious choice of normalization (Brio & Wu, 1988; Roe & Balsara, 1996). To this end, solving the eigenvalue problem $M_{q}R_{q}=R_{q}\Lambda_{q}$ (or similarly $L_{q}M_{q}=\Lambda_{q}L_{q}$), wherein $R_{q}$ ($L_{q}$) is the right (left) eigenmatrix, we find the eigenvalues to be

where $a_{q}=B_{q}/\sqrt{\mu_{0}\rho}$ is the familiar $q$-directed Alfvén speed, and $c_{\text{s},q}$ and $c_{\text{f},q}$ are the slow and fast speeds, respectively, given by

and $c=\sqrt{\gamma p/\rho}$ is the (isotropic) adiabatic sound speed. The corresponding right and left eigenmatrices ($R_{q}$ and $L_{q}$) are provided in Appendix A.

Each eigenvalue in Equation 4 pertains to a unique mode of the system of Equation 1, and prescribes the group velocity of disturbances along the direction $q$. Therefore, given the convective nature of the first two eigenvalues (as implied by their respective group velocities being equal to the plasma velocity), they must correspond to (zero-frequency) pseudo-modes. This means that, so long as the plasma behaves precisely according to the set of Equation 1, they do not propagate any energy, and are merely there to ensure that the zero magnetic field divergence requirement, as well as the plasma properties defining the entropy (which is constant under our polytropic assumption) are identically preserved and convectively transported throughout the system. For reasons that will become more evident from the mathematics, we shall term these to be the field divergence and entropy pseudo-modes. On the other hand, as signified by the propagation velocities $\pm a_{q}$, $\pm c_{\text{s},q}$, and $\pm c_{\text{f},q}$ shifted relative to the bulk velocity $v_{q}$, the remaining six eigenvalues belong to the Alfvén, slow, and fast modes, respectively, with the $+$ and $-$ signs indicating the direction of propagation on the corresponding characteristic curve. Note that the characteristic curves for the Alfvén waves are simply the magnetic field lines, while the fast and slow waves can, in general, propagate obliquely.

We are now in a position to derive the mode-decomposed expression for the total energy $E_{\rm{tot}}$. For a non-gravitational ideal plasma, $E_{\rm{tot}}$ is made out of kinetic (kin), internal (int), and magnetic (mag) parts given by

Differentiating subsection 2.2 with respect to time yields

On the other hand, using the relation $M_{q}=R_{q}\Lambda_{q}L_{q}$, we may rewrite Equation 2 in the following format

where $\boldsymbol{\mathcal{L}}_{q}(\boldsymbol{x},\boldsymbol{\nabla},t)=L_{q}\Lambda_{q}\partial_{q}{\boldsymbol{P}}$ is the vector of characteristic derivatives in the direction $q$. The $\boldsymbol{\nabla}$ argument means that $\boldsymbol{\mathcal{L}}_{q}$ contain spatial derivatives, while it is devoid of time derivatives. Therefore, Equation 8 provides an important connection between the local rate of change of the state vector and the mode-decomposed quantity $R_{q}\boldsymbol{\mathcal{L}}_{q}$ that solely depends upon local spatial variations. It is the second term that encapsulates the physics of the system, and hence prescribes appropriate local transformations to be utilized in order to break down the total energy.

Equation 8 provides a method to decompose the rate of change of each of the state variables in terms of each of the eight modes of the system. Thus, substituting the time derivatives obtained from Equation 8 into the RHS of Equation 7, we may rewrite the rate of change of the total energy in terms of an 8-sum composed of terms representing the contribution of each mode:

in which $\boldsymbol{c}_{n}$ are the local vectors of proportionality formed by combining the coefficients of the time derivative terms in Equation 7 with elements of the $R_{q}$ matrices, and $\overline{\mathcal{L}}_{n}$ are ($3\times 8$) matrices whose columns are made out of the $\boldsymbol{\mathcal{L}}_{q}$ vectors, with $n$ representing the mode numbers according to: ($n=1,2$) the field divergence (div) and entropy (ent) pseudo-modes, ($n=3,4$) the forward ($+$) and reverse ($-$) Alfvén (A) modes, ($n=5,6$) the forward and reverse slow (s) modes, and ($n=7,8$) the forward and reverse fast (f) modes, respectively, where forward and reverse refer to the two directions along each characteristic curve. The second line of subsection 2.2 shows the mode-expanded summation in the above-mentioned order, with this summation being over the directions $q$ instead.

The individual terms in the summation above, henceforth referred to as the eigenenergy time derivatives, are given by

in which $\boldsymbol{a}$ is the Alfvén speed vector, and the subscript $\perp$ denotes 2-vectors comprised of the components along the axes transverse to the direction $q$ according to the familiar Cartesian right-hand permutations (i.e., if $q=y$, then $\boldsymbol{q}_{\perp}=\left(z,x\right)$). We emphasize that Equation 10 are exact, nonlinear, and analytical everywhere in an ideal plasma. Note the entirely magnetic nature of Equation 10a and 10c associated with the field-divergence pseudo-mode and the Alfvén mode, and the fully acoustic nature of Equation 10b corresponding to the entropy pseudo-mode, as one would expect. Moreover, when summed over $q$ and time-integrated, Equation 10a provides the energy version of the constraint on the magnetic field to ensure that it remains divergence free, while Equation 10b provides such a constraint on our assumed thermodynamical nature of the gas to remain polytropic. Of course, for numerical simulations, such constraints are not always precisely met at all times and locations and there will be deviations, due to resistivity, shock heating, etc.. As we will later see, this will cause some amounts of energy to be carried by the pseudo-modes, despite our analytically deduced statement that they should not contribute to the overall energy transport under the strictly exact ideal-MHD framework.

Note the relationship between the $\pm$ signs in Equation 10c through to 2.2, where the first two $\pm$ signs on the RHS of Equation 10c are redundant, but purposely placed for clarity and consistency with the two subsequent equations. In all these cases, the reverse ($-$, $n=3,5,7$) and forward ($+$, $n=4,6,8$) expressions are simply connected by flipping the sign of the entire equation in addition to flipping the sign of all the velocity components and their derivatives on the RHS. This further provides insight on the differences between the forward and reverse modes, indicating that they indeed describe the direction of the flow of energy along the characteristics.

In light of Equation 10, it is noteworthy that subsection 2.2 can also be recast as follows

where $\boldsymbol{\lambda}_{q}$ and $\boldsymbol{F}_{q}$ are 8-vectors comprised, respectively, out of the eight eigenvalues in each direction and their corresponding coefficients found in Equation 10. The entries of $\boldsymbol{F}_{q}$ will all then share the same dimension of energy density per unit length, or simply force density.

The assumption that the solutions are known a priori means that the quantities in subsection 2.2 are also known and quantifiable everywhere in the plasma. Thus, we postulate that the total energy itself can be broken down in terms of the individual energies carried by each ideal-MHD mode, derived by performing time integration of each of the terms in that equation according to

where $E_{\rm{rec}}$ (given by the terms in the summation’s brackets) is the net amount of pure propagating wave-energy of the dynamic ideal plasma as recovered by the decomposition method, and $f$ is the function of integration corresponding to the static background energy. subsection 2.2, together with Equation 10 form the core of the eigenenergy decomposition method (EEDM).

As a first sanity test of the decomposition method formulated here, we refer the reader to Appendix B, where we consider pure magnetoacoustic waves, and illustrate that the method correctly attributes no energy to the Alfvén branch.

subsection 2.2 gives an expression for the total energy in terms of the component wave energies. Analytically, this should match the total energy as calculated from primitive code variables (Eq. 2.2). However, this may not be the case in a numerical implementation such as the one described below (due to numerical differentiation and integration in the decomposition method, see below and Appendix E). In order to assess the fidelity of the decomposition we write the energy obtained from primitive code variables (Eq. 2.2) as $E_{\text{tot}}=E_{\text{orig}}+E_{0}$, with $E_{0}$ being the energy of the static initial (or background) state, and $E_{\text{orig}}$ being the (“original”) energy of the fluctuating (wave) component – obtained by subtracting the total energy from the initial total energy, using primitive code variables. When presenting our results in the subsequent sections, we will test the numerical fidelity of the decomposition method by providing comparisons between the net wave energies pre- and post-decomposition given by $E_{\rm orig}$ and $E_{\rm rec}$, respectively.

Due to the mentioned lack of 3D nonlinear solutions, and in order to provide more insight into the method, we will have to resort to direct numerical simulations that are representative of sufficiently simple conditions to be able to interpret and substantiate the results. In the two following sections, we will discuss two ideal-MHD simulation setups along with their results. These are made of a uniform plasma superimposed by (I) a vertical magnetic field and excited by a circularly polarized driver aimed at primarily generating upward travelling Alfvén waves, and for comparison with previous work, (II) the 3D structured magnetic field of Galsgaard et al. (2003) containing a null point at the center and excited by the circular driver described therein, again aimed at maximising the generation of Alfvén waves. Due to the core differences between the two simulation setups, we will fully describe each case and discuss its results individually in the two separate sections that follow.

The mentioned drivers are introduced as boundary conditions (BCs), and once again we remind the reader that these do not constitute exact ideal-MHD solutions (but rather numerical ones), and hence some discrepancy is expected with the theoretical exposition above (explored below). This means that Alfvén waves will not be the only product of such drivers. Note that the purpose of the simulations introduced here is to verify the capacity of the decomposition method to identify wave modes, rather than to model any specific physical process in the solar atmosphere or elsewhere.

The numerical package used throughout this work is Lare3d (Arber et al., 2001), which is second-order accurate in space and time. The code is built on a staggered grid, and at each step advances the solutions to the dimensionless Lagrangian 3D MHD equations and then conservatively remaps them back onto the original grid using volume averaging techniques. So far we introduced and formulated the method using the dimensional MHD equations for the sake of more clarity. However, in view of the dimensionless nature of the Lare3d solver used for our simulations, for the remainder of this study we shall work with non-dimensional units. This is simply done by setting $\mu_{0}=1$ throughout our equations and treating all the other variables as non-dimensional.

For the initial state we consider a uniform plasma threaded by a uniform magnetic field of non-dimensional strength $B_{0}=0.1$ in the $z$-direction. In simulation Ia, the uniform density and temperature in the initial state are chosen as $\rho_{0}=0.0198$, and $T_{0}=2.2748$, giving a plasma-$\beta$ of $9.03$. In Ib, $\rho_{0}=2.4614\times 10^{-5}$, $T_{0}=111.21$, and $\beta=0.54$. The remaining initial conditions on the pressure and internal energy are then easily found using the ideal gas equation. Lastly, to ensure the system is initially at rest, we set $\boldsymbol{v}_{0}=0$. We note that these parameter combinations are chosen to match two different heights from a plane-stratified model of the solar atmosphere (to be considered in a future paper).

Given the symmetry of the model, we set the horizontal BCs to periodic, while choosing fully symmetric boundaries at the top plane. For the bottom boundary plane, we wish to set up simple torsional Alfvén conditions that incompressively excite travelling waves at $z=0$. The goal is to primarily generate travelling Alfvén waves. However, as there are no general exact nonlinear solutions to the ideal-MHD Equation 1 that would guarantee pure Alfvén waves (except for the semi-trivial case discussed in Appendix D), we may seek approximate linear solutions of this kind, with the caveat that when used in the nonlinear regime, they will excite other waves as well.

Turkmani & Torkelsson (2004) present 1.5D WKB solutions describing the approximate linear behaviour of “circularly polarized” torsional Alfvén waves in a gravitationally stratified isothermal plasma superimposed with a uniform vertical magnetic field $\boldsymbol{B}=B_{0}\hat{\bf z}$. For the purpose of this study under the uniform plasma assumption, we shall do away with the stratification and modify their solutions accordingly as follows

where $a_{0}=B_{0}/\sqrt{\rho_{0}}$ is the initial vertical component of the Alfvén speed, $k=2\pi/5$ is the wavenumber, $\omega$ is the Alfvén frequency, $\varphi=\pi/4$ sets the phase, and the subscript $\perp$ represents the directions perpendicular to the direction of propagation, i.e., the $z$-axis. To configure the driver to produce travelling waves moving initially at the Alfvén speed, we may localise the amplitude of the field fluctuations in space and time by choosing

in which $A$ is a constant (the maximum amplitude, set to $2\times 10^{-4}$ in simulation Ia and $10^{-2}$ in simulation Ib), $z_{0}=-4.1556$ is the initial location of the pulse for a smooth introduction of the disturbances into the main domain, and $\sigma=0.5$ sets the width of the pulse. Finally, setting $v_{z}=0$ completes setting up our desired torsional Alfvénic driver in the ghost cells. It should be noted that the described driver does not satisfy the $z$-component of the full momentum equation in the main simulation domain, which means that a component $v_{z}$ is generated in the first internal grid cell to balance this deficiency—see Appendix C for a further discussion.

Before proceeding with presenting the simulation results, we note that the difficulties in finding an exact solution to the nonlinear MHD equations arise from the space and time localization of the Alfvén wave. In Appendix D we show that a monochromatic Alfvén wave, constructed by taking the limit of Equation 14 as $\sigma\to\infty$, does produce a semi-trivial solution of the MHD equations. This is equivalent to setting $B_{\perp}$ to a constant. The eigenenergy decomposition method can be analytically verified in that case, thus providing another sanity check on the validity of the method.

We shall now individually present the results of the first family of simulations (Ia and Ib). The general formats in which the figures are presented are standardized to remain the same across this and the next section, for better tractability and readability, unless stated otherwise.

The MHS model atmosphere of simulation II is comprised of a uniform plasma, superimposed with a 3D structured potential magnetic field including a null at the center of the box at $x=y=0$ and $z=1$, given by

where $B_{0}=0.1$. Note that the $z$ component of the Alfvén speed is zero across the fan at $z=1$, indicating that pure Alfvén waves cannot transmit through this plane. The thermodynamics of the system are prescribed by a uniform density and temperature $\rho_{0}=25$ and $T_{0}=0.1$ (the pressure and internal energy then follow from the ideal gas law). This sets the initial sound speed at $c_{0}\approx 0.4$. Furthermore, given the nature of the magnetic field, the model atmosphere will have only one equipartition contour enclosing the null point, where the sound and Alfvén speeds coincide. It has a radius of 0.106 in the $y=0$ plane. Inside this contour $\beta>1$ and is formally infinite at the null (numerically, $\beta_{\text{max}}\approx 1230$ at cell centers), while outside of it $\beta<1$, with a minimum of $\beta_{\text{min}}\approx 0.0125$. The magnetic field increasingly dominates the dynamics for decreasing $\beta$. Outside the equipartition layer, fast (slow) waves are magnetically (acoustically) dominated, and vice versa inside. The non-dimensionalization factors for the Lare simulation are given by $L_{N}=150$ km, $\rho_{N}=2.007\times 10^{-11}$ kg.m^-3, $T_{N}=10^{6}$ K.

The BCs are set to fully symmetric, i.e., $\partial_{q}{v_{q}}=0$, $\partial_{q}{B_{q}}=0$, and $\partial_{q}{Q}=0$, where $Q$ is any of the thermodynamical variables, on all faces of the simulation box, except for the bottom boundary plane where we implement another version of the torsional Alfvénic driver suitable for handling the now variable magnetic field. To formulate the driver, this time we start off by prescribing a 2D incompressive flow velocity profile, which is a slightly modified version to that of Galsgaard et al. (2003), given by

in which $t_{0}=0.581$ is the peak time of the Alfvénic pulse, $\sigma_{R}=0.1$ sets the width of the pulse in the $xy$ plane, $\sigma_{\text{P}}=0.07$ (in combination with the Alfvén speed) sets the width of the propagating pulse along the $z$-axis, and $v_{0}=1$ is the amplitude of the driver. In a similar fashion as before, the desired BCs on the magnetic field are then set by

This ensures that the field is initially line-tied, and incompressive perturbations are introduced in keeping with the velocity field according to the WKB solutions of Turkmani & Torkelsson (2004). Finally, the BCs on the thermodynamics of the system are fixed to their respective initial values.

It is worth emphasising that these BCs do not strictly satisfy the full nonlinear MHD equations, which will result in the type of boundary layer common in MHD simulations with driven boundaries. This has a (relatively minor) impact on how we interpret the energy decomposition results, described further below, and in Appendix E.

Figure 7 depicts the directionally summed (over both spatial and characteristic directions) branches of the physical and pseudo energies of simulation II, again calculated by numerically integrating the eigenenergy time derivatives over time. In this figure, panels (a) through to (f) to the left of the vertical separator line exhibit the net amounts of the slow ($E_{\text{s}}$), field-divergence ($E_{\text{div}}$), Alfvén ($E_{\text{A}}$), entropy ($E_{\text{ent}}$), fast ($E_{\text{f}}$), and the recovered ($E_{\text{rec}}$) energies and pseudo-energies, respectively. The panels on the right of the separator represent the same quantities as before (as per Figure 2), and are calculated directly from the single simulation time step. We chose this more compact $q$- and characteristic direction-summed representation of the results to capture the most important aspects of each of the branches, as there are cancellations between the different directional components of each branch due to the azimuthal symmetry of the simulation. Note how all the energy branches lie entirely in the vicinity of the driver at the bottom boundary, implying that the rotational driver does not induce any flow of energy in the horizontal directions: the energy flux is predominantly upwards along the field lines. The only added feature in this figure is the vertical dashed line in panel (j), which marks the $x$-coordinate where we slice the data further to inspect the 1D behaviour of the energies, as shown in Figure 8.

Contrary to simulations Ia and Ib, panel (c) of Figure 7 shows significant amplitudes associated with the Alfvén energy. Note how the Alfvén energy spreads out across the fan plane at $z=1$, but cannot transmit through. This is the expected behaviour since $a_{z}=0$ across this plane, meaning that it acts as a shield against the upward propagation of Alfvén waves. However, upon interacting with the equipartition contour surrounding the null, the Alfvén waves can partially convert to magnetoacoustic waves which in turn find their way into the upper half above the fan (see the features above the fan in panels a and e). Furthermore, note how the two horn-like features in panel (b) are aligned with the field lines, as one would expect from pure Alfvén waves. The footpoints of these features are rooted in the vicinity of the maxima of the spatial Gaussian function in Equation 16, around the $x^{2}+y^{2}-\sigma_{R}^{2}=0$ circle. In view of the incompressibility of Alfvén waves and judging by the zero Alfvén-energy in the white regions of this panel, it appears that the driver must have triggered some thermal (and magnetic) pressure gradients in the volume adjacent to the spine ($z$-axis). This means that some amounts of magnetoacoustic energy are likely to have been generated and channeled upwards through the ‘tube’ formed in the middle of the Alfvén energy features, which is correctly seen to be the case in panel (e) belonging to the fast energy.

In the low $\beta$ region (away from the equipartition contour) below the fan, the Alfvén and fast modes are near-degenerate, yet have distinct variations to the MHD state vector $\boldsymbol{P}$ (compare, e.g., the left and right eigenvectors for the Alfvén and fast modes as given in Appendix A). The near-degeneracy is indicated by the similar group velocities (blue and green dots), as well as the coincidence of the horn-like features in panels (c) and (e), which further unveils that up to a certain height below the equipartition contour, the fast energy is also propagating along the field lines. These two modes slightly alter the background state as they are introduced by the driver and propagate through the system. In the EEDM, these appear as a decrease in $E_{\text{f}}$ and increase in $E_{\text{A}}$ (both relative to the background) below the fan surface. The EEDM accurately captures the properties of two simultaneously propagating wave systems, as given by the fact that the energy in the Alfvén mode cannot and does not cross the fan surface at $z=1$, while the fast mode can and, indeed, does.

As the pulse nears the fan, all the energy branches show some contribution to the disturbance that follows the field lines, spreading out against the fan plane from underneath. At the same time, as shown in panel (e), a substantial amount of the energy begins to wrap around the null and transmit through to the upper region (as observed by, e.g., McLaughlin & Hood, 2004; Galsgaard et al., 2003). This fast mode is in part generated by the boundary driving, and in part by mode conversion at the equipartition contour shown in panel (h). This is exactly what was previously observed by Galsgaard et al. (2003). Examining panel (e) in more detail, the fast energy can be broken down into three parts: (i) the vertical positive-energy orange beam extending along the spine up to the null, (ii) the purple negative-energy features in the driven sheath, and (iii) the V-shaped structure above the null. As touched on before, (i) is the result of thermal pressure as well as magnetic pressure fluctuations (triggered, at least in part, by the induced vertical velocity described in Appendix C) inside the cavity confined by the Alfvén energy (in panel c). Note that this feature is positive in value and does not coincide with any other features from other energy branches (except for the faint vertical blue beam in panel b, which is negligible), and therefore appears to be a pure fast mode confined near the spine. Thus, it is at least partially responsible for the transmission of fast energy into the upper half in the form of the V-shape above the fan. The purple features (ii), are negative in value, which simply means that this part of the fast wave is extracting energy from the background due to the nonlinear nature of the evolution changing the background state. Notice how (ii) can be further broken down to the lower dimmer purple part extending from the bottom boundary up to $z\approx 0.25$, and the upper darker part. While the darker portion is the actual propagating pulse itself, the dimmer part marks its wake, where nonlinear background disturbances are excited by the passing of the pulse through the plasma. Such nonlinear disturbances lead to numerical heating, and hence generate the (very small) entropy pseudo-energy in panel (d). Due to its coincidence and similarity (both in the dimmer and darker regions) with $E_{\text{div}}$ (b), it appears that the mentioned extracted fast energy is somehow linked to the observed field-divergence pseudo-energy. We speculate that perhaps the field-divergence pseudo-energy, which is the result of some numerical pitfalls (see Appendix E), is essentially leaking from the fast energy. We will have a closer look at this speculation shortly. Lastly, (iii) is the net result of the transmission of (i), as well as mode conversions from $E_{\text{A}}$ and/or $E_{\text{div}}$.

The wrapping of the fast mode around the null is mainly due to the variable fast speed radially away from the spine, where at any given height, it is at its minimum at the location of the spine ($x=y=0$), and increases hyperbolically with $x$ and $y$. Consequently, the fast waves, either directly transmitted or converted from Alfvén waves, that are near the spine will refract in towards the null and criss-cross above it in higher energy concentrations forming the maximum of the fast energy right above the null (the sharp bright orange vertex of the V-shape), while the waves farther from the spine move at higher speed and refract in to locations higher above the fan (the wings of the V-shape). One can take on a WKB approach and think about the refraction pattern in terms of rays intersecting in higher densities right above the null and rarefying in the wings.

There is a rather minimal energy contribution from slow waves (panel a). A faint signal does propagate upward from the lower boundary at the background slow wave speed (black/red dot), providing further validation of the decomposition. In addition, there is a second set of features evident in $E_{\text{s}}$ that appear around the fan surface after the Alfvén and fast modes reach that height. This is also seen in the red curve of Figure 8, inside the equipartition contour. In the time animations this second feature does not propagate all the way from the boundary—but rather appears following the arrival of the fast/Alfvén pulse in the vicinity of the equipartition surface—indicating that this feature in the slow mode energy is generated by mode conversion of the cospatial Alfvén and fast mode pulses. These slow modes were not identified in Galsgaard et al. (2003), but appear analogous to the converted slow modes exiting a null point identified in 2D by Tarr et al. (2017) and in 3D by Thurgood et al. (2017), after some effort; they are identified trivially by the present decomposition scheme.

One important caveat to the above results is illustrated in Figure 7(b), which shows that significant energy appears in the field-divergence pseudo-mode, and (to a much smaller extent) the entropy pseudo-mode in panel (d). We remind the reader that the very presence of any non-zero pseudo-energies ultimately boils down to the detection of non-adiabaticity (departure from the polytropic gas assumption) and/or a non-zero field divergence as computed by the EEDM. Thus, to correctly interpret the results, one must take heed in identifying any sources of differences between how the simulating code and the decomposition scheme numerically evaluate these two quantities, i.e., how the spatial partial derivatives are calculated. This is expounded in detail in Appendix E. But in short, in order to numerically construct the terms in Equation 10, we have first interpolated all physical variables (i.e., the plasma state vector $\boldsymbol{{P}}$) to the same grid locations (due to the staggered stencil used by Lare), and then evaluated the required spatial derivatives using a cubic B-spline scheme. This means that there is a mismatch between the derivative values used in the decomposition scheme and those same derivatives represented in Lare’s native numerical scheme. This is particularly illustrated by the non-zero value of $\mbox{\boldmath$\nabla$}{\mathbf{\cdot}}{{\boldsymbol{\mathit{B}}}}$ in Figure 7, which is in contrast to the native zero value (to machine precision) of $\mbox{\boldmath$\nabla$}{\mathbf{\cdot}}{{\boldsymbol{\mathit{B}}}}$ in Lare. We stress that the high-order differentiation methods were found to be required to accurately recover the original wave-energy (Figure 8). However, this comes at the cost of a leaky partitioning of energy between the modes in certain circumstances that remain to be explored fully. This illustrates the care that should be taken in implementing the decomposition scheme and interpreting its results. Additionally, recall that the quantities shown to the left of the vertical separator in the figure are integrated over time, whereas those in the rightmost column are computed instantaneously for the given snapshot. This implies that even low amounts of $\mbox{\boldmath$\nabla$}{\mathbf{\cdot}}{{\boldsymbol{\mathit{B}}}}$ or non-adiabaticity in the decomposition scheme in individual snapshots can potentially generate high values of pseudo-mode energies, due the time integration accumulating effects from all the previous snapshots.

Figure 8 depicts 1D plots of the various energy-decomposed branches of simulation II along a line in $z$ at $x=0.107$, $y=0.003$, which is signified by the vertical dashed line in Figure 7(j), and $t=1.335$. The layout of the panels follows from Figure 3. However, for more clarity and in view of the mentioned similarities between the $E_{\text{f}}$ and $E_{\text{div}}$ profiles, we have now added the dashed brown curve which represents $E_{\text{f}}+E_{\text{div}}$. Note the rather small mismatch between the original and recovered energies in the bottom panel over $0\leqslant z\lessapprox 0.25$. This is discussed in detail in Appendix E. But in short, it stems from discontinuities introduced by the driver (seen as the cusp point near $z=0$ in the blue curve belonging to $E_{\text{orig}}$), owing to the lack of complete nonlinear ideal-MHD solutions discussed before. As a result, the energy allocations performed by the method are slightly off in this region. Nevertheless, the actual pulse is located above this height, where the original energy recovery is precise, making the decomposition scheme valid inside the propagating wave.

Notice the near anti-symmetry in trend (not magnitude) of the blue ($E_{\text{f}}$) and the magenta ($E_{\text{div}}$) curves in the left side of the upper frame up to the first equipartition level. This is more clearly observed in the dashed brown curve plotting their net effect $E_{\text{f}}$ + $E_{\text{div}}$, which exhibits a similar trend to that of the Alfvén wave. This provides further evidence on the speculation touched on above that the detected $E_{\text{div}}$ might have been leaked from the fast branch, and that had the field-divergence been numerically zero as measured by the EEDM, the actual fast energy branch would have been the brown dashed curve belonging to $E_{\text{f}}$ + $E_{\text{div}}$. However, further investigation is required to substantiate this conjecture. Overall, comparing the two panels from $z\approx 0.25$ up to the left vertical dashed line, we can see that the total positive propagating wave-energy supplied by the driver at the particular point $x=0.107$ and $y=0.003$ is mainly being carried by the Alfvén and slow modes. However, as the wave approaches the equipartition contour, mode conversion switches the wave-energy to be predominantly carried by the fast branch. Moreover, in the time evolution of this plot, it is seen in our animations that as the fast branch hits the upper equipartition level (right vertical dashed line), it seems to start propagating faster than the background fast speed (blue dot) at first glance. However, this is only due to fast rays from elsewhere having refracted and found their way into this location.

Figure 9 shows the time-distance ($t$-$z$) plots of the $q$- and characteristic direction-summed energy components decomposed by the method, taken at the same location as marked by the vertical dashed line in Figure 7(j), i.e., at $x=0.107$ and $y=0.003$. We remind the reader that the local slope of the distinctive filamentary features are given by the projection of the local speed of the corresponding wave branch onto the spatial cut. For instance, notice the constant speed of propagation of the slow branch in panel (a), as indicated by the straight filaments at the top and bottom of the panel, due to the nearly constant slow wave speed along this line, except near the equipartition layer. This is in turn due to the constant background sound speed on the one hand, and the acoustically dominated nature of the slow waves away from the null on the other hand. In contrast, the fast mode speed increases away from the null along this cut, causing a telltale bending of the energy density curve in panel (e). The more smeared out features in all panels indicate refracted waves originally launched from other locations crossing paths with the vertical dashed line in Figure 7(j). The Alfvén energy depicted in panel (c) is seen to be continually blocked by the fan plane at $z=1$, as expected. The fast energy in panel (e) is comprised of (i) an upward propagating positive part seen as the orange filament, and (ii) the purple regions, which manifest the near anti-symmetry with $E_{\text{div}}$ (of panel b) discussed above. Finally, the discontinuous slow energy features above the equipartition contour in panel (a) imply that the resultant slow wave must have been produced via mode transformations from other branches.