A&A 369, 1122-1139 (2001)
DOI: 10.1051/0004-6361:20010171

Numerical simulations of linear magnetohydrodynamic waves in two-dimensional force-free magnetic fields

I. Arregui - R. Oliver - J. L. Ballester

Departament de Física, Universitat de les Illes Balears, 07071 Palma de Mallorca, Spain

Received 4 August 2000 / Accepted 26 January 2001

Abstract
High resolution observations of the solar corona made with instruments onboard the SOHO and TRACE spacecrafts have provided new evidence for the presence of oscillations in a variety of coronal magnetic structures. Most of these observations have been interpreted in terms of linear standing or propagating magnetohydrodynamic (MHD) waves, but the theoretical understanding is far from complete since analytical solutions to the linearised MHD wave equations can only be found for very simple magnetic configurations. Taking into account that the solar corona is basically structured by force-free magnetic fields, our purpose in this paper is to present the derivation of the linear MHD wave equations for a two-dimensional force-free magnetic field configuration having longitudinal invariance, as well as to introduce a numerical code to solve the resulting system of coupled partial differential equations. The accuracy of the code has been checked by numerically solving two cases for which analytical or simple numerical solutions exist. To our knowledge, this is the only two-dimensional code developed to study the normal MHD modes of oscillation of a general force-free field with longitudinal invariance.

Key words: magnetohydrodynamics - methods: numerical - Sun: oscillations - Sun: magnetic fields

  
1 Introduction

Multiwavelength observations of the solar corona, obtained by both ground-based telescopes and instruments onboard spacecrafts, have revealed the high complexity of its structuring by magnetic fields. These observations show a variety of coronal structures such as loops and arcades, where the magnetic field is closed, or plumes and coronal holes, where the magnetic field is open.

Coronal seismology is becoming an important field of research within solar coronal physics and firm evidence about the existence of oscillations in different coronal structures has been provided by instruments onboard SOHO and TRACE spacecrafts. For instance, DeForest & Gurman (1998) have reported observations, made with EUV on board SOHO, of quasi-periodic perturbations in the brightness of FeIX and FeX emission lines in coronal plumes. The perturbations propagate outward at 75-150 km s^-1 with periods of 10-15 min. These compressive perturbations were interpreted as slow magnetoacoustic waves by Ofman et al. (1999, 2000a). EUV brightenings in active region coronal loops have been observed with TRACE by Nightingale et al. (1999), who have also interpreted them as compressional waves. Robbrecht et al. (1999) using simultaneously EIT, onboard SOHO, and TRACE, have pointed out the existence of weak perturbations in loops, originating at the footpoints, propagating along them with an apparent speed of 130 km s^-1. They suggest that these disturbances are slow magnetoacoustic waves propagating in coronal loops. Aschwanden et al. (1999), using TRACE, detected spatial oscillations of coronal loops with a transverse amplitude of about 4000 km and a mean period of 280 s. They suggest that the oscillation appears to be a standing wave with fixed nodes at the footpoints and that a fast kink mode seems to provide the best agreement with observations. Moreover, the damping of the oscillation has been explained by dissipation through viscous or resistive processes (Nakariakov et al. 1999; Roberts 2000) or by means of a mechanism based on the topology of magnetic field lines (Schrijver & Brown 2000). De Moortel et al. (2000) have also reported TRACE observations of coronal loop oscillations having periods of 180-420 s and outward propagating velocities of 70-165 kms^-1. They have explained these observations in terms of slow magnetoacoustic waves propagating along the loops. Ofman et al. (1997, 2000b), using UVCS onboard SOHO, have detected quasi-periodic variations in the polarised brightness in polar coronal holes. Again, these observations have been explained in terms of density fluctuations caused by compressional waves propagating in polar coronal holes.

Most of the above observations, and many others reported in previous works, can be interpreted in terms of linear MHD waves and, from this point of view, investigations of adiabatic standing MHD waves in coronal magnetic configurations, such as coronal magnetic arcades, have been undertaken (e.g. Oliver et al. 1993, 1996; Terradas et al. 1999). Existing theoretical studies indicate that analytical solutions to the linearised wave equations can only be found for very simple configurations and that in most cases it is necessary to resort to numerical techniques. For this reason, we have created a numerical code to investigate the MHD modes of a general two-dimensional force-free magnetic configuration with longitudinal invariance. There are no restrictions on the equilibrium structure other than longitudinal invariance and the zero-$\beta$ assumption, which is considered reasonable when modelling coronal structures. By setting $\beta=0$ the slow mode is removed and only the fast and Alfvén modes remain, and although this may appear as a restriction it will allow us to better understand the interplay between these two modes. Our future plans are to exclude the cold plasma assumption and to modify the code to include plasma pressure terms. To our knowledge, the present code is the only one freely available to study the (normal and continuum) modes of two-dimensional configurations.

In this paper, we first present the derivation of the linear MHD wave equations, expressed in terms of the velocity components, for a two-dimensional force-free coronal magnetic configuration (Sect. 2) and next introduce the numerical code to solve the resulting system of coupled partial differential equations (Sect. 3). The user need only supply analytical expressions for the Cartesian magnetic field components and for their first- and second-order derivatives with respect to the Cartesian coordinates. Hence, equilibrium configurations such as twisted cylindrical magnetic fields with invariance along the z-axis can also be investigated. In addition, the user has the option of solving the wave equations using non-Cartesian coordinates by providing analytical expressions for the new coordinates in terms of the Cartesian ones together with their first- and second-order derivatives. The goodness of the numerical code has been checked by numerically computing the normal modes of two configurations (a straight, uniform magnetic field and a potential arcade) for which analytical or simple numerical solutions exist (Sect. 4). Finally, in Sect. 5 conclusions are drawn.

  
2 Basic equations

  
2.1 Equilibrium configuration

We assume that our coronal structure may be modelled as a force-free magnetic field satisfying

$$% \frac{1}{\mu}\left(\nabla\times{\vec B}\right)\times{\vec B}=0,$$ (1)

complemented with the divergence-free field condition

$$% \nabla\cdot{\vec B}=0.$$ (2)

Since the Lorentz force vanishes in the equilibrium, the plasma current, ${\vec j}=\nabla\times{\vec B}/\mu$, is parallel to the equilibrium field, ${\vec B}$, and the magnetic forces (tension and gradient of magnetic pressure) are mutually balanced.

Let us consider a two-dimensional equilibrium in which the plasma and field variables do not vary along the y-axis of a Cartesian frame x, y, z. Then, the equilibrium magnetic field can be written as

$$% {\vec B}=\nabla A\times\mbox{$\hat{{\vec e}}_y$ }+B_{y}\mb... ...ox{$\hat{{\vec e}}_{\rm p}$ }+B_{y}\mbox{$\hat{{\vec e}}_y$ },$$ (3)

with a component in the longitudinal ( $\mbox{$\hat{{\vec e}}_y$ }$) direction and a component in the poloidal ( $\mbox{$\hat{{\vec e}}_{\rm p}$ }$) direction, contained in the xz-plane and described by a flux function $A{\mbox{$\left(x,z\right)$ }}$. The force-free condition (1) dictates that $\mbox{$B_{y}$ }(x,z)$ is a strict function of $A{\mbox{$\left(x,z\right)$ }}$ and that the two are related by

$$% \mbox{$\nabla^2 A$ }+\mbox{$B_{y}$ }\frac{{\rm d} \mbox{$B_{y}$ }}{{\rm d}A}=0,$$ (4)

which is the well-known Grad-Shafranov equation for this type of equilibrium and determines the flux function, $A{\mbox{$\left(x,z\right)$ }}$, once the function $\mbox{$B_{y}$ }(A)$ has been prescribed. Analytical solutions to Eq. (4) are known for particular forms of $\mbox{$B_{y}$ }(A)$, e.g. $\mbox{$B_{y}$ }(A)=constant$ (potential field), $\mbox{$B_{y}$ }(A)=cA$ (constant-$\alpha$ field) and $\mbox{$B_{y}$ }(A)=cA^{\frac{1}{2}}$ (constant-current field).

Since we are neglecting gas pressure and gravity terms in the force-balance equation, the density does not appear in Eq. (4) and one can impose the equilibrium density or Alfvén speed profiles arbitrarily. In our general two-dimensional equilibrium the density will be given by a function of the form $\rho=\rho{\mbox{$\left(x,z\right)$ }}$, which defines the Alfvén speed by means of $\mbox{$v^2_{\rm A}$ }{\mbox{$\left(x,z\right)$ }}=B^2/\mu\rho$.

2.2 Linear and adiabatic wave equations

In this section the equations of linear, ideal MHD waves are derived for our general, two-dimensional force-free equilibrium. Let us start with the ideal, cold plasma MHD equations for adiabatic changes of state,

   

$$\frac{{\rm D}\rho}{{\rm D}t}+\rho\nabla\cdot{\vec v}=0,$$ (5)

$$\rho\frac{{\rm D}{\vec v}}{{\rm D}t}=\frac{1}{\mu}\left(\nabla\times{\vec B}\right)\times{\vec B},$$ (6)

$$\frac{\partial{\vec B}}{\partial t}=\nabla\times\left({\vec v}\times{\vec B}\right),$$ (7)

with

$$% \frac{\rm D}{{\rm D} t}=\frac{\partial}{\partial t}+{\vec v}\cdot\nabla.$$ (8)

Equations (5)-(7) are, respectively, the continuity equation, the momentum equation for a cold, nonviscous plasma and the induction equation.

If we now consider small displacements about the equilibrium state, and assuming time variations of the form $\exp(i\omega t)$, with $\omega$ the frequency, Eqs. (5)-(7) can be written in the following form

   

$$% i\omega\rho_{1} -\rho\nabla\cdot{\vec v}-{\vec v}\cdot\nabla\rho,$$ = (9)

$$i\omega\rho{\vec v} \frac{1}{\mu}\left(\nabla\times{\vec B}\right)\times{\vec B}_{1}+ \frac{1}{\mu}\left(\nabla\times{\vec B}_{1}\right)\times{\vec B},$$ = (10)

$$i\omega{\vec B}_{1} \nabla\times\left({\vec v}\times{\vec B}\right),$$ = (11)

with ${\vec v}$ the plasma velocity, $\rho$ and ${\vec B}$ the equilibrium density and magnetic field and $\rho_{1}$ and ${\vec B}_{1}$ the density and magnetic field perturbations. A single expression for ${\vec v}$ can be derived by eliminating ${\vec B}_{1}$ in Eq. (10) with the help of Eq. (11),

 

$$% -\rho\omega^2{\vec v} \frac{1}{\mu}\left(\nabla\times{\vec B}\right)\times\left\{\nabla... ...[\nabla\times \left({\vec v}\times{\vec B}\right)\right]\right\}\times{\vec B}.$$ (12)

Once ${\vec v}$ is obtained from this equation, the perturbations $\rho_{1}$ and ${\vec B}_{1}$ can be calculated from Eqs. (9) and (11) and some properties of the modes can be investigated. Further information on the nature of waves comes from the magnetic pressure,

$$% p_{1m}=\frac{{\vec B}\cdot{\vec B}_{1}}{\mu}$$ (13)

and from the perturbation to the Lorentz force (i.e. the vector on the right-hand side of Eq. (10)), which is responsible for driving the oscillatory plasma motions. This force is the sum of the gradient of the perturbation to the magnetic pressure ( $\nabla p_{\rm 1m}$) and the perturbation to the magnetic tension,

$$% \mbox{${\bf {\cal T}}_1$ }=\frac{1}{\mu}\mbox{$\left({\vec... ...{1}+\frac{1}{\mu}\left({\vec B}_{1}\cdot\nabla\right){\vec B}.$$ (14)

Equation (12) is a vector formula that actually corresponds to three second-order partial differential equations for the three velocity components. Here, the different vectors (and not just ${\vec v}$) have been decomposed in their components by projecting them onto the directions defined at each point by the unit vectors, $\mbox{$\hat{{\vec e}}_n$ }=\mbox{$\nabla A$ }/\vert\mbox{$\nabla A$ }\vert$ (normal to magnetic surfaces), $\mbox{$\hat{{\vec e}}_{\Vert}$ }={\vec B}/B$ (parallel to field lines) and $\mbox{$\hat{{\vec e}}_\perp$ }=\mbox{$\hat{{\vec e}}_{\Vert}$ }\times\mbox{$\hat{{\vec e}}_n$ }$ (tangent to magnetic surfaces and perpendicular to field lines). A representation of these three basis vectors for a particular magnetic configuration can be seen in Fig. 1. Then, any vector, and in particular the velocity perturbation, can be written as

$$% {\vec v}=\mbox{$v_n$ }\mbox{$\hat{{\vec e}}_n$ }+\mbox{$v_\... ..._\perp$ }+\mbox{$v_{\Vert}$ }\mbox{$\hat{{\vec e}}_{\Vert}$ }.$$ (15)

Since the equilibrium quantities do not depend on the y-coordinate, we can Fourier analyse the perturbed quantities with respect to this coordinate and make them proportional to $\exp(ik_y y)$. After some manipulations (see Appendix A for a detailed derivation), the parallel component of Eq. (12) yields $\mbox{$v_{\Vert}$ }=0$ and the normal and perpendicular components can be cast as (note the inclusion of the factor $i\omega/\left(\mbox{$v^2_{\rm A}$ }\rho\right)$ on both sides of these expressions),

 

$$% -\frac{\omega^2}{\mbox{$v^2_{\rm A}$ }}\mbox{$v_n$ } \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cd... ... p}$ }}\right]\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$ =

$$+\left(C_{1}- \frac{\mbox{$B^2_{y}$ }}{B^2}k_y^2\right)\mbox{$v_n... ...ht)\right]\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_\perp$ }$$

$$-\left[\frac{1}{B}\frac{{\rm d}\mbox{$B_{y}$ }}{{\rm d}A}+\frac{\... ...box{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }} \right]ik_y\right\}\mbox{$v_\perp$ },$$ (16)

 

$$% -\frac{\omega^2}{\mbox{$v^2_{\rm A}$ }}\mbox{$v_\perp$ } -\frac{\mbox{$B_{y}$ }}{B\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\vec... ...}$ }}{B_p^4}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_n$ }$$ =

$$+\left(C_{3}+\frac{B^2\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}... ...t\!\nabla\right)$ }\mbox{$v_\perp$ }+\left(C_{4}-k_y^2\right)\mbox{$v_\perp$ },$$ (17)

where expressions for the coefficients C₁-C₄ and the quantities $\mbox{$B_{\rm s}$ }$, $\mbox{$B_{\rm ps}$ }$, etc. can be found in Appendix B. It can be appreciated that the derivatives in the xz-plane have been expressed in terms of the field-related operators $\nabla A\cdot\nabla$and ${\vec B}_{\rm p}\cdot\nabla$, which are simply the spatial derivatives across magnetic surfaces and along the poloidal magnetic field, respectively. Thus, the two non-zero components of the vector wave Eq. (12) reduce to a set of two second-order coupled partial differential equations (describing coupled fast and Alfvén modes) for the components of the perturbed velocity in the directions normal and perpendicular to the equilibrium magnetic field. As the equilibrium configuration is force-free, the gas pressure gradient plays no role and the slow mode is absent ( $\mbox{$v_{\Vert}$ }=0$).

	Figure 1: Three dimensional view of magnetic field lines (solid) in a force-free coronal arcade with uniform shear. Dashed lines are the projections of field lines onto the xy- and xz-planes. The field related vectors $\mbox{$\hat{{\vec e}}_n$ }$, $\mbox{$\hat{{\vec e}}_{\Vert}$ }$ and $\mbox{$\hat{{\vec e}}_\perp$ }$ are plotted at a particular point
Open with DEXTER

Once the velocity of the plasma motions has been calculated a task that is performed numerically as explained in Sect. 3, the perturbations to the other physical quantities and forces can be computed. These quantities are given in Appendix C.

The derivation of our main formule are now complete, but we can check to detect the possible presence of errors in the derivation and/or in the transcription of equations to the numerical code. First, take Eq. (10) and multiply it by $i\omega/\left(\mbox{$v^2_{\rm A}$ }\rho\right)$, such as was done before in the process of deriving Eqs. (16) and (17). After expressing the perturbed Lorentz force as the sum of the perturbed pressure gradient and magnetic tension we have

$$% -\frac{\omega^2}{\mbox{$v^2_{\rm A}$ }}{\vec v}=-i\frac{\om... ...{\omega}{\mbox{$v^2_{\rm A}$ }\rho}\mbox{${\bf {\cal T}}_1$ }.$$ (18)

By considering the normal component of this expression, one can confirm that the right-hand side of Eq. (16) coincides with the sum of Eqs. (C.8) and (C.20). Similarly, for the perpendicular component, the right-hand side of Eq. (17) coincides with the sum of Eqs. (C.9) and (C.21). And finally, the sum of Eqs. (C.10) and (C.22) is zero, so $\mbox{$v_{\Vert}$ }=0$ is recovered.

  
3 Numerical method

This section is devoted to describing the numerical procedure used in solving Eqs. (16) and (17).

3.1 Generalised $\psi$, $\chi$ coordinates

The derivatives of the velocity components in Eqs. (16) and (17) have been written as combinations of the operators $\nabla A\cdot\nabla$ and ${\vec B}_{\rm p}\cdot\nabla$. For the numerical integration of the wave equations it will be necessary to write these operators explicitly in terms of partial derivatives in the directions defined by two spatial coordinates, for which the symbols $\psi$ and $\chi$ are used, which are functions of x and z and map the part of the xz-plane in which the numerical solution is to be computed. Examples of how to choose these two coordinates as functions of x and z are given in Sect. 4.

Now, the operators

$$% \nabla A\cdot\nabla=B_z\frac{\partial}{\partial x}-B_x\frac{\partial}{\partial z}$$ (19)

and

$$% {\vec B}_{\rm p}\cdot\nabla=B_x\frac{\partial}{\partial x}+B_z\frac{\partial}{\partial z}$$ (20)

can be transformed, with the help of the chain rule, into operators with derivatives of $\psi$ and $\chi$ instead of x and z (see Appendix D).

The possibility of using coordinates other than x and z can be of interest in a number of cases. One can, for example, choose flux coordinates (linked to the shape of magnetic surfaces) to investigate the Alfvén continuum, while selecting Cartesian coordinates may be used to study the discrete fast and Alfvén modes of the system. On the other hand, cylindrical coordinates can also be imposed for an equilibrium with longitudinal invariance such as the one assumed here. It will suffice to interchange the names for the y- and z-axes and to express r and $\theta$ as functions of the poloidal coordinates x and z.

3.2 Finite difference discretisation of wave equations

One can now insert the equations for the operators $\nabla A\cdot\nabla$, ${\vec B}_{\rm p}\cdot\nabla$, etc. from Appendix D into Eqs. (16) and (17), which results in the following two partial differential equations with non-constant and complex coefficients,

 

$$a_{11}\frac{\mbox{$\mbox{$\partial^2$ }\mbox{$v_n$ }$ }}{\mbox{$\... ...l\psi$ }}+ e_{11}\frac{\mbox{$\partial\mbox{$v_n$ }$ }}{\mbox{$\partial\chi$ }}$$

$$+f_{11}\mbox{$v_n$ }+a_{12}\frac{\mbox{$\mbox{$\partial^2$ }\mbox... ...^2$ }}+d_{12}\frac{\mbox{$\partial\mbox{$v_\perp$ }$ }}{\mbox{$\partial\psi$ }}$$

$$+e_{12}\frac{\mbox{$\partial\mbox{$v_\perp$ }$ }}{\mbox{$\partial\chi$ }}+f_{12}\mbox{$v_\perp$ }=\mbox{$\omega^2$ }\mbox{$v_n$ }$$ (21)

and

 

$$a_{21}\frac{\mbox{$\mbox{$\partial^2$ }\mbox{$v_n$ }$ }}{\mbox{$\... ...l\psi$ }}+ e_{21}\frac{\mbox{$\partial\mbox{$v_n$ }$ }}{\mbox{$\partial\chi$ }}$$

$$+f_{21}\mbox{$v_n$ }+a_{22}\frac{\mbox{$\mbox{$\partial^2$ }\mbox... ...^2$ }}+d_{22}\frac{\mbox{$\partial\mbox{$v_\perp$ }$ }}{\mbox{$\partial\psi$ }}$$

$$+e_{22}\frac{\mbox{$\partial\mbox{$v_\perp$ }$ }}{\mbox{$\partial\chi$ }}+f_{22}\mbox{$v_\perp$ }=\mbox{$\omega^2$ }\mbox{$v_\perp$ }.$$ (22)

Expressions for the coefficients a₁₁, b₁₁, ... are given in Appendix D.

We now construct a finite difference replacement of Eqs. (21) and (22). The region to be considered is a rectangle with sides parallel to the $\psi$- and $\chi$-axes, covered by a mesh of $N_\psi\times N_\chi$ points in which approximations to the functions $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ are to be computed. At the point $(\mbox{$\psi_i$ },\mbox{$\chi_j$ })$ the first and second derivatives of these functions are approximated in the usual manner (Mitchell & Griffiths 1980; Ames 1992),

  

$$% \left.\frac{\partial\mbox{$v_n$ }}{\mbox{$\partial\psi$ }}\right\vert _{i,j} \frac{\mbox{$v_n^{({i+1},{j})}$ } -\mbox{$v_n^{({i-1},{j})}$ }}{2\mbox{$h_\psi$ }}+O(\mbox{$h_\psi^2$ }),$$ = (23)

$$% \left.\frac{\partial\mbox{$v_n$ }}{\mbox{$\partial\chi$ }}\right\vert _{i,j} \frac{\mbox{$v_n^{({i},{j+1})}$ } -\mbox{$v_n^{({i},{j-1})}$ }}{2\mbox{$h_\chi$ }}+O(\mbox{$h_\chi^2$ }),$$ = (24)

$$\left.\frac{\mbox{$\partial^2$ }\mbox{$v_n$ }}{\mbox{$\partial\psi^2$ }}\right\vert _{i,j} \frac{\mbox{$v_n^{({i+1},{j})}$ } -2\mbox{$v_n^{({i},{j})}$ }+\mbox{$v_n^{({i-1},{j})}$ }}{\mbox{$h_\psi^2$ }}+O(\mbox{$h_\psi^2$ }),$$ = (25)

$$\left.\frac{\mbox{$\partial^2$ }\mbox{$v_n$ }}{\mbox{$\partial\chi^2$ }}\right\vert _{i,j} \frac{\mbox{$v_n^{({i},{j+1})}$ } -2\mbox{$v_n^{({i},{j})}$ }+\mbox{$v_n^{({i},{j-1})}$ }}{\mbox{$h_\chi^2$ }}+O(\mbox{$h_\chi^2$ }),$$ = (26)

and similarly for $\mbox{$v_\perp$ }$. In these expressions, $\mbox{$v_n^{({i},{j})}$ }$ is the value of $\mbox{$v_n$ }$ at the point $(\mbox{$\psi_i$ },\mbox{$\chi_j$ })$ and $\mbox{$h_\psi$ }$, $\mbox{$h_\chi$ }$ are the uniform grid spacings in the $\psi$- and $\chi$-directions. For the mixed derivative one has (see Mitchell & Griffiths 1980)

 

$$% \left.\frac{\mbox{$\partial^2$ }\mbox{$v_n$ }}{\mbox{$\partial\psi\partial\chi$ }}\right\vert _{i,j} \frac{1}{\mbox{$h_\psi$ }\mbox{$h_\chi$ }} \left\{\left(\alpha_2-... ...ht)\left[\mbox{$v_n^{({i+1},{j})}$ }+\mbox{$v_n^{({i},{j+1})}$ } \right.\right.$$ =

$$+\left.\mbox{$v_n^{({i-1},{j})}$ }+\mbox{$v_n^{({i},{j-1})}$ }-2\mbox{$v_n^{({i},{j})}$ }\right]$$

$$+\alpha_1\left[\mbox{$v_n^{({i+1},{j+1})}$ }+\mbox{$v_n^{({i-1},{j-1})}$ }\right]$$

$$-\left.\alpha_2\left[\mbox{$v_n^{({i-1},{j+1})}$ }+\mbox{$v_n^{({i+1},{j-1})}$ }\right]\right\}$$

$$+ O(\max\{\mbox{$h_\psi^2$ },\mbox{$h_\chi^2$ }\}).$$ (27)

The parameters $\alpha_1$ and $\alpha_2$ in this expression can be chosen so as to minimise the error of the difference approximation. This results in the formula

$$% \alpha_1+\alpha_2=\frac{1}{2}\cdot$$ (28)

Equation (27) as it stands makes use of the function whose derivative is being approximated at nine points. Setting $\alpha_1=1/2$, $\alpha_2=0$ or $\alpha_1=0$, $\alpha_2=1/2$ has the effect of reducing the number of points to just seven, while taking $\alpha_1=\alpha_2=1/4$ results in only four points being used. This freedom in the selection of the quantities $\alpha_1$, $\alpha_2$ can be used to improve the stability and convergence of the numerical solution. It is worth noticing that they can both be changed from one point to another and/or have different values for the approximation of the mixed derivatives of $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$.

Using expressions (23)-(27) to discretise Eqs. (21) and (22) in the $\left(N_\psi-2\right)\times\left(N_\chi-2\right)$ interior points yields $2\left(N_\psi-2\right)\left(N_\chi-2\right)$ equations for the unknowns at these points and the $2\left(N_\psi+N_\chi-2\right)$ unknowns at the points on the boundary of the system. The same discretised equations can be assumed to apply at the boundary points, but then some boundary conditions on $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ are needed to eliminate the unknowns that are introduced by the use of points outside the region of integration. This issue will be discussed next.

3.3 Boundary conditions

At present, only homogeneous boundary conditions have been implemented in our numerical code. Thus, for each of the two velocity components it is possible to choose between Dirichlet ( $\mbox{$v_n$ }=0$ or $\mbox{$v_\perp$ }=0$) or Neumann ( $\frac{\partial\mbox{$v_n$ }}{\partial n}=0$ or $\frac{\partial\mbox{$v_\perp$ }}{\partial n}=0$) conditions. Here $\frac{\partial}{\partial n}$ represents the derivative normal to the boundary, which means that $\frac{\partial}{\partial n}\equiv\frac{\partial}{\partial\psi}$ on a $\chi=constant$ boundary and $\frac{\partial}{\partial n}\equiv\frac{\partial}{\partial\chi}$ on a $\psi=constant$ boundary. The first of these boundary conditions is useful when modelling a boundary that behaves as a rigid, perfectly conducting wall, while the second one can be applicable to exploit spatial symmetry in the solutions, as will be shown in the next section.

There are four possible combinations of boundary conditions, namely

   

$$% \mbox{$v_n$ } \mbox{$v_\perp$ }=0,$$ = (29)

$$% \frac{\partial\mbox{$v_n$ }}{\partial n} \mbox{$v_\perp$ }=0,$$ = (30)

$$% \mbox{$v_n$ } \frac{\partial\mbox{$v_\perp$ }}{\partial n}=0$$ = (31)

and

$$% \frac{\partial\mbox{$v_n$ }}{\partial n}=\frac{\partial\mbox{$v_\perp$ }}{\partial n}=0.$$ (32)

From the numerical viewpoint, the second and third combinations are equivalent and so only the implementation of Eqs. (29), (30) and (32) are explained next. The first case is the simplest to treat and one only needs to consider the finite difference equivalent to Eqs. (21) and (22) at interior points while substituting $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ by zero for the values of these functions on the boundary. This does not introduce any extra unknowns, i.e. there is no need to use any points exterior to the system. Moreover, the finite difference version of Eqs. (21) and (22) must not be considered on the boundary, so the number of algebraic equations is reduced in $N_\psi$ (for a $\chi=constant$ boundary) or $N_\chi$ (for a $\psi=constant$ boundary).

The cases given by Eqs. (30), (31) and (32) are a bit more difficult to handle since the use of the difference equations at boundary points forces one to introduce fictitious (or "ghost'') points exterior to the region of integration. Equation (32) poses little problem since it implies that $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ on the ghost points can be substituted by their counterparts directly across the boundary. Thus, on the $\psi=\psi_1$ boundary Eq. (32) is equivalent to $\frac{\partial\mbox{$v_n$ }}{\partial\chi}=\frac{\partial\mbox{$v_\perp$ }}{\partial\chi}=0$, which by virtue of Eq. (24) yields $\mbox{$v_n^{({0},{j})}$ }=\mbox{$v_n^{({2},{j})}$ }$ and $\mbox{$v_\perp^{({0},{j})}$ }=\mbox{$v_\perp^{({2},{j})}$ }$, with the index i=0 indicating fictitious points and $j=1,2,\ldots,N_\chi$.

Finally, dealing with boundary conditions of the form (30) is still more complex. Now $\mbox{$v_\perp$ }$ is known on the boundary and Eq. (22) needs not be used. However, Eq. (21) discretised on the boundary contains the values of the two velocity components on the fictitious points and while $\frac{\partial\mbox{$v_n$ }}{\partial n}$ helps eliminate the values of $\mbox{$v_n$ }$ outside the computational domain in favour of values inside the domain, it is not possible to do the same for $\mbox{$v_\perp$ }$. Our approach has been to assume that $\mbox{$v_\perp$ }$ is antisymmetric about the boundary, so the values of this function on the ghost mesh points are equal to minus the interior values just across the boundary. Hence, on the $\psi=\psi_1$ boundary Eq. (30) is interpreted as $\mbox{$v_n^{({0},{j})}$ }=\mbox{$v_n^{({2},{j})}$ }$ and $\mbox{$v_\perp^{({0},{j})}$ }=-\mbox{$v_\perp^{({2},{j})}$ }$, with $j=1,2,\ldots,N_\chi$.

3.4 Solving the complex eigenvalue problem

After discarding the boundary points on which Eqs. (21) and (22) do not need to be discretised, we are left with a set of algebraic expressions for the unknowns $\mbox{$v_n^{({i},{j})}$ }$ and $\mbox{$v_\perp^{(i,j)}$ }$ with an additional parameter that must be determined: the frequency $\omega$. This set of equations constitute an eigenvalue problem of the form

$$% {\vec A \vec U}=\lambda {\vec U},$$ (33)

where $\lambda$ is defined as $\omega^2$, $\vec U$ is the vector of unknown velocity components (i.e., the $\mbox{$v_n^{({i},{j})}$ }$ and $\mbox{$v_\perp^{(i,j)}$ }$ ordered in some manner) and $\vec A$ is the matrix of coefficients constructed from the finite difference counterparts of Eqs. (21)-(22).

To set up $\vec U$ one must keep in mind that each point on the two-dimensional grid can be labelled with two indices, i and j, related to its coordinates $\psi_i$ and $\chi_j$. It is then necessary to introduce a way of labelling all points with a single index, k. To do so, the mesh points are ordered consecutively, starting from $(\psi_1,\chi_1)$, next $(\psi_2,\chi_1)$, etc., until $(\psi_{N_\psi},\chi_1)$. The next point has coordinates $(\psi_1,\chi_2)$, then $(\psi_2,\chi_2)$, and so on until the point $(\psi_{N_\psi},\chi_{N_\chi}$) is reached; see Fig.  2 for an example with $N_\psi =5$ and $N_\chi =6$.

	Figure 2: Mesh points for a finite difference computation with $N_\psi =5$ and $N_\chi =6$. The label k used in defining the vector $\vec U$ is printed next to each point on the grid
Open with DEXTER

Finally, the unknowns are put together in the vector $\vec U$ constructed as follows

$$% {\vec U}=\left(\mbox{$v_n^{1}$ },\mbox{$v_\perp^{1}$ },\mb... ...k}$ },\ldots,\mbox{$v_n^{N}$ },\mbox{$v_\perp^{N}$ }\right)^T,$$ (34)

where $\mbox{$v_n^{k}$ }$ and $\mbox{$v_\perp^{k}$ }$, with $k=1,2,\ldots,N$, stand for the approximate values of the eigenfunctions $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ at the kth grid point. Here $N=N_\psi\times N_\chi$ is the total number of grid points. From the previous discussion of boundary conditions it is obvious that the arrangement given in Eq. (34) is valid for the case in which Neumann conditions are imposed on the two velocity components on all four boundaries. Nevertheless, whenever Dirichlet conditions are imposed on an unknown on a given boundary, the value of that unknown is not included in the vector $\vec U$.

To solve the eigenvalue problem (33) it is important to keep in mind that $\vec A$ is a large, sparse, non-symmetric matrix and that common methods used for calculating eigenvalues and/or eigenvectors of a matrix are computationally expensive. These methods become impractical when dealing with very large matrices, as in the present case, in which the size of $\vec A$ is of the order of $(2N_\psi N_\chi)\times(2N_\psi N_\chi)$, i.e. $5000\times 5000$ for the somewhat moderate values $N_{\psi }=N_{\chi }=50$.

The approach chosen here has been to select numerical methods that exploit the sparseness of $\vec A$. First, one can use a coarse grid to compute selected eigenvectors, e.g. those whose moduli lie between two values, $\lambda_{\rm min}$ and $\lambda_{\rm max}$. This task has been accomplished with the NAG subroutine F02BDF, which first reduces the complex matrix to Hessenberg form and then makes use of the LR algorithm and inverse iteration to obtain the desired eigensolutions. It must be emphasised that to keep the order of $\vec A$ reasonably small the mesh must be rather coarse (we have used values up to $N_\psi=90$ and $N_\chi=31$), since all elements (both zero and non-zero) of the matrix are stored and this results in large memory requirements. Moreover, the computational time scales as n³, where $n\simeq 2N_\psi N_\chi$ as we mentioned before.

At this stage some information has been determined regarding the spectrum of the system, so it is now possible to refine the solutions by using a finer grid together with the frequency values coming from the previous step. For the purpose of obtaining single eigenvalues and their corresponding eigensolution, the inverse iteration has been used. Details of this technique can be found in Kerner (1989) and Press et al. (1992); we have followed the second of these references but in some aspects our implementation is different from theirs and in addition there seems to be an error in their Eq. (11.7.10), we therefore describe our procedure.

Recall that the aim is to obtain an approximation, as accurate as possible, of an eigensolution of the eigenvalue problem (33). We start from an initial approximation, $\lambda_0$, of the desired eigenvalue and construct the matrix ${\vec A}-\lambda_0{\vec I}$, with ${\vec I}$ the identity matrix. The idea now is to express $\lambda$ satisfying Eq. (33) as

$$% \lambda=\lambda_0+\Delta\lambda,$$ (35)

where $\Delta\lambda$ is determined using an iterative process. At the ith step we solve the algebraic system of equations

$$% \left({\vec A}-\lambda_0{\vec I}\right)\cdot{\vec y}={\vec b}_i,$$ (36)

where ${\vec b}_i$ is the solution $\vec y$ to this equation in the (i-1)th step, normalised so that ${\vec b}_i\cdot{\vec b}_i=1$. In the first iteration (i=1) one can simply take ${\vec b}_0$ as a random normalised vector, so there is no need for a priori knowledge of the eigenvector being searched. After Eq. (36) has been solved, one computes the new vector ${\vec b}_{i+1}$ by normalising $\vec y$, i.e.

$$% {\vec b}_{i+1}=\frac{\vec y}{\vert\vec y\vert}\cdot$$ (37)

This vector will then become the right-hand side of Eq.  (36) in the next iteration. Now an estimate of $\Delta\lambda$ is calculated from

$$% \Delta\lambda_{i+1}=\frac{{\vec b}_i\cdot{\vec b}_{i+1}}{\vert\vec y\vert}\cdot$$ (38)

At the end of the iterative process the eigensolution is simply

$$% \lambda\simeq \lambda_0+\Delta\lambda_{i+1}$$ (39)

and

$$% {\vec U}\simeq {\vec b}_{i+1}.$$ (40)

The iterations are stopped when either of the following criteria are fulfilled,

$$% \vert\Delta\lambda_{i+1}-\Delta\lambda_i\vert < \epsilon_1,$$ (41)

$$% \vert{\vec b}_{i+1}-{\vec b}_i\vert < \epsilon_2,$$ (42)

where $\epsilon_1$ and $\epsilon_2$ are two small constants (e.g., 10^-12 and 10^-10, respectively). In some cases, after a few iterations there is no progress and the quantities $\vert\Delta\lambda_{i+1}-\Delta\lambda_i\vert$ and $\vert{\vec b}_{i+1}-{\vec b}_i\vert$ do not decrease rapidly enough. It is then better to restart the iterations with a new guess for $\lambda$ equal to $\lambda_0+\Delta\lambda_{i+1}$.

The advantage of the inverse iteration method is that it is much faster than other methods which calculate all or selected eigensolutions of $\vec A$. Nevertheless, one still must solve a large algebraic problem (Eq. (36)) a number of times. It is here that the sparseness of $\vec A$ can be exploited, since there are some algorithms that efficiently solve this kind of problem. One of these algorithms is provided as the package ME48, which is part of the Harwell Sparse Matrix Library. There is one subroutine in this package that first selects permutations ${\vec P}$ and ${\vec Q}$ so that the matrix ${\bf PAQ}$ is in block triangular form and can thus be easily factorised. Next, a second subroutine performs the factorisation of each diagonal block and, finally, a third subroutine solves the system by block solution through the block triangular form. By far the most time consuming of these tasks is the matrix factorisation, which must be carried out only once at the beginning of the iterative process (i.e. for i=1). The details of the factorisation are stored and in subsequent steps one can use them to solve Eq. (36) at a small computational cost.

3.5 Computation of perturbed physical variables

After solving the eigenvalue problem, Eq. (33), we have an approximation to the spatial distribution of the functions $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ from which the perturbed quantities (density, magnetic field components and magnetic pressure) and perturbed forces (components of magnetic pressure gradient and magnetic tension) can be obtained. We use Eqs. (C.1)-(C.10) and (C.20)-(C.22) and again write the directional derivatives as in Appendix D. Next, the derivatives of $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ in the resulting expressions are approximated by the finite difference formulas (23)-(27). We have found it very useful to compute these physical quantities for they can shed light on how the oscillatory motions are driven and, therefore, on the nature of the modes.

3.6 Practical use of the numerical code

Despite the complexity of the coefficients in Eqs. (21) and (22), which contain many derivatives of equilibrium quantities and of the generalised coordinates $\psi$ and $\chi$, the numerical code is rather simple to use. The end user only needs to define a few of subroutines, providing details on the equilibrium (in particular, B_x, B_y, B_z and $v_{\rm A}^2$ as functions of position) and also the first- and second-order derivatives of the Cartesian magnetic field components with respect to x and z. For the computation of the density perturbation, the first-order derivatives of $\rho$ with respect to x and z must also be given. Moreover, if generalised coordinates are to be used instead of the x- and z-coordinates, then the definition of $\psi$ and $\chi$ in terms of x and z must be provided together with their first- and second-order derivatives.

  
4 Numerical results

In this section we analyse numerically the normal modes of oscillation of two magnetic configurations: a uniform, straight magnetic field and a potential arcade. These two cases have analytical or simple numerical solutions so that the goodness of the numerical code can be checked.

4.1 Uniform straight magnetic field

We first consider a cold plasma in a horizontal, straight and uniform magnetic field ${\vec B}=B\mbox{$\hat{{\vec e}}_x$ }$. The density is given by $\rho=\rho\left(z\right)$ so that the Alfvén speed is $\mbox{$v^2_{\rm A}$ }\left(z\right)=B^2/\mu\rho$. Also, the plasma is assumed to be contained in a box limited by rigid plates at $x=\pm L$ and z=0, H, which corresponds to Dirichlet boundary conditions for $\mbox{$v_n$ }$ and $\mbox{$v_\perp$ }$ on the edges of the computational domain.

The wave Eqs. (16) and (17) reduce to

 

$$-\frac{\omega^2}{\mbox{$v^2_{\rm A}$ }}\mbox{$v_n$ }=\partial_{x}^2\mbox{$v_n$ }+\partial_{z}^2\mbox{$v_n$ },$$ (43)

 

$$\left(\partial_{x}^2-\frac{\omega^2}{\mbox{$v^2_{\rm A}$ }}\right)\mbox{$v_\perp$ }=0.$$ (44)

Equations (43) and (44), having analytical solutions, describe uncoupled fast and Alfvén waves. Fast modes are characterised by motions in the normal direction ( $\mbox{$v_n$ }\neq 0$, $\mbox{$v_\perp$ }=0$) while Alfvén modes have motions in the perpendicular direction ( $\mbox{$v_n$ }=0$, $\mbox{$v_\perp$ }\neq 0$).

4.1.1 Alfvén continuum modes

Equation (44) describes Alfvén continuum modes. The continuous part of the MHD wave spectrum is well documented in the plasma physics literature (see e.g. Appert 1974; Goedbloed 1975, 1983) and distribution theory yields a solution to Eq. (44) in terms of the delta function $\delta\left(z-z_{\rm A}\right)$, where $z=z_{\rm A}$ is the singular surface, multiplied by a function of x. After imposing the vanishing of $\mbox{$v_\perp$ }$ at $x=\pm L$ we obtain

$$% \mbox{$v_\perp$ }=C \sin\left(k_{x}x\right)\delta\left(z-z_{\rm A}\right),$$ (45)

with C an arbitrary constant, $k_{x}=n\pi/2L$ and n an integer. In addition, the frequency is given by

$$% \omega^2\equiv\omega_{\rm A}^2=k_{x}^2\mbox{$v^2_{\rm A}$ }\left(z_{\rm A}\right).$$ (46)

The collection of improper eigenvalues ${\omega_{\rm A}^2}$ for the interval $0\le z \le H/L$ constitutes the continuous Alfvén spectrum.

Figure 3: a) Surface plot of the perpendicular velocity component for the Alfvén fundamental continuum mode in a uniform, straight magnetic field configuration ($\delta =6$, H/L=1). This solution, with $\omega L/\mbox{$v_{\rm A}$ }_0=4.33$, has been computed numerically with $N_{\psi }=N_{\chi }=50$ grid points and boundary conditions (29). b) Spatial profile of $\mbox{$v_\perp$ }$ on the singular surface $z_{\rm A}=0.5L$. The solid line is the x-dependent part of the analytical solution (Eq. (45)) and the circles are the finite difference solution $\mbox{$v_\perp$ }(x,z_{\rm A})$ to Eqs. (16) and (17)

Open with DEXTER

We have next assumed an Alfvén speed profile of the form

$$% \mbox{$v^2_{\rm A}$ }\left(z\right)=\mbox{$v^2_{\rm A}$ }_{0} {\rm e}^{-\left(2-\delta\right)z/L},$$ (47)

$\delta$ being a parameter that determines the rate of change of the Alfvén speed with height from the value $v_{\rm A0}$ at z=0. Now, the continuum Alfvén mode frequency is given by

$$% \omega_{\rm A}=k_{x}\mbox{$v_{\rm A}$ }_0{\rm e}^{-\left(2-\delta\right)z/2L}.$$ (48)

Next, Eqs. (16) and (17) have been solved numerically for $\delta =6$, H/L=1 with the aim of recovering the analytical eigenfunction (Eq. (45)) and eigenfrequency (Eq. (48)). Cartesian coordinates have been chosen by setting $\psi$ and $\chi$ as follows

  

$$% \psi z/L, {\mbox{\hspace{2cm}}} 0\leq\psi\leq H/L,$$ = (49)

$$% \chi x/L, {\mbox{\hspace{2cm}}} -1\leq\chi\leq 1.$$ = (50)

A mesh with 50 points in the two spatial directions ( $N_{\psi }=N_{\chi }=50$) has been used and boundary conditions given by Eq. (29) have been imposed.

Continuum solutions, like the one shown in Fig. 3a, are characterised by $\mbox{$v_n$ }\simeq 0$ (to an extremely high accuracy) in the whole computational domain and $\mbox{$v_\perp$ }\simeq 0$ everywhere except on a single magnetic surface. Therefore, our numerical code recovers the delta function in the analytical solution (45) by renormalising the numerical eigenfunctions in such a way that finite values of $\mbox{$v_\perp$ }$ outside the singular surface and of $\mbox{$v_n$ }$ everywhere become zero and infinite values of $\mbox{$v_\perp$ }$ on the singular surface become finite. To confirm the coincidence between the numerical and analytical solutions we have plotted them together in Fig. 3b, that shows a perfect agreement.

Finally, the Alfvén continuous spectrum has been computed numerically and compared with the analytical formula (48). The results, displayed in Fig. 4, again show an excellent agreement.

	Figure 4: Plot of Alfvén continuum frequency, $\omega$, versus the position of the singular magnetic surface, $z_{\rm A}$, for the fundamental mode in a uniform, straight magnetic field configuration with $\delta =6$, H/L=1. The solid line represents the analytical solution, given by Eq. (48), and the circles are the numerical solution computed with $N_{\psi }=N_{\chi }=50$
Open with DEXTER

  
4.1.2 Fast modes

Upon assuming that the normal velocity in Eq. (43) can be written as a separable function of x and z, we have

$$% \mbox{$v_n$ }=X\left(x\right)Z\left(z\right).$$ (51)

Now, Eq. (43) reduces to two ordinary differential equations for X and Z. The first one can be readily solved and, after imposing the boundary condition $X\left(-L\right)=X\left(L\right)=0$, results in

$$% X\left(x\right)=C \sin\left(k_{x} x\right),$$ (52)

where C is an arbitrary constant and $k_{x}=n_x\pi/2L$ with $n_x=1,2,\ldots$ Here n_x is the number of extrema (i.e. maxima and minima) of $\mbox{$v_n$ }$ in the x-direction.

For the vertical dependence of the normal velocity we obtain a second-order ordinary differential equation,

$$% \frac{{\rm d}^2 Z}{{\rm d}z^2}+\left[\frac{\omega^2}{\mbox{$v^2_{\rm A}$ }\left(z\right)}-k_{x}^2\right] Z=0.$$ (53)

To solve this equation we consider the new independent variable (valid for $\delta\neq 2$)

  

$$u=\alpha \ {\rm e}^{-\left(\delta-2\right)z/2L},$$ (54)

$$\mbox{with}$$

$$\alpha=\frac{2}{\vert\delta-2\vert}\frac{\omega L}{\mbox{$v_{\rm A}$ }_0}\cdot$$ (55)

$$\mbox{Then Eq.~(\ref{vertical1}) becomes }$$

$$u\frac{\rm d}{{\rm d}u}\left(u\frac{{\rm d}Z}{{\rm d}u}\right)+\left(u^2-m^2\right)Z=0,$$ (56)

where

$$m=\frac{2k_{x}L}{\vert\delta-2\vert}\cdot$$ (57)

The solution to Eq. (56) can be written as

 

$$Z\left(u\right)=C_{1}\ J_{m}\left(u\right)+C_{2}\ Y_{m}\left(u\right),$$ (58)

$$\mbox{or}$$

$$Z\left(z\right)=C_{1}\ J_{m}\left(\alpha \ {\rm e}^{-\left(\delta-2\right)z/2L}\right)$$

$$+ C_{2}\ Y_{m}\left(\alpha \ {\rm e}^{-\left(\delta-2\right)z/2L}\right).$$ (59)

The dispersion relation, which follows after imposing Z(0)=Z(H)=0, is thus

$$J_{m}\left(\alpha\right)Y_{m}\left(u_{\rm H}\right) -J_{m}\left(u_{\rm H}\right)Y_{m}\left(\alpha\right)=0,$$ (60)

where

$$u_{\rm H}=\alpha \ {\rm e}^{-\left(\delta-2\right)H/2L}$$ (61)

is simply the value of u at the boundary z=H.

Figure 5: a) Surface plot of the normal velocity component for the fundamental (i.e., nx=1, nz=1) fast mode in a straight magnetic field with $\delta =4$, H/L=4. This solution, with $\omega L/\mbox{$v_{\rm A}$ }_0=4.5817$, has been computed using a grid of $61\times 61$ points in the xz-plane and rigid plate boundary conditions (Eq. (29)). b) Vertical dependence of $\mbox{$v_n$ }$ for the same fast mode. Solid line: analytical solution Eq. (59). Circles: numerical solution to Eqs. (16) and (17)

Open with DEXTER

An analytical solution can be derived by first computing $\omega$ from the transcendental expression (60) once $\delta$, H/L and n_x are imposed. The frequency of the mode can next be substituted into Eq. (55) to compute $\alpha$. This parameter is finally inserted into Eq. (59) and the spatial structure of the eigenfunction $\mbox{$v_n$ }$ is completely determined. There is a whole array of eigensolutions, each with its own frequency, characterised by the number n_x and n_z of extrema of $\mbox{$v_n$ }$ in the x- and z-directions.

A numerical solution has been computed for a $\delta =4$, H/L=4 configuration (see Fig. 5a) using $\psi$ and $\chi$ given by Eqs. (49) and (50). This mode has $\mbox{$v_n$ }\neq 0$ and $\mbox{$v_\perp$ }=0$ (within machine precision), such as corresponds to a fast mode, and a close inspection reveals $\mbox{$v_n$ }$ is a separable function of x and z. Since $\mbox{$v_n$ }$ displays a single extrema in the z-direction we label this as a n_z=1 mode. A comparison between the z-dependence of this numerical solution and the analytical one is shown in Fig. 5b, which displays the fine performance of the numerical code.

The last check we perform with the straight magnetic field configuration consists in comparing the numerical frequency for various fast modes with the solutions calculated from Eq. (60); see Fig. 6. As expected, the two sets of frequencies agree quite well, the larger difference appearing for modes with higher spatial structure (that is, with larger values of n_x and/or n_z) for which a finer grid is required to achieve an accuracy similar to that of smoother modes.

Figure 6: Comparison of fast mode frequencies in a uniform, straight magnetic field configuration with $\delta =4$, H/L=4. Diamonds represent $\omega$ computed by solving Eq. (60). Filled circles represent $\omega$ computed numerically by solving Eqs. (16) and (17) with a mesh of $N_{\psi }=61$, $N_{\chi }=61$ points. The number of extrema of the eigenfunction $\mbox{$v_n$ }$ in the x-direction is nx. For a fixed nx and with increasing $\omega$, the number of extrema of $\mbox{$v_n$ }$ in the vertical direction is $n_z=1,2, \ldots$

Open with DEXTER

4.2 Coronal potential arcade

Another particular solution of the force-free Eq. (1) is a potential structure such as a coronal arcade of magnetic field lines (see Fig. 7) with $\mbox{$B_{y}$ }=0$ and a flux function of the form

$$% A\left(x,z\right)=\frac{B_{0}}{k}\cos kx \, \ {\rm e}^{-kz},$$ (62)

where $\Lambda_{B}\equiv k^{-1}=2L/\pi$ is the magnetic scale height. If a density profile of the form $\rho=\rho_{0}\ {\rm e}^{-z/\Lambda}$ is assumed, then the Alfvén speed of this potential coronal arcade is given by

$$% \mbox{$v_{\rm A}$ }\left(z\right)=\mbox{$v_{\rm A}$ }_{0}\ {\rm e}^{-\left(2-\delta\right)k z/2L},$$ (63)

with $\mbox{$v_{\rm A}$ }_{0}$ the Alfvén speed at z=0 and $\delta$ the ratio of the magnetic scale height to the density scale height ( $\delta=\Lambda_{\rm B}/\Lambda$).

	Figure 7: Magnetic field lines in a coronal potential arcade with a flux function given by Eq. (62). The field related vectors $\mbox{$\hat{{\vec e}}_n$ }$ and $\mbox{$\hat{{\vec e}}_{\Vert}$ }$, both contained in the xz-plane, have been drawn at a particular point
Open with DEXTER

The MHD modes of oscillation of this potential arcade were investigated by Terradas et al. (1999) and by Oliver et al. (1993), who wrote the linear ideal MHD wave equations for this zero-$\beta$ configuration as follows

$$% \mbox{$v^2_{\rm A}$ }\nabla^2\left(\mbox{$v_n$ }\, {\rm e}^{-kz}\right)+ \omega^2\mbox{$v_n$ }\, {\rm e}^{-kz}=0$$ (64)

and

$$% -\rho\omega^2\mbox{$v_\perp$ }=\frac{1}{\mu}\mbox{$\left({\vec B}\cdot\nabla\right)$ }^2\mbox{$v_\perp$ }.$$ (65)

4.2.1 Alfvén continuum modes

Equation (65) describes the Alfvén mode, characterised by plasma displacements in the perpendicular direction, i.e. along the y-axis. The ${\vec B}\cdot\nabla$ derivative makes this mode highly anisotropic and motions are confined to a particular magnetic surface (described by A=constant). Hence, solutions to Eq. (65) are of the form

$$% \mbox{$v_\perp$ }=\tilde v_\perp(x)\delta(\psi-\psi_{\rm A}),$$ (66)

where $\psi=\cos kx\,{\rm e}^{-kz}$ is a new spatial coordinate proportional to the flux function and therefore constant on a given magnetic surface. The delta function in this expression accounts for the non-square integrable behaviour of Alfvén continuum modes, whereas $\tilde v_\perp$ gives the variation of the perpendicular velocity component on the singular surface $\psi=\psi_{\rm A}$.

Figure 8: a) Surface plot of the perpendicular velocity component for the Alfvén fundamental continuum mode in a potential arcade with $\delta =6$. This numerical solution has been calculated using $N_{\chi }=N_{\psi }=100$. Plasma motions are confined to the singular surface $\psi _{\rm A}\simeq 0.7070$, which corresponds to a footpoint position x0=0.5L. The mode frequency is $\omega L/\mbox{$v_{\rm A}$ }_0=5.3384$. b) Spatial profile of $\mbox{$v_\perp$ }$ on the singular surface between x=0 and x=x0. The coordinate $\chi$ has been eliminated in favour of x by means of Eq. (72). The solid line is the solution to Eq. (68) and the circles are the finite difference solution to Eqs. (16) and (17). Symmetry of $\mbox{$v_\perp$ }$ about x=0 has been used when computing the numerical solution by selecting the appropriate boundary conditions

Open with DEXTER

	Figure 9: Plot of the frequency, $\omega$, of the Alfvén fundamental mode versus the footpoint position, x0, in a potential arcade with the same parameter values as in Fig. 8. The solid line represents the solution to Eq. (68) and the circles the two-dimensional numerical solution obtained with $N_\psi =N_\chi =100$
Open with DEXTER

In order to obtain $\tilde v_\perp$, Oliver et al. (1993) considered a particular field line whose footpoints have coordinates $x=\pm x_0$, z=0. The directional derivative ${\vec B}\cdot\nabla$ along this field line can be written as

$$% {\vec B}\cdot\nabla=\mbox{$B_{x}$ }\frac{\partial}{\partial... ...{$B_{x}$ }\frac{\rm d}{{\rm d}x}=kA(x_0)\frac{\rm d}{{\rm d}x}$$ (67)

with $A(x_0)=B_{0}/k\cos kx_0$. After some manipulations Eq. (65) reduces to

$$% \frac{{\rm d}^2\tilde v_\perp}{{\rm d}x^2} +\frac{\omega^2... ...{\cos kx}\right]^{\delta} \cos^{-2} kx_{0}\,\,\tilde v_\perp=0$$ (68)

(cf. Oliver et al. 1993 for further details). Assuming that the photosphere (z=0) behaves as a rigid and perfectly conducting medium to coronal perturbations, the boundary condition for $\tilde v_\perp$ is

$$% \tilde v_\perp\left(\pm x_{0}\right)=0.$$ (69)

For the case $\delta=0$ Eq. (68) becomes an ordinary differential equation with constant coefficients and analytical solutions can be easily obtained. For $\delta\neq 0$, however, Oliver et al. (1993) failed to find an analytical solution and had to integrate Eq. (68) numerically.

For the numerical computations we have considered the case $\delta =6$ and have chosen generalised coordinates that facilitate the recovery of continuum modes. By selecting

  

$$% \psi \cos kx\,{\rm e}^{-kz},{\mbox{\hspace{2cm}}} 0\le\psi\le 1,$$ = (70)

$$% \chi \frac{\sin kx\,{\rm e}^{-kz}}{\left(1-\psi^2\right)^{1/2}},{\mbox{\hspace{1.5cm}}} -1\le\chi\le +1,$$ = (71)

field lines become straight lines with $\psi=constant$and, for each value of $\psi$, $\chi$ varies between -1 at the left footpoint (x=-x₀, z=0) and +1 at the right footpoint (x=x₀, z=0).

Numerical solutions have been found by solving Eqs. (16) and (17) in a mesh with $N_\psi=100$, $N_\chi=100$ points. The eigenfunction plotted in Fig. 8a possesses, as expected, $\mbox{$v_n$ }=0$ on the whole system and $\mbox{$v_\perp$ }=0$ except on a singular surface $\psi=\psi_{\rm A}$. Using Eqs. (70) and (71) one can express $\chi$ in terms of x on the singular surface,

$$% x=\frac{1}{k}\arctan\frac{\chi\left(1-\psi_{\rm A}^2\right)^{1/2}}{\psi_{\rm A}}\cdot$$ (72)

Now, $\mbox{$v_\perp$ }$ can be drawn as a function of x on $\psi=\psi_{\rm A}$ and a comparison with the solution to Eq. (68) can be made (see Fig. 8b). Both the eigenfunction and the eigenvalue ( $\omega L/\mbox{$v_{\rm A}$ }_0=5.3384$) are in good agreement.

Changing the position of the field line footpoint, x₀, produces a variation of the mode frequency that generates the Alfvén continuous spectrum. Figure 9 shows the frequency of the fundamental Alfvén continuum mode versus the footpoint position, with the numerical values fitting well the results of Oliver et al. (1993).

Figure 10: a) Surface plot of $\mbox{$v_n$ }$ in a potential arcade with $\delta =4$, H/L=4 for a mode with nx=2, nz=1 and $\omega L/\mbox{$v_{\rm A}$ }_0=8.06168$. The solution has been obtained using a grid of $60\times 170$ points in the xz-plane. b) Vertical dependence of $\mbox{$v_n$ }$. Solid line: analytical solution from Oliver et al. (1993). Circles: numerical solution to Eqs. (16) and (17)

Open with DEXTER

4.2.2 Fast modes

In order to obtain analytical solutions to Eq. (64), Oliver et al. (1993) expressed the variable $\mbox{$v_n$ }\,{\rm e}^{-kz}$ as a separable function in the x- and z-coordinates,

$$% \mbox{$v_n$ }\left(x,z\right)\,{\rm e}^{-kz}=X\left(x\right)Z\left(z\right),$$ (73)

which can be inserted in Eq. (64) to give

$$% X\left(x\right)=C \sin\left(k_{x} x\right),$$ (74)

where C is an arbitrary constant and $k_{x}=n_x\pi/2L$ with $n_x=1,2,\ldots$ Again n_x is the number of maxima and minima of $\mbox{$v_n$ }$ in the horizontal direction.

For the vertical dependence a second-order differential equation similar to the one for the straight field (Eq. (53)) is obtained,

$$% \frac{{\rm d}^2Z}{{\rm d}z^2}+ \left[\frac{\omega^2}{\mbox{... ...{0}^2}\,{\rm e}^{-k\left(\delta-2\right)z} -k_{x}^2\right]Z=0.$$ (75)

Following the procedure used to solve Eq. (53), the solution to this equation can be written in terms of the Bessel functions J_m and Y_m. Then, imposing Z(0)=Z(H)=0 a dispersion relation, similar to Eq. (60) is derived. The mechanism for obtaining the frequency and spatial profile of the normal velocity component is the same of Sect. 4.1.2.

Numerical solutions have been computed for a potential arcade configuration with $\delta =4$ and H/L=4. Eqs. (16) and (17) are now solved in a Cartesian mesh with $-1\leq x/L \leq +1$ and $0\leq z/L \leq H/L$ and $60\times 170$ points in the x- and z-directions. Figure 10a shows a surface plot of $\mbox{$v_n$ }$ for the n_x=2, n_z=1 fast mode (which, as all other fast modes in this configuration, have $\mbox{$v_\perp$ }=0$). As can be seen, the spatial structure of the eigenfunction appears to be a separable function of the x and z variables, as was assumed by Oliver et al. (1993) in obtaining analytical solutions to the wave Eq. (73).

Figure 11: Comparison of fast mode frequencies in a potential arcade with $\delta =4$, H/L=4. Diamonds represent $\omega$ computed as in Oliver et al. (1993) by solving Eq. (64). Circles represent $\omega$ computed numerically by solving Eqs. (16) and (17) with a mesh of Nx=60, Nz=170 points. The number of extrema of the eigenfunction $\mbox{$v_n$ }$ in the x-direction is nx. For a given nx and with increasing $\omega$, the number of extrema of $\mbox{$v_n$ }$ in the vertical direction is $n_z=1,2, \ldots$

Open with DEXTER

Figures 10b and 11 show the perfect agreement between the numerical and analytical solutions in both the spatial structure of the eigenfunctions as well as in the values of the wave frequencies. The agreement between the two sets of solutions in Fig. 11 is better for modes with low frequency and with few spatial oscillations in the xz-domain.

  
5 Summary

In this paper, the linearised MHD wave equations for a two-dimensional force-free magnetic field, having longitudinal invariance, have been derived and a numerical code to solve the resulting system of coupled partial differential equations has been introduced. Two simple magnetic configurations (a straight magnetic field and a potential arcade), having analytical or simple numerical solutions, have been used to test the goodness and accuracy of the code, and the comparison of the results points out an excellent agreement. The code also proves to be a good tool for the investigation of Alfvén continuum modes in two-dimensional equilibria, which are obtained by selecting flux coordinates instead of Cartesian coordinates. Of course, the code can be used to study linear waves in more complex force-free fields and our aim is to use it to investigate linear MHD waves in sheared coronal magnetic arcades. In these structures, the longitudinal component of the magnetic field as well as the longitudinal wavenumber will produce the coupling of fast and Alfvén modes, changing substantially the spatial structure of the perturbed physical quantities in such a way that modes can possess both fast and Alfvén mode properties. Also, in the case of a curved magnetic field, curvature will introduce extra effects worth being studied.

Despite the complexity of Eqs. (16) and (17), the numerical code is simple to use since the various quantities coming into the wave equations (like the coefficients presented in Appendices B and D) are computed from subroutines provided by the user in which the equilibrium variables and their (first- and/or second-order) derivatives are specified. In addition, the possibility of using different coordinate systems, apart from the Cartesian ones, allows the oscillatory modes of different equilibrium configurations to be explored. Finally, the output of the code consists not only of the normal and perpendicular velocity components, but also the perturbed density and the three components of the perturbed magnetic tension and perturbed magnetic pressure gradient. All this information will be helpful when investigating the physical nature of coupled fast and Alfvén modes.

Acknowledgements

I. A. acknowledges DGES for a fellowship. The authors wish to acknowledge the financial support received from DGES under project PB96-0092.

References

 

• Ames, W. F. 1992, Numerical Methods for Partial Differential Equations, 3rd ed. (New York, Academic Press) In the text  
• Appert, K., Gruber, R., & Vaclavik, J. 1974, Phys. Fluids, 17, 1471 In the text  
• Aschwanden, M. J., Fletcher, L., Schrijver, C. J., & Alexander, D. 1999, ApJ, 520, 880 In the text NASA ADS  
• DeForest, C. E., & Gurman, J. B. 1998, ApJ, 501, L217 In the text NASA ADS  
• De Moortel, I., Ireland, J., & Walsh, R. W. 2000, A&A, 355, 23 In the text NASA ADS  
• Goedbloed, J. P. 1975, Phys. Fluids, 18, 1258 In the text NASA ADS  
• Goedbloed, J. P. 1983, Lecture Notes on Ideal Magnetohydro-dynamics, Rijnhuizen Report, 83 In the text  
• Kerner, W. 1989, J. Comp. Phys., 85, 1 In the text  
• Mitchell, A. R., & Griffiths, D. F. 1980, The Finite Difference Method in Partial Differential Equations (New York, John Wiley & Sons) In the text  
• Nakariakov, V. M., Ofman, L., DeLuca, E. E., Roberts, B., & Davila, J. M. 1999, Science, 285, 862 In the text  
• Nightingale, R. W., Aschwanden, M. J., & Hurlburt, N. E. 1999, Sol. Phys., 190, 249 In the text NASA ADS  
• Ofman, L., Nakariakov, V. M., & DeForest, C. E. 1999, ApJ, 514, 441 In the text NASA ADS  
• Ofman, L., Nakariakov, V. M., & Sehgal, N. 2000a, ApJ, 533, 1071 In the text NASA ADS  
• Ofman, L., Romoli, M., Poletto, G., Noci, G., & Kohl, J. L. 1997, ApJ, 491, L111 In the text NASA ADS  
• Ofman, L., Romoli, M., Poletto, G., Noci, G., & Kohl, J. L. 2000b, ApJ, 529, 592 In the text NASA ADS  
• Oliver, R., Ballester, J. L., Hood, A. W., & Priest, E. R. 1993, A&A, 273, 647 In the text NASA ADS  
• Oliver, R., Hood, A. W., & Priest, E. R. 1996, ApJ., 461, 424 In the text NASA ADS  
• Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in FORTRAN. The Art of Scientific Computing, 2nd edition (Cambridge University Press, Cambridge, USA), 487 In the text  
• Robbrecht, E., Berghmans, D., & Poedts, S. 1999, in SOHO 8, Plasma dynamics and diagnostics in the solar transition region and corona, ed. J. C. Vial, & B. Kaldeich-Schürman, European Space Agency, SP-446, 575 In the text  
• Roberts, B. 2000, Sol. Phys., 193, 139 In the text NASA ADS  
• Schrijver, C. J., & Brown, D. S. 2000, ApJ, 537, L69 In the text NASA ADS  
• Terradas, J., Oliver, R., & Ballester, J. L. 1999, ApJ, 517, 488 In the text NASA ADS

 

  
Appandix A: Derivation of the wave equations

In this appendix, a detailed derivation of the wave Eqs. (16) and (17) is presented. To make progress with Eq. (12) some expressions are needed. First, the curl of the equilibrium magnetic field is given by

$$% \nabla\times{\vec B}= B\frac{{\rm d}\mbox{$B_{y}$ }}{{\rm d}A}\mbox{$\hat{{\vec e}}_{\Vert}$ }.$$ (A.1)

The following vector identity is also useful

$$% \nabla\times(f{\vec a})=\nabla f\times {\vec a}+f\left(\nabla\times{\vec a}\right),$$ (A.2)

together with the formula for the gradient of a function f(x,y,z) and the divergence of a vector function ${\vec a}=a_n\, \mbox{$\hat{{\vec e}}_n$ }+a_\perp \, \mbox{$\hat{{\vec e}}_\perp$ }+a_\Vert\, \mbox{$\hat{{\vec e}}_{\Vert}$ }$

 

$$% \nabla f=\frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabl... ...ec B}_{\rm p}\!\cdot\!\nabla\right)$ }f\right)\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ (A.3)

 

$$% \nabla\cdot{\vec a} \frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot... ...}}{\mbox{$B^3_{\rm p}$ }} a_n + \frac{\mbox{$B_{\rm p}$ }}{B}\partial_y a_\perp$$ =

$$\frac{\mbox{$B_{y}$ }}{B\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}... ...eft(B\mbox{$B_{\rm s}$ }+\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right) a_\perp$$ -

$$\frac{\mbox{$B_{y}$ }}{B}\partial_y a_\Vert+\frac{1}{B}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }a_\Vert-\mbox{$B_{\rm s}$ }a_\Vert,$$ + (A.4)

where expressions for the derivatives of B and $\mbox{$B_{\rm p}$ }$ along and across magnetic field lines (namely $\mbox{$B_{\rm s}$ }$, $\mbox{$B_{\rm a}$ }$, $\mbox{$B_{\rm ps}$ }$, ...) can be found in Appendix B.

Next, the curl and divergence of the basis vectors can be determined,

  

$$% \nabla\times\mbox{$\hat{{\vec e}}_n$ } -\frac{\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }}{B}\,\mbox{$\hat{{... ...frac{\mbox{$B_{y}$ }\mbox{$B_{\rm ps}$ }}{B}\,\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (A.5)

$$\nabla\times\mbox{$\hat{{\vec e}}_\perp$ } \frac{1}{B}\left(BB_{\rm s}-\mbox{$B_{\rm p}$ }B_{\rm ps}\right)\... ...B\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }}\,\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (A.6)

$$\nabla\times\mbox{$\hat{{\vec e}}_{\Vert}$ } \frac{\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{\mbox{$B_{\rm p}$ }}\,\... ...$ } +\frac{{\rm d}\mbox{$B_{y}$ }}{{\rm d}A}\,\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (A.7)

$$% \nabla\cdot\mbox{$\hat{{\vec e}}_n$ } -\frac{B^3\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }},$$ = (A.8)

$$% \nabla\cdot\mbox{$\hat{{\vec e}}_\perp$ } \frac{\mbox{$B_{y}$ }}{B\mbox{$B_{\rm p}$ }}\left(B\mbox{$B_{\rm s}$ }+\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right)$$ = (A.9)

and

$$% \nabla\cdot\mbox{$\hat{{\vec e}}_{\Vert}$ }=-\mbox{$B_{\rm s}$ }.$$ (A.10)

Now, with the help of the above formulas, one can show that the normal and perpendicular components of Eq. (12) reduce to Eqs. (16) and (17).

  
Appendix B: Formulas for quantities in the wave Eqs. (16) and (17)

The following short-hand notation for the derivatives of the equilibrium (total and poloidal) magnetic field strength are used

  

$$% \mbox{$B_{\rm a}$ } \frac{1}{B^2}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }B,$$ = (B.1)

$$\mbox{$B_{\rm s}$ } \frac{1}{B^2}\mbox{$\left({\vec B}\cdot\nabla\right)$ }B,$$ = (B.2)

$$\mbox{$B_{\rm pa}$ } \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$B_{\rm p}$ },$$ = (B.3)

$$\mbox{$B_{\rm ps}$ } \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$B_{\rm p}$ },$$ = (B.4)

$$\mbox{$B_{\rm aa}$ } \frac{1}{B^3}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }^2 B,$$ = (B.5)

$$\mbox{$B_{\rm ss}$ } \frac{1}{B^3}\mbox{$\left({\vec B}\cdot\nabla\right)$ }^2 B,$$ = (B.6)

$$\mbox{$B_{\rm sa}$ } \frac{1}{B^3}\mbox{$\left({\vec B}\cdot\nabla\right)$ }\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }B,$$ = (B.7)

$$\mbox{$B_{\rm as}$ } \frac{1}{B^3}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$\left({\vec B}\cdot\nabla\right)$ }B,$$ = (B.8)

$$\mbox{$B_{\rm paa}$ } \frac{1}{\mbox{$B^3_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }^2 \mbox{$B_{\rm p}$ },$$ = (B.9)

$$\mbox{$B_{\rm pss}$ } \frac{1}{\mbox{$B^3_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }^2 \mbox{$B_{\rm p}$ },$$ = (B.10)

$$\mbox{$B_{\rm psa}$ } \frac{1}{\mbox{$B^3_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cd... ...ght)$ }\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$B_{\rm p}$ },$$ = (B.11)

$$% \mbox{$B_{\rm pas}$ } \frac{1}{\mbox{$B^3_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cd... ...t)$ }\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$B_{\rm p}$ }.$$ = (B.12)

The constants C₁ to C₄ are given by

$$\frac{\mbox{$B^2_{\rm p}$ }}{B^2}\left(\mbox{$B_{\rm pss}$ }-2\mb... ...g_1 +\frac{1}{B}g_{\rm 2a}-\frac{\mbox{$B_{\rm a}$ }}{\mbox{$B_{\rm p}$ }} g_2,$$ C₁ = (B.13)

$$-\frac{1}{B}g_{\rm 3a}-\frac{\mbox{$B_{\rm a}$ }}{\mbox{$B_{\rm p}$ }}g_3,$$ C₂ = (B.14)

$$\frac{\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }}{B}\frac{{\rm d}\mbox{$B_{y}$ }}{{\rm d}A}-\frac{\mbox{$B_{\rm p}$ }}{B^2} g_{\rm 1s}$$ C₃ =

$$+\frac{1}{B^2}\left(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox... ..._{\rm 2s} +\frac{\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B\mbox{$B_{\rm p}$ }} g_2,$$ (B.15)

$$\mbox{$B_{\rm ss}$ }-\mbox{$B^2_{\rm s}$ }-\frac{\mbox{$B^2_{\rm p}$ }}{B^2}\left(\mbox{$B_{\rm pss}$ }-\mbox{$B^2_{\rm ps}$ }\right)$$ C₄ =

$$-\frac{1}{B^2}\left(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right)^2$$

$$+\frac{\mbox{$B_{y}$ }}{B^2} g_{\rm 3s}-\frac{\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B\mbox{$B_{\rm p}$ }} g_3,$$ (B.16)

where the following short-hand notation has been used

$$B\frac{{\rm d}\mbox{$B_{y}$ }}{{\rm d}A}-\frac{2\mbox{$B_{y}$ }B^2\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }},$$ g₁ = (B.17)

$$\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)B^2\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }},$$ g₂ = (B.18)

$$\frac{\mbox{$B_{y}$ }}{\mbox{$B_{\rm p}$ }}\left(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right),$$ g₃ = (B.19)

$$% g_{\rm 1s}=\frac{B}{\mbox{$B_{\rm p}$ }}\frac{{\rm d}\mbox{$B_{... ...eft(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right)\right],$$ (B.20)

$$% g_{\rm 2a}= \frac{\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }}{\mbo... ...box{$B_{\rm p}$ }\mbox{$B_{y}$ }\frac{{\rm d}\mbox{$B_{y}$ }}{{\rm d}A}\right),$$ (B.21)

$$% g_{\rm 2s} \frac{\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }}{\mbox{$B^2_{\rm p}... ...ht] + \frac{B^4\mbox{$B_{\rm a}$ }\mbox{$B_{\rm ps}$ }}{\mbox{$B^3_{\rm p}$ }},$$ (B.22)

$$% g_{\rm 3a}= \frac{\mbox{$B_{y}$ }}{\mbox{$B^2_{\rm p}$ }}\left[... ...rm d}A}-\frac{\mbox{$B_{y}$ }\mbox{$B_{\rm pa}$ }}{\mbox{$B_{\rm p}$ }}\right),$$ (B.23)

$$% g_{\rm 3s}= \frac{\mbox{$B_{y}$ }}{\mbox{$B^2_{\rm p}$ }}\left[... ...}$ }}\left(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right).$$ (B.24)

  
Appendix C: Expressions for the perturbed variables and forces

From the perturbed induction Eq. (11), the normal, perpendicular and parallel components of the perturbed magnetic field are

$$% i\omega \frac{B_{1n}}{B}=\frac{1}{B}\mbox{$\left({\vec B}_... ...box{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }}{B}\right)\mbox{$v_n$ },$$ (C.1)

$$% i\omega \frac{B_{1\perp}}{B} \frac{1}{B}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }... ...\mbox{$B_{y}$ }B\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}\right)\mbox{$v_n$ }$$ = (C.2)

and

$$% i\omega \frac{B_{1\Vert}}{B} \frac{\mbox{$B_{y}$ }}{B\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}... ...B^2_{y}$ }\right)B\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }}\mbox{$v_n$ }\cdot$$ = (C.3)

In addition, Eqs. (9) and (13) for the perturbed density and magnetic pressure together with Eqs. (A.3), (A.4) and (A.8)-(A.10) yield

$$% i\omega\frac{\rho_{1}}{\rho} -\frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdo... ...rm s}$ }+\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right)\right]\mbox{$v_\perp$ }$$ = (C.4)

and

 

$$% i\omega \frac{p_{\rm 1m}}{p_{\rm m}} -\frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdo... ..._{\rm ps}$ }\right) -\frac{\mbox{$B_{\rm p}$ }}{B}ik_y\right]\mbox{$v_\perp$ }.$$ = (C.5)

where $p_{\rm m}$ is the equilibrium magnetic pressure and $\rho_{\rm a}$ and $\rho_{\rm s}$are

$$% \rho_{\rm a} \frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\rho,$$ = (C.6)

$$\rho_{\rm s} \frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\rho.$$ = (C.7)

Now, Eqs. (C.5) and (A.3) can be used to determine the three components of the magnetic pressure gradient force (symbolised here as $\left[\nabla p_{\rm 1m}\right]_{n}$, $\left[\nabla p_{\rm 1m}\right]_{\perp}$ and $\left[\nabla p_{\rm 1m}\right]_{\Vert}$),

 

$$% -i\frac{\omega}{\mbox{$v^2_{\rm A}$ }\rho}\left[\nabla p_{\rm 1m}\right]_{n} \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }^2\mbox{$v_n$ }$$ =

$$+\left[\frac{1}{\mbox{$B_{\rm p}$ }}\left(B\mbox{$B_{\rm a}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm pa}$ }\right)\right.$$

$$+\left.\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)... ... p}$ }}\right]\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$

$$+ h_{1}\mbox{$v_n$ }-\frac{\mbox{$B_{y}$ }}{B}\frac{1}{\mbox{$B^2... ...ight)$ }\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_\perp$ }$$

$$+\left[\frac{\mbox{$B_{\rm p}$ }}{B}ik_y-\frac{\mbox{$B_{y}$ }}{B... ...\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_\perp$ }$$

$$-\left[\frac{\mbox{$B_{\rm p}$ }}{B}\frac{{\rm d}\mbox{$B_{y}$ }}... ...m p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_\perp$ }$$

$$+ \left[\frac{1}{B}\left(B\mbox{$B_{\rm a}$ }+\mbox{$B_{\rm p}$ }\mbox{$B_{\rm pa}$ }\right)ik_y+h_{2}\right]\mbox{$v_\perp$ },$$ (C.8)

 

$$% -i\frac{\omega}{\mbox{$v^2_{\rm A}$ }\rho}\left[\nabla p_{\rm 1m}\right]_{\perp} -\frac{\mbox{$B_{y}$ }} {B\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\ve... ...abla\right)$ }\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$ =

$$-\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)\mbox{... ..._{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_n$ }$$

$$+ \left[\frac{1}{B}\ ik_y-\frac{\mbox{$B_{y}$ }}{B\mbox{$B^2_{\rm... ...\right)\right]\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$

$$+ \left[\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}ik_y + h_{3}\right]\mbox{$v_n$ }$$

$$+\frac{\mbox{$B^2_{y}$ }}{B^2\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }^2\mbox{$v_\perp$ }$$

$$+ 2 \left[ -\frac{\mbox{$B_{\rm p}$ }\mbox{$B_{y}$ }}{B^2}ik_y\right.$$

$$+\left.\frac{\mbox{$B^2_{y}$ }}{B^2\mbox{$B_{\rm p}$ }}\left(B\mb... ...m p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_\perp$ }$$

$$+\left(-2\frac{\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B}ik_y+h_{4}-\frac{\mbox{$B^2_{\rm p}$ }}{B^2}k_y^2\right)\mbox{$v_\perp$ },$$ (C.9)

 

$$% -i\frac{\omega}{\mbox{$v^2_{\rm A}$ }\rho}\left[\nabla p_{\rm 1m}\right]_{\Vert} \frac{1} {B\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cd... ...abla\right)$ }\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$ =

$$+\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)\mbox{... ..._{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_n$ }$$

$$+ \frac{1}{B}\left[\mbox{$B_{y}$ }ik_y\right.$$

$$+\left.\left(2B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\... ...$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$

$$+ \left[\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right... ...rm a}$ }\mbox{$B_{y}$ }}{\mbox{$B^3_{\rm p}$ }}ik_y + h_{5}\right]\mbox{$v_n$ }$$

$$-\frac{\mbox{$B_{y}$ }}{B^2\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }^2\mbox{$v_\perp$ }$$

$$+\left[ \frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)} {B^2}ik_y\right.$$

$$-\left.2\frac{\mbox{$B_{y}$ }}{B^2}\left(B\mbox{$B_{\rm s}$ }-\mb... ...m p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_\perp$ }$$

$$-\left[-\frac{\mbox{$B_{\rm p}$ }\mbox{$B_{y}$ }}{B^2}ik_y+\frac{... ...eft(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right) \right.$$

$$-\left.\frac{\mbox{$B_{\rm p}$ }}{B^2} \left(B\mbox{$B_{\rm s}$ }... ...mbox{$B_{\rm ps}$ }\right)\right]ik_y\mbox{$v_\perp$ } +h_{6}\mbox{$v_\perp$ }.$$ (C.10)

Finally, the calculation of the perturbed magnetic tension from Eq. (14) requires the derivatives of $\mbox{$\hat{{\vec e}}_n$ }$, $\mbox{$\hat{{\vec e}}_\perp$ }$ and $\mbox{$\hat{{\vec e}}_{\Vert}$ }$ in the directions of these three vectors. These are

  

$$% \mbox{$\left(\mbox{$\hat{{\vec e}}_n$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_n$ } -\frac{B_{y}\mbox{$B_{\rm ps}$ }}{B}\mbox{$\hat{{\vec e}}_\perp$ ... ...ac{\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }}{B}\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (C.11)

$$\mbox{$\left(\mbox{$\hat{{\vec e}}_n$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_\perp$ } \frac{\mbox{$B_{y}$ }\mbox{$B_{\rm ps}$ }}{B}\mbox{$\hat{{\vec e}... ...x{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}\right)\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (C.12)

$$\mbox{$\left(\mbox{$\hat{{\vec e}}_n$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_{\Vert}$ } -\frac{\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }}{B}\mbox{$\hat{{\v... ...box{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}\right)\mbox{$\hat{{\vec e}}_\perp$ },$$ = (C.13)

$$\mbox{$\left(\mbox{$\hat{{\vec e}}_\perp$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_n$ } -\frac{\mbox{$B^2_{y}$ }B\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$... ... }B\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (C.14)

$$\mbox{$\left(\mbox{$\hat{{\vec e}}_\perp$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_\perp$ } \frac{\mbox{$B^2_{y}$ }B\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ ... ...\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right)\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (C.15)

$$\mbox{$\left(\mbox{$\hat{{\vec e}}_\perp$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_{\Vert}$ } -\frac{\mbox{$B_{y}$ }B\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }... ...}-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right)\mbox{$\hat{{\vec e}}_\perp$ },$$ = (C.16)

$$\mbox{$\left(\mbox{$\hat{{\vec e}}_{\Vert}$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_n$ } \frac{\mbox{$B_{y}$ }B\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}... ...rac{B\mbox{$B_{\rm a}$ }}{\mbox{$B_{\rm p}$ }}\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (C.17)

$$\mbox{$\left(\mbox{$\hat{{\vec e}}_{\Vert}$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_\perp$ } -\frac{\mbox{$B_{y}$ }B\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }... ...y}$ }\mbox{$B_{\rm s}$ }}{\mbox{$B_{\rm p}$ }}\mbox{$\hat{{\vec e}}_{\Vert}$ },$$ = (C.18)

$$% \mbox{$\left(\mbox{$\hat{{\vec e}}_{\Vert}$ }\cdot\nabla\right)$ }\mbox{$\hat{{\vec e}}_{\Vert}$ } \frac{B\mbox{$B_{\rm a}$ }}{\mbox{$B_{\rm p}$ }}\mbox{$\hat{{\vec... ...{y}\mbox{$B_{\rm s}$ }}{\mbox{$B_{\rm p}$ }}\mbox{$\hat{{\vec e}}_\perp$ }\cdot$$ = (C.19)

Then, the three components of the perturbed magnetic tension force are

 

$$% i\frac{\omega}{\mbox{$v^2_{\rm A}$ }\rho}\mbox{${\bf {\cal T}}_1$ }_{n} \frac{1}{B^2}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$... ...^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$ =

$$+ \left(h_{7}-\frac{\mbox{$B^2_{y}$ }}{B^2}k_{y}^2\right)\mbox{$v... ...ac{2B^2\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}ik_y\right)\mbox{$v_\perp$ },$$ (C.20)

 

$$% i\frac{\omega}{\mbox{$v^2_{\rm A}$ }\rho}\mbox{${\bf {\cal T}}_1$ }_{\perp} \frac{2\mbox{$B_{y}$ }\mbox{$B_{\rm a}$ }} {\mbox{$B^2_{\rm p}$ }... ...^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$ =

$$+\left(\frac{2\mbox{$B^2_{y}$ }\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }} ik_y+h_{9}\right)\mbox{$v_n$ }$$

$$+\frac{1}{B^2}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }^2\mbox{$v_\perp$ }$$

$$+\left(\frac{2\mbox{$B_{y}$ }}{B^2}ik_y-\frac{2\mbox{$B^2_{y}$ }\... ... \right)\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_\perp$ }$$

$$+ \left(\frac{2\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B}ik_y+h_{10}-\frac{\mbox{$B^2_{y}$ }}{B^2} k_y^2\right)\mbox{$v_\perp$ },$$ (C.21)

 

$$% i\frac{\omega}{\mbox{$v^2_{\rm A}$ }\rho}\mbox{${\bf {\cal T}}_1$ }_{\Vert} -\frac{1} {B\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\c... ...abla\right)$ }\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$ =

$$-\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)\mbox{... ..._{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_n$ }$$

$$- \frac{1}{B}\left[\mbox{$B_{y}$ }ik_y\right.$$

$$+\left.\left(2B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\... ...$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$v_n$ }$$

$$- \left[\frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right... ...m a}$ }\mbox{$B_{y}$ }}{\mbox{$B^3_{\rm p}$ }}ik_y + h_{11}\right]\mbox{$v_n$ }$$

$$+\frac{\mbox{$B_{y}$ }}{B^2\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }^2\mbox{$v_\perp$ }$$

$$-\left[ \frac{\left(\mbox{$B^2_{\rm p}$ }-\mbox{$B^2_{y}$ }\right)} {B^2}ik_y\right.$$

$$-\left.2\frac{\mbox{$B_{y}$ }}{B^2}\left(B\mbox{$B_{\rm s}$ }-\mb... ...m p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$v_\perp$ }$$

$$+ \left[-\frac{\mbox{$B_{\rm p}$ }\mbox{$B_{y}$ }}{B^2}ik_y+\frac... ...eft(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }\right) \right.$$

$$-\left.\frac{\mbox{$B_{\rm p}$ }}{B^2} \left(B\mbox{$B_{\rm s}$ }... ...box{$B_{\rm ps}$ }\right)\right]ik_y\mbox{$v_\perp$ } +h_{12}\mbox{$v_\perp$ }.$$ (C.22)

The coefficients h₁ to h₁₂ are given by

$$\frac{1}{B}g_{\rm 2a}+\frac{\mbox{$B_{\rm a}$ }}{\mbox{$B_{\rm p}$ }} g_2,$$ h₁ = (C.23)

$$-\frac{1}{B}g_{\rm 3a}-\frac{\mbox{$B_{\rm a}$ }}{\mbox{$B_{\rm p}$ }}g_3,$$ h₂ = (C.24)

$$-\frac{\mbox{$B_{y}$ }}{B^2} g_{\rm 2s} -\frac{\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B\mbox{$B_{\rm p}$ }} g_2,$$ h₃ = (C.25)

$$\frac{\mbox{$B_{y}$ }}{B^2} g_{\rm 3s}+\frac{\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B\mbox{$B_{\rm p}$ }} g_3,$$ h₄ = (C.26)

$$\frac{\mbox{$B_{\rm p}$ }}{B^2} g_{\rm 2s}+\frac{\mbox{$B_{\rm s}$ }}{B} g_2,$$ h₅ = (C.27)

$$-\frac{\mbox{$B_{\rm p}$ }}{B^2} g_{\rm 3s}-\frac{\mbox{$B_{\rm s}$ }}{B}g_3,$$ h₆ = (C.28)

$$\frac{\mbox{$B^2_{\rm p}$ }}{B^2}\left(\mbox{$B_{\rm pss}$ }-2\mb... ...\mbox{$B^2_{\rm p}$ }} g_1-\frac{2\mbox{$B_{\rm a}$ }}{\mbox{$B_{\rm p}$ }}g_2,$$ h₇ = (C.29)

h₈ = 0, (C.30)

$$\frac{\mbox{$B_{\rm p}$ }\mbox{$B_{\rm ps}$ }}{B}\frac{{\rm d}\mbox{$B_{y}$ }}{{\rm d}A}-\frac{\mbox{$B_{\rm p}$ }}{B^2} g_{\rm 1s}$$ h₉ =

$$+\frac{1}{B^2}\left(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox... ..._{\rm 2s} +\frac{\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B\mbox{$B_{\rm p}$ }} g_2,$$ (C.31)

$$\frac{1}{B^2}\left[B^2\left(\mbox{$B_{\rm ss}$ }-\mbox{$B^2_{\rm ... ...$B^2_{\rm p}$ }\left(\mbox{$B_{\rm pss}$ }-\mbox{$B^2_{\rm ps}$ }\right)\right]$$ h₁₀ =

$$-\frac{1}{B^2}\left(B\mbox{$B_{\rm s}$ }-\mbox{$B_{\rm p}$ }\mbox... ...\right)^2 -\frac{2\mbox{$B_{y}$ }\mbox{$B_{\rm s}$ }}{B\mbox{$B_{\rm p}$ }}g_3,$$ (C.32)

$$-\frac{\mbox{$B_{\rm p}$ }}{B^2}g_{\rm 2s}-\frac{\mbox{$B_{\rm s}$ }}{B}g_2,$$ h₁₁ = (C.33)

$$\frac{\mbox{$B_{\rm p}$ }}{B^2}g_{\rm 3s} +\frac{\mbox{$B_{\rm s}$ }}{\mbox{$B_{\rm p}$ }}g_3.$$ h₁₂ = (C.34)

  
Appendix D: Formulas related to the numerical solution of the wave equations

The following expressions are the directional derivatives written in terms of $\partial /\mbox{$\partial\psi$ }$ and $\partial /\mbox{$\partial\chi$ }$, with $\psi$ and $\chi$a pair of independent coordinates.

$$% \frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ } \mbox{$\psi_{\rm a}$ }\frac{\partial}{\mbox{$\partial\psi$ }}+\mbox{$\chi_{\rm a}$ }\frac{\partial}{\mbox{$\partial\chi$ }},$$ = (D.1)

$$\frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ } \mbox{$\psi_{\rm s}$ }\frac{\partial}{\mbox{$\partial\psi$ }}+\mbox{$\chi_{\rm s}$ }\frac{\partial}{\mbox{$\partial\chi$ }},$$ = (D.2)

$$% \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\... ...partial\psi$ }}+\mbox{$\chi_{\rm aa}$ }\frac{\partial}{\mbox{$\partial\chi$ }},$$ (D.3)

$$% \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\... ...artial\psi$ }} +\mbox{$\chi_{\rm as}$ }\frac{\partial}{\mbox{$\partial\chi$ }},$$ (D.4)

$$% \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\... ...artial\psi$ }} +\mbox{$\chi_{\rm sa}$ }\frac{\partial}{\mbox{$\partial\chi$ }},$$ (D.5)

 

$$% \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\... ...ial\psi$ }}+\mbox{$\chi_{\rm ss}$ }\frac{\partial}{\mbox{$\partial\chi$ }}\cdot$$ (D.6)

The quantities $\mbox{$\psi_{\rm a}$ }$, $\mbox{$\psi_{\rm s}$ }$, etc. have been defined as

$$% \mbox{$\psi_{\rm a}$ } \frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\psi,$$ = (D.7)

$$\mbox{$\psi_{\rm s}$ } \frac{1}{\mbox{$B_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\psi,$$ = (D.8)

$$\mbox{$\psi_{\rm aa}$ } \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }^2\psi,$$ = (D.9)

$$\mbox{$\psi_{\rm as}$ } \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\psi,$$ = (D.10)

$$\mbox{$\psi_{\rm sa}$ } \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }\mbox{$\left(\mbox{$\nabla A$ }\cdot\nabla\right)$ }\psi,$$ = (D.11)

$$\mbox{$\psi_{\rm ss}$ } \frac{1}{\mbox{$B^2_{\rm p}$ }}\mbox{$\left({\vec B}_{\rm p}\!\cdot\!\nabla\right)$ }^2\psi$$ = (D.12)

and similarly for the derivatives of $\chi$.

Finally, the expressions for the coefficients in Eqs. (21) and (22) are given by

$$% a_{11}=-\mbox{$v^2_{\rm A}$ }\left(\mbox{$\psi_{\rm a}^2$ ... ...ac{\mbox{$B^2_{\rm p}$ }}{B^2}\mbox{$\psi_{\rm s}^2$ }\right),$$ (D.13)

$$% b_{11}=-2\mbox{$v^2_{\rm A}$ }\left(\mbox{$\psi_{\rm a}$ }\... ...}$ }}{B^2}\mbox{$\psi_{\rm s}$ }\mbox{$\chi_{\rm s}$ }\right),$$ (D.14)

$$% c_{11}=-\mbox{$v^2_{\rm A}$ }\left(\mbox{$\chi_{\rm a}^2$ }... ...ac{\mbox{$B^2_{\rm p}$ }}{B^2}\mbox{$\chi_{\rm s}^2$ }\right),$$ (D.15)

$$d_{11}=-\mbox{$v^2_{\rm A}$ }\left\{\mbox{$\psi_{\rm aa}$ }+\frac... ...box{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }}\right]\mbox{$\psi_{\rm a}$ }\right\},$$ (D.16)

$$e_{11}=-\mbox{$v^2_{\rm A}$ }\left\{\mbox{$\chi_{\rm aa}$ }+\frac... ...ox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }} \right]\mbox{$\chi_{\rm a}$ }\right\},$$ (D.17)

$$% f_{11}=-\mbox{$v^2_{\rm A}$ }\left(C_{1}-\frac{\mbox{$B^2_{y}$ }}{B^2}k_{y}^2\right),$$ (D.18)

$$% a_{12}=\mbox{$v^2_{\rm A}$ }\frac{\mbox{$B_{y}$ }}{B}\mbox{$\psi_{\rm a}$ }\mbox{$\psi_{\rm s}$ },$$ (D.19)

$$% b_{12}=\mbox{$v^2_{\rm A}$ }\frac{\mbox{$B_{y}$ }}{B}\left(... ...\rm s}$ }+\mbox{$\psi_{\rm s}$ }\mbox{$\chi_{\rm a}$ }\right),$$ (D.20)

$$% c_{12}=\mbox{$v^2_{\rm A}$ }\frac{\mbox{$B_{y}$ }}{B}\mbox{$\chi_{\rm a}$ }\mbox{$\chi_{\rm s}$ },$$ (D.21)

$$d_{12}=-\mbox{$v^2_{\rm A}$ }\left\{-\frac{\mbox{$B_{y}$ }}{B}\mb... ...{$B_{\rm p}$ }\mbox{$B_{\rm pa}$ }\right)\right]\mbox{$\psi_{\rm s}$ }\right\},$$ (D.22)

$$e_{12}=-\mbox{$v^2_{\rm A}$ }\left\{-\frac{\mbox{$B_{y}$ }}{B}\mb... ...{$B_{\rm p}$ }\mbox{$B_{\rm pa}$ }\right)\right]\mbox{$\chi_{\rm s}$ }\right\},$$ (D.23)

$$f_{12}=-\mbox{$v^2_{\rm A}$ }\left\{C_{2}-\left[\frac{1}{B}\left(... ...x{$B^2_{y}$ }\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }}\right]i k_{y}\right\},$$ (D.24)

$$% a_{21}=\mbox{$v^2_{\rm A}$ }\frac{\mbox{$B_{y}$ }}{B}\mbox{$\psi_{\rm a}$ }\mbox{$\psi_{\rm s}$ }=a_{12},$$ (D.25)

$$% b_{21}=\mbox{$v^2_{\rm A}$ }\frac{\mbox{$B_{y}$ }}{B}\left(... ... }+\mbox{$\psi_{\rm s}$ }\mbox{$\chi_{\rm a}$ }\right)=b_{12},$$ (D.26)

$$% c_{21}=\mbox{$v^2_{\rm A}$ }\frac{\mbox{$B_{y}$ }}{B}\mbox{$\chi_{\rm a}$ }\mbox{$\chi_{\rm s}$ }=c_{12},$$ (D.27)

$$d_{21}=-\mbox{$v^2_{\rm A}$ }\left[-\frac{\mbox{$B_{y}$ }}{B}\mbo... ...$ }B^2\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }}\mbox{$\psi_{\rm s}$ }\right],$$ (D.28)

$$-\mbox{$v^2_{\rm A}$ }\left[-\frac{\mbox{$B_{y}$ }}{B}\mbox{$\chi... ...$ }B^2\mbox{$B_{\rm a}$ }}{\mbox{$B^3_{\rm p}$ }}\mbox{$\chi_{\rm s}$ }\right],$$ (D.29)

$$% f_{21}=-\mbox{$v^2_{\rm A}$ }\left(C_{3}+\frac{B^2\mbox{$B_{\rm a}$ }}{\mbox{$B^2_{\rm p}$ }} ik_{y}\right),$$ (D.30)

$$% a_{22}=-\mbox{$v^2_{\rm A}$ }\mbox{$\psi_{\rm s}^2$ },$$ (D.31)

$$% b_{22}=-2\mbox{$v^2_{\rm A}$ }\mbox{$\psi_{\rm s}$ }\mbox{$\chi_{\rm s}$ },$$ (D.32)

$$% c_{22}=-\mbox{$v^2_{\rm A}$ }\mbox{$\chi_{\rm s}^2$ },$$ (D.33)

$$% d_{22}=-\mbox{$v^2_{\rm A}$ }\left(\mbox{$\psi_{\rm ss}$ }-... ...{y}$ }\mbox{$B_{\rm ps}$ }}{B^2}\mbox{$\psi_{\rm s}$ }\right),$$ (D.34)

$$% e_{22}=-\mbox{$v^2_{\rm A}$ }\left(\mbox{$\chi_{\rm ss}$ }-... ...{y}$ }\mbox{$B_{\rm ps}$ }}{B^2}\mbox{$\chi_{\rm s}$ }\right),$$ (D.35)

$$% f_{22}=-\mbox{$v^2_{\rm A}$ }\left(C_{4}-k_{y}^2\right).$$ (D.36)